Matlab code to evaluate the performace of differencing schemes on solution accuracy for a second order derivative of a function

Objectives:

• Derive the 4th order approximations of the second order derivative using:

1. Central difference (CD)

2. Skewed right sided difference (S-RSD)

3. Skewed left sided difference (S-LSD)

• Write a Matlab code to evaluate the second order derivative of the analytical function exp(x)*cos(x) and compare it with the 3 numerical approximations derived.
• Compare the errors given by 4th order approximations and 1st order approximations i.e. Forward Differencing schemes (FDS) and Backward Differencing Scheme (BDS).
• Describe the advantages of skewed differencing scheme over Central differencing scheme.

Introduction:

The taylor series expansion is given by the following expression:

`f(x+ trianglex) = f(x) + trianglex f\'(x) + ((trianglex)^2)/(2!) f\'\'(x) + ((trianglex)^3)/(3!) f\'\'\'(x) + ((trianglex)^4)/(4!) f\'\'\'\'(x) + ((trianglex)^5)/(5!) f\'\'\'\'\'(x)`

Since, we need to do the 4th order approximation, information from 5 nodes on either side of the zero node is required.

Using the taylor series expansion, the following terms are estimated:

`f(i-5) = f(i) + (-5 trianglex) f\'(i) + ((-5trianglex)^2)/(2!) f\'\'(i) + ((-5trianglex)^3)/(3!) f\'\'\'(i) + ((-5trianglex)^4)/(4!) f\'\'\'\'(i) + ((-5trianglex)^5)/(5!) f\'\'\'\'\'(i)`

`f(i-4) = f(i) + (-4trianglex) f\'(i) + ((-4trianglex)^2)/(2!) f\'\'(i) + ((-4trianglex)^3)/(3!) f\'\'\'(i) + ((-4trianglex)^4)/(4!) f\'\'\'\'(i) + ((-4trianglex)^5)/(5!) f\'\'\'\'\'(i)`

`f(i-3) = f(i) + (-3trianglex) f\'(i) + ((-3trianglex)^2)/(2!) f\'\'(i) + ((-3trianglex)^3)/(3!) f\'\'\'(i) + ((-3trianglex)^4)/(4!) f\'\'\'\'(i) + ((-3trianglex)^5)/(5!) f\'\'\'\'\'(i)`

`f(i-2) = f(i) + (-2trianglex) f\'(i) + ((-2trianglex)^2)/(2!) f\'\'(i) + ((-2trianglex)^3)/(3!) f\'\'\'(i) + ((-2trianglex)^4)/(4!) f\'\'\'\'(i) + ((-2trianglex)^5)/(5!) f\'\'\'\'\'(i)`

`f(i-1) = f(i) + (-trianglex) f\'(i) + ((-trianglex)^2)/(2!) f\'\'(i) + ((-trianglex)^3)/(3!) f\'\'\'(i) + ((-trianglex)^4)/(4!) f\'\'\'\'(i) + ((-trianglex)^5)/(5!) f\'\'\'\'\'(i)`

`f(i) = f(i)`

`f(i+1) = f(i) + (trianglex) f\'(i) + ((trianglex)^2)/(2!) f\'\'(i) + ((trianglex)^3)/(3!) f\'\'\'(i) + ((trianglex)^4)/(4!) f\'\'\'\'(i) + ((trianglex)^5)/(5!) f\'\'\'\'\'(i)`

`f(i+2) = f(i) + (2trianglex) f\'(i) + ((2trianglex)^2)/(2!) f\'\'(i) + ((2trianglex)^3)/(3!) f\'\'\'(i) + ((2trianglex)^4)/(4!) f\'\'\'\'(i) + ((2trianglex)^5)/(5!) f\'\'\'\'\'(i)`

`f(i+3) = f(i) + (3trianglex) f\'(i) + ((3trianglex)^2)/(2!) f\'\'(i) + ((3trianglex)^3)/(3!) f\'\'\'(i) + ((3trianglex)^4)/(4!) f\'\'\'\'(i) + ((3trianglex)^5)/(5!) f\'\'\'\'\'(i)`

`f(i+4) = f(i) + (4trianglex) f\'(i) + ((4trianglex)^2)/(2!) f\'\'(i) + ((4trianglex)^3)/(3!) f\'\'\'(i) + ((4trianglex)^4)/(4!) f\'\'\'\'(i) + ((4trianglex)^5)/(5!) f\'\'\'\'\'(i)`

`f(i+5) = f(i) + (5trianglex) f\'(i) + ((5trianglex)^2)/(2!) f\'\'(i) + ((5trianglex)^3)/(3!) f\'\'\'(i) + ((5trianglex)^4)/(4!) f\'\'\'\'(i) + ((5trianglex)^5)/(5!) f\'\'\'\'\'(i)`

Based on the order of approximation, we can determine the number of terms based on each sceheme by :

Central differencing scheme (CDS):

Number of terms = Order of approximation + Order of accuracy - 1

Skewed Differencing Scheme (SDS):

Number of terms = Order of approximation + Order of accuracy

• CDS

Assume `(del^2f)/(delx^2) = a.f(i-2) + b.f(i-1) + c.f(i) + d.f(i+1) + e.f(i+2)`

Table 1. Central Differencing Scheme

By solving the equations, the values of coefficients are:

a = -0.0833, b = 1.3333, c = -2.5000, d = 1.3333, e = -0.0833

• Skewed Right-Sided Differencing Scheme

Assume `(del^2f)/(delx^2) = a.f(i) + b.f(i+1) + c.f(i+2) + d.f(i+3) + e.f(i+4) + f.f(i+5)`

Table 2. Skewed Right-Sided Differencing Scheme

By solving the equations, the values of coefficients are:

a = 3.75, b = -12.8333, c = 17.8333, d = -13, e = 5.0833, f = -0.8333

• Skewed Left-Sided Differencing Scheme

Assume `(del^2f)/(delx^2) = a.f(i) + b.f(i-1) + c.f(i-2) + d.f(i-3) + e.f(i-4) + f.f(i-5)`

Table 3. Skewed Left-Sided Differencing Scheme

By solving the equations, the values of coefficients are:

a = 3.75, b = -12.8333, c = 17.8333, d = -13, e = 5.0833, f = -0.8333

Matlab code:

• Estimation of coefficients

clear all
close all
clc

% Central Differencing Scheme (CDS)
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];
Coefficient_CDS = inv(A) * B;
Coefficient_CDS = Coefficient_CDS\';

% Skewed Right - Sided Differencing Scheme
C = [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];
D = [0; 0; 1; 0; 0; 0];
Coefficient_Skewed_RDS = inv(C) * D;
Coefficient_Skewed_RDS = Coefficient_Skewed_RDS\';

% Skewed Left - Sided Differencing Scheme
E = [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];
F = [0; 0; 1; 0; 0; 0];
Coefficient_Skewed_LDS = inv(E) * F;
Coefficient_Skewed_LDS = Coefficient_Skewed_LDS\';

• Function for Central Differencing Scheme

% Function for central differencing scheme
function out = central_diff_approximation(x, dx)
% f(x) = exp(x)*cos(x)
% INPUT
% syms x
% a = diff(exp(x)*cos(x));
% b = diff(a)
% OUTPUT
% b = -2*exp(x)*sin(x)
% So, the analytical derivative, f\'\'(x) = -2*exp(x)*sin(x)

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

% Numerical derivative
% Fourth order central differencing approximation
a = -0.0833;
b = 1.3333;
c = -2.5000;
d = 1.3333;
e = -0.0833;

%central_differencing = (a*f(i-2) + b*f(i-1)+ c*f(i) +d*f(i+1) + e*f(i+2))/dx^2
central_differencing = ((a*exp(x-(2*dx))*cos(x-(2*dx))) + (b*exp(x-dx)*cos(x-dx))+ (c*exp(x)*cos(x)) + (d*exp(x+dx)*cos(x+dx)) + (e*exp(x+(2*dx))*cos(x+(2*dx))))/(dx^2);
%output
out = abs(central_differencing - analytical_derivative);
end

• Function for Skewed right sided differencing scheme

% Function for Skewed right sided differencing scheme

function out = skewed_rightsided_approximation(x, dx)
%f(x) = exp(x)*cos(x)
% INPUT
% syms x
% a = diff(exp(x)*cos(x));
% b = diff(a)
% OUTPUT
% b = -2*exp(x)*sin(x)
% So, the analytical derivative, f\'\'(x) = -2*exp(x)*sin(x)

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

% Numerical derivative
% Fourth order Skewed right sided differencing approximation
a = 3.75000;
b = -12.83333;
c = 17.83333;
d = -13.0000;
e = 5.08333;
f = -0.83333;

%skewed_rightsided_differencing = ( a*f(i) + b*f(i+1)+ c*f(i+2) +d*f(i+3) + e*f(i+4) + f*f(i+5))/dx^2
skewed_rightsided_differencing = ((a*exp(x)*cos(x)) + (b*exp(x+dx)*cos(x+dx))+ (c*exp(x+(2*dx))*cos(x+(2*dx))) + (d*exp(x+(3*dx))*cos(x+(3*dx))) + (e*exp(x+(4*dx))*cos(x+(4*dx))) + (f*exp(x+(5*dx))*cos(x+(5*dx))))/(dx^2);
%output
out = abs(skewed_rightsided_differencing - analytical_derivative);
end

• Function for Skewed left sided differencing scheme

% Function for Skewed left sided differencing scheme

function out = skewed_leftsided_approximation(x, dx)
%f(x) = exp(x)*cos(x)
% INPUT
% syms x
% a = diff(exp(x)*cos(x));
% b = diff(a)
% OUTPUT
% b = -2*exp(x)*sin(x)
% So, the analytical derivative, f\'\'(x) = -2*exp(x)*sin(x)

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

% Numerical derivative
% Fourth order Skewed left sided differencing approximation
a = 3.75000;
b = -12.83333;
c = 17.83333;
d = -13.0000;
e = 5.08333;
f = -0.83333;

%skewed_leftsided_differencing = ( a*f(i) + b*f(i-1)+ c*f(i-2) +d*f(i-3) + e*f(i-4) + f*f(i-5))/dx^2
skewed_leftsided_differencing = ((a*exp(x)*cos(x)) + (b*exp(x-dx)*cos(x-dx))+ (c*exp(x-(2*dx))*cos(x-(2*dx))) + (d*exp(x-(3*dx))*cos(x-(3*dx))) + (e*exp(x-(4*dx))*cos(x-(4*dx))) + (f*exp(x-(5*dx))*cos(x-(5*dx))))/(dx^2);
%output
out = abs(skewed_leftsided_differencing - analytical_derivative);
end

• Function for forward differencing scheme

% Function for Forward differencing scheme

function out = forward_approximation(x, dx)
%f(x) = exp(x)*cos(x)
% INPUT
% syms x
% a = diff(exp(x)*cos(x));
% b = diff(a)
% OUTPUT
% b = -2*exp(x)*sin(x)
% So, the analytical derivative, f\'\'(x) = -2*exp(x)*sin(x)

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

% Numerical derivative
% Forward_differencing = (f(x+dx) - f(x))/dx
forward_differencing = ((exp(x+dx)*cos(x+dx)) - (exp(x)*cos(x)))/dx;
%output
out = abs(forward_differencing - analytical_derivative);
end

• Function for backward differencing scheme

% Function for Backward differencing scheme

function out = backward_approximation(x, dx)
%f(x) = exp(x)*cos(x)
% INPUT
% syms x
% a = diff(exp(x)*cos(x));
% b = diff(a)
% OUTPUT
% b = -2*exp(x)*sin(x)
% So, the analytical derivative, f\'\'(x) = -2*exp(x)*sin(x)

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

% Numerical derivative
% Backward_differencing = (f(x-dx) - f(x))/dx
backward_differencing = ((exp(x-dx)*cos(x-dx)) - (exp(x)*cos(x)))/dx;
%output
out = abs(backward_differencing - analytical_derivative);
end

• Evaluation of error for different differencing schemes

clear all
close all
clc

% Values for x and dx
x = pi/4;
dx = linspace(pi/4,pi/2000,100);

for i = 1: length(dx)
    Central_error(i) = central_diff_approximation(x, dx(i)); % Error by central differencing
    Skewed_rightsided_error(i) = skewed_rightsided_approximation(x, dx(i)); % Error by skewed right sided differencing
    Skewed_leftsided_error(i) = skewed_leftsided_approximation(x,dx(i)); % Error by skewed left sided differencing
    FDS_error(i) = forward_approximation(x,dx(i)); % Error by forward differencing
    BDS_error(i) = backward_approximation(x,dx(i)); % Error by backward differencing
end

%plotting the dx values against errors for each type of schemes
figure(1)
loglog(dx, Central_error,\'m\',\'linewidth\',3)
hold on
loglog(dx, Skewed_rightsided_error, \'r\',\'linewidth\',3)
hold on
loglog(dx, Skewed_leftsided_error,\'g\',\'linewidth\',3)
hold on
loglog(dx, FDS_error, \'--\',\'linewidth\',3)
hold on
loglog(dx, BDS_error,\'--\',\'linewidth\', 3)
grid
xlabel(\'Grid spacing (dx)\')
ylabel(\'Absolute error\')
legend(\'Central error\', \'Right skewed error\', \'Left skewed error\', \'FDS error\', \'BDS error\')

Results:

Conclusions:

1. The fourth order differencing approaches are very much more accurate than the first order difference approaches like FDS and BDS.
2. While the grid spacing (dx) decreases, fourth order approximations get closer to the exact value in an accelerated rate. Although, the FDS and BDS shows no significant change in the errors. Hence, it denotes that higher order approximations are significantly more accurate than first order approximations.
3. For obtaining high accuracy in computed solutions, smaller values of dx is recommeded.
4. Comparison of the fourth order approximations indicate that the skewed difference scheme (right sided and left sided) gives the lower errors than the central difference scheme. With the decrease in grid spacing, the errors of skewed difference schemes are also reducing. But for the central difference scheme, after some point, the error value is a constant even after further reduction in dx value.
5. When we consider a scenario where no nodes are available on either sides of the node which we consider to discretize, we cannot use central differece scheme because it requires informations from nodes of either sides. In this type of scenarios, skewed difference schemes are useful because they require information only from any one of the sides of the node under consideration.