Week 3 - Taylor table method and Matlab code

INTRODUCTION

The Taylor table method was used to derive an expression for central-difference, skew-right and skew-left schemes which are fourth-order accurate.

DERIVATION OF SCHEMES

For a fourth-order central difference scheme the number of nodes to be considered = p+q-1. Where p is the order of the derivative and q is the order of approximation. The second-order derivative and a fourth-order scheme gives a value of p = 2 and q =4. Therefore the number of nodes required for a fourth-order central difference scheme of second derivative is 5

The fourth-order central difference scheme can be written as :

`\frac{d^2f(x)}{dx^2} = af(x-2dx) + bf(x -dx) + cf(x) + df(x+dx) +ef(x+2dx)`

The Taylor table for the above scheme is :

	`f(x)`	`f'(x)dx`	`f''(x)(dx)^2`	`f'''(x)(dx)^3`	`f^{iv}(x)(dx)^4`	`f^{v}(x)(dx)^5`
`af(x-2dx)`	`a`	`-2a`	`\frac{4a}{2}`	`\frac{-8c}{6}`	`\frac{16c}{24}`	`\frac{-32c}{120}`
`bf(x-dx)`	`b`	`-b`	`b/2`	`-b/6`	`b/24`	`-b/120`
`cf(x)`	`c`	0	0	0	0	0
`df(x+dx)`	`d`	`d`	`\frac{d}{2}`	`\frac{d}{6}`	`\frac{d}{24}`	`\frac{d}{120}`
`ef(x+2dx)`	`e`	`2e`	`\frac(4e)(2)`	`\frac{8e}{6}`	`\frac{16e}{24}`	`\frac{32e}{120}`
`Sigma_{n=0)^{n=n}f^n(x)dx^n`	0	0	1	0	0	0

The coefficients can be easily calculated since each of the columns represents a linear algebraic equation. This can be solved using a simple Matlab code

A =[1,1,1,1,1;
    -2,-1,0,1,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)]

B =[0;0;1;0;0]
A_inv = inv(A)
X = A_inv*B

The output of the above code is shown below:

For any skewed scheme the number of nodes to be considered = p+q. Where p is the order of the derivative and q is the order of approximation. The second-order derivative and a fourth-order scheme give a value of p = 2 and q =4. Therefore the number of nodes required for a fourth-order skewed difference scheme of the second derivative is 6.

The fourth-order skew-right scheme can be written as :

`\frac{d^2f(x)}{dx^2} = af(x) + bf(x +dx) + cf(x+2dx) + df(x+3dx) +ef(x+4dx) +gf(x+5dx)`

The Taylor table for this scheme can be represented as follows:

	`f(x)`	`f'(x)dx`	`f''(x)(dx)^2`	`f'''(x)(dx)^3`	`f^{iv}(x)(dx)^4`	`f^{v}(x)(dx)^5`
`af(x)`	`a`	0	0	0	0	0
`bf(x+dx)`	`b`	`b`	`b/2`	`b/6`	`b/24`	`b/120`
`cf(x+2dx)`	`c`	`2c`	`\frac{4c}{2}`	`\frac{8c}{6}`	`\frac{16c}{24}`	`\frac{32c}{120}`
`df(x+3dx)`	`d`	`3d`	`\frac{9d}{2}`	`\frac{27d}{6}`	`\frac{81d}{24}`	`\frac{243d}{120}`
`ef(x+4dx)`	`e`	`4e`	`\frac(16e)(2)`	`\frac{64e}{6}`	`\frac{256e}{24}`	`\frac{1024e}{120}`
`gf(x+5dx)`	`g`	`5g`	`\frac{25g}{2}`	`\frac{125g}{6}`	`\frac{625g}{24}`	`\frac{3125g}{120}`
Sigma_{n=0)^{n=n}f^n(x)dx^n	0	0	1	0	0	0

The coefficients can be found out as in the previous case, using Matlab. The code and the output is shown below:

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]

A_inv = inv(A)

X = A_inv*B

The output of the above code isshown below:

The fourth-order skew left scheme can be written as:

`\frac{d^2f(x)}{dx^2} = af(x) + bf(x -dx) + cf(x-2dx) + df(x-3dx) +ef(x-4dx) +gf(x-5dx)`

The Taylor table that can be used for deriving the above scheme is shown below:

	`f(x)`	`f'(x)dx`	`f''(x)(dx)^2`	`f'''(x)(dx)^3`	`f^{iv}(x)(dx)^4`	`f^{v}(x)(dx)^5`
`af(x)`	`a`	0	0	0	0	0
`bf(x-dx)`	`b`	`-b`	`b/2`	`-b/6`	`b/24`	`-b/120`
`cf(x-2dx)`	`c`	`-2c`	`\frac{4c}{2}`	`\frac{-8c}{6}`	`\frac{16c}{24}`	`\frac{-32c}{120}`
`df(x-3dx)`	`d`	`-3d`	`\frac{9d}{2}`	`\frac{-27d}{6}`	`\frac{81d}{24}`	`\frac{-243d}{120}`
`ef(x-4dx)`	`e`	`-4e`	`\frac(16e)(2)`	`\frac{-64e}{6}`	`\frac{256e}{24}`	`\frac{-1024e}{120}`
`gf(x-5dx)`	`g`	`-5g`	`\frac{25g}{2}`	`\frac{-125g}{6}`	`\frac{625g}{24}`	`\frac{-3125g}{120}`
Sigma_{n=0)^{n=n}f^n(x)dx^n	0	0	1	0	0	0

The coefficients are calculated using Matlab as in the previous case.

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]

A_inv = inv(A)

X = A_inv*B

The output of the above code is

RESULTS AND CODE

The code is written for implementing fourth-order accurate central difference, skew-right and skew-left schemes in Matlab.

clear all 
close all
clc


x = pi/3;
dx = linspace(pi/20,pi/2000,100)


for i=1:length(dx)
    
    error_central(i) = fourth_order_central(x,dx(i))
    error_skew_right(i) = fourth_order_skewright(x,dx(i))
    error_skew_left(i) = fourth_order_skewleft(x,dx(i))
end
figure(1)

loglog(dx,error_central,'b','linewidth',1);
hold on;

loglog(dx,error_skew_right,'r','linewidth',1);
hold on;

loglog(dx,error_skew_left,'g','linewidth',1);
hold off;
grid on
title('Absolute error in different schemes')
xlabel('log(dx)')
ylabel('log(error)')
legend('Fourth order central','Fourth order skew right','Fourth order skew left','Location','northwest')

The second-order derivative of the given function is calculated at x = `pi/3` using different schemes which are defined as separate functions.

The function calculating the analytical function is shown below:

function f = expcos(x)
f = cos(x)*exp(x);
end

The function calculating the second derivative of the above function is shown below:

function analy_der = dexpcos(x)

analy_der = -2*exp(x)*sin(x)

end

The fourth-order central difference scheme is implemented using the function fourth_order_central(x,dx) which is shown below:

function central_error = fourth_order_central(x,dx)

central_difference = -expcos(x-2*dx) + (16*expcos(x-dx)) - (30*expcos(x)) + (16*expcos(x+dx)) - expcos(x+2*dx);

central_difference = central_difference/(12*dx*dx);

central_error = abs(dexpcos(x)-central_difference);

end

The function defined in the above snippet of code returns the absolute error between the analytical value and numerical value calculated using the fourth-order central difference scheme.

The fourth-order skew left scheme is shown below:

function skew_left_error = fourth_order_skewleft(x,dx)

skew_left_difference = (45*expcos(x)) - (154*expcos(x-dx));

skew_left_difference = skew_left_difference + (214*expcos(x-2*dx)) - (156*expcos(x-3*dx));

skew_left_difference = skew_left_difference + (61*expcos(x-4*dx)) - (10*expcos(x-5*dx));

skew_left_difference = skew_left_difference/(12*dx*dx);

skew_left_error = abs(dexpcos(x)-skew_left_difference);

end

The function in the above snippet of code calculates the absolute error between the fourth-order skew-left scheme and the analytical derivative.

The function of the fourth-order skew right scheme is shown below:

function skew_right_error = fourth_order_skewright(x,dx)

skew_right_difference = (45*expcos(x)) - (154*expcos(x+dx));

skew_right_difference = skew_right_difference + (214*expcos(x+2*dx)) - (156*expcos(x+3*dx));

skew_right_difference = skew_right_difference + (61*expcos(x+4*dx)) - (10*expcos(x+5*dx));

skew_right_difference = skew_right_difference/(12*dx*dx);

skew_right_error = abs(dexpcos(x)-skew_right_difference);

end

The function in the above snippet of code calculates the absolute error between the fourth-order skew-right scheme and the analytical derivative.

The output of the program is shown below:

From the output of the program, it can be seen that the fourth-order central difference scheme is more accurate in most of the cases. At first, both skewed right and left schemes had the same accuracy but as the value of dx is increased the error in the skew right scheme is more than that in the skew left scheme.

A central difference scheme can be used only when there are points on both the left and right side of the given location. While evaluating properties at the boundary, when there are nodes only to the left or right the skew schemes can be used.

It is also interesting to note that the skewed schemes can be used when the property value at any given point is only determined by the properties at the left or right of the given point as in the problems where the direction of propagation of information is important.

CONCLUSION

From the above study, the derivation of the numerical scheme using Taylor Table is implemented.

It is seen that the central difference scheme happened to be more accurate than the skew left and skew right scheme.

All the different schemes have applications in different domains in solving problems in computational fluid mechanics.