Taylor table method and Matlab code
- Central difference
Numerical stencil 2nd order derivative for 4th order accuracy by central difference:-
(x-2Δx),(x-Δx),x,(x+Δx),(x+2Δx)
Now computing(2x)/(x2):-
(2x)/(x2)=a.f(x-2Δx)+b.f(x-Δx)+c.f(x)+d.f(x+Δx)+e.f(x+2Δx)
Now making the taylor table:-
|
(2x)/(x2) |
f(x) |
f’(x).Δx |
f’’(x).Δx2 |
f’’’(x).Δx3 |
fiv(x).Δx4 |
|
a.f(x-2Δx) |
a |
-2a |
2a |
-8/6.a |
16/24.a |
|
b.f(x-Δx) |
b |
-b |
b/2 |
-⅙.8 |
1/24.b |
|
c.f(x) |
c |
0 |
0 |
0 |
0 |
|
d.f(x+Δx) |
d |
d |
d/2 |
1/6.d |
1/24.d |
|
e.f(x+2Δx) |
e |
2e |
2e |
8/6.e |
16/24.e |
|
Σ |
0 |
0 |
1 |
0 |
0 |
Now using matlab to find out the values of each coefficient:-
a= -0.083333, b=1.333333, c=-2.500000, d=1.333333, e=-0.083333
- Skewed right-sided difference
Numerical stencil 2nd order derivative for 4th order accuracy by skewed right sided difference:-
x,(x+Δx),(x+2Δx),(x+3Δx),(x+4Δx)
Now computing(2x)/(x2):-
(2x)/(x2)=a.f(x)+b.f(x+Δx)+c.f(x+2Δx)+d.f(x+3Δx)+e.f(x+4Δx)
Now making the taylor table:-
|
(2x)/(x2) |
f(x) |
f’(x).Δx |
f’’(x).Δx2 |
f’’’(x).Δx3 |
fiv(x).Δx4 |
fv(x).Δx5 |
|
a.f(x) |
a |
0 |
0 |
0 |
0 |
0 |
|
b.f(x+Δx) |
b |
b |
b/2 |
1/6.b |
1/24.b |
1/120.b |
|
c.f(x+2Δx) |
c |
2.c |
4/2.c |
8/6.c |
16/24.c |
32/120.c |
|
d.f(x+3Δx) |
d |
3.d |
9/2.d |
27/6.d |
81/24.d |
243/120.d |
|
e.f(x+4Δx) |
e |
4.e |
16/2.e |
64/6.e |
256/24.e |
1024/120.e |
|
f.f(x+5Δx) |
f |
5.f |
25/2.f |
125/6.f |
625/24.f |
3125/120.f |
|
Σ |
0 |
0 |
1 |
0 |
0 |
0 |
Now using matlab to find out the values of each coefficient:-
a=3.75000,b=-12.83333, c=17.83333, d=-13.00000, e=5.08333, f=-0.83333
- Skewed left-sided difference:-
Numerical stencil 2nd order derivative for 4th order accuracy by skewed left sided difference:-
x,(x-Δx),(x-2Δx),(x-3Δx),(x-4Δx),
Now computing(2x)/(x2):-
(2x)/(x2)=a.f(x)+b.f(x-Δx)+c.f(x-2Δx)2+d.f(x-3Δx)3+e.f(x-4Δx)4+f.(x-5Δx)5
Now making the taylor table:-
|
(2x)/(x2) |
f(x) |
f’(x).Δx |
f’’(x).Δx2 |
f’’’(x).Δx3 |
fiv(x).Δx4 |
fv(x).Δx5 |
|
a.f(x) |
a |
0 |
0 |
0 |
0 |
0 |
|
b.f(x+Δx) |
b |
-b |
b |
-1/6.b |
1/24.b |
-1/120.b |
|
c.f(x+2Δx) |
c |
-2.c |
4/2.c |
-8/6.c |
16/24.c |
-32/120.c |
|
d.f(x+3Δx) |
d |
-3.d |
9/2.d |
-27/6.d |
81/24.d |
-243/120.d |
|
e.f(x+4Δx) |
e |
-4.e |
16/2.e |
-64/6.e |
256/24.e |
-1024/120.e |
|
f.f(x+5Δx) |
f |
-5.f |
25/2.f |
-125/6.f |
625/24.f |
-3125/120.f |
|
Σ |
0 |
0 |
1 |
0 |
0 |
0 |
Now using matlab to find out the values of each coefficient:-
a=3.75000,b=-12.83333, c=17.83333, d=-13.00000, e=5.08333, f=-0.83333
Main Code:-
clear all
close all
clc
%4th order approximation of 2nd order derivative
%value of x
x=pi;
%analytical_derivative
% analytical finction = exp(x)*cos(x)
% 2nd order derivative= -2*sin(x)*exp(x)
analytical_derivative=(-2*(sin(x))*(exp(x)));
%numerical_derivative
%value of dx
dx=linspace(pi/2,pi/4000000000,2000);
for i=1:length(dx)
%2nd order central differencing
central_differencing(i)=((-0.083333)*exp(x-2*(dx(i)))*cos(x-2*(dx(i)))+(1.333333)*exp(x-(dx(i)))*cos(x-(dx(i)))+(-2.500000)*exp(x)*cos(x)+(1.333333)*exp(x+(dx(i)))*cos(x+(dx(i)))+(-0.083333)*exp(x+2*(dx(i)))*cos(x+2*(dx(i))))/((dx(i))^2);
%error in 2nd order central differencing
error_central_diff(i)=abs(central_differencing(i)-analytical_derivative);
%2nd order skewed right sided differencing
skewed_right_sided(i)=(((3.75)*exp(x)*cos(x))+((-12.833)*exp(x+dx(i))*cos(x+dx(i)))+((17.833)*exp(x+(2*dx(i)))*cos(x+(2*dx(i))))+((-13)*exp(x+(3*dx(i)))*cos(x+(3*dx(i))))+((5.0833)*exp(x+(4*dx(i)))*cos(x+(4*dx(i))))+((-0.833)*exp(x+(5*dx(i)))*cos(x+(5*dx(i)))))/((dx(i))^2);
%error in 2nd order skewed right sided differencing
error_skewed_right(i)=abs(skewed_right_sided(i)-analytical_derivative);
%2nd order skewed left sided differencing
skewed_left_sided(i)=(((-0.833)*exp(x-(5*dx(i)))*cos(x-(5*dx(i))))+((5.0833)*exp(x-(4*dx(i)))*cos(x-(4*dx(i))))+((-13)*exp(x-(3*dx(i)))*cos(x-(3*dx(i))))+((17.833)*exp(x-(2*dx(i)))*cos(x-(2*dx(i))))+((-12.833)*exp(x-dx(i))*cos(x-dx(i)))+((3.75)*exp(x)*cos(x)))/((dx(i))^2);
%error in 2nd order skewed left sided differencing
error_skewed_left(i)=abs(skewed_left_sided(i)-analytical_derivative);
end
loglog(dx,error_central_diff,'r')
hold on
loglog(dx,error_skewed_right,'g')
hold on
loglog(dx,error_skewed_left,'b')
Result:-
red line= central scheme
green line=Skewed right-sided difference
blue line=Skewed left-sided difference
Observation:-
1)Central differencing scheme seems more stable than the other two methods and error is in fixed proportions with dx.
2)for certain range of values, the error given by central scheme is smaller then the other two schemes.
3)the sharp dips in the curve is due to unstablility introduced in the system.