Week 3.5 - Deriving 4th order approximation of a 2nd order derivative using Taylor Table method
Tayler's Series Expansion
The Taylor series of a function is an infinite sum of terms that are expressed in terms of the function's derivatives at a single point ,
f(x+Δx) = f(x) + f'(x).Δx + f''(x).(Δx)^2/2! + f'''(x).(Δx)^3/3! + ..........
Now finding the value of f'(x) , f''(x) in 4rh order approximation
Deriving the 4th order approximations of the second-order derivative
1. Central Difference
The number of Nodes required for 4th order approximation of a second-order derivative can be given by the equation
N=P+Q">N=P+Q
N=P+Q">N= 2+4-1 = 5
where
N">N is the number of nodes
Q">Q is the order of derivative
P">P is the order of approximation
So, here we need 5 nodes in order to find the 2nd order derivative with the 4th order of approximation.
Consider the value of ∂^2f/(∂x)^2 with the 4th order of approximation is given by the equation
∂^2f/(∂x)^2 = af(x-2Δx) + bf(x-Δx) + cf(x) + df(x+Δx) + ef(x+2Δx)
Now, using taylor series to find the values of discrete poits with repect to point x
af(x-2Δx) = af(x) - af'(x).Δx + af''(x).(Δx)^2/2! - af'''(x).(Δx)^3/3! + af''''(x).(Δx)^4/4! - af'''''(x).(Δx)^5/5! +...........
bf(x-Δx) = bf(x) - bf'(x).Δx + bf''(x).(Δx)^2/2! - bf'''(x).(Δx)^3/3! + bf''''(x).(Δx)^4/4! - bf'''''(x).(Δx)^5/5! + .........
cf(x) = cf(x)
df(x+Δx) = df(x) + df'(x).Δx + df''(x).(Δx)^2/2! + df'''(x).(Δx)^3/3! + df''''(x).(Δx)^4/4! + df'''''(x).(Δx)^5/5! + .........
ef(x+2Δx) = ef(x) + ef'(x).Δx + ef''(x).(Δx)^2/2! + ef'''(x).(Δx)^3/3! + ef''''(x).(Δx)^4/4! + ef'''''(x).(Δx)^5/5! +...........
adding the LHS and RHS of the above equations , we get
af(x-2Δx) + bf(x-Δx) + cf(x) + df(x+Δx) + ef(x+2Δx) = f(x)[a+b+c+d+e] - f'(x).Δx[2a+b+d+2e] + f''(x).(Δx)^2[4a/2+b+d+4e/2] - f'''(x).(Δx)^3[8a/6+b+d+8e/6] + f''''(x).(Δx)^4[16a/24+b+d+16e/24] - f'''''(x).(Δx)^5[32a/120+b+d+32e/120] + ..........
In the above equations, we want only the coefficients of f''(x).Δx^2 and we assume the coefficient to be one and for the remaining terms we assume the coefficients to be zero.
1.1 TAYLER'S TABLE
|
Assumed function ∂^2f/(∂x)^2/(Δx)^2 |
f(x) | f'(x).Δx | f''(x).(Δx)^2/2! |
f'''(x). (Δx)^3/3! |
f''''(x).(Δx)^4/4! | f'''''(x).(Δx)^5/5! | |
| af(x-2Δx) | a | -2a | 4a | -8a | 16a | -32a | |
| bf(x-Δx) | b | -b | b | -b | b | -b | |
| cf(x) | c | 0 | 0 | 0 | 0 | 0 | |
| df(x+Δx) | d | d | d | d | d | d | |
| ef(x+Δx) | e | 2e | 4e | 8e | 16e | 32e | |
| sum of co-efficients | 0 | 0 | 1 | 0 | 0 | ? |
By coupling above co-efficients
a+b+c+d+e = 0
-2a-b+d+2e = 0
(4a/2)+(b/2)+(d/2)+(4e/2) = 1
-(8a/6)-(b/6)+(d/6)+(8e/6) = 0
(16a/24)+(b/24)+(d/24)(16e/24) = 0
-(32a/120)-(b/120)+(d/120)+(32e/120) = ?
Above algebric equation written into matrix i.e
A = [1 1 1 1 1 ;-2 -1 0 1 2 ;4/2 1/2 0 1/2 4/2 ;-8/6 -1/6 0 1/6 8/6 ;16/24 1/24 0 1/24 16/24]
B = [0;0;1;0;0]
X = AB
A = 1.0000 1.0000 1.0000 1.0000 1.0000
-2.0000 -1.0000 0 1.0000 2.0000
2.0000 0.5000 0 0.5000 2.0000
-1.3333 -0.1667 0 0.1667 1.3333
0.6667 0.0417 0 0.0417 0.6667
B = 0
0
1
0
0
X= -0.0833
1.3333
-2.5000
1.3333
-0.0833
Now we can write the equation using the Taylors Table as
∂^2f/(∂x)^2 = ( af(x-2Δx) + bf(x-Δx) + cf(x) + df(x+Δx) + ef(x+2Δx) )/(Δx)^2
2. Skewed right-sided difference
The number of Nodes required for 4th order approximation of a second-order derivative can be given by the equation
N=P+Q">N=P+Q
N=P+Q">N=2+4 = 6
where
N">N is the number of nodes
Q">Q is the order of derivative
P">P is the order of approximation
So, here we need 6 nodes in order to find the 2nd order derivative with the 4th order of approximation.
Consider the value of ∂^2f/(∂x)^2 with the 4th order of approximation is given by the equation
1.2 TAYLOR'S TABLE
Tayer's Table Expansion
|
Assumed function ∂^2f/(∂x)^2/(Δx)^2 |
f(x) | f'(x).Δx | f''(x).(Δx)^2/2! |
f'''(x). (Δx)^3/3! |
f''''(x).(Δx)^4/4! | f'''''(x).(Δx)^5/5! | f'''''(x).(Δx)^6/6! | |
| af(x) | a | 0 | 0 | 0 | 0 |
0 |
0 | |
| bf(x+Δx) | b | b | b | b | b | b | b | |
| cf(x+2Δx) | c | 2c | 4c | 8c | 16c | 32c | 64c | |
| df(x+3Δx) | d | 3d | 9d | 27d | 81d | 243d | 729d | |
| ef(x+4Δx) | e | 4e | 16e | 64e | 256e | 1024e | 4096e | |
| gf(x+5Δx) | g | 5g | 25g | 125g | 625g | 3125g | 15625g | |
| sum of co-efficients | 0 | 0 | 1 | 0 | 0 | 0 | ? |
By coupling the sum of co-efficients
a+b+c+d+e+g=0
b+2c+3d+4e+5g=0
b+4c/2+9d/2+16e/2+25g/2=1
b+8c/6+27d/6+64e/6+125g/6=0
b+16c/24+81d/24+256e/24+625g/24=0
b+32c/120+243d/120+1024e/120+3125g/120=0
b+64c/720+729d/720+4096e/720+15625g/720=?
Above algebric equation written into matrix i.e
A = [1 1 1 1 1 1 ;0 1 2 3 4 5 ;0 1/2 4/2 9/2 16/2 25/2 ;0 1/6 8/6 27/6 64/6 125/6 ;0 1/24 16/24 81/24 256/24 625/24 ;0 1/120 32/120 243/120 1024/120 3125/120]
B = [0;0;1;0;0;0]
X=AB
A = 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000
0 1.0000 2.0000 3.0000 4.0000 5.0000
0 0.5000 2.0000 4.5000 8.0000 12.5000
0 0.1667 1.3333 4.5000 10.6667 20.8333
0 0.0417 0.6667 3.3750 10.6667 26.0417
0 0.0083 0.2667 2.0250 8.5333 26.0417
B = 0
0
1
0
0
0
X = 3.7500
-12.8333
17.8333
-13.0000
5.0833
-0.8333
Now we can write the equation using the Taylors Table as
∂^2f/(∂x)^2 = ( af(x)+bf(x+Δx)+cf(x+2Δx)+df(x+3Δx)+ef(x+4Δx)+gf(x+5Δx) )/(Δx)^2
3. Skewed left-sided difference
The number of Nodes required for 4th order approximation of a second-order derivative can be given by the equation
N=P+Q">N=P+Q
N=P+Q">N=2+4 = 6
where
N">N is the number of nodes
Q">Q is the order of derivative
P">P is the order of approximation
So, here we need 6 nodes in order to find the 2nd order derivative with the 4th order of approximation,
Consider the value of ∂^2f/(∂x)^2 with the 4th order of approximation is given by the equation
1.3 TAYLOR'S TABLE
Tayer's Table Expansion
|
Assumed function ∂^2f/(∂x)^2/(Δx)^2 |
f(x) | f'(x).Δx | f''(x).(Δx)^2/2! |
f'''(x). (Δx)^3/3! |
f''''(x).(Δx)^4/4! | f'''''(x).(Δx)^5/5! | f'''''(x).(Δx)^6/6! | |
| af(x) | a | 0 | 0 | 0 | 0 |
0 |
0 | |
| bf(x-Δx) | b | -b | b | -b | b | -b | b | |
| cf(x-2Δx) | c | -2c | 4c | -8c | 16c | -32c | 64c | |
| df(x-3Δx) | d | -3d | 9d | -27d | 81d | -243d | 729d | |
| ef(x-4Δx) | e | -4e | 16e | -64e | 256e | -1024e | 4096e | |
| gf(x-5Δx) | g | -5g | 25g | -125g | 625g | -3125g | 15625g | |
| sum of co-efficients | 0 | 0 | 1 | 0 | 0 | 0 | ? |
By coupling the sum of co-efficients
a+b+c+d+e+g=0
-b-2c-3d-4e-5g=0
b+4c/2+9d/2+16e/2+25g/2=1
-b-8c/6-27d/6-64e/6-125g/6=0
b+16c/24+81d/24+256e/24+625g/24=0
-b-32c/120-243d/120-1024e/120-3125g/120=0
b+64c/720+729d/720+4096e/720+15625g/720=?
Above algebric equation written into matrix i.e
A = [1 1 1 1 1 1 ;0 -1 -2 -3 -4 -5 ;0 1/2 4/2 9/2 16/2 25/2 ;0 -1/6 -8/6 -27/6 -64/6 -125/6 ;0 1/24 16/24 81/24 256/24 625/24 ;0 -1/120 -32/120 -243/120 -1024/120 -3125/120]
B = [0;0;1;0;0;0]
X=AB
A = 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000
0 -1.0000 -2.0000 -3.0000 -4.0000 -5.0000
0 0.5000 2.0000 4.5000 8.0000 12.5000
0 -0.1667 -1.3333 -4.5000 -10.6667 -20.8333
0 0.0417 0.6667 3.3750 10.6667 26.0417
0 -0.0083 -0.2667 -2.0250 -8.5333 -26.0417
B = 0
0
1
0
0
0
X = 3.7500
-12.8333
17.8333
-13.0000
5.0833
-0.8333
Now we can write the equation using the Taylors Table as
∂^2f/(∂x)^2 = ( af(x)+bf(x-Δx)+cf(x-2Δx)+df(x-3Δx)+ef(x-4Δx)+gf(x-5Δx) )/(Δx)^2
clear all
close all
clc
%analytical fucntion=(e^x)cos(x);
%analytical derivative
%f(x)=(e^x)cos(x)
%f'(x)=(e^x)cos(x)-(e^x)sin(x)=e^x[cosx-sinx]
%f''(x)=(e^x)cos(x)-(e^x)sin(x)-[(e^x)sin(x)+(e^x)cos(x)]
%f''(x)=-2*e^x*sin(x)
%Given inputs
x=pi/3; %derivative to be computed at x=pi/3
dx=linspace(pi/4000,pi/40); %grid spacing of dx(Assumed)
%analytical derivative
analytical_derivative=-2*exp(x)*sin(x);
%Skewed right-sided(4th order approimation)
%Solving for the coefficients
A=[1 1 1 1 1 1;0 1 2 3 4 5;0 1/2 4/2 9/2 16/2 25/2;0 1/6 8/6 27/6 64/6 125/6;0 1/24 16/24 81/24 256/24 625/24;0 1/120 32/120 243/120 1024/120 3125/120];
B=[0;0;1;0;0;0];
C_r=inv(A)*B;
%Skewed left-sided(4th order approimation)
%Solving for the coefficients
A=[1 1 1 1 1 1;0 -1 -2 -3 -4 -5;0 1/2 4/2 9/2 16/2 25/2;0 -1/6 -8/6 -27/6 -64/6 -125/6;0 1/24 16/24 81/24 256/24 625/24;0 -1/120 -32/120 -243/120 -1024/120 -3125/120];
B=[0;0;1;0;0;0];
C_l=inv(A)*B;
%central differencing(4th order approimation)
%Solving for the coefficients
A=[1, 1, 1, 1, 1;-2, -1, 0, 1, 2;4/2, 1/2, 0, 1/2, 4/2;-8/6, -1/6, 0, 1/6, 8/6;16/24, 1/24, 0, 1/24, 16/24];
B=[0;0;1;0;0];
C_c=inv(A)*B;
for i=1:length(dx)
%numerical derivative
%Skewed right-sided(4th order approimation)
%f''(x)=af(x)+bf(x+Δx)+cf(x+2Δx)+df(x+3Δx)+ef(x+4Δx)+gf(x+5Δx)
right_difference=(C_r(1)*(x)+C_r(2)*(x+dx(i))+C_r(3)*(x+2*dx(i))+C_r(4)*(x+3*dx(i))+C_r(5)*(x+4*dx(i))+C_r(6)*(x+5*dx(i)))/(dx(i)^2);
%Skewed left-sided(4th order approimation)
%f''(x)=af(x)+bf(x-Δx)+cf(x-2Δx)+df(x-3Δx)+ef(x-4Δx)+gf(x-5Δx)
left_difference=(C_l(1)*(x)+C_l(2)*(x-dx(i))+C_l(3)*(x-2*dx(i))+C_l(4)*(x-3*dx(i))+C_l(5)*(x-4*dx(i))+C_l(6)*(x-5*dx(i)))/(dx(i)^2);
%central differencing(4th order approimation)
%f''(x)=af(x-2Δx)+bf(x-Δx)+cf(x)+df(x+Δx)+ef(x+2Δx)
central_difference=(C_c(1)*(x-2*dx(i))+C_c(2)*(x-1*dx(i))+C_c(3)*(x)+C_c(4)*(x+dx(i))+C_c(5)*(x+2*dx(i)))/(dx(i)^2);
%absolute error at x = pi/3
error_r(i)=abs(analytical_derivative-right_difference);
error_l(i)=abs(analytical_derivative-left_difference);
error_c(i)=abs(analytical_derivative-central_difference);
error=[error_r;error_l;error_c];
end
%Plotting the Absolute Error Curves
figure(1)
loglog(dx,error_r,'color','r')
grid on
hold on
loglog(dx,error_l,'color','b')
loglog(dx,error_c,'color','g')
ylabel('Absolute Error (Log Scale)')
xlabel('Grid Spacing (dx in Log Scale)')
legend('Right Differencing','Left Differencing','Central Differencing')
title('Absolute Error Comparison')
The above graph shows that the central scheme is most accurate compared to the skewed differencing. The curves also show that when the value of dx is very small or large the approximating function becomes unstable. This instability is caused due to the combination of various errors such as roundoff error, residual error and truncation error.
coclusion: A skewed scheme is useful when we have to predict using the boundary condition as we have only one direction to go, then the skewed scheme will be useful and we cannot use the central differencing scheme for such an application because it required information from both directions of the point. (i.e in fig. shows that in condition of skewes direction points are unstable with bounary condition and for the central differencing it has stable condition from both the direction.)