Week 3.5 - Deriving 4th order approximation of a 2nd order derivative using Taylor Table method

Week 3.5

Deriving 4th order approximation of a 2nd order derivative using Taylor Table method

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

1. Central difference
2. Skewed right-sided difference
3. Skewed left-sided difference

After deriving these schemes, prove that your skewed schemes are fourth-order accurate i.e the thumb rule given in HINT (ii).

Once you have done this, 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 that you have derived.

Provide a plot that compares the absolute error between the above-mentioned schemes

In your own words, describe why a skewed scheme is useful? What can a skewed scheme do that a CD scheme cannot?

Aim : 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 that we have derived using

1. Central difference
2. Skewed right-sided difference
3. Skewed left-sided difference

Approximations of the second-order derivative:

CFD is the analysis of physical systems like fluid flow and heat transfer by performing computer-based simulations. We take up a system, find the underlying governing differential equations, and determine suitable numerical methods to approximate them. A numerical approximation is necessary because a problem in differential equations can rarely be solved analytically. Taylor table method is another numerical approximation method to compute the derivatives of the function. In the Taylor table method, we find numerical schemes to meet the target accuracy.

As order of approximation get higher we get higher accuracy. Hence we are going to derive 3 methods of approximation for given function.

For each type of difference method we have to find the No.of nodes in a stencil to form an equation.

from hint we get that two equations.

P=N-Q+1

P=N-Q

dx - represents the difference of two consecutive nodes

Here, function f(x) =

To calculate 4th order approximation of 2nd order derivative of function f(x);

Here, P = 4 and Q = 2

Taylor series general equation:

`f(x+ⅆx)=f(x)+(∂f(x))/∂x.dx+(∂^2 f(x))/〖∂x〗^2 .((ⅆ〖x)〗^2)/2!+⋯+(∂^n f(x))/〖∂x〗^n .((ⅆ〖x)〗^n)/n!+⋯⋅`

Governing Equations:

1. Central Difference Method :

For central schemes; no of nodes = P + Q -1 = 4 + 2 – 1 = 5

N = 5

Nodes considered: [ (i -2) , (i -1), (i), (i +1), (i+2) ]

Preparing Taylor table from above equations

Taylor Table for central difference method

Converting above Taylor table into matrix form

using Matlab Solving above equation we get,

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

2. Skewed right-sided difference:

No. of nodes = p+q = 4+2 = 6

p = order of derivative , q = order of approximation

Derivation: To derive fourth order approximation of second order derivative, we start with below approximation,

Preparing Taylor table from above equations

using Matlab Solving above equation we get,

a = 3.7500 , b = -12.8333, c = 17.8333, d =-13.0000,

e = 5.0833, g = -0.8333

3.Skewed Left-sided difference :

No. of nodes = p+q = 4+2 = 6

p = order of derivative

q = order of approximation

Derivation: To derive fourth order approximation of second order derivative, we start with below approximation,

using Matlab Solving above equation we get,

a = 3.7500 , b = -12.8333, c = 17.8333, d =-13.0000,

e = 5.0833, g = -0.8333

MATLAB CODE:

Programming Language: MATLAB

Function for central difference

function error = central_difference(analytical_derivative,x,dx)
%central differencing=
%a*f(x-2*dx)+b*f(x-dx)+c*f(x)+d*f(x+dx)+e*f(x+2*dx)
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]
X = inv(A)*B
f1=exp(x-2*dx)*cos(x-2*dx);
f2=exp(x-dx)*cos(x-dx);
f3=exp(x)*cos(x);
f4=exp(x+dx)*cos(x+dx);
f5=exp(x+2*dx)*cos(x+2*dx);
central_diff=((X(1)*f1)+(X(2)*f2)+(X(3)*f3)+(X(4)*f4)+(X(5)*f5))/(dx^2);
error=abs(analytical_derivative-central_diff);
end

Function for Skewed Right-sided difference

function error = skewed_right_difference(analytical_derivative,x,dx)
%skewed right sided difference=
%(a*f(x)+b*f(x+dx)+c*f(x+2*dx)+d*f(x+3*dx)+e*f(x+4*dx)+g*f(x+5*dx))/(dx^2)
C = [1 1 1 1 1 1;0 1 2 3 4 5;0 1/2 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]
Y = inv(C)*D
f1=exp(x)*cos(x);
f2=exp(x+dx)*cos(x+dx);
f3=exp(x+2*dx)*cos(x+2*dx);
f4=exp(x+3*dx)*cos(x+3*dx);
f5=exp(x+4*dx)*cos(x+4*dx);
f6=exp(x+5*dx)*cos(x+5*dx);
skewed_right_diff=((Y(1)*f1)+(Y(2)*f2)+(Y(3)*f3)+(Y(4)*f4)+(Y(5)*f5)+(Y(6)*f6))/(dx^2);
error=abs(analytical_derivative-skewed_right_diff);
end

Function for Skewed Left-sided difference

function error = skewed_left_difference(analytical_derivative,x,dx)
%skewed left sided difference=
%(a*f(x-5*dx)+b*f(x-4*dx)+c*f(x-3*dx)+d*f(x-2*dx)+e*f(x-dx)+g*f(x))/(dx^2)
E = [1 1 1 1 1 1;-5 -4 -3 -2 -1 0;25/2 16/2 9/2 4/2 1/2 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]
F = [0;0;1;0;0;0]
Z = inv(E)*F
f1=exp(x-5*dx)*cos(x-5*dx);
f2=exp(x-4*dx)*cos(x-4*dx);
f3=exp(x-3*dx)*cos(x-3*dx);
f4=exp(x-2*dx)*cos(x-2*dx);
f5=exp(x-dx)*cos(x-dx);
f6=exp(x)*cos(x);
skewed_left_diff=((Z(1)*f1)+(Z(2)*f2)+(Z(3)*f3)+(Z(4)*f4)+(Z(5)*f5)+(Z(6)*f6))/(dx^2);
error=abs(analytical_derivative-skewed_left_diff);
end

Main Program

Discretization_Main

clear all
close all
clc
x=pi/3;
f=exp(x)*cos(x);
dx=linspace(pi/4,pi/400,20);
analytical_derivative = -2*exp(x)*sin(x);
for i=1:length(dx)
    central_error(i) = central_difference(analytical_derivative,x,dx(i));
    skewed_right_error(i) = skewed_right_difference(analytical_derivative,x,dx(i));
    skewed_left_error(i) = skewed_left_difference(analytical_derivative,x,dx(i));
end
%plotting log graph
figure(1)
loglog(dx,central_error,'color','b','marker','square')
hold on
loglog(dx,skewed_right_error,'color','k','marker','O')
loglog(dx,skewed_left_error,'color','r','marker','diamond')
xlabel('dx(grid size)')
ylabel('Absolute error')
title('Consistency of different scheme')
legend('central difference','skewed right sided difference','skewed left sided difference')

Output :

• After analyzing the graph between ‘error’ vs. ‘dx’ for central, skewed right-sided, and skewed left-sided difference 4th order error; we can easily notice that central difference approximation produces more accurate solutions; hence produces lesser errors as compared to right-skewed and left-skewed difference schemes.

Why a skewed scheme is useful? What can a skewed scheme do that a CD scheme cannot do?

• Since we can see that the central differencing scheme gives us more immaculate solutions or lesser errors as compared to the skewed difference scheme. But it is not always viable to use the central differencing to evaluate solutions.
• Suppose, if we want to compute numerical approximations for a function at the boundary or corner node; we cannot use the central difference scheme at this particular situation because it requires data from both sides of the node or boundary. So, in this situation skewed difference method plays a major role; as in the skewed scheme data from one side can be used to compute numerical approximation at a point.
• Skewed schemes are used when we are unknown about the data from both the sides of a point. Because we can not use central scheme at that situation.

Conclusion:

• From the above plot, we can observe that the order of truncation error is almost same for both skewed right and left difference. for a given grid size(dx), truncation error of central difference is lesser compared to both skewed left and right sided scheme therefore central differencing scheme is comparatively better numerical approximation.
• Since, the central difference numerical scheme takes information from both right and left grid points. Central difference scheme cannot be calculated or may not provide appropriate solution in some flow region such as boundary, shock wave, Central difference fails in boundary regions because it requires information from points beneath the boundary which is out of the computational domain. In such cases, skewed scheme may be preferred.
• Central differencing scheme has very less error compared to other schemes. Hence, it is better to use central differencing scheme because it uses information from both the sides of the point of consideration to give the approximate and accurate results.
• As left and right sided skew methods uses the data only from one side of point of interest it gives less accurate results.