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

Aim:

Deriving 4th order approximation of a 2nd order derivative using the Taylor Table method.

Objective:

To derive the following 4th order approximation of the second order derivative using

Central difference

Skewed right-sided difference

Skewed left-sided difference

The first step is to find the number of nodes required for deriving each difference method using the given formula.

After deriving all the different methods, the other step is to compare all the derived difference methods with the exact solution and find the error range for each different way using MATLAB.

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

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

1. Central difference:

Finding the no of nodes:

The order of accuracy P = Number of nodes in the stencil N - order of derivative Q + 1 (symmetric/central schemes)

Than the Number of nodes N = The order of accuracy P + Order of derivative Q - 1

N = 4+2-1

N = 5

Derivation:

To derive fourth-order approximations for second-order derivatives, we start with these approximations, Here the nodes selected are symmetrical i.e. the data is taken from both sides of the point of focus. Using Taylor’s method, we get

`(d^2f)/dx^2 = af(i-2) + bf(i-1) + cf(i) + df(i+1) + ef(i+2)`

Taylor series expansion for the above terms:

`af(i-2) = af(i) - (2dxf'(i))/1 + (4adx^2f''(i))/2 - (4adx^3f'''(i))/3 + (2adx^4f''''(i))/3`

`bf(i-1) = bf(i) - (bdxf'(i))/1 + (bdx^2f''(i))/2 - (bdx^3f'''(i))/6 + (bdx^4f''''(i))/24`

`cf(i) =cf(i)`

`df(i+1) = df(i) + (d*dxf'(i))/(1!) + (d*dx^2f''(i))/(2!) + (d*dx^3f'''(i))/(3!) + (d*dx^4f''''(i))/(4!)`

`ef(i+2) = ef(i) + (2edxf'(i))/(1!) + (4edx^2f''(i))/(2!) + (8edx^3f'''(i))/(3!) + (16edx^4f''''(i))/(4!)`

Representing in the form of a table:

Calculate the values of a, b, c, d, and e, then you get the required approximation. These values can be found by using the Taylor table. Here there are two nodes taken from both sides of the point of focus i.e. f(i).

`dx^2f''(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-efficient	0	0	1	0	0

`A = [[1,1,1,1,1],[-2,-1,0,0,2],[2,1/2,0,1/2,2],[4/3,1/6,0,1/6,4/3],[2/3,1/24,0,1/24,2/3]]`

`X = [[a],[b],[c],[d],[e]]`

`B = [[0],[0],[1],[0],[0]]`

These coefficients will be further used in code to get the value of the derivative.

2. Skewed right-sided difference

The order of accuracy P = Number of nodes in the stencil N - order of derivative Q for one-sided/skewed schemes.

Here N = 6

`(d^2f)/dx^2 = af(i) + bf(i+1) + cf(i+2) + df(i+3) + ef(i+4) + ff(i+5)`

When the data point of the left nodes is not available (at initial boundary conditions) then we have to rely on the available data set.

Representing in the form of a table:

Calculate the values of a, b, c, d, and e, then you get the required approximation. These values can be found by using the Taylor table.

	`f(x)`	`dxf'(x)`	`dx^2f''(x)`	`dx^3f'''(x)`	`dx^4f''''(x)`	`dx^5f'''''(x)`
`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`
`ff(i+5)`	f	5f	`(25f)/2`	`(125f)/6`	`(625f)/24`	`(3125f)/120`
Sum of co-efficient	0	0	1	0	0	0

`A = [[1,1,1,1,1,1],[0,1,2,3,4,5],[0,1/2,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]]`

`X = [[a],[b],[c],[d],[e],[f]]`

`B = [[0],[0],[1],[0],[0],[0]]`

These coefficients will be further used in code to get the value of the derivative.

3. Skewed left-sided difference:

The order of accuracy P = Number of nodes in the stencil N - order of derivative Q for one-sided/skewed schemes.

Here N = 6.

`(d^2f)/dx^2 = af(i-5) + bf(i-4) + cf(i-3) + df(i-2) + ef(i-1) + ff(i)`

When the data point of the right nodes is not available (at boundary conditions) then we have to rely on the available data set.

Representing in the form of a table.

Calculate the values of a, b, c, d, and e, then you get the required approximation. These values can be found by using the Taylor table.

	`f(x)`	`dxf'(x)`	`dx^2f''(x)`	`dx^3f'''(x)`	`dx^4f''''(x)`	`dx^5f'''''(x)`
`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)/12`
`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`
`ff(i)`	f	0	0	0	0	0
Sum of coefficient	0	0	1	0	0	0

`A = [[1,1,1,1,1,1],[-5,-4,-3,-2,-1,0],[25/2,16/2,9/2,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]]`

`X = [[a],[b],[c],[d],[e],[f]]`

`B = [[0],[0],[1],[0],[0],[0]]`

These coefficients will be further used in code to get the value of the derivative.

Function code for central differencing:

function Error_central_differencing = fourth_order_central_differencing(x,dx);

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

%analytical 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_central_differencing=abs(central_differencing-analytical_derivative);
end

Function code for Right-sided Skew difference:

function Error_right_sided_skew = fourth_order_right_sided_skew(x,dx);
%analytic function = exp(x)*cos(x)

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

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

%numerical derivative

%right sided skew
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_right_sided_skew=abs(right_sided_skew-analytical_derivative);
end

Function code for left-sided Skew difference:

function Error_left_sided_skew = fourth_order_left_sided_skew(x,dx);
%analytic function = exp(x)*cos(x)

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

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

%numerical derivative

%Left sided skew
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_left_sided_skew=abs(left_sided_skew-analytical_derivative);
end

Matlab Program main program:

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,'b')
hold on
loglog(dx,right_sided_skew,'r')
hold on
loglog(dx,left_sided_skew,'g')
xlabel('range of dx')
ylabel('truncation error')
title('range of dx vs error')
legend('central differencing','skewed rigth sided difference','skewed left sided difference','location','northwest');

Screenshot of the main program:

Steps involved:

1. We assign the value of x and the range of dx
2. We define an anonymous function of analytical function and its derivative.
3. Use the matrices from the above-discussed section and solved it using the inversion method i.e.

AX=B ==> X=A^-1 B

And stored the coefficients.

4. Assigned the rows & columns to the appropriate matrices.
5. Loop the order of approximation expression in a for loop under the range of dx. Here absolute error is also calculated from the derivative value calculated analytically.
6. The calculated values are plotted to do a better comparison.

Output:

By seeing the plots that how the error is varying as a factor of given dx for various fourth order schemes, The error is decreasing as the grid size is decreasing & also central difference scheme is more accurate than the skewed right side and skewed left side scheme. The higher the order of the method the lower the error is going to be for the same grid size. But in reality when the geometry becomes very complex and when the flow becomes very complex you cannot use the higher-order method always, and in order for the simulation to run the commercial codes always go to first-order schemes in order to get a solution.

Conclusion:

As seen from the plot the central differencing scheme has very less error when compared to the other schemes. Therefore, the central differencing scheme is better than other schemes because it uses information from both sides from the point of consideration to giving appropriate and accurate values. As the left & right skewed difference methods use data only from one side of the point of interest it gives less approximation and accurate values about that point.

The system doesn't contain sufficient information on both sides of the point of interest in such cases we cannot use a central differencing scheme so the skewed scheme is useful in such cases. When we want to calculate the value of a function at the boundary we cannot use central difference schemes because it requires the value at the nodes on both sides whereas in skewed schemes only one-sided values are used. This is where skewed schemes come into play. This is what a skewed scheme can do that a central difference scheme cannot, but when it comes to accuracy point of view as we have seen above central difference is more accurate than skewed schemes for the same order of error.