Deriving and comparing the fourth order approximation of the 2nd order derivative

The following 4th order approximations of the second-order derivative. 

1. Central difference - we will consider f\'\'(x) as:

  f\'\'(x) = a*f(i-2) + b*f(i-1) + c*f(i) + d*f(i+1) + e*f(i+2)           -------(eq.1)

now, by using Taylor\'s method we will get a system of five simultaneous equation of a,b,c,d,e which is solved by matrix inversion and multiplication method, done in MATLAB for which the respective code is:

Here, A is the coefficient matrix, B is the matrix containing right-sided element ,and X is the solution matrix.

A =

    1.0000    1.0000    1.0000    1.0000    1.0000
   -2.0000   -1.0000         0    1.0000    2.0000
    2.0000    0.5000         0    0.5000    2.0000
   -1.3333   -0.1667         0    0.1667    1.3333
    0.6667    0.0417         0    0.0417    0.6667

B =

     0
     0
     1
     0
     0

X = inv(A) * B

X =

   -0.0833
    1.3333
   -2.5000
    1.3333
   -0.0833

Thus the value of coefficients are :

a = -0.0833

b =  1.3333

c = -2.5000

d =  1.3333

e = -0.0833

substituting the values in (eq.1) the equation becomes zero, so

Rewriting the order term and dividend the equation by we get

Now the truncation error term denotes that the order of the approximation is 4th order.

2. Skewed right-sided difference - Our scheme for f\"\"(x) will be:

  f\'\'(x) = a*f(i) + b*f(i+1) + c*f(i+2) + d*f(i+3) + e*f(i+4)           -------(eq.2)

now, by using Taylor\'s method we will get a system of five simultaneous equation of a,b,c,d,e which is solved by matrix inversion and multiplication method, done in MATLAB for which the respective code is:

Here, A is the coefficient matrix, B is the matrix containing right-sided element, and X is the solution matrix.

>> A = [1, 1, 1, 1, 1;0, 1, 2, 3, 4;0, 0.5, 2, 4.5, 8;0, 1/6, 8/6, 27/6, 64/6;0, 1/24, 16/24, 81/24,256/24]

A =

    1.0000    1.0000    1.0000    1.0000    1.0000
         0    1.0000    2.0000    3.0000    4.0000
         0    0.5000    2.0000    4.5000    8.0000
         0    0.1667    1.3333    4.5000   10.6667
         0    0.0417    0.6667    3.3750   10.6667

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

B =

     0
     0
     1
     0
     0

>> X = inv(A) * B

X =

    2.9167
   -8.6667
    9.5000
   -4.6667
    0.9167

Thus the value of coefficient are:

a =  2.9167

b = -8.6667

c =  9.5000

d = -4.6667

e =  0.9167

thus (eq.2) will become :

`(del^2f)/(delx^2) = (2.9167*f(i) - 8.6667*f(i+1) + 9.5*f(i+2) - 4.6667*f(i+3) + 0.9167f(i+4))/(Deltax^2)`

3. Skewed left-sided difference -Our scheme for f\"\"(x) will be:

  f\'\'(x) = a*f(i) + b*f(i-1) + c*f(i-2) + d*f(i-3) + e*f(i-4)           -------(eq.2)

now, by using Taylor\'s method we will get a system of five simultaneous equation of a,b,c,d,e which is solved by matrix inversion and multiplication method, done in MATLAB for which the respective code is:

Here, A is the coefficient matrix, B is the matrix containing right-sided element, and X is the solution matrix.

>> A = [1, 1, 1, 1, 1;0, -1, -2, -3, -4;0, 0.5, 2, 4.5, 8;0, -1/6, -8/6, -27/6, -64/6;0, 1/24, 16/24, 81/24,256/24]

A =

    1.0000    1.0000    1.0000    1.0000    1.0000
         0   -1.0000   -2.0000   -3.0000   -4.0000
         0    0.5000    2.0000    4.5000    8.0000
         0   -0.1667   -1.3333   -4.5000  -10.6667
         0    0.0417    0.6667    3.3750   10.6667

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

B =

     0
     0
     1
     0
     0

>> X = inv(A) * B

X =

    2.9167
   -8.6667
    9.5000
   -4.6667
    0.9167

Thus the value of coefficient are:

a =  2.9167

b = -8.6667

c =  9.5000

d = -4.6667

e =  0.9167

thus (eq.2) will become :

`(del^2f)/(delx^2) = (2.9167*f(i) - 8.6667*f(i-1) + 9.5*f(i-2) - 4.6667*f(i-3) + 0.9167f(i-4))/(Deltax^2)`

Now we are going to compare the absolute Value of errors obtained by Central difference, Skewed right-sided difference and, Skewed left-sided difference for 4th order approximations of the second-order derivative. 

By using the following MATLAB code we can obtain our results:

%% calculating the value of function
function y = value_fxn(x)
y = exp(x) * cos(x);
end

%% Function to calculate truncation error for the central difference

function out  = central_diff_error(x, dx)
%% Analytical Value
exact_value = -2*exp(x)*sin(x);
%% Central difference value for (i-2,i-1,i,i+1,i+2)
central = (-1*value_fxn(x-2*dx) + 16*value_fxn(x-dx) - 30*value_fxn(x) + 16*value_fxn(x+dx) - 1*value_fxn(x+2*dx))/(12*(dx^2));
%% absolute truncuation error
out = abs(central - exact_value);
end

%% Function to calculate truncation error for the skewed right-sided difference

function out  = skewed_right_diff_error(x, dx)
%% Analytical Value
exact_value = -2*exp(x)*sin(x);
%% for skewed right sided differnce for (i,i+1,i+2,i+3,i+4)
skewed_right = (35*value_fxn(x) - 104*value_fxn(x+dx) + 114*value_fxn(x+2*dx) - 56*value_fxn(x+3*dx) + 11*value_fxn(x+4*dx))/(12*(dx^2));
%% calculating absolute error
out = abs(skewed_right - exact_value);
end


%% Function to calculate truncation error for the skewed left-sided difference

function out  = skewed_left_diff_error(x, dx)
%% Analytical Value
exact_value = -2*exp(x)*sin(x);
%% for skewed left sided differnce for (i,i+1,i+2,i+3,i+4)
skewed_left = (35*value_fxn(x) - 104*value_fxn(x-dx) + 114*value_fxn(x-2*dx) - 56*value_fxn(x-3*dx) + 11*value_fxn(x-4*dx))/(12*(dx^2));
%% calculating absolute truncuation error
out = abs(skewed_left - exact_value);
end


%% Main run file
close all;
clear all;
clc;

%% Assuming x = pi/4
x = pi/4;
%% now defining step size range
dx = linspace(pi/10,pi/1000,100);

%% calculating error
for i =1 : length(dx)
    error_central(i) = central_diff_error(x,dx(i));
    error_right(i) = skewed_right_diff_error(x,dx(i));
    error_left(i) = skewed_left_diff_error(x,dx(i));
end
%% ploting the required graph
loglog(dx,error_central,\'k\',\'linewidth\',1.5);
hold on;
loglog(dx,error_right,\'r-*\',\'linewidth\',1.5);
hold on;
loglog(dx,error_left,\'b--\',\'linewidth\',1.5);

xlabel(\'step size\');
ylabel(\'error due to truncation\');
title(\'Comparison of accuracy for differend order of scheme for second order derivative\');
legend(\'Central difference order\',\'Skewed right sided order\',\'Skewed left sided error\');

This is the log plot obtained from the above MATLAB function which is shown below:

We can following conclusion from the above figure are:

1. It is observed that the central difference scheme always gives the most accurate or less truncation error at any step size when compared to skewed right-sided and skewed left-sided.

2. For a relatively high value of step size, we can observe that error in skewed right-sided is less that skewed left-sided scheme.

3. Although we observed that the central difference scheme is most accurate between all three but in certain cases when complexity increases central difference scheme may result in a larger error than skewed right-sided and skewed left-sided scheme.