Taylor table method and Matlab code

Aim:

1. To derive the fourth-order approximation of second-order derivative using central difference, skewed right-sided, and skewed left-sided difference scheme.
2. To prove that the skewed schemes are fourth-order accurate
3. To write a program in MATLAB/Octave to evaluate the second-order derivative of the analytical function: exp(x)*cos(x) and compare it with the above-mentioned approximation methods

Analytical solution:

f(x) =`e^x cos⁡x`

By using the product rule, f’(x), and f’’(x) can be found

f’(x) =`-e^x.sin⁡(x)+cos⁡(x).(e^x)`

f’’(x) =`-[e^x.cos⁡x+e^x.sin⁡x]+[-e^x.sin⁡x+e^x.cos⁡x]`

                =`-2e^x.sin⁡x`

Number of nodes calculation:

1. For skewed difference schemes, n = p + q

Where p = order of derivative

          q = order of approximation.

Hence, n = 2+4 = 6

2. For central difference schemes, n = p+q-1

Hence, n = 2+4-1 = 5

Taylor table for the skewed right-sided difference scheme:

The total number of nodes required to solve this problem is 6

 From this we get six unknown variables and six equations:

a+b+c+d+e+g = 0

b+2c+3d+4e+5g = 0

1/2b+4/2c+9/2d+16/2e+25/2g =  1

1/6b+8/6c+27/6d+64/6e+125/6g = 0

1/24b+16/24c+81/24d+256/24e+625/24g = 0

1/120b+32/120c+243/120d+1024/120e+3125/120g = 0

Solving in matlab/octave, we get the value of the variables:

a =3.75, b = -12.833, c = 17.833, d = -13, e = 5.0833, g = -0.833

Hence,

`(d^2 u)/(dt^2 )= (af(i)+bf(i+1)+cf(i+2)+df(i+3)+ef(i+4)+gf(i+5))/(Δx^2 )`

`= (3.75f(i)-12.833f(i+1)+17.833f(i+2)-13f(i+3)+5.0833f(i+4)-0.833f(i+5))/(Δx^2 )`

Taylor table for the skewed left-sided difference scheme:

The total number of nodes required to solve this problem is 6

 From this we get six unknown variables and six equations:

a+b+c+d+e+g = 0

-5a-4b-3c-2d-e = 0

25/2a+16/2b+9/2c+4/2d+1/2e = 1

-125/6a-64/6b-27/6c-8/6d-1/6e = 0

625/24a+256/24b+81/24c+16/24d+1/24e = 0

-3125/120a-1024/120b-243/120c-32/120d-1/120e = 0

Solving in matlab/octave, we get the value of the variables:

a = -0.833, b = 5.0833, c = -13, d = 17.833, e = 12.833, g = 3.75

Hence,

 `(d^2 u)/(dt^2 )= (af(i-5)+bf(i-4)+cf(i-3)+df(i-2)+ef(i-1)+gf(i))/(Δx^2 )`

`= (-0.833f(i-5)+5.0833(f-4)-13f(i-3)+17.833f(i-2)+12.833f(i-1)+3.75f(i))/(Δx^2 )`

Taylor table for Central difference scheme:

The total number of nodes required to solve this problem is 5

 From this we get five unknown variables and five equations:

a+b+c+d+e = 0

-2a-b+d+2e = 0

4/2a+1/2b+1/2d+4/2e = 1

-8/6a-1/6b+1/6d+8/6e = 0

16/24a+1/24b+1/24d+16/24e = 0

Solving in matlab/octave, we get the value of the variables:

a = 0.28571, b = -0.14286, c = -0.28571, d = -0.14286, e = 0.28571

here, a = e, and b = d

Hence,

`(d^2 u)/(dt^2 )= (af(i-2)+bf(i-1)+cf(i)+df(i+1)+ef(i+2))/(Δx^2 )`

`= (0.28571f(i-2)-0.14286f(i-1)-0.28571f(i)-0.14286f(i+1)+0.28571f(i+2))/(Δx^2 )`

Proof skewed schemes are fourth-order accurate:

1) Skewed right scheme

WKT,

`Δx^2 (d^2 u)/(dt^2 )= af(i)+bf(i+1)+cf(i+2)+df(i+3)+ef(i+4)+gf(i+5)`

From the above Taylor table, we get,

`Δx^2 (d^2 u)/(dt^2 )= (a+b+c+d+e+g)f(i) + (b+2c+3d+4e+5g)f’(i)Δx + (b/2+2c+9/2d+8e+25/2g)f’’(i)Δx^2 + (1/6b+8/6c+27/6d+64/6e+125/6g)f3Δx^3 + (1/24b+16/24c+81/24d+256/24e+625/24g)f4Δx^4 + (1/120b+32/120c+243/120d+1024/120e+3125/120g)f5Δx^5 + O(Δx^6)`

From our assumption from above Taylor table to get the value of unknowns, we assumed that the sum of terms of orders 0,1,3,4,5 is zero. Only sum of terms of order2, and sum of terms of order six and above are non-zero. Hence,

`Δx^2 (d^2 u)/(dt^2 )= (b/2+2c+9/2d+8e+25/2g)f’’(i)Δx^2 + O(Δx^6)`

Dividing throughout by ,`Δx^2`

`(d^2 u)/(dt^2 )= (b/2+2c+9/2d+8e+25/2g)f’’(i) + O(Δx^4)`

O(`Δx^4`)  proves that the right skewed scheme is fourth order accurate.

2) Skewed left scheme

WKT,

`Δx^2 (d^2 u)/(dt^2 )= af(i-5)+bf(i-4)+cf(i-3)+df(i-2)+ef(i-1)+gf(i)`

From the above Taylor table, we get,

`Δx^2 (d^2 u)/(dt^2 )= (a+b+c+d+e+g)f(i) - (5a+4b+3c+2d+e)f’(i)Δx + (25a/2+16b/2+9c/2+4/2d+0.5e)f’’(i)Δx^2 - (125/6a+64/6b+27/6c+8/6d+1/6e)f3Δx^3 + (625/24a+256/24b+81/24c+16/24d+1/24e)f4Δx^4- (3125/120a+1024/120b+243/120c+32/120d+1/120e)f5Δx^5 + O(Δx^6)`

From our assumption from above Taylor table to get the value of unknowns, we assumed that the sum of terms of orders 0,1,3,4,5 is zero. Only sum of terms of order2, and sum of terms of order six and above are non-zero. Hence,

`Δx^2 (d^2 u)/(dt^2 )= (25/2a+16/2b+9/2c+2d+0.5e)f’’(i)Δx^2 + O(Δx^6)`

Dividing throughout by ,`Δx^2`

`(d^2 u)/(dt^2 )= (25/2a+16/2b+9/2c+2d+0.5e)f’’(i) + O(Δx^4)`

O(`Δx^4`)  proves that the right-skewed scheme is fourth-order accurate.

MATLAB/Octave code:

1) Central difference scheme

function out = central_difference(x,dx)

% Function = exp(x)*cos(x);

% Second_order_analytical_derivative;
Second_order_analytical_derivative = -2*sin(x)*exp(x);

% Numerical_derivative;

% To find coefficients
Ac = [1 1 1 1 1;
-2 -1 0 1 2;
2 0.5 0 0.5 2;
-8/6 -1/6 0 1/6 8/6;
16/24 1/24 0 1/24 16/24];

Bl = [0;0;1;0;0];

Xc = AcBl

central_difference = (Xc(1)*exp(x-2*dx)*cos(x-2*dx)+Xc(2)*exp(x-dx)*cos(x-dx)+Xc(3)*exp(x)*cos(x)+Xc(4)*exp(x+dx)*cos(x+dx)+Xc(5)*exp(x+2*dx)*cos(x+2*dx))/dx^2;

out = abs(central_difference-Second_order_analytical_derivative);

end

This code is used to compute coefficients in the central difference scheme. The result of the coefficients is stored in the array Xc. These coefficients are used to find central difference approximation by substituting them in the equation we derived earlier. To reduce round-up error, the coefficient terms are called directly from Xc array, rather than substituting the value manually. Once the central difference is found, the absolute difference between central approximation and analytical solution is found.

2) Skewed right-sided difference scheme

function out = skewed_right_side(x,dx)
% f(x) = exp(x)*cos(x);

% Second_order_analytical_derivative;
second_order_analytical_derivative=-2*sin(x)*exp(x);

% Numerical_derivative;
% finding coefficients

Ar = [1 1 1 1 1 1;
0 1 2 3 4 5;
0 0.5 2 4.5 8 12.5;
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];

Xr = ArB

skewed_right_side = (Xr(1)*exp(x)*cos(x)+Xr(2)*exp(x+dx)*cos(x+dx)+Xr(3)*exp(x+2*dx)*cos(x+2*dx)+Xr(4)*exp(x+3*dx)*cos(x+3*dx)+Xr(5)*exp(x+4*dx)*cos(x+4*dx)+Xr(6)*exp(x+5*dx)*cos(x+5*dx))/(dx^2);

out = abs(skewed_right_side-second_order_analytical_derivative);

end

This code is used to compute coefficients in the skewed right difference scheme. The result of the coefficients is stored in the array Xr. These coefficients are used to find the right-skewed difference approximation by substituting them in the equation we derived earlier. To reduce round-up error, the coefficient terms are called directly from Xc array, rather than substituting them manually. Once skewed right approximation is found, the absolute difference between skewed right approximation and analytical solution is found.

3) Skewed left-sided difference scheme

function out = skewed_left_side(x,dx)
% f(x) = exp(x)*cos(x);

% Second_order_analytical_derivative;
second_order_analytical_derivative = - 2*sin(x)*exp(x);

% Numerical_derivative;

% finding coefficients

Al = [1 1 1 1 1 1;
-5 -4 -3 -2 -1 0;
25/2 16/2 9/2 2 0.5 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];

B = [0;0;1;0;0;0];

Xl = AlB

skewed_left_side = (Xl(1)*exp(x-5*dx)*cos(x-5*dx)+Xl(2)*exp(x-4*dx)*cos(x-4*dx)+Xl(3)*exp(x-3*dx)*cos(x-3*dx)+Xl(4)*exp(x-2*dx)*cos(x-2*dx)+Xl(5)*exp(x-dx)*cos(x-dx)+Xl(6)*exp(x)*cos(x))/dx^2;

out = abs(skewed_left_side-second_order_analytical_derivative);

end

This code is used to compute coefficients in a skewed left difference scheme. The result of the coefficients is stored in the array Xl. These coefficients are used to find left-skewed difference approximation by substituting them in the equation we derived earlier. To reduce round-up error, the coefficient terms are called directly from the Xl array, rather than substituting them manually. Once skewed left approximation is found, the absolute difference between skewed left approximation and analytical solution is found.

Main Code:

clear all
close all
clc

% Input parameters
x = pi/3;
dx = linspace(pi/4,pi/4000,100);

% f(x)=exp(x)*cos(x);

% Second_order_analytical_derivative;
second_order_analytical_derivative = -2*sin(x)*exp(x);

% For loop;

for i = 1:length(dx)

central_difference_error(i) = central_difference(x,dx(i));

skewed_right_side_error(i) = skewed_right_side(x,dx(i));

skewed_left_side_error(i) = skewed_left_side(x,dx(i));

end

% Plotting;

 figure(1)

 loglog(dx,central_difference_error,'r')
 hold on

 loglog(dx,skewed_right_side_error,'g')
 hold on

 loglog(dx,skewed_left_side_error,'b')

 % Labeling;
 xlabel('dx')
 ylabel('error')
 legend('central difference error','skewed right side error','skewed left side error')

In this program, the value of x is given and the range of values of dx is defined. The above programs are called within this program and finally a graph is plotted which gives us a clear idea of the magnitude of error with the increase/decrease in value of dx. Also we can find which approximation has the least magnitude of error for a particular value of dx.

Results:

The main code was run for 100 iterations, with x = pi/3 and dx ranging from pi/4 to pi/4000, and the following solution(plot) was obtained,

By reading the graph we can say that

a) The magnitude of error of the central difference scheme is comparatively low than the skewed difference schemes fo most values of dx

b) The magnitude of the error is almost the same for both the skewed difference schemes, except for a particular value of dx, the skewed left difference error becomes lower than that of the central difference scheme. so if we want to solve the problem for that particular value of dx, then a skewed left-sided difference scheme may be used.

Now, the same code was run for 1000 iterations, with x = pi/3 and dx ranging from pi/4 to pi/40000, and the following plot was obtained,

Again, the same code was run for 10000 iterations, with x = pi/3 and dx ranging from pi/4 to pi/400000, and the following plot was obtained,

Now, when one sees the first plot only, s/he may conclude that the magnitude of error reduces with the value of dx. But this is not true. As the value of dx decreases from 1, the magnitude of error also decreases, up to around dx ~ 1e-2. When the value of dx further decreases the magnitude of error also increases, as shown in plot 2. This is due to the introduction of an error term called 'round-off error' with the decrease in the value of dx. 

From the third plot, we can see that the magnitude of error at around dx~ 1e-3 is inconsistent. There is a brief region around dx~1e-3 that is unstable. And with the further reduction in dx, the magnitude of the error increases. Now for any CFD simulations, the unsteady nature of error magnitude is undesirable, so it is advisable to keep the value of dx within the steady region, where the magnitude of error can be predicted.

Use of skewed scheme:

Skewed difference schemes can be used to calculate the magnitude of error at the boundary since they require the nodes only on one side of the boundary. On the other hand, the central difference scheme cannot be applied on the boundary as it requires nodes on either side of the boundary.

Conclusion:

After completing this assignment, we understood the following concepts:

• Even though the central difference scheme and both the skewed difference schemes have the same order of approximation, the magnitude of error obtained by the central difference scheme is comparatively lower than skewed schemes.
• The number of nodes required for computing the central difference scheme is one less than that of skewed difference schemes for second-order derivatives.
• Even though the central difference scheme has a lower magnitude of error compared to the skewed difference scheme, it cannot be used to calculate the error at boundaries as it requires nodes on either side. Also, skewed difference schemes may show a lower magnitude of the error when it is applied for flow problems
• For all the approximation schemes, the magnitude of error remains consistent for only up to a particular value of dx, further reduction in dx makes the error magnitude destabilize, and hence it cannot be predicted in advance, so the value is usually fixed within the steady zone, where the magnitude of error can be predicted, in this example from 1 to 1e-2.