Deriving 4th order approximation of a 2nd order derivative by differencing schemes and evaluating them using MATLAB

Objective

> Derive the following 4th order approximations of the second-order derivative.

• Central difference
• Skewed right-sided difference
• Skewed left-sided difference

> To prove that the 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 derived.

> Generate a plot that compares the absolute error between the above-mentioned schemes.

> Explain why a skewed scheme is useful? What can a skewed scheme do that a CD scheme cannot?

4th Order approximations of the second-order derivative

Central difference

The number of nodes in the stencil can be calculated by the formula

N = P + Q - 1

Where,

N - Number of nodes

P - Order of approximations

Q - Order of derivative

Note: 1 is subtracted to derive the central scheme

To derive 4th order approximation of the second-order derivative, the number of nodes required for a central scheme is

N = 4 + 2 -1

N = 5

The hypothesis equation for central scheme with 5 nodes will be

`(del^2x)/(dely^2) = af(i-2)+bf(i-1)+cf(i)+df(i+1)+ef(i+2)`

Expanding all the terms in taylor series we get

`af(i-2) = af(i)-a(2Deltax)/(1!)f'(i)+a(2Deltax)^2/(2!)f''(i)-a(2Deltax)^3/(3!)f'''(i)+a(2Deltax)^4/(4!)f''''(i)-a(2Deltax)^5/(5!)f'''''(i)+a(2Deltax)^6/(6!)f''''''(i)+...`

`bf(i-1) = bf(i)-b(Deltax)/(1!)f'(i)+b(Deltax)^2/(2!)f''(i)-b(Deltax)^3/(3!)f'''(i)+b(Deltax)^4/(4!)f''''(i)-b(Deltax)^5/(5!)f'''''(i)+b(Deltax)^6/(6!)f''''''(i)-...`

`cf(i) = cf(i)`

`df(i+1) = df(i)+d(Deltax)/(1!)f'(i)+d(Deltax)^2/(2!)f''(i)+d(Deltax)^3/(3!)f'''(i)+d(Deltax)^4/(4!)f''''(i)+d(Deltax)^5/(5!)f'''''(i)+d(Deltax)^6/(6!)f''''''(i)+...`

`ef(i+2) = ef(i)+e(2Deltax)/(1!)f'(i)+e(2Deltax)^2/(2!)f''(i)+e(2Deltax)^3/(3!)f'''(i)+e(2Deltax)^4/(4!)f''''(i)+e(2Deltax)^5/(5!)f'''''(i)+e(2Deltax)^6/(6!)f''''''(i)+...`

Taylor table representation

`Terms->`	`f(i)`	`Deltaxf'(i)`	`Deltax^2f''(i)`	`Deltax^3f'''(i)`	`Deltax^4f''''(i)`	`Deltax^5f'''''(i)`
`af(i-2)`	`a`	`-2a`	`(4a)/2`	`-(8a)/6`	`(16a)/24`	`-(32a)/120`
`bf(i-1)`	`b`	`-b`	`b/2`	`-b/6`	`b/24`	`-b/120`
`cf(i)`	`c`	`0`	`0`	`0`	`0`	`0`
`df(i+1)`	`d`	`d`	`d/2`	`d/6`	`d/24`	`d/120`
`ef(i+2)`	`e`	`2e`	`(4e)/2`	`(8e)/6`	`(16e)/24`	`(32e)/120`
`Sum->`	`0`	`0`	`1`	`0`	`0`	`?`

Matrix representation

`[[1,1,1,1,1],[-2,-1,0,1,2],[4/2,1/2,0,1/2,4/2],[-8/6,-1/6,0,1/6,8/6],[16/24,1/24,0,1/24,16/24],[-32/120,-1/120,0,1/120,32/120]][[a],[b],[c],[d],[e]] = [[0],[0],[1],[0],[0]]`

MATLAB is used to solve the above matrix equation

The solved variable "X" has the values of a,b,c,d,e which are substituted in the hypothesis equation

`(del^2x)/(dely^2) = (-0.833f(i-2)+1.3333f(i-1)-2.5f(i)+1.3333f(i+1)-0.0833f(i+2))/(Deltax^2)`

Skewed Right-side difference

The number of nodes in the stencil can be calculated by the formula

N = P + Q

Where,

N - Number of nodes

P - Order of approximations

Q - Order of derivative

To derive 4th order approximation of the second-order derivative, the number of nodes required for a skewed right side scheme is

N = 4 + 2

N = 6

The hypothesis equation for skewed right side scheme with 6 nodes will be

`(del^2x)/(dely^2) = af(i)+bf(i+1)+cf(i+2)+df(i+3)+ef(i+4)+gf(i+5)`

Taylor table representation

`Terms->`	`f(i)`	`Deltaxf'(i)`	`Deltax^2f''(i)`	`Deltax^3f'''(i)`	`Deltax^4f''''(i)`	`Deltax^5f'''''(i)`	`Deltax^6f'''''"(i)`
`af(i)`	`a`	`0`	`0`	`0`	`0`	`0`	`0`
`bf(i+1)`	`b`	`b`	`b/2`	`b/6`	`b/24`	`b/120`	`b/720`
`cf(i+2)`	`c`	`2c`	`(4c)/2`	`(8c)/6`	`(16c)/24`	`(32c)/120`	`(64c)/720`
`df(i+3)`	`d`	`3d`	`(9d)/2`	`(27d)/6`	`(81d)/24`	`(243d)/120`	`(729d)/720`
`ef(i+4)`	`e`	`4e`	`(16e)/2`	`(64e)/6`	`(256e)/24`	`(1024e)/120`	`(4096e)/720`
`gf(i+5)`	`g`	`5g`	`(25g)/2`	`(125g)/6`	`(625g)/24`	`(3125g)/120`	`(15625g)/720`
`Sum->`	`0`	`0`	`1`	`0`	`0`	`0`	`?`

Matrix representation

`[[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]][[a],[b],[c],[d],[e],[g]] = [[0],[0],[1],[0],[0],[0]]`

MATLAB is used to solve the above matrix equation

The solved variable "X" has the values of a,b,c,d,e which are substituted in the hypothesis equation

`(del^2x)/(dely^2) = (3.75f(i)-12.8333f(i+1)+17.8333f(i+2)-13f(i+3)+5f(i+4)-0.8333f(i+5))/(Deltax^2)`

Skewed Left-side difference

The number of nodes in the stencil can be calculated by the formula

N = P + Q

Where,

N - Number of nodes

P - Order of approximations

Q - Order of derivative

To derive 4th order approximation of the second-order derivative, the number of nodes required for a skewed right side scheme is

N = 4 + 2

N = 6

The hypothesis equation for skewed right side scheme with 6 nodes will be

`(del^2x)/(dely^2) = af(i)+bf(i-1)+cf(i-2)+df(i-3)+ef(i-4)+gf(i-5)`

Taylor table representation

`Terms->`	`f(i)`	`Deltaxf'(i)`	`Deltax^2f''(i)`	`Deltax^3f'''(i)`	`Deltax^4f''''(i)`	`Deltax^5f'''''(i)`	`Deltax^6f'''''"(i)`
`af(i)`	`a`	`0`	`0`	`0`	`0`	`0`	`0`
`bf(i-1)`	`b`	`-b`	`b/2`	`-b/6`	`b/24`	`-b/120`	`b/720`
`cf(i-2)`	`c`	`-2c`	`(4c)/2`	`(-8c)/6`	`(16c)/24`	`(-32c)/120`	`(64c)/720`
`df(i-3)`	`d`	`-3d`	`(9d)/2`	`(-27d)/6`	`(81d)/24`	`(-243d)/120`	`(729d)/720`
`ef(i-4)`	`e`	`-4e`	`(16e)/2`	`(-64e)/6`	`(256e)/24`	`(-1024e)/120`	`(4096e)/720`
`gf(i-5)`	`g`	`-5g`	`(25g)/2`	`(-125g)/6`	`(625g)/24`	`(-3125g)/120`	`(15625g)/720`
`Sum->`	`0`	`0`	`1`	`0`	`0`	`0`	`?`

Matrix representation

`[[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]][[a],[b],[c],[d],[e],[g]] = [[0],[0],[1],[0],[0],[0]]`

MATLAB is used to solve the above matrix equation

The solved variable "X" has the values of a,b,c,d,e which are substituted in the hypothesis equation

`(del^2x)/(dely^2) = (3.75f(i)-12.8333f(i+1)+17.8333f(i+2)-13f(i+3)+5f(i+4)-0.8333f(i+5))/(Deltax^2)`

Matlab program to evaluate the second-order derivative of the analytical function exp(x)*cos(x) and to compare it with the 3 numerical approximations derived.

clear all
close all
clc

% Inputs
x = pi/3;
dx = linspace(pi/700,pi/4,100);

% Analytical_Function
% f(x) = exp(x)*cos(x)
% f''(x) = -2*exp(x)*sin(x)
a_d = -2*exp(x)*sin(x);

% Creating anonymous function
f =@(x) exp(x)*cos(x);

% Deriving the numerical schemes
% Central difference
A = [1,1,1,1,1;-2,-1,0,1,2;4/2,1/2,0,1/2,4/2;-8/6,-1/6,0,1/6,8/6;16/24,1/24,0,1/24,16/24];
B = [0;0;1;0;0];
Xc= linsolve(A,B);

% Skewed Right side difference
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];
Xr= linsolve(A,B);

% Skewed Left side difference
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];
Xl= linsolve(A,B);

% Numerical derivative
for i = 1:length(dx)
    % Central Difference     
    Nd_Cd(i) = (Xc(1)*f(x-2*dx(i))+Xc(2)*f(x-1*dx(i))+Xc(3)*f(x)+Xc(4)*f(x+1*dx(i))+Xc(5)*f(x+2*dx(i)))/(dx(i)^2);
    % Skewed right side difference
    Nd_Sr(i) = (Xr(1)*f(x)+Xr(2)*f(x+1*dx(i))+Xr(3)*f(x+2*dx(i))+Xr(4)*f(x+3*dx(i))+Xr(5)*f(x+4*dx(i))+Xr(6)*f(x+5*dx(i)))/(dx(i)^2);
    % Skewed left side difference
    Nd_Sl(i) = (Xl(1)*f(x)+Xl(2)*f(x-1*dx(i))+Xl(3)*f(x-2*dx(i))+Xl(4)*f(x-3*dx(i))+Xl(5)*f(x-4*dx(i))+Xl(6)*f(x-5*dx(i)))/(dx(i)^2);
    
    % Absolute error value calculation 
    % Error values of Central difference    
    E_Cd(i) = abs(Nd_Cd(i) - a_d);
    % Error values of skewed right side difference     
    E_Sr(i) = abs(Nd_Sr(i) - a_d);
    % Error values of skewed left side difference    
    E_Sl(i) = abs(Nd_Sl(i) - a_d);
end

% Plot
figure(1)
loglog(dx,E_Cd,'b','LineWidth',1.5)
hold on
loglog(dx,E_Sr,'k','LineWidth',1.5)
loglog(dx,E_Sl,'r','LineWidth',1.5)
title('Comparison of differencing schemes')
xlabel('Range of dx')
ylabel('Error values of different schemes')
legend('Central difference','Skewed right side difference','Skewed left side difference')

Program explanation

> At first, the inputs x and dx are given.

> Then the second derivative of the analytical function is stored to a variable which can be used for error calculation later in the program.

> An anonymous function for the function `exp(x)*cos(x)` is created.

> The numerical schemes "Central differencing, Skewed right side differencing, Skewed left side differencing" are derived using the pre-work done before programming.

> A loop for the length of dx is generated in which, the function is numerically solved using the 3 derived schemes.

> The anonymous function is called to solve the numerical derivatives which are then stored in their respective arrays.

> The absolute error values of the three schemes are also calculated in the same loop using the value of the analytical derivative.

> The log log plot is made between the error values of these three schemes and the range of dx where the solution is stable.

Result

> The range of dx (pi/4 - pi/700) is where the solution is stable comparatively.

> From the above plot, it can be seen that the central difference scheme is more stable and have the lowest error.

> The instability in the skewed left side difference scheme might be formed due to the truncation error or by the value of dx.

> Hence, the errors of the three schemes for the function `exp(x)*cos(x)` are compared.

Skewed schemes

> Although skewed schemes are less accurate in the above case they hold the advantage to calculate the approximation in any order.

> If boundary conditions are given then, skewed schemes come in handy compared to the CD scheme.

> For convection and diffusion equations, the CD scheme produces unstable solutions.

> Skewed schemes can be used for backwards differencing(left side) as well as forward differencing(right side) whereas, the central scheme can only be used for an internal solution.

> Skewed schemes require information from only one direction but, CD schemes require information from both sides.