DERIVING 4TH ORDER APPROXIMATION FOR SECOND-ORDER DERIVATIVE USING TAYLOR TABLE METHOD

Aim: To Derive 4th order approximations for second-order derivative.

Theory:

the taylor series for function is written as,

`f(x+deltax)= f(x)+(del f)/(delx).(deltax)+(del^2f)/(del x^2).((deltax^2)/2) +`........(1)

the abvoe equation is valid only if :

• if the number of terms are infinite and the series converges.
• `deltax` tends to zero`

If u(i,j) is the x component of velocity at point u(i,j ); u(i+1,j) is x component of velocity at pount(i+1,j)

Then from equation 1:

u(i+1,j)=u(i,j)+`(delu)/(delx)`(i,j).`deltax`+`(delta^2u)/(deltax^2)`(i,j).`((delta x^2)/2)`+.........(2)

solving for `(delu)/(delx)(i,j)` we get;

`(del u)/(del x)(i,j) =(u(i+1 ,j)-u(i,j))/(del x ) -(del^2u)/(del x^2)(i,j).((del x)/(2 ))+`.................(3)

Then term in the RHS of equation 3 containing u(i+j) - u(i,j) is the finite-difference. The other term that is neglected is the truncation error.

Since in the majority of the CFD, a first-order approximation is not sufficient. WE turn our heads to central difference schemes.

`(del^2 u)/(del x)(x^2 ) =(u(i+j ,j) -2.u(i,j) +u(i-j ,j))/(delta x^2) +`errors ....................(4)

Taylor table is a method that comes up with schemes that meet the targeted accuracy.

The formula for several nodes is as follow:

• for the central differences scheme, the number of nodes = p+q-1, where p is the order of derivative and q is the order of approximation.
• for the skewed schemes, be it right-sided or left-sided, the number of nodes is p+q, where p is the order of derivative q is the order of approximation.

Since the order is second and the order of approximation is 4, the number of nodes in the central difference scheme is 5, and the number of nodes is skewed right-sided, and skewed left-sided schemes are 6.

CENTRAL DIFFERENCING SCHEME:

We take up a hypothetical equation like:

`((del^2u)/(del x^2))x^2 = a.f(i-2) +b.f(i-1 ) + c.f(i) + d.f(i+1) + e.f(i+2)` ................(5)

Here the a,b,c,d,e are coefficients, The values of these coefficients are need to be found out.

The mode of Operandi will be like this:

• Firstly we expand every term of the above equation according to the Taylor series.
• Then we write the value obtained in a table format. The columns are the stencils.
• Then we equations with unknown to find out the values.

From Equation 5:

a.f(i - 2) = a.f(i) +a.`(delu)/(delx).(-2.(delta x)/1)` + a.`((del^2u)/(delx^2)).((-2.deltax)^2/4)` + a.`(del^3u)/(del x^3).((-2.deltax)^3/6) +a.(del^4u)/(del x^4).((-2.deltax^4)/24)`.........(6)

b.f(i-1) = b.f(i) + `b.(delu)/(delx).((-deltax)/1) + b.(del^2u)/(delx^2).((-deltax)^2/2) + b.(del^3u)/(delx^3).(((-deltax)^3)/6)`+`b.(del^4u)/(delx^4).((-deltax)^4/24)`..............(7)

c.f(i) = c.f(i).......................(8)

d. f(i+1) = d.f(i) +`d.(delu)/(delx).deltax + d.(del^2x)/(del x^2).((deltax)^2/2)+ d.(del^3u)/(delx^3).(((deltax)^3)/6) + d.(del^4u)/(delx^4).(((deltax)^4)/24)`..............(9)

e.f(i+2) = e.f(i) + `e.(delu)/(delx).(2.deltax) + e.(del^2u)/(delx^2).((2.deltax)^2/2) + e.(del^3u)/(delx^3).((2.deltax)^3/6 )+ e.(del^4u)/(delx^4).((2.deltax)^4/24)`..................(10)

Rearranging, we get:

`(del^2u)/(delx^2)` = f(i)(a+b+c+d+e) + `(delu)/(delx).deltax`.(-2.a + (-b) + d +2.e) +`(del^2u)/(delx^2).(deltax)^2`.(2a`+b/2 + d/2 +2e`) + `(del^3u)/(delx^3).(deltax^3).(((-8a)/6)+(-b/6)+d/6+(8e)/6)` + `(del^4u)/(delx^4).(deltax^4).((16a)/24+b/24+d/24+(16e)/24)`.......................(11)

Now writting equation 11 in tabular form, we get the Taylor Table as follow:

f(i) `(delu)/(delx).deltax` `(del^2u)/(delx^2).(deltax^2)` `(del^3u)/(delx^3).(deltax^3)` `(del^4u)/(delx^4).(deltax^4)`

a.f(i-2) a -2a 2a -`(8a)/6` `(16a)/24`

b.f(i-1) b -b `b/2` `-(b/6)` `b/24`

c.f(i) c 0 0 0 0

d.f(i+1) d d d/2` `b/6` `d/24`

e.f(i+2) e 2e 2e `(8e)/6` `(16e)/24`

Aprroximating and adding up the coefficient with the logic of the terms f(i),`(delu)/(delx), (del^3u)/(delx^3),(del^4u)/(delx^4)`will be approximated to 0 and the coefficients of `(del^2u)/(delx^2)`will be approximated to 1 , we get the following set of equation:

a+b+c+d+e =0

-2a-b+d+2e=0

2a+b/2+d/2+2e=1

-(8/6)*a-b/6+d/6+(8e/6)=0

(16a/24)+b/24+d/24+(16e/24)=0

Here we have 5 eqautions and 5 unknown

To solve this we we will write this in matrix form:

`[[1,1,1,1,1], [-2,-1,0,1,2],[2,1/2,0,1/2,2],[-8/6,-1/6,0,1/6,8/6],[16/24,1/24,0,1/24,16/24]]*[(a),(b),(c),(d)]=[(0),(0),(1),(0),(0)]`

Solving thr above matrix we get a=(-1/12), b=(4/3), c=(-5/2), d=(4/3), e=(-1/12)

So by inserting the above values of the coefficients in the main equation:

`(del^2u)/(del x^2) =(-1/12).f(i-2) +(4/3).f(i-1)+(-5/2).f(i)+(4/3).f(i+1)+(-1/12).f.(i+2)`..................(**)

Nownfor skewed right sided and skewed left sided schemes the same above way is followed to create the Taylor table. Herenwe directly input the table by not repeating the procedure:

SKEWED RIGHT SIDED DIFFERENCE:

The hypothetical equation is:

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

Expanding each of the above terms in the order of taylor series ; we obtain the Taylor table as :

f(i) `(delu)/(delx).deltax` `(del^2u)/(delx^2).(deltax^2)` `(del^3u)/(delx^3).(deltax^3)` `(del^4u)/(delx^4).(deltax^4)` `(del^5u)/(delx^5).(deltax^5)`

a.f(i) a 0 0 0 0 0

b.f(i+1) b b b/2 b/6 b/24 b/120

c.f(i+2) c 2c 2c 8/6*c 16/24*c 32/12*c

d.f(i+3) d 3d 9/2*d 27/6*d 81/24*d 243/120*d

e.f(i+4) e 4e 8*e 64/6*e 256/24*e 1024/120*e

g.f(i+5) g 5g 25/2*g 125/6*g 625/24*g 3125/120*g

Solving the above equation by applying the same logic that the term containing second order derivative is approximated to 1 and all others are 0, we get the coefficients as:

a=`15/4, b=(-77/6), c=(107/6), d=(-13), e=(61/12), g=(-5/6)`

Inserting the values in the parent equation we get:

`(del^2u)/(delx^2) = (15/4).f(i) +(-77/6).f(i+1) + (107/6).f(i+2) +(-13).f(i+3) + (61/12).f(i+4) + (-5/6).f(i+5)`

SKEWED LEFT-SIDED DIFFERENCING:

In this scheme we will consider the discretization points at the left side. The hypothetical equation for this scheme is written as :

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

By expanding each and every terms of the RHS according to the taylor series we get:

f(i) `(delu)/(delx).deltax` `(del^2u)/(delx^2).(deltax^2)` `(del^3u)/(delx^3).(deltax^3)` `(del^4u)/(delx^4).(deltax^4)` `(del^5u)/(delx^5).(deltax^5)`

a.f(i) a 0 0 0 0 0

b.f(i-1) b -b `b/2` `-b/6` `-b/24` `-b/120`

c.f(i-2) c -2c 2c `-8/6c` `16/24c` `-32/120c`

d.f(i-3) d -3d `9/2d` `-27/6d` `81/24d` `-243/120d`

e.f(i-4) e -4e 8e `-64/6e` `256/24e` `-1024/120e`

g.f(i-5) g -5g `25/2g` `-125/6g` `325/24g` `-3125/120g`

In matrix form the equations stand as:

`[[1,1,1,1,1,1],[0,-1,-2,-3,-4,-5],[0,1/2,2,9/2,8,25/2],[0,-1/6,-8/6,-27/6,-24/6,-125/6],[0,1/24,16/24,81/24,256/24,625/24],[0,-1/12,-32/120,-243/120,-1024/120,-3125/120]]*[(a),(b),(c),(d),(e),(g)]` =`[(0),(0),(1),(0),(0),(0)]`

Solving it by matrix inversion method in matlab, we get the coefficients as:

a= 3.75 b= -12.8333 c=17.8333 d=13.000 e=5.083 g=-0.8333

Inserting the values of coefficients:

`(del^2u)/(delx^2)=`3.75.f(i) + (-12.833).f(i-1) + 17.8333.f(i-2) + 13.00.f(i-3) + 5.083.f(i-4) + (-0.8333).f(i-5)

CODE EXPLANATION:

We aim to compare the errors calculated through the three different schemes.

• The function code is written for the three different schemes, the use of a function code is that we can call that value again and again in the main body. Care is taken for the syntax of a function. The coefficient values are provided that are obtained through the matrix inversion method by solving the linear equation above. The analytical derivative value is calculated. Since it is of second order, the second derivative is calculated. The error is calculated by using the abs command and subtracting the analytical derivative function from the three different schemes.
• In the main body, there are basically three parts; the First one is assigning the value of X, the range of dx.
• The linspace command is used to provide input for the dx as we intended to create an array that will store the dx values between two limits. dx is basically the grid size or the mesh size.
• The second one is creating a loop.
• The main idea behind the loop is that we want to calculate the errors for the three different schemes throughout the range of dx values. i is the loop counter variable, it keeps on updating the values of dx.
• The third one is creating a plot. A log log plot is created, proper labels are given using the x label and y label.

Function Codes:

function left_error = left_skewed_difference(x,dx)
  
  % analytical deivative
  analytical_derivative = -2*exp(x)*sin(x);
  
  a= (15/4);
  b= (-77/6);
  c= (107/6);
  d= (-13);
  e= (61/12);
  g= (-5/6);
  
  left_scheme= (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)+g*exp(x-5*dx)*cos(x-5*dx))/(dx^2);
  left_error= abs(left_scheme-analytical_derivative);
  end

function right_error = right_skewed_difference(x,dx)
  
  % analytical deivative
  analytical_derivative = -2*exp(x)*sin(x);
  
  a= (15/4);
  b= (-77/6);
  c= (107/6);
  d= (-13);
  e= (61/12);
  g= (-5/6);
  
  right_scheme= (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)+g*exp(x+5*dx)*cos(x+5*dx))/(dx^2);
  right_error= abs(right_scheme-analytical_derivative);
  end

function central_error = central_skewed_difference(x,dx)
  
  % analytical deivative
  analytical_derivative = -2*exp(x)*sin(x);
  
  a= (-1/12);
  b= (4/3);
  c= (-5/2);
  d= (4/3);
  e= (-1/12);

  
  central_difference_scheme= (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);
  central_error= abs(central_difference_scheme-analytical_derivative);
  end

MATLAB CODE:

clear all
close all
clc

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

x= pi/3;
dx= linspace(pi/400,pi/4,50);

for i= 1:length(dx)
  
     central_difference_scheme_error(i) = central_skewed_difference(x,dx(i));
     right_skewed_scheme_error(i) = right_skewed_difference(x,dx(i));
     left_skewed_scheme_error(i) = left_skewed_difference(x,dx(i));
end

%plotting the error graphs

figure(1)
loglog(dx,right_skewed_scheme_error,'color','b','marker','diamond')
hold on
loglog(dx,left_skewed_scheme_error, 'color','g','marker','*')
hold on 
loglog(dx,central_difference_scheme_error, 'color','r','marker','.')
grid on
xlabel('grid size')
ylabel('error')
legend('skewed right sided ','skewed left sided', 'central difference schmene')

OUTPUT:

CONCLUSION:

For the same value of dx, the central difference scheme gives a minimum error than the other two schemes.

• For the right-sided skewed scheme, we take nodes right of the point. For the left-sided skewed scheme, we take nodes left of the point.
• In CFD, the first-order accuracy is not sufficient. We require second-order accuracy.
• The main demerit of higher-order accuracy is that more grid points are required, so more computing time.
• Higher-order accuracy gives a more qualitative solution.
• The skewed scheme schemes find their importance in the discretization in the boundaries. If we tend to apply a central difference scheme at the boundary, we will need to have information about the parameter at a point below the boundary that we don't have. The reflection method was used earlier which is of little help. Skewed schemes find their importance here.