Taylor table method and matlab code for a fourth order approximation of second order derivative

Aim: A)To derive the 4th order approximations of the second order derivative

1. Central difference
2. Skewed right sided difference
3. skewed left handed difference

B) To write a Matlab Program to evaluate the second order derivative of analytical function exp(x)*cos(x) and compare it with numerical approximations, and also to compare the absolute errors between the second order schemes.

C) second order schemes comparison with FDS and BDS.

1. Fourth order Central difference of the second derivative:

frac{∂^2 f}{∂x^2 } = 4 th order accurate approximation using i -2, i -1, i, i +1, i +2

Here, the stencil is symmetric, it is a symmetric 2nd order scheme.

Using Taylor series method,

(∂^2 f)/(∂x^2 ) = af(i-2) + bf(i-1) +cf(i) + df(i+1)+ ef(i+2) ......(1)

Writing the tabular form, we have

 f(i) Δx fI(i) Δx2 fII(i) Δx3 fIII(i) Δx4 fIV(i) Δx5 fV(i) Δx6 fVI(i) af(i-2) a -2a (4a)/2 -(8a)/6 (16a)/24 -(32a)/120 (64a)/720 bf(i-1) b -b b/2 -b/6 b/24 -b/120 b/720 cfi c 0 0 0 0 0 0 df(i+1) d d d/2 d/6 d/24 d/120 d/720 ef(i+2) e 2e (4e)/2 (8e)/6 (16e)/24 (32e)/120 (64e)/720 0 0 1 0 0 ? ?

On further simplification, we get

[[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]][(a),(b),(c),(d),(e),(f)] = [(0),(0),(1),(0),(0)]

This is in the form of ''AX = B'' . Rewriting them in equation form,

a + b + c + d + e =0           ...........(2)

-2a - b + 0 + d + 2e = 0       ...........(3)

2a + b/2 + 0 + d/2 + 2e = 1  ...........(4)

-4/3 - b/6 + 0 + d/6 + 4/3 =0 ...........(5)

2a/3 + b/24 + 0 + d/24 + 2a/3 = 0 .....(6)

we have 5 equations and 5 unknowns, finally we get

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

substituting a,b,c,d,e in equation (i) so we get

(∂^2 f)/(∂x^2 ) approx (-0.0833f(i-2) + 1.3333f(i-1) - 2.5f(i) + 1.3333f(i+1)+ 0.08333f(i+2))/(Δx^2)

2. Fourth order skewed right sided difference of the second derivative:

frac{∂^2 f}{∂x^2 } = 4 th order accurate approximation using i , i +1, i, i +2, i +3,

i +4, i+5

Here, the stencil is not symmetric, but it is a skewed right sided difference of 2nd order scheme.

Using Taylor series method,

(∂^2 f)/(∂x^2 ) = af(i) + bf(i+1) +cf(i +2)+ df(i+3)+ ef(i+4)+gf(i+5)

Writing the tabular form, we have

 f(i) Δx fI(i) Δx2 fII(i) Δx3 fIII(i) Δx4 fIV(i) Δx5 fV(i) Δx6 fVI(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 2c (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 (32e)/3 (256e)/24 (1024e)/120 (4096e)/720 gf(i+5) g 5g (25g)/2 (125g)/6 (625g)/24 (3125g)/120 (15625g)/720 0 0 1 0 0 ? ?

On further simplification, we get

[[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,32/120,243/120,1024/120,3125/120]][(a),(b),(c),(d),(e),(g)] = [(0),(0),(1),(0),(0),(0)]

This is in the form of ''AX = B'' . Rewriting them in equation form,

a + b + c + d + e +g=0           ...........(7)

0 + b + 2c + 3d + 4e + 5g= 0       ...........(8)

0 + b/2 + 2c + 9d/2 + 8e + 25g/2 = 1  ...........(9)

0 + b/6 + 4c/3 + 9d/2 + 32e/3 + 125g/6 =0 ...........(10)

0 + b/24 + 2c/3 + 27d/8 + 32e/3 + 625g/24 = 0 .....(11)

0 + b/120 + 32c/120+ 243d/120 + 1024e/120 + 3125g/120= 0.......(12)

we have 6 equations and 6 unknowns, finally we get

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

substituting a,b,c,d,e in equation (i) so we get

(∂^2 f)/(∂x^2 ) approx (1.3396f(i) - 2.9892f(i+1) + 1.9301f(i+2) - 0.2508f(i+3) - 0.0295f(i+4))/(Δx^2)

3. Fourth order skewed left sided difference of the second derivative:

frac{∂^2 f}{∂x^2} = 4 th order accurate approximation using i, i -1, i-2, i -3, i-4 and

i-5

Here, the stencil is not symmetric, but it is a skewed left sided difference of 2nd order scheme.

Using Taylor series method,

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

Writing the tabular form, we have

 f(i) Δx fI(i) Δx2 fII(i) Δx3 fIII(i) Δx4 fIV(i) Δx5 fV(i) Δx6 fVI(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 2c (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 (32e)/3 (256e)/24 (1024e)/120 (4096e)/720 gf(i+5) g 5g (25g)/2 (125g)/6 (625g)/24 (3125g)/120 (15625g)/720 0 0 1 0 0 ? ?

On further simplification, we get

[[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,-32/120,-243/120,-1024/120,-3125/120]][(a),(b),(c),(d),(e),(g)] = [(0),(0),(1),(0),(0),(0)]

This is in the form of ''AX = B'' . Rewriting them in equation form,

a + b + c + d + e + g=0            ...........(13)

0- b - 2c - 3d - 4e - 5g = 0       ..........(14)

0+ b/2 + 2c + {9d}/2 + 8e + {25g}/2= 1  ...........(15)

0 - b/6 - {4c}/3 - {9d}/2 - {32e}/3 - {125g}/6  =0 ...........(16)

0+ b/24 + {2c}/3 + {27d}/8 + {32e}/3 + {625g}/24 = 0 .....(17)

0 - b/120 - {32c}/120- {243d}/120 - {1024e}/120 - {3125g}/120 = 0.....(18))

we have 5 equations and 5 unknowns, finally we get

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

Calling a  function of fourth order second derivative Central differencing :

function out =fourth_order_second_derivative_central_differencing(x,dx)

% analytical function = exp(x)*cos(x)
% Exact second derivative = -2*exp(x)*sin(x)

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

% fourth order second derivative central differencing = (a*f(i-2) + b*f(i-1) + c*f(i) + d*f(i+1) + e*f(i+2))/(dx^2);
% The constants calculated using 5 equations are:  a = -0.0833, b = 1.3333, c = -2.5, d = 1.3333, e = -0.0833
fourth_order_second_derivative_central_differencing = ((-0.0833*exp(x-(2*dx))*cos(x-(2*dx))) + (1.3333*exp(x-dx)*cos(x-dx)) + (-2.5*exp(x)*cos(x)) + (1.3333*exp(x+dx)*cos(x+dx)) + (-0.0833*exp(x+(2*dx))*cos(x+(2*dx))))/(dx^2)

% fourth order second derivation approximation(central differencing)
out = abs(fourth_order_second_derivative_central_differencing-exact_second_derivative)

end

Calling a function of fourth order second derivative Right side Skeweddifferencing:

function out =fourth_order_second_derivative_skewed_right_side_differencing(x,dx)

% analytical function = exp(x)*cos(x)
% Exact second derivative = -2*exp(x)*sin(x)

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

% numerical derivative
% skewed_right_sided_differencing = (a*f(i) + b*f(i+1) + c*f(i+2) + d*f(i+3) + e*f(i+4) + g*f(i+5)/(dx^2);
% The constants calculated using 5 equations are:  a = 3.75, b = -12.8333, c = 17.8333, d = -13.0000, e = 5.0833, g = -0.8333
fourth_order_second_derivative_skewed_right_side_differencing = ((3.75*exp(x)*cos(x)) - (12.8333*exp(x+dx)*cos(x+dx)) + (17.8333*exp(x+(2*dx))*cos(x+(2*dx))) - (13*exp(x+(3*dx))*cos(x+(3*dx))) + (5.0833*exp(x+(4*dx))*cos(x+(4*dx))) - (0.83333*exp(x+(5*dx))*cos(x+(5*dx))))/(dx^2)

% fourth order second derivation approximation(right sided differencing)
out = abs(fourth_order_second_derivative_skewed_right_side_differencing-exact_second_derivative)

end

Calling a function of fourth order second derivative Left side Skewed differencing:

function out =fourth_order_second_derivative_skewed_left_side_differencing(x,dx)

% analytical function = exp(x)*cos(x)
% Exact second derivative = -2*exp(x)*sin(x)

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

% numerical derivative
% skewed_left_sided_differencing = (a*f(i) + b*f(i-1) + c*f(i-2) + d*f(i-3) + e*f(i-4) + gf(i-5)/(dx^2);
% The constants calculated using 5 equations are:  a = 3.75, b = -12.8333, c = 17.8333, d = -13, e = 5.0833, g = -0.8333
fourth_order_second_derivative_skewed_left_side_differencing = ((3.75*exp(x)*cos(x))-(12.8333*exp(x-dx)*cos(x-dx))+(17.8333*exp(x-(2*dx))*cos(x-(2*dx)))-(13*exp(x-(3*dx))*cos(x-(3*dx)))+(5.0833*exp(x-(4*dx))*cos(x-(4*dx)))-(0.8333*exp(x-(5*dx))*cos(x-(5*dx))))/(dx^2)

% fourth order second derivation approximation(left sided differencing)
out = abs(fourth_order_second_derivative_skewed_left_side_differencing-exact_second_derivative)

end

Calling a function of forward Differencing:

function out = forward_differencing(x,dx)

% analytical function = exp(x)*cos(x)
% Exact second derivative = -2*exp(x)*sin(x)

exact_second_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

% First order approximation
out = abs(forward_differencing - exact_second_derivative)

end

Calling a function of Backward Differencing:

function out = backward_differencing(x,dx)

% analytical function = exp(x)*cos(x)
% Exact second derivative = -2*exp(x)*sin(x)

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

% numerical derivative
% forward differencing = (f(x+dx)-f(x))/dx

backward_differencing = ((exp(x)*cos(x)) - (exp(x-dx)*cos(x-dx)))/dx

% First order approximation
out = abs(backward_differencing - exact_second_derivative)

end

The Main Matlab Program:

clear all
close all
clc

% analytical function = exp(x)*cos(x)
% f"(x) = -2*exp(x)*sin(x)

% Assuming the input value of x
x = pi/3;
% Defining the grid size,dx
dx = linspace(pi/10,pi/100,10)

% second_order_derivative = -2*exp(x)*sin(x);
for i = 1:length(dx)

fourth_order_second_derivative_central_differencing_error(i) = fourth_order_second_derivative_central_differencing(x,dx(i))
fourth_order_second_derivative_skewed_right_side_differencing_error(i) = fourth_order_second_derivative_skewed_right_side_differencing(x,dx(i))
fourth_order_second_derivative_skewed_left_side_differencing_error(i) = fourth_order_second_derivative_skewed_left_side_differencing(x,dx(i))
forward_differencing_error(i) = forward_differencing(x,dx(i))
backward_differencing_error(i) = backward_differencing(x,dx(i))
end

figure(1)
loglog(dx,fourth_order_second_derivative_central_differencing_error)
hold on
loglog(dx, fourth_order_second_derivative_skewed_right_side_differencing_error,'r')
hold on
loglog(dx, fourth_order_second_derivative_skewed_left_side_differencing_error,'g')
xlabel('Grid size (dx)')
ylabel('Error')
legend('central differencing','Skewed Right sided differencing','skewed Left sided differencing')
title('Comparison of errors for different fourth order second derivative schemes')

figure(2)
loglog(dx,fourth_order_second_derivative_central_differencing_error)
hold on
loglog(dx, fourth_order_second_derivative_skewed_right_side_differencing_error,'r')
hold on
loglog(dx, fourth_order_second_derivative_skewed_left_side_differencing_error,'g')
hold on
loglog(dx, forward_differencing_error,'y')
hold on
loglog(dx, backward_differencing_error,'black')
xlabel('Grid size (dx)')
ylabel('Error')
legend('central differencing','Skewed Right sided differencing','skewed Left sided differencing','forward differencing','backward differencing')
title('Comparison of errors for different fourth order second derivative schemes VS Forward Differencing VS Backward Differencing')`

Plot 1:

Plot 2:

The absolute error for skewed schemes and cenral scheme with many neighbouring points results is  approximately a curve with respect to grid size dx. This is due to reason that the influence of many neighbouring points and a very small value of dx, generates large no.of airthmetic operations, which results in round off error. where for the same range of values dx, the error obtained by FDS and BDS are not curved, due to less no. of airthmetic operations.

Central differencing scheme error is less, but skewed schemes are mostly used because it can compute a single or higher order derivative at any corner node, in this case is no more  central scheme is useful. If scheme is not symmetric in those cases only skewed schemes are helpful, which means no. of node points on right and left are not same .

Parametric Study of Gate Valve Shravankumar Nagapuri · 2019-10-22 03:14:10

Objectives: To perform a parametric study on the gate valve simulation. To obtain the mass flow rates at the outlet for 5 design points. To show the cut section view for different design points, and also to show the gate disc lift and fluid volume extraction. To s Read more

Symmetry vs Wedge vs HP equation Shravankumar Nagapuri · 2019-10-09 22:50:52

Objectives: To Write a Matlab program that can generate the computational blockMesh file  automatically for a wedge angle of 3 degrees and compare the results obtained for Symmetry BC and Wedge BC. To Write a Matlab program that takes an angle as input and gene Read more

Combustion simulation on the combustor model Shravankumar Nagapuri · 2019-10-05 16:23:22

1Q. Briefly explain about the possible types of combustion simulations in FLUENT. A:  Combustion or burning is a high temperature exothermic chemical reaction between a fuel and an oxidant accompanied by the production of heat, light and unburnt gases in the fo Read more

Cyclone separator challenge Shravankumar Nagapuri · 2019-10-01 17:50:02

Cyclone Separator: Principle and Working: A high speed rotating (air)flow is established within a cylindrical or conical container called a cyclone. Air flows in a spiral pattern, beginning at the top (wide end) of the cyclone and ending at the bottom (narrow) end bef Read more

Gearbox sloshing effect Shravankumar Nagapuri · 2019-09-28 15:37:58

Objective: To analyse the flow pattern of the fluid inside the gearbox, for two different clearances of the same geometry and each geometry flow is analyzed by two different fluids (Water and Oil) as lubricants. Gearbox sloshing effect:  Slosh refers Read more

Simulation of Flow through a pipe in OpenFoam Shravankumar Nagapuri · 2019-09-26 19:10:15

Objectives: To write a program in Matlab that can generate the computational mesh automatically for any wedge angle and Grading scheme. To calculate, length of the pipe using the entry length formula for laminar flow through a pipe To show that entry length is suf Read more

Conjugate Heat Transfer CHT Analysis on a graphics card Shravankumar Nagapuri · 2019-09-20 16:08:56

Objectives: To find out the maximum temperature attained by the processor. To Prove that the simulation has achieved convergence with appropriate images and plots. To Find out the heat transfer coefficient at appropriate areas of the model. To Identify potential h Read more

Finite Volume Method FVM Literature Review Shravankumar Nagapuri · 2019-09-17 13:14:28

Finite volume Method: The finite volume method (FVM) is a method for representing and evaluating partial differential equations in the form of algebraic equations. Similar to the finite difference method (FDM) or finite element method (FEM), values are calculated at Read more

Rayleigh Taylor instability Shravankumar Nagapuri · 2019-09-14 17:48:30

Q1. What are some practical CFD models that have been based on the mathematical analysis of Rayleigh Taylor waves? In your own words, explain how these mathematical models have been adapted for CFD calculations  Kelvin–Helmholtz instability: The Kelvin Read more

BlockMesh Drill down challenge Shravankumar Nagapuri · 2019-09-11 18:17:33

Objectives: To Generate the BlockMesh file for the given Geometry using OpenFOAM and use the icoFoam solver to simulate the flow through a backward-facing step. To create multiple meshes and will be comparing the results obtained from each mesh. To show the velocit Read more