Week 3.5 - Deriving 4th order approximation of a 2nd order derivative using Taylor Table method

Theory :

During solving a CFD problem we may come across partial derivative equation which consists 2^nd order derivative. Since due to complexity of the equation we may have to chose numerical approach rather than analytical. But when dealing with numerical approach we try to make the error margin as low as possible.

As the order of approximation get higher, we get higher accuracy. Hence, we are going to derive the three methods of approximation for an analytical function. The methods are following.

1. CENTRAL DIFFERENCING SCHEME:

NO OF NODES:

The number of nodes `= p+q-1`

Where, `p` is order of derivative.

`q` is order of approximations.

The no of nodes `= 2+4-1 = >5.`

DERIVATION:

To derive a fourth order approximation for the second order derivative,we start with these approximation,

`(∂^2u)/(∂x^2)=(af(i−2)+bf(i−1)+cf(i)+df(i+1)+ef(i+2))/(dx^2)`..........(1)

taylor series derivation for the above terms,

`af(i−2)=af(i)−(2af'(i)dx)/(1!)+(4af''(i)dx^2)/(2!)+(8af'''(i)dx^3)/(3!)+(16af''''(i)dx^4)/(4!)`

Simillarly expand all the terms in the equation (1) and are represented in the form of table using taylor table method.

`dx2f''(i)`	`f(i)`	`f'(i)dx`	`f''(i)dx^2`	`f'''(i)dx^3`	`f''''(i)dx^4`
`af(i-2)`	`a`	`-2a`	`2a`	`−(4a)/3`	`(2a)/3`
`bf(i-1)`	`b`	`-b`	`b/2`	`−b/6`	`b/24`
`cf(i)`	`c`	`0`	`0`	`0`	`0`
`df(i+1)`	`d`	`d`	`d/2`	`d/6`	`d/24`
`ef(i+2)`	`e`	`2e`	`2e`	`(4e)/3`	`(2e)/3`
sum of co_efficients	`0`	`0`	`1`	`0`	`0`

the linear equation obtained from the above table are solved using matrix inversion method,

`A.X = B`

`[[1,1,1,1,1],[-2,-1,0,1,2],[2,1/2,0,1/2,2],[-4/3,-1/6,0,1/6,1/4],[2/3,1/24,0,1/24,2/3]]*[[a],[b],[c],[d],[e]]=[[0],[0],[1],[0],[0]]`

`X = text(inv)(A)*B`

values of the unknown are,

substituting these coefficient in equation(1) we get,

`(∂^2u)/(∂x^2)=((−0.0833)f(i−2)+(1.3333)f(i−1)+(−2.5000)f(i)+(1.3333)f(i+1)+(−0.0833)f(i+2))/(dx^2)`

2. SKEWED RIGHT HAND DIFFERENCE SCHEME :

NO OF NODES :

The number of nodes `= p+q`

Where, `p` is order of derivative.

`q` is order of approximations.

The no of nodes `= 2+4 = >6.`

DERIVATION :

To derive a fourth order approximation for the second order derivative,we start with these approximation,

`(∂2u)/(∂x^2)=(af(i)+bf(i+1)+cf(i+2)+df(i+3)+ef(i+4)+gf(i+5))/(dx^2)`..........(2)

expand all the terms in the equation (2) using taylor series expansion and are represented in the form of table using taylor table method.

`dx2f''(i)`	`f(i)`	`f'(i)dx`	`f''(i)dx^2`	`f'''(i)dx^3`	`f''''(i)dx^4`	`f'''''(i)dx^5`
`af(i)`	`a`	`0`	`0`	`0`	`0`	`0`
`bf(i+1)`	`b`	`b`	`b/2`	`b/6`	`b/24`	`b/120`
`cf(i+2)`	`c`	`2c`	`(4c)/2`	`(8c)/6`	`(16c)/24`	`(32c)/120`
`df(i+3)`	`d`	`3d`	`(9d)/2`	`(27d)/6`	`(81d)/24`	`(243d)/120`
`ef(i+4)`	`e`	`4e`	`(16e)/2`	`(64e)/6`	`(256e)/24`	`(1024e)/120`
`gf(i+5)`	`g`	`5g`	`(25g)/2`	`(125g)/6`	`(625g)/24`	`(3125g)/120`
sum of co_efficients	`0`	`0`	`1`	`0`	`0`	`0`

the linear equation obtained from the above table are solved using matrix inversion method,

`A.X = B`

`[[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]]*[[a],[b],[c],[d],[e],[g]]=[[0],[0],[1],[0],[0],[0]]`

`X = text(inv)(A)*B`

we get the coefficients as,

substituting these coefficient in equation(2) we get,

`(∂^2u)/(∂x^2)=((3.7500)f(i)+(−12.8333)f(i+1)+(17.8333)f(i+2)+(−13.0000)f(i+3)+(5.0833)f(i+4)+(−0.8333)f(i+5))/(dx^2)`

3. SKEWED LEFT HAND DIFFERENCE SCHEME :

NO OF NODES :

The number of nodes `= p+q`

Where, `p` is order of derivative.

`q` is order of approximations.

The no of nodes `= 2+4 = >6.`

DERIVATION :

To derive a fourth order approximation for the second order derivative,we start with these approximation,

`(∂^2u)/(∂x^2)=(af(i−5)+bf(i−4)+cf(i−3)+df(i−2)+ef(i−1)+gf(i))/(dx^2)`..........(3)

expand all the terms in the equation (3) using taylor series expansion and are represented in the form of table using taylor table method.

`dx2f''(i)`	`f(i)`	`f'(i)dx`	`f''(i)dx^2`	`f'''(i)dx^3`	`f''''(i)dx^4`	`f'''''(i)dx^5`
`af(i-5)`	`a`	`-5a`	`(25a)/2`	`−(125a)/6`	`(625a)/24`	`−(3125a)/120`
`bf(i-4)`	`b`	`-4b`	`(16b)/2`	`−(64b)/6`	`(256b)/24`	`−(1024b)/120`
`cf(i-3)`	`c`	`-3c`	`(9c)/2`	`−(27c)/6`	`(81c)/24`	`−(243c)/120`
`df(i-2)`	`d`	`-2d`	`(4d)/2`	`−(8d)/6`	`(16d)/24`	`−(32d)/120`
`ef(i-1)`	`e`	`-e`	`e/2`	`−e/6`	`e/24`	`−e/120`
`gf(i)`	`g`	`0`	`0`	`0`	`0`	`0`
sum of co_efficients	`0`	`0`	`1`	`0`	`0`	`0`

the linear equation obtained from the above table are solved using matrix inversion method,

`A.X = B`

`[[1,1,1,1,1,1],[-5,-4,-3,-2,-1,0],[25/2,16/2,9/2,4/2,1/2,0],[-125/6,-64/6,-27/6,-8/6,-1/6,0],[625/24,256/24,81/24,16/24,1/24,0],[-3125/120,-1024/120,-243/120,-32/120,-1/120,0]].[[a],[b],[c],[d],[e],[g]]=[[0],[0],[1],[0],[0],[0]]`

`X = text(inv)(A)*B`

we get the coefficients as,

substituting these coefficient in equation(3) we get

`(∂^2u)/(∂x^2)=((−0.8333)f(i−5)+(5.0833)f(i−4)+(−13.0000)f(i−3)+(17.8333)f(i−2)+(−12.8333)f(i−1)+(3.7500)f(i))/(dx^2)`

Main Code :

clear all
close all
clc

%value of x
x= pi/3;

%range of dx
dx= linspace(pi/4,pi/4000,20);

%for loop
for i=1:length(dx)
    
    central_differencing(i) = fourth_order_central_differencing(x,dx(i));
    
    right_sided_skew(i) = fourth_order_right_sided_skew(x,dx(i));
    
    left_sided_skew(i) = fourth_order_left_sided_skew(x,dx(i));
    
end

%plotting
figure(1)
loglog(dx,central_differencing,'r','LineWidth',2)
hold on
loglog(dx,right_sided_skew,'g','LineWidth',2)
hold on
loglog(dx,left_sided_skew,'b','LineWidth',2)
xlabel('range of dx')
ylabel('truncation error')
title('range of dx vs error')
legend('central differencing','skewed right sided difference','skewed left sided difference','location','northwest')

%% Functions 
% Function code for Central Differencing
function Error_central_differencing = fourth_order_central_differencing(x,dx)

% analytic function = exp(x)*cos(x)

% analytic derivative
%f''(x) = -2*exp(x)*sin(x)

analytical_derivative = -2*exp(x)*sin(x);

%numerical derivative

%central differencing
central_differencing = (((-0.0833)*((exp(x-(2*dx)))*(cos(x-(2*dx)))))+((1.3333)*((exp(x-dx))*(cos(x-dx))))-((2.5000)*((exp(x))*(cos(x))))+((1.3333)*((exp(x+dx))*(cos(x+dx))))-((0.0833)*((exp(x+(2*dx)))*(cos(x+(2*dx))))))/(dx^2);

%error
Error_central_differencing = abs(central_differencing- analytical_derivative);
end

% Function code for Skewed Right Sided Differencing
function Error_right_sided_skew = fourth_order_right_sided_skew(x,dx)

% analytic function = exp(x)*cos(x)

% analytic derivative
%f''(x) = -2*exp(x)*sin(x)

analytical_derivative = -2*exp(x)*sin(x);

%numerical derivative

%central differencing
right_sided_skew = (((3.7500)*((exp(x))*(cos(x))))-((12.8333)*((exp(x+dx))*(cos(x+dx))))+((17.8333)*((exp(x+(2*dx)))*(cos(x+(2*dx)))))-((13.0000)*((exp(x+(3*dx)))*(cos(x+(3*dx)))))+((5.0833)*((exp(x+(4*dx)))*(cos(x+(4*dx)))))-((0.8333)*((exp(x+(5*dx)))*(cos(x+(5*dx))))))/(dx^2);

%error
Error_right_sided_skew = abs(right_sided_skew- analytical_derivative);
end

% Function code for Skewed Left Sided Differencing
function Error_left_sided_skew = fourth_order_left_sided_skew(x,dx)

% analytic function = exp(x)*cos(x)

% analytic derivative
%f''(x) = -2*exp(x)*sin(x)

analytical_derivative = -2*exp(x)*sin(x);

%numerical derivative

%central differencing
left_sided_skew = (((-0.8333)*((exp(x-(5*dx)))*(cos(x-(5*dx)))))+((5.0833)*((exp(x-(4*dx)))*(cos(x-(4*dx)))))-((13.0000)*((exp(x-(3*dx)))*(cos(x-(3*dx)))))+((17.8333)*((exp(x-(2*dx)))*(cos(x-(2*dx)))))-((12.8333)*((exp(x-dx))*(cos(x-dx))))+((3.7500)*((exp(x))*(cos(x)))))/(dx^2);

%error
Error_left_sided_skew = abs(left_sided_skew- analytical_derivative);
end

Result :

Conclusion :

As seen from the plot, the central differencing scheme has very less error compared to other schemes. Therefore, central differencing is better than other schemes because it uses information from both the sides from the point of consideration to give approximative and accurate values.

As the left and right skewed difference methods uses data only from one side of point of interest it gives less approximative and accurate values about the point.