Week 3 - Taylor table method and Matlab code
Aim:
Derive the Central difference, Skewed right-sided difference, Skewed left-sided difference scheme for 4th order approximation of the second-order derivative(f'') and we have to prove the Skewed difference has 4th order accuracy.
Then, write a program in Matlab to evaluate the second-order derivative of the analytical function exp(*cos(x) and compare it with the 3 numerical approximations that you have derived and plot a graph that compares the absolute error between the above-mentioned schemes. Explain the code and analyse the plot.x)
Finally, explain the necessity of Skewed scheme approximation.
Main Code:
%Computation of Higher order approximation using symmetric and skewed scheme
%The given function f(x) = exp(x)*cos(x)
close all
clear all
clc
%Assuming suitable value for x and for the range of dx
x = pi/3;
dx = linspace (0.01,1,60);
%Declaring the size of the variable
fourthorder_central_sol = zeros (1,length(dx));
fourthorder_sright_sol = zeros (1,length(dx));
fourthorder_sleft_sol = zeros (1,length(dx));
secondorder_forward_sol = zeros(1,length(dx));
secondorder_backward_sol = zeros(1,length(dx));
%loop to compute the various numericaal solutions for each vale of dx
for i = 1:length(dx)
%f(x-5*dx)
f5l = exp(x-5*dx(i))*cos(x-5*dx(i));
%f(x-4*dx)
f4l = exp(x-4*dx(i))*cos(x-4*dx(i));
%f(x-3*dx)
f3l = exp(x-3*dx(i))*cos(x-3*dx(i));
%f(x-2*dx)
f2l = exp(x-2*dx(i))*cos(x-2*dx(i));
%f(x-dx)
f1l = exp(x-dx(i))*cos(x-dx(i));
%f(x)
fm = exp(x)*cos(x);
%f(x+dx)
f1r = exp(x+dx(i))*cos(x+dx(i));
%f(x+2*dx)
f2r = exp(x+2*dx(i))*cos(x+2*dx(i));
%f(x+3*dx)
f3r = exp(x+3*dx(i))*cos(x+3*dx(i));
%f(x+4*dx)
f4r = exp(x+4*dx(i))*cos(x+4*dx(i));
%f(x=5*dx)
f5r = exp(x+5*dx(i))*cos(x+5*dx(i));
%Computing f(x) by calling their respective function for numerical solution
fourthorder_central_sol(i) = forder_central_diff(x,f2l, f1l, fm, f1r, f2r, dx(i));
fourthorder_sright_sol(i) = forder_eskewedright_diff (x,fm,f1r,f2r,f3r,f4r,f5r,dx(i));
fourthorder_sleft_sol(i) = forder_skewedleft_diff (x,f5l,f4l,f3l,f2l,f1l,fm,dx(i));
secondorder_forward_sol(i) = secondorder_forward_diff (x,fm,f1r,f2r,dx(i));
secondorder_backward_sol(i) = secondorder_backward_diff (x,f2l,f1l,fm,dx(i));
end
% Ploting the graph between dx and various scheme
figure(1)
loglog (dx,fourthorder_central_sol,'m','linewidth',1.5)
hold on
loglog (dx,fourthorder_sright_sol,'b','linewidth',1.5)
hold on
loglog (dx,fourthorder_sleft_sol,'g','linewidth',1.5)
hold on
loglog (dx,secondorder_forward_sol,'+r','linewidth',1.5)
hold on
loglog (dx,secondorder_backward_sol,'--c','linewidth',1.5)
title('Fourth order approximation for second order derivative')
grid on
xlabel(' "Value of dx"')
ylabel('Value of Error for x = pi/3')
legend('central difference error','Skewed Right difference error','Skewed Left difference error', 'Forward differencce','Backward difference')
Derivation of fourth order approximation for the central difference scheme of second order derivative:
First we have to find the number of points we need to consider in the Numerical stencil, for that we can use the formula N = A + O - 1 for symmetric or central difference,
where N = Number of points, A = Order of Approximation, O= Order of derivative
N = 4 + 2 -1 = 5, SO we need to take five points in the Numerical Stencil, symmetric from f(x).
Which gives as f''(x) = a*f(x-2dx) + b*f(x-dx) + c*f(x) + d*f(x+dx) + e*f(x+2dx), all we have to find the co-efficient.
For that we have to Taylor expand the above equation and using Taylor table we can find the coefficient and the desire formula for calculating f''(x).
| f''(x) | f(x) | f'(x)(dx) | f''(x)(dx^2/2!) | f'''(x)(dx^3/3!) | f""(x)(dx^4/4!) | f""'(x)(dx^5/5!) | f"""(x)(dx^6/6!) |
| a*f(x-2*dx) | a | -2a | 4a | -8a | 16a | -32a | 64a |
|
b*f(x-dx) |
b | -b | b | -b | b | -b | b |
| c*f(x) | c | 0 | 0 | 0 | 0 | 0 | 0 |
| d*f(x+dx) | d | d | d | d | d | d | d |
| e*f(x+2*dx) | e | 2e | 4e | 8e | 16e | 32e | 64e |
| ∑row | 0 | 0 | 1 | 0 | 0 | ? | ? |
From the above attached image we can see the Taylor expansion for all the necessary therms and written in the form of Taylor Table . From the Taylor table we can say the sum of coefficient for the f'' (x)(dx)^2 is 1 and every other co-efficient with-in the number of points we took in the stencil is zero, So we can come up with 5 equations with 5 unknown. The equation is written below
equ_1 = a+b+c+d+e = 0
equ_2 = -2*a-b+d+2*e = 0
equ_3 = (4*a+b+d+4*e)/2=1
equ_4 = -8*a -b+0+d+8*e= 0
equ_5 = 16*a+b+0+d+16*e = 0
When we got this equation, basically, we can find the coefficient by Inverse matrix method or simply solve the equation in the Matlab. The are several ways to solve using different functions like
"X= linsolve (A,B);" or "X=SOLVE ([equ_1, equ_2,equ_3,equ_4,equ_5], [a,b,c,d,e]);" or "X=inv(A)*B;".
[A] = [ 1 1 1 1 1 ;
-2 -1 0 1 2 ;
4 1 0 1 4 ;
-8 -1 0 1 8 ;
16 1 0 1 16 ];
[B] = [0;
0;
2;
0;
0];
X= linsolve (A,B)
[a: [-0.0833 :
b: 1.3333 :
X = c: = -2.5000 :
d: 1.3333 :
e]: -0.0833 ]:
Therefore, the fourth order right skewed approximation for the second order derivative is
(d^2f/dx^2) = (-0.0833*f(x-2*dx) + 1.3333*f(x-dx) - 2.500*f(x) +1.333*f(x+dx) - 0.0833*f(x+2*dx))/(dx^2)
Function to calculate Fourth order Central approximation:
% Function to calculate Central difference fourth order approximation
function out = forder_central_diff (x,f2l,f1l,fm,f1r,f2r,dx)
%given analytical function exp(x)*cos(x)
%analytical_derivative_first_order f'(x)= exp(x)*cos(x)- exp(x)*sin(x)
%analytical derivative second order f''(x) = -2*exp(x)*sin(x)
analytical_solution = -2*exp(x)*sin(x);
%For the numerical derivative. We need to slove taylor table to find the co-efficient for each derivative.
%From taylor's table we can get the below equation
%equ_1 = a + b + c + d + e = 0;
%equ_2 = -2*a - b + d + 2*e = 0;
%equ_3 = 4*a + b + d + 4*e = 2;
%equ_4 = -8*a -b + 0+ d + 8*e = 0;
%equ_5 = 16*a + b + 0 + d + 16*e = 0;
%We can represent the equation in the form of a matrix AX = B
%Matrix A
[A] = [ 1 1 1 1 1 ;
-2 -1 0 1 2 ;
4 1 0 1 4 ;
-8 -1 0 1 8 ;
16 1 0 1 16 ];
% Matrix B
[B] = [0;
0;
2;
0;
0];
%Solving the equation for the constant
X = linsolve(A,B);
%Extracting the values of co-efficient from the solution
a = X(1);
b = X(2);
c = X(3);
d = X(4);
e = X(5);
%Central diffrence sceme fourth order approximation
fourthorder_central_app = (a*f2l + b*f1l + c*fm + d*f1r + e*f2r)/(dx^2);
out = (analytical_solution - fourthorder_central_app);
endHow do we know that the above scheme is fourth order:
So, from the way we took the discriminate points in the stencil itself shows that the order of accuracy is fourth order (i.e. N = A + O for symmetric points).
However, we can prove this by Expanding the Taylor's table with considering the higher order terms.
We can clearly see that the summation of terms for f(x),f''(x),f""(x) are 0, So considering only the remaining terms
d^2f/dx^2 = (4a+b+d+4e) f"(x)(dx^2/2!) + (-32a - b +d +32e)f""'(x)(dx^5/5!) + (64a +b+b+64a)f"""(x)(dx^6/6!)+.....
We know that the co_efficient of f"(x) = 1 and we can clearly see the co-efficient a=e and b=d, So the f""'(x) becomes zero (0).
Substituting the the values of d^2f/dx^2 = a*f(x-2dx) + b*f(x-dx) + c*f(x) + d*f(x+dx) + e*f(x+2dx) from the Taylor's table
a*f(x-2dx) + b*f(x-dx) + c*f(x) + d*f(x+dx) + e*f(x+2dx) = (4a+b+d+4e) f"(x)(dx^2/2!) + (64a +b+b+64a)f"""(x)(dx^6/6!) +.....
to get f"(x) on the right side, we have to divide the both side by dx^2
a*f(x-2dx) + b*f(x-dx) + c*f(x) + d*f(x+dx) + e*f(x+2dx)/dx^2 = (4a+b+d+4e) f"(x)(1/2!) + (64a +b+b+64a)f"""(x)(dx^4/6!) +......
In the above equation we can see that the first dx term has the power of 4 , So, if we truncate all those terms from dx^4, we get
(-0.0833*f(x-2*dx) + 1.3333*f(x-dx) - 2.500*f(x) +1.333*f(x+dx) - 0.0833*f(x+2*dx))/(dx^2) = (4a+b+d+4e) f"(x)(1/2!) + error in Order of (dx^4) or O(dx^4)
Hence, we proved that the scheme we come up with have an error of order (dx^4)
Derivation of fourth order approximation for the Skewed right side difference scheme of second order derivative:
First we have to find the number of points we need to consider in the Numerical stencil, for that we can use the formula N = A + O for Skewed scheme,
where N = Number of points, A = Order of Approximation, O= Order of derivative
N = 4 + 2 = 6, SO we need to take Six points in the Numerical Stencil, towards right from f(x).
Which gives as f''(x) = af(x) + bf(x+dx) + cf(x+2dx)+ df(x+3dx) + ef(x+4dx) + gf(x+5dx), all we have to find the co-efficient.
For that we have to Taylor expand the above equation and using Taylor table we can find the coefficient and desire formula for calculating f''(x).
| f''(x) | f(x) | f'(x)(dx) | f''(x)(dx^2/2!) | f'''(x)(dx^3/3!) | f''''(x)(dx^4/4!) | f""'(x)(dx^5/5!) | f"""(x)(dx^6/6!) |
| a*f(x) | a | 0 | 0 | 0 | 0 | 0 | 0 |
| b*f(x+dx) | b | b | b | b | b | b | b |
| c*f(x+2*dx) | c | 2c | 4c | 8c | 16c | 32c | 64c |
| d*f(x+3*dx) | d | 3d | 9d | 27d | 81d | 243d | 729d |
| e*f(x+4*dx) | e | 4e | 16e | 64e | 256e | 1024e | 4096e |
| g*f(x+5*dx) | g | 5g | 25g | 125g | 625g | 3125g | 15625g |
| ∑rows | 0 | 0 | 1 | 0 | 0 | 0 | ? |
From the above attached image we can see the Taylor expansion for all the necessary therms and written in the form of Taylor Table . From the Taylor table we can say the sum of coefficient for the f'' (x)(dx)^2 is 1 and every other co-efficient with-in the number of points we took in the stencil is zero, So we can come up with 6 equations with 6 unknown. The equation is written below
equ_1 = a+b+c+d+e+g = 0
equ_2 = 0+b+2*c+3*d+4*e+5*g = 0
equ_3 = (0+b+4*c+9*d+16*e+25*g)/2 = 0
equ_4 = 0+b+8*c+27*d+64*e+125*h = 0
equ_5 = 0+b+16*c+81*d+256*e+625*g = 0
equ_6 = 0+b+32*c+243*d+1024*e+3125*g = 0
When we got this equation, basically, we can find the coefficient by Inverse matrix method or simply solve the equation in the Matlab. The are several ways to solve using different functions like
"X= linsolve (A,B);" or "X=SOLVE ([equ_1, equ_2,equ_3,equ_4,equ_5], [a,b,c,d,e]);" or "X=inv(A)*B;".
[A] = [ 1 1 1 1 1 1 ;
0 1 2 3 4 5 ;
0 1 4 9 16 25 ;
0 1 8 27 64 125 ;
0 1 16 81 256 625 ;
0 1 32 243 1024 3125 ]
[B] = [0;
0;
2;
0;
0;
0];
X= linsolve (A,B);
[a: [3.7500 :
b: -12.8333 :
X = c: = 17.8333 :
d: -13.000:
e: 5.0833:
g]: -0.8333 ]:
Therefore, the fourth order right skewed approximation for the second order derivative is
(d^2f/dx^2) = (3.750*f(x) -12.8333*f(x+dx) + 17.8333*f(x+2dx) - 13.00*f(x+3dx) + 5.0833*f(x+4dx) -0.8333*f(x+5dx))/(dx^2)
Function to Calculate the fourth order approximation skewed-right scheme:
%Function to calculate Skewed right difference fourth order approximation
function out = forder_eskewedright_diff (x,fm,f1r,f2r,f3r,f4r,f5r,dx)
%given analytical function exp(x)*cos(x)
%analytical_derivative_first_order f'(x)= exp(x)*cos(x)- exp(x)*sin(x)
%analytical derivative second order f''(x) = -2*exp(x)*sin(x)
analytical_solution = -2*exp(x)*sin(x);
%For the numerical derivative. We need to slove taylor table to find the co-efficient for each derivative.
%From taylor's table we can get the below equation
%equ_1 = a + b + c + d + e + g == 0;
%equ_2 = 0 + b + 2*c + 3*d + 4*e + 5*g == 0;
%equ_3 = 0 + b + 4*c + 9*d + 16*e + 25*g == 2;
%equ_4 = 0 + b + 8*c + 27*d + 64*e + 125*h == 0;
%equ_5 = 0 + b + 16*c + 81*d + 256*e + 625*g == 0;
%equ_6 = 0 + b + 32*c + 243*d + 1024*e + 3125*g == 0;
%We can represent the equation in the form of a matrix AX = B
%Matrix A
[A] = [ 1 1 1 1 1 1 ;
0 1 2 3 4 5 ;
0 1 4 9 16 25 ;
0 1 8 27 64 125 ;
0 1 16 81 256 625 ;
0 1 32 243 1024 3125 ] ;
% Matrix B
[B] = [ 0; 0; 2; 0; 0; 0 ];
% Solving the equation for the constant
X = linsolve(A,B);
% Extracting the values of co-efficient from the solution
a = X(1);
b = X(2);
c = X(3);
d = X(4);
e = X(5);
g = X(6);
% Skewed Right diffrence scheme fourth order approximation
fourthorder_sright_app = (a*fm + b*f1r + c*f2r + d*f3r + e*f4r + g*f5r)/(dx^2);
out = abs(analytical_solution - fourthorder_sright_app);
endHow do we know that the scheme is fourth order:
So, from the way we took the discriminate points in the stencil itself shows that the order of accuracy is fourth order (i.e. N = A + O -1 for SKEWED points).
However, we can prove this by EXPANDING the Taylor's table with considering the higher order terms.
d^2f/dx^2 = (b +4c +9d +16e +25g)f"(x)(dx^2/2!) + (4b+64c+729d+4096e+15625g) f"""(x)(dx^6/6!)+.....
sub d^2f/dx^2 = a*f(x) + b*f(x+dx) + c*f(x+2dx)+ d*f(x+3dx) + e*f(x+4dx) + g*f(x+5dx)
a*f(x) + b*f(x+dx) + c*f(x+2dx)+ d*f(x+3dx) + e*f(x+4dx) + g*f(x+5dx) = (b +4c +9d +16e +25g)f"(x)(dx^2/2!) + (4b+64c+729d+4096e+15625g) f"""(x)(dx^6/6!)+.....
divide both side by(dx^2)
(a*f(x) + b*f(x+dx) + c*f(x+2dx)+ d*f(x+3dx) + e*f(x+4dx) + g*f(x+5dx))/(dx^2) = (b +4c +9d +16e +25g)f"(x)(1/2!) + (4b+64c+729d+4096e+15625g) f"""(x)(dx^4/6!)+.....
In the above equation we can see that the first dx term has the power of 4 , So, if we truncate all those terms from dx^4, we get
(3.750*f(x) -12.8333*f(x+dx) + 17.8333*f(x+2dx) - 13.00*f(x+3dx) + 5.0833*f(x+4dx) -0.8333*f(x+5dx))/(dx^2) = (b +4c +9d +16e +25g)f"(x)(1/2!) + Error in order of (dx^4)
Hence, we proved that the scheme we come up with, have an error of order (dx^4)
Derivation of fourth order approximation for the Skewed left side difference scheme of second order derivative:
First we have to find the number of points we need to consider in the Numerical stencil, for that we can use the formula N = A + O for Skewed scheme,
where N = Number of points, A = Order of Approximation, O= Order of derivative
N = 4 + 2 = 6, SO we need to take Six points in the Numerical Stencil, towards left from f(x).
Which gives as f''(x) = a*f(x-5dx) + b*f(x-4dx) + c*f(x-3dx) + d*f(x-2dx) + e*f(x-dx) + g*f(x), all we have to find the co-efficient.
For that we have to Taylor expand the above equation and using Taylor table we can find the coefficient and desire formula for calculating f''(x).
| f''(x) | f(x) | f'(x)(dx) | f''(x)(dx^2/2!) | f'''(x)(dx^3/3!) | f''''(x)(dx^4/4!) | f""'(x)(dx^5/5!) | f''''''(x)(dx^6/6!) |
| a*f(x-5*dx) | a | -5a | 25a | -125a | 625a | -3125a | 15625a |
| b*f(x-4*dx) | b | -4b | 16b | -64b | 256b | -1024b | 4096b |
| c*f(x-3*dx) | c | -3c | 9c | -27c | 81c | -243c | 729c |
| d*f(x-2*dx) | d | -2d | 4d | -8d | 16d | -32d | 64d |
| e*f(x-dx) | e | -e | e | -e | e | -e | e |
| g*f(x) | g | 0 | 0 | 0 | 0 | 0 | 0 |
| ∑row | 0 | 0 | 1 | 0 | 0 | 0 | ? |
From the above attached image we can see the Taylor expansion for all the necessary therms and written in the form of Taylor Table . From the Taylor table we can say the sum of coefficient for the f'' (x)(dx)^2 is 1 and every other co-efficient with-in the number of points we took in the stencil is zero, So we can come up with 6 equations with 6 unknown. The equation is written below
equ_1 = a+b+c+d+e+g = 0
equ_2 = -5a-4b-3c-2d-e-0 = 0
equ_3 = (25a+16b+9c+4d+e+0)/2 = 1
equ_4 = -125a-64b-27c-8d-e-0= 0
equ_5 = 625a+256b+81c+16d+e+0 = 0
equ_6 = -3125a-1024b-243c-32d-e-0= 0
When we got this equation, basically, we can find the coefficient by Inverse matrix method or simply solve the equation in the Matlab. The are several ways to solve using different functions like
"X= linsolve (A,B);" or "X=SOLVE ([equ_1, equ_2,equ_3,equ_4,equ_5], [a,b,c,d,e]);" or "X=inv(A)*B;".
[A] = [ 1 1 1 1 1 1 ;
-5 -4 -3 -2 -1 0 ;
25 16 9 4 1 0 ;
-125 -64 -27 -8 -1 0 ;
625 256 81 16 1 0 ;
-3125 -1025 -243 -32 -1 0 ] ;
[B] = [0;
0;
2;
0;
0;
0];
X= linsolve (A,B)
[a; [-0.8775 ;
b; 5.3043;
X = c; = -13.4420;
d; 18.2754;
e; -13.0543;
g]; 3.7942 ];
Therefore, the fourth order right skewed approximation for the second order derivative is
(d^2f/dx^2) = (-0.8775*f(x-5*dx) +5.3043*f(x-4*dx) -13.4420*f(x-3*dx) +18.2754*f(x-2*dx) -13.0543*f(x-dx) +3.7942*f(x))/(dx^2)
Function to Calculate the fourth order approximation skewed-left scheme:
%Function to calculate Skewed left difference fourth order approximation
function out = forder_skewedleft_diff (x,f5l,f4l,f3l,f2l,f1l,fm,dx)
%given analytical function exp(x)*cos(x)
%analytical_derivative_first_order f'(x)= exp(x)*cos(x)- exp(x)*sin(x)
%analytical derivative second order f''(x) = -2*exp(x)*sin(x)
analytical_solution = -2*exp(x)*sin(x);
%For the numerical derivative. We need to slove taylor table to find the co-efficient for each derivative.
%From taylor's table we can get the below equation
%equ_1 = a + b + c + d + e + g == 0;
%equ_2 = -5*a - 4*b - 3*c - 2*d - e - 0 == 0;
%equ_3 = 25*a + 16*b + 9*c + 4*d + e + 0 == 2;
%equ_4 = -125*a - 64*b - 27*c - 8*d - e - 0 == 0;
%equ_5 = 625*a + 256*b + 81*c + 16*d + e + 0 == 0;
%equ_6 = -3125*g - 1024*e - 243*d - 32*c - b - 0 == 0;
%We can represent the equation in the form of a matrix AX = B
%Matrix A
[A] = [ 1 1 1 1 1 1 ;
-5 -4 -3 -2 -1 0 ;
25 16 9 4 1 0 ;
-125 -64 -27 -8 -1 0 ;
625 256 81 16 1 0 ;
-3125 -1025 -243 -32 -1 0 ] ;
% Matrix B
[B] = [0;
0;
2;
0;
0;
0];
% Solving the equation for the constant
X = linsolve(A,B);
% Extracting the values of co-efficient from the solution
a = X(1);
b = X(2);
c = X(3);
d = X(4);
e = X(5);
g = X(6);
%Skewed left difference scheme fourth Order approximation
fourthorder_sleft_app = (a*f5l+b*f4l+c*f3l+d*f2l+e*f1l+g*fm)/(dx^2);
out = abs(analytical_solution - fourthorder_sleft_app);
endHow do we know that the scheme is fourth order:
So, from the way we took the discriminate points in the stencil itself shows that the order of accuracy is fourth order (i.e. N = A + O -1 for SKEWED points).
However, we can prove this by EXPANDING the Taylor's table with considering the higher order terms.
d^2f/dx^2 = (25a+16b +9c +4d +e)f"(x)(dx^2/2!) + (15625a+ 4096b+ 729c+64d+e) f"""(x)(dx^6/6!) +.....
sub d^2f/dx^2 = a*f(x-5dx) + b*f(x-4dx) + c*f(x-3dx) + d*f(x-2dx) + e*f(x-dx) + g*f(x)
a*f(x-5dx) + b*f(x-4dx) + c*f(x-3dx) + d*f(x-2dx) + e*f(x-dx) + g*f(x) = (25a+16b +9c +4d +e)f"(x)(dx^2/2!) + (15625a+ 4096b+ 729c+64d+e) f"""(x)(dx^6/6!) +.....
divide both side by(dx^2)
(a*f(x-5dx) + b*f(x-4dx) + c*f(x-3dx) + d*f(x-2dx) + e*f(x-dx) + g*f(x))/(dx^2) = (25a+16b +9c +4d +e)f"(x)(1/2!) + (15625a+ 4096b+ 729c+64d+e) f"""(x)(dx^4/6!) +.....
In the above equation we can see that the first dx term has the power of 4 , So, if we truncate all those terms from dx^4, we get
(-0.8775*f(x-5*dx) +5.3043*f(x-4*dx) -13.4420*f(x-3*dx) +18.2754*f(x-2*dx) -13.0543*f(x-dx) +3.7942*f(x))/(dx^2)= (25a+16b +9c +4d +e)f"(x)(1/2!) + Error in Order of (dx^4)
Hence, we proved that the scheme we come up with have an error of order (dx^4)
Derivation of second order approximation forward difference scheme of second order derivative:
First we have to find the number of points we need to consider in the Numerical stencil, for that we can use the formula N = A + O -1 for Skewed scheme,
where N = Number of points, A = Order of Approximation, O= Order of derivative
N = 2 + 2 -1 = 3, SO we need to take 3 points in the Numerical Stencil, from f(x) forward.
Which gives as f''(x) = a*f(x) + b*f(x+dx) + c*f(x+2*dx), all we have to find the co-efficient.
For that we have to Taylor expand the above equation and using Taylor table we can find the coefficient and desire formula for calculating f''(x).
| f''(x) | f(x) | f'(x)(dx) | f''(x)(dx^2/2!) |
| a*f(x) | a | 0 | 0 |
| b*f(x+dx) | b | b | b |
| c*f(x+2*dx) | c | 2c | 4c |
| ∑row | 0 | 0 | 1 |
From the above attached image we can see the Taylor expansion for all the necessary therms and written in the form of Taylor Table . From the Taylor table we can say the sum of coefficient for the f'' (x)(dx)^2 is 1 and every other co-efficient with-in the number of points we took in the stencil is zero, So we can come up with 6 equations with 6 unknown. The equation is written below
equ 1 :- a + b + c = 0
equ 2 :- b + 2*c = 0
equ 3 :- b + 4*c = 2
From equ (2), we can say b = -2c sub in equ (3), we get c = 1 sub b and c in equ (1) we get a = 1
So, a = 1, b = -2, c = 1
Therefore, the second order forward approximation for the second order derivative is
(d^2f/dx^2) = 1*f(x) - 2*f(x+dx) + 1*f(x+2*dx) /(dx^2)
Function to Calculate the second order approximation forward difference:
%Function to calculate Second order forward difference
function out = secondorder_forward_diff (x,fm,f1r,f2r,dx)
%given analytical function exp(x)*cos(x)
%analytical_derivative_first_order f'(x)= exp(x)*cos(x)- exp(x)*sin(x)
%analytical derivative second order f''(x) = -2*exp(x)*sin(x)
analytical_solution = -2*exp(x)*sin(x);
%Forward difference approximation
secondorder_forward_app = (fm - 2*f1r + f2r)/(dx^2);
out = abs(analytical_solution - secondorder_forward_app);
endDerivation of second order approximation backward difference scheme of second order derivative:
First we have to find the number of points we need to consider in the Numerical stencil, for that we can use the formula N = A + O -1 for Skewed scheme,
where N = Number of points, A = Order of Approximation, O= Order of derivative
N = 2 + 2 -1 = 3, SO we need to take 3 points in the Numerical Stencil, from f(x) backward.
Which gives as f''(x) = a*f(x-2*dx) + b*f(x-dx) + c*f(x), all we have to find the co-efficient.
For that we have to Taylor expand the above equation and using Taylor table we can find the coefficient and desire formula for calculating f''(x).
From the above attached image we can see the Taylor expansion for all the necessary therms and written in the form of Taylor Table . From the Taylor table we can say the sum of coefficient for the f'' (x)(dx)^2 is 1 and every other co-efficient with-in the number of points we took in the stencil is zero, So we can come up with 6 equations with 6 unknown. The equation is written below
| f''(x) | f(x) | f'(x)(dx) | f''(x)(dx^2/2!) |
| a*f(x-2*dx) | a | -2a | 4a |
| b*f(x-dx) | b | -b | b |
| c*f(x) | c | 0 | 0 |
| ∑row | 0 | 0 | 1 |
equ 1 :- a + b + c = 0
equ 2 :- -2a-b = 0
equ 3 :- (dx^2) ( 4a+ b) = 2
From equ (2), we can say b = -2a sub in equ (3), we get a= 1 sub b and c in equ (1) we get c = 1
So, a = 1, b = -2, c = 1
Therefore, the second order fbackward approximation for the second order derivative is
(d^2f/dx^2) = a*f(x-2*dx) + b*f(x-dx) + c*f(x)
Function to Calculate the second order approximation forward difference:
%Function to calculate Second order backward difference
function out = secondorder_backward_diff (x,f2l,f1l,fm,dx)
%given analytical function exp(x)*cos(x)
%analytical_derivative_first_order f'(x)= exp(x)*cos(x)- exp(x)*sin(x)
%analytical derivative second order f''(x) = -2*exp(x)*sin(x)
analytical_solution = -2*exp(x)*sin(x);
%Backward difference approximation
secondorder_backward_app = (f2l - 2*f1l + fm)/(dx^2);
out = abs (analytical_solution - secondorder_backward_app);
endGraph: Output
From the above comparison of different scheme for fourth order approximation second order derivative, we can come to a conclusion that the central difference approximation yields better accuracy than the Skewed right or Skewed left difference scheme . Since, forward and backward approximation have second order accuracy, both are less accurate and have higher error that the fourth order accuracy schemes.
Note:
So, the question should be like, does the error become arbitrarily small if we reduce the value of dx? The answer is NO.
Because of the Numerical truncation error and important thing is that every machine (computer) have a Round off error around (er) ≈10^-16 for double precision. So, if we wanna store any arbitrarily small number in the memory it might be only represented in 10^-16 of the number, what we actually want to represent. Which means, If we choose to store A + (er/3), our computer cant tell the difference between A or (A+er/2) or (A=er/3). Simply, all the three numbers A or [(A+er/2) or (A=er/3)] are the same for our computer cause it doesn't have that much precision to differentiate those numbers. Thus, there is limit that how small we can take dx, where after this point the error starts to grow. So, the key take away is that value of error at any point is the net sum of the numerical error and the round off error. The round off error increases with the order of approximation. Therefore, the higher order scheme is not always appreciable. So, there is a trade off between this Taylor's series truncation error and machine round-off error while choosing the value of dx.
Below you can find the screenshot of the output that has some round-off error: Here, I took x= (2*pi) and dx = (pi/4,pi/4000,60), So the error starts to grow after certain point.
Importance of skewed scheme. What a skewed scheme do that a CD scheme cannot do:
Generally, Central difference approximation is better compare to the forward difference approximation and backward difference approximation. Because we have seen previously that the Forward and Backward difference approximation systematically over predicted and under predicted respectively. And, it's always a good to go when it come to Central difference scheme.
But, there are limitations, that's where Skewed difference approximation come handy.
So, Basically its not possible to compute the derivative f'(x) at the boundaries of the given x data. Because, every function is built in a boundary characterized by domain and anywhere beyond the domain the function can not exist. Since, Central difference demand data both side from the boundary (if z is the boundary of the function CD need f(z+dx) thats why it is impossible to calculate the central difference scheme. Whereas, we can use Forward or Skewed right scheme to get higher order accuracy and at initial point and Backward difference or Skewed Left scheme at the endpoint of the function.
If we think a bit technical, Practically, the Central Difference scheme can be powerless, when we want to calculate the function f'(t) with respect to t in Real Time. We cant use Central difference, because only at every t sec the new data will be updated and we can't predict f(t+dt) the future in reality to calculate the central difference approximation. An illustration, I would like to think about is an asteroid falling, If we have the position or velocity f(t) of the asteroid we can't take the derivative f'(t) of the function using central difference scheme to predict the position or the velocity of the asteroid.