To evaluate higher order approximations of a function by Finite Differencing Central and Skewed schemes using Taylor Table Method on MATLAB

1.Project Aim: To evaluate higher order approximations of the function by Finite Differencing (Central and Skewed schemes) using Taylor Table Method on MATLAB

2.Objectives of the Project :

• Derive the following 4th order approximations of the second-order derivative using: Central difference, Skewed right-sided difference and Skewed left-sided difference
• To show that skewed schemes are fourth-order accurate.
• Write a program in Matlab to evaluate the second-order derivative of the analytical function exp(x)*cos(x) and compare it with the 3 numerical approximations.
• To provide a plot that compares the absolute error between the above-mentioned schemes.
• To describe why a skewed scheme is useful? What can a skewed scheme do that a CD scheme cannot do?

3.Procedure :

3.1. Introduction:

The Taylor Table method is used to get the higher order approximations for a function by expanding it as per Taylor\'s series and then tabulating it. Suppose if \'n\' is the order of derivative that is to be approximated and \'p\' is the order of approximation, then the difference equation involving the skewed approximation will have (n+p) terms. for, this particular case, we wish to approximate 2nd order derivative (n=2) using 4th order approximations (p=4), which means that we would require 6 terms for the respective skewed approximation.

Consideration of the grid determines whether the grid is symmetrical or skewed. The figure below explains these types of 1D grids.

Let us try to understand the meaning of 4th order symmetric and skewed numerical stencils from the above figure. \'i\' is the point at which numerical derivative is to be calculated with other being its neighboring points.

• Fourth order symmetric scheme (Central Differencing Scheme) considers 2 points each at right and left side of the point \'i\'.
• Fourth order right skewed scheme (Forward Differencing Scheme) considers 4 points on the right side of the point \'i\'.
• Fourth order left skewed scheme (Backward differencing Scheme) considers 4 points on the left side of the point \'i\'.

3.2. Numerical discretization:

The above mentioned schemes are written in terms of linear equation of the function values at the neighboring points having constant coefficients as shown below. These coefficients are to be calculated using the Taylor table method.

• Central Differencing Scheme

`f\'\'(x)=(a_1f(x-2h)+b_1f(x-h)+c_1f(x)+d_1f(x+h)+e_1f(x+2h))/h^2`

• Left sided Skewed Differencing Scheme

`f\'\'(x)=(a_2f(x)+b_2f(x-h)+c_2f(x-2h)+d_2f(x-3h)+e_2f(x-4h)+g_2f(x-5h))/h^2`

• Right sided Skewed Differencing Scheme

`f\'\'(x)=(a_3f(x)+b_3f(x+h)+c_3f(x+2h)+d_3f(x+3h)+e_3f(x+4h)+g_3f(x+5h))/h^2`

The Taylor Table for the three schemes are as shown

• Central Differencing Scheme:

	`f(x)`	`Deltaxf^((1))(x)`	`Deltax^2f^((2))(x)`	`Deltax^3f^((3))(x)`	`Deltax^4f^((4))(x)`	`Deltax^5f^((5))(x)`	`Deltax^6f^((6))(x)`
`a_1f(x-2h)`	`a_1`	`-2a_1`	`2a_1`	`-(4a_1)/3`	`(2a_1)/3`	`-(4a_1)/15`	`(4a_1)/45`
`b_1f(x-h)`	`b_1`	`-b_1`	`b_1/2`	`-b_1/6`	`b_1/24`	`-b_1/120`	`b_1/720`
`c_1f(x)`	`c_1`	0	0	0	0	0	0
`d_1f(x+h)`	`d_1`	`d_1`	`d_1/2`	`d_1/6`	`d_1/24`	`d_1/120`	`d_1/720`
`e_1f(x+2h)`	`e_1`	`2e_1`	`2e_1`	`(4e_1)/3`	`(2e_1) / 3`	`(4e_1 )/ 15`	`(4e _1)/ 45`
Total	0	0	1	0	0	?	?

These coefficients can be solved in Matlab(in codes below) and the results are as shown

>> [a1 b1 c1 d1 e1]

ans =

   -0.0833    1.3333   -2.5000    1.3333   -0.0833

Central difference approximation can be given as:

`f\'\'(x)approx(-0.0833f(x-2h)+1.3333f(x-h)-2.50f(x)+1.3333f(x+h)-0.0833f(x+2h))/(h^2)+O(h^4)`

• Left-sided skewed scheme

These coefficients can be solved in Matlab(in codes below) and the results are as shown

>> [a2 b2 c2 d2 e2 g2]

ans =

  Columns 1 through 5

    3.7500  -12.8333   17.8333  -13.0000    5.0833

  Column 6

   -0.8333

Left sided skewed approximation can be given as :

`f\'\'(x)approx(3.75f(x)-12.8333f(x-h)+17.8333f(x-2h)-13.0f(x-3h)+5.0833f(x-4h)-0.833f(x-5h))/(h^2)+O(h^4)`

• Right-sided skewed scheme

These coefficients can be solved in Matlab(in codes below) and the results are as shown

>> [a3 b3 c3 d3 e3 g3]

ans =

  Columns 1 through 4

    3.7500  -12.8333   17.8333  -13.0000

  Columns 5 through 6

    5.0833   -0.8333

Right sided skewed approximation can be given as

`f\'\'(x)approx(3.75f(x)-12.8333f(x+h)+17.8333f(x+2h)-13.0f(x+3h)+5.0833f(x+4h)-0.833f(x+5h))/(h^2)+O(h^4)`

3.3 MALAB Codes

code for the main function

clear all
close all
clc
 % Calculation of the Coefficients for all the three schemes 
 
 % 1. Fourth order CDS
 % f\'\'(i) = a1*f(i-2) + b1*f(i-1) + c1*f(i) + d1*f(i+1) + e1*f(i+2);
A1 = [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];

B1 = [0;0;1;0;0];
X1 = A1\\B1;
a1 = X1(1);b1 = X1(2);c1 = X1(3);d1 = X1(4);e1 = X1(5);
 
 % 2. Fourth order left-sided skewed
 % f\'\'(i) = a2*f(i) + b2*f(i-1) + c2*f(i-2) + d2*f(i-3) + e2*f(i-4) + g2f(i-5);
A2 = [ 1        1         1        1        1       1;
       0        -1       -2       -3       -4      -5;   
       0        1/2       2        9/2      8       25/2;
       0        -1/6     -4/3     -9/2     -32/3   -125/6 ;
       0        1/24      2/3      27/8     32/3    625/24;
       0        -1/120   -4/15    -81/40   -128/15 -625/24];       

B2 = [0;0;1;0;0;0];
X2 = A2\\B2;
a2 = X2(1);b2 = X2(2);c2 = X2(3);d2 = X2(4);e2 = X2(5);g2 = X2(6);

 % 2. Fourth order right-sided skewed
 % f\'\'(i) = a3*f(i) + b3*f(i+1) + c3*f(i+2) + d3*f(i+3) + e3*f(i+4) + g3f(i+5);
A3 = [ 1        1         1        1        1       1;
       0        1         2        3        4       5;   
       0        1/2       2        9/2      8       25/2;
       0        1/6      4/3      9/2      32/3     125/6 ;
       0        1/24      2/3      27/8     32/3    625/24;
       0        1/120    4/15     81/40    128/15   625/24];       

B3 = [0;0;1;0;0;0];
X3 = A3\\B3;
a3 = X3(1);b3 = X3(2);c3 = X3(3);d3 = X3(4);e3 = X3(5);g3 = X3(6);

% Given function f(x) = exp(x) * cos(x)
% Analytical second order derivative f\'\'(x) = -2* exp(x)*sin(x)
x = pi/3;
dx = linspace(pi/400,pi/4,1000);
m1 = zeros(length(dx));
analytical_derivative = -2* exp(x)*sin(x);
% 
 for i = 1:length(dx)
     CDS_fourth_order_approx(i) = CDS(X1,x,dx(i));
     BDS_fourth_order_approx(i) = BDS(X2,x,dx(i));
     FDS_fourth_order_approx(i) = FDS(X3,x,dx(i));
      error_CDS(i) = abs(CDS_fourth_order_approx(i)- analytical_derivative);
      error_FDS(i) = abs(FDS_fourth_order_approx(i)- analytical_derivative);
      error_BDS(i) = abs(BDS_fourth_order_approx(i)- analytical_derivative);
 end

 % Plotting the errors
 figure(1)
 loglog(dx,error_CDS,\'LineWidth\',2 )
 hold on
 loglog(dx,error_FDS,\'LineWidth\',2 )
 hold on
 loglog(dx,error_BDS,\'LineWidth\',2 )
 grid on
 ax = gca; ax.FontSize = 14;
 lgd = legend (\'fourth order CDS\',\'fourth order Right Skewed\',\'fourth order Left Skewed\',\'location\',\'southeast\');
 lgd.FontSize = 14;
 xlabel(\'log(dx)\',\'FontSize\',14);
 ylabel(\'log(error)\',\'FontSize\',14);
 title(\'Error comparison of different schemes\',\'FontSize\',18);

code for CDS function

function CDS_fourth_approx = CDS(X,x,dx)

f_i_minus_2 = exp (x - 2*dx)* cos(x - 2*dx);
f_i_minus_1 = exp (x - dx)* cos(x - dx);
f_i = exp (x)* cos(x);
f_i_plus_1 = exp (x + dx)* cos(x + dx);
f_i_plus_2 = exp (x + 2*dx)* cos(x + 2*dx);


CDS_fourth_approx = (X(1)*f_i_minus_2 + X(2)*f_i_minus_1 + X(3)*f_i + X(4)*f_i_plus_1 + X(5)*f_i_plus_2)/dx^2;

code for BDS (left sided skewed) function

function BDS_fourth_approx = BDS(X,x,dx)



f_i_minus_5 = exp (x - 5*dx)* cos(x - 5*dx);
f_i_minus_4 = exp (x - 4*dx)* cos(x - 4*dx);
f_i_minus_3 = exp (x - 3*dx)* cos(x - 3*dx);
f_i_minus_2 = exp (x - 2*dx)* cos(x - 2*dx);
f_i_minus_1 = exp (x - dx)* cos(x - dx);
f_i = exp (x)* cos(x);

BDS_fourth_approx = (X(6)*f_i_minus_5 + X(5)*f_i_minus_4 + X(4)*f_i_minus_3 + X(3)*f_i_minus_2 + X(2)*f_i_minus_1 + X(1)*f_i)/dx^2 ;

code for FDS function

function FDS_fourth_approx = FDS(X,x,dx)

f_i = exp (x)* cos(x);
f_i_plus_1 = exp (x + dx)* cos(x + dx);
f_i_plus_2 = exp (x + 2*dx)* cos(x + 2*dx);
f_i_plus_3 = exp (x + 3*dx)* cos(x + 3*dx);
f_i_plus_4 = exp (x + 4*dx)* cos(x + 4*dx);
f_i_plus_5 = exp (x + 5*dx)* cos(x + 5*dx);


FDS_fourth_approx = (X(1)*f_i + X(2)*f_i_plus_1 + X(3)*f_i_plus_2 + X(4)*f_i_plus_3 + X(5)*f_i_plus_4 + X(6)*f_i_plus_5)/dx^2;

4. Results:

Figure1. Comparision for the absolute errors for the given scheme at the same grid size indicated in codes

As shown from the figure, we can say that the fourth order CDS is more consistent and accurate than the skewed schemes. This can be very well explained by the fact that it uses the points very close to the point \'i\' (i.e. i-2, i-1, i, i+1, i+2) for calculations whereas skewed schemes use points that are far away from the point \'i\'. Moreover the constant slope of CDS explains the fact about consistency of the CDS at higher values of dx.

Advantages of the skewed schemes over Symmetric schemes:

Though, the error of skewed schemes is greater than that of symmetric scheme, they have some advantages associated with them. In case of boundary nodes, where the symmetric schemes cannot be used, skewed schemes can be used by making sure that the grid is fine enough to reduce the error and to prevent the skewed schemes from going in consistent.