BlogFor BusinessTalk to a counsellor
CoursesIIT ProgramsPlacementsLog in
BlogFor Business
Talk to a counsellor
Challenge

Week 3 - Taylor table method and Matlab code

SShivshankar Nadarajan6 Aug 2020

Aim:

To derive the following 4th order approximations of the second-order derivative using Taylor table

1. Central difference

2. Skewed right-sided difference

3. Skewed left-sided difference

and prove that your skewed schemes are fourth-order accurate.

Theory:

The discretization of the geometric domain is a must for any type of numerical solution. In finite difference discretization method, the grid is locally structured according to the local coordinate system.

In this method, nodes are defined on the basis of set of indices (i,j) and according to the method of differencing neighbouring nodes are either increasing or decreasing.

For calculating the neighbouring nodes either in increasing or decreasing trend Taylor series is used.

Taylors series for any function will be : `f(x+Deltax)=f(x) + (Deltax)/(1!)*f'(x) + (Deltax)^2/(2!)*f''(x) + (Deltax)^3/(3!)*f'''(x) +.....`

where `Deltax=` distance between two stencils.

The number of required nodes to solve the problem depends on P (order of derivative) and Q (order of approximation)

The method of difference is based on their use of neighbouring nodes

Right Skewed differencing:- only the forward nodes are used. If i is the target node than i+1, i+2, i+3....will be used for right-skewed differencing (forward differencing).

Left Skewed differencing:- only the backwards nodes are used. If i is the target node than i-1, i-2, i-3...will be used for left-skewed differencing (backward differencing).

Central differencing:- Both the forward and backward node is used. If i is the target node than .....i-1, i, i+1,... will be used for central differencing.

Derivation:

CENTRAL DIFFERENCE

  • Here the number of nodes will be P+Q-1 i.e for 4th order approximation of the second order derivative it will be 4+2-1=5, thus 5 nodes will be required to derive the solution.
  • The 5 nodes for central difference scheme will be: i-2, i-1, i, i+1, i+2.
  • `Deltax^2.(d^2u)/dx^2= a.f(i-2)+b.f(i-1)+c.f(i)+d.f(i+1)+e.f(i+2)`
  • Taylor table for the following is given as,

f(i)

`Deltax*f(i)`

`Deltax^2*f(i)`

`Deltax^3*f(i)`

`Deltax^4*f(i)`

a.f(i-2)

a

-2a

`(4a)/2`

`(-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`

`(d)/6`

`(d)/24`

e.f(i+2)

e

2e

`(4e)/2`

`(8e)/6`

`(16e)/24`

`Deltax^2.f''(x)`

0

0

1

0

0

Thus there are 5 unknown and 5 equations

`a+b+c+d+e=0`

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

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

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

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

Converting them into matrix form,

`[[1,1,1,1,1],[-2,-1,0,1,2], [4,1,0,1,4],[-8,-1,0,1,8],[16,1,0,1,16]]*[[a],[b],[c],[d],[e]]=[[0],[0],[1],[0],[0]]`

close all
clear all
clc
 syms a b c d e 
 eqn1= a+b+c+d+e==0;
 eqn2= -2*a-b+d+2*e==0;
 eqn3= (4*a)/2+b/2+d/2+(4*e)/2==1;
 eqn4= -(8*a)/6-(b)/6+(d)/6+(8*e)/6==0;
 eqn5= (16*a)/24+(b)/24+(d)/24+(16*e)/24==0;
 
 [A,B] = equationsToMatrix([eqn1, eqn2, eqn3, eqn4, eqn5], [a, b, c, d, e])
 
 X = linsolve(A,B)

Output

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
 
 
X =
 
 -1/12
   4/3
  -5/2
   4/3
 -1/12

Thus the new equation can be written by substituing the values of unknown in equation,

  • `Deltax^2.(d^2u)/dx^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)`
  • `(d^2u)/dx^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))/((Deltax)^2)+O*(Deltax)^4`

SKEWED RIGHT-SIDED DIFFERENCE

  • Here the number of nodes will be P+Q i.e for 4th order approximation of the second-order derivative it will be 4+2=6, thus 6 nodes will be required to derive the solution.
  • The 6 nodes for skewed right sided difference(forward differencing scheme) will be: i, i+1, i+2, i+3, i+4, i+5.
  • thus the equation will be
  • `Deltax^2.(d^2u)/dx^2= a.f(i)+b.f(i+1)+c.f(i+2)+d.f(i+3)+e.f(i+4)+g.f(i+5)`
  • Taylor table for the following is given as,

f(i)

`Deltax*f(i)`

`Deltax^2*f(i)`

`Deltax^3*f(i)`

`Deltax^4*f(i)`

`Deltax^5*f(i)`

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

`(4c)/2`

`(8c)/6`

`(16c)/24`

`(32c)/120`

d.f(i+3)

d

3d

`(9d)/2`

`(27d)/6`

`(81d)/24`

`(243d)/120`

e.f(i+4)

e

4e

`(16e)/2`

`(64e)/6`

`(256e)/24`

`(1024e)/120`

g.f(i+5)

g

5g

`(25g)/2`

`(125g)/6`

`(625g)/24`

`(3125g)/120`

`Deltax^2.f''(x)`

0

0

1

0

0

0

Thus there are 6 equation and 6 unkown, therefore equating them we get

`a+b+c+d+e+g=0`

`b+2c+3d+4e+5g=0`

`b/2+(4c)/2+(9d)/2+(16e)/2+(25g)/2=1`

`b/6+(8c)/6+(27d)/6+(64e)/6+(125g)/6=0`

`b/24+(16c)/24+(81d)/24+(256e)/24+(625g)/24=0`

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

Converting them into matrix form,

`[[1,1,1,1,1,1],[0,1,2,3,4,5], [0,1,4,9,16,25],[0,1,8,27,64,125],[0,1,16,81,256,625],[0,1,32,243,1024,3125]]*[[a],[b],[c],[d],[e],[g]]=[[0],[0],[1],[0],[0],[0]]`

close all
clear all
clc
 syms a b c d e g
 eqn1= a+b+c+d+e+g==0;
 eqn2= b+2*c+3*d+4*e+5*g==0;
 eqn3= b/2+(4*c)/2+(9*d)/2+(16*e)/2+(25*g)/2==1;
 eqn4= b/6+(8*c)/6+(27*d)/6+(64*e)/6+(125*g)/6==0;
 eqn5= b/24+(16*c)/24+(81*d)/24+(256*e)/24+(625*g)/24==0;
 eqn6= b/120+(32*c)/120+(243*d)/120+(1024*e)/120+(3125*g)/120==0;

 [A,B] = equationsToMatrix([eqn1, eqn2, eqn3, eqn4, eqn5, eqn6], [a, b, c, d, e, g])

 X = linsolve(A,B)

Output

A =
 
[ 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, 4/15, 81/40, 128/15, 625/24]
 
 
B =
 
 0
 0
 1
 0
 0
 0
 
 
X =
 
  15/4
 -77/6
 107/6
   -13
 61/12
  -5/6

Thus the new equation can be written by substituing the values of unknown in equation,

`Deltax^2.(d^2u)/dx^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)`

`(d^2u)/dx^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))/((Deltax)^2)+ o*(Deltax)^4`

SKEWED LEFT-SIDED DIFFERENCE

  • Here the number of nodes will be P+Q i.e for 4th order approximation of the second order derivative it will be 4+2=6, thus 6 nodes will be required to derive the solution.
  • The 6 nodes for skewed right sided difference(forward differencing scheme) will be: i, i-1, i-2, i-3, i-4, i-5.
  • `Deltax^2.(d^2u)/dx^2= a.f(i)+b.f(i-1)+c.f(i-2)+d.f(i-3)+e.f(i-4)+g.f(i-5)`
  • Taylor table for the following is given as,

f(i)

`Deltax*f(i)`

`Deltax^2*f(i)`

`Deltax^3*f(i)`

`Deltax^4*f(i)`

`Deltax^5*f(i)`

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

`(4c)/2`

`-(8c)/6`

`(16c)/24`

`-(32c)/120`

d.f(i-3)

d

-3d

`(9d)/2`

`-(27d)/6`

`(81d)/24`

`-(243d)/120`

e.f(i-4)

e

-4e

`(16e)/2`

`-(64e)/6`

`(256e)/24`

`-(1024e)/120`

g.f(i-5)

g

-5g

`(25g)/2`

`-(125g)/6`

`(625g)/24`

`-(3125g)/120`

`Deltax^2.f''(x)`

0

0

1

0

0

0

Thus there are 6 equation and 6 unkown, therefore equating them we get

`a+b+c+d+e+g=0`

`-b-2c-3d-4e-5g=0`

`b/2+(4c)/2+(9d)/2+(16e)/2+(25g)/2=1`

`-b/6-(8c)/6-(27d)/6-(64e)/6-(125g)/6=0`

`b/24+(16c)/24+(81d)/24+(256e)/24+(625g)/24=0`

`-b/120-(32c)/120-(243d)/120-(1024e)/120-(3125g)/120=0`

Converting them into matrix form,

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

close all
clear all
clc
 syms a b c d e g
 eqn1= a+b+c+d+e+g==0;
 eqn2= -b-2*c-3*d-4*e-5*g==0;
 eqn3= b/2+(4*c)/2+(9*d)/2+(16*e)/2+(25*g)/2==1;
 eqn4= -b/6-(8*c)/6-(27*d)/6-(64*e)/6-(125*g)/6==0;
 eqn5= b/24+(16*c)/24+(81*d)/24+(256*e)/24+(625*g)/24==0;
 eqn6= -b/120-(32*c)/120-(243*d)/120-(1024*e)/120-(3125*g)/120==0;
 
 [A,B] = equationsToMatrix([eqn1, eqn2, eqn3, eqn4, eqn5, eqn6], [a, b, c, d, e, g])
 
 X = linsolve(A,B)

Output:

A =
 
[ 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, -4/15, -81/40, -128/15, -625/24]
 
 
B =
 
 0
 0
 1
 0
 0
 0
 
 
X =
 
  15/4
 -77/6
 107/6
   -13
 61/12
  -5/6

Thus the new equation can be written by substituing the values of unknown in equation,

`Deltax^2.(d^2u)/dx^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)`

`(d^2u)/dx^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))/((Deltax)^2)+ o*(Deltax)^4`

MATLAB CODE

Main function:

close all
clear all
clc

x=pi/3;
dx= linspace(pi/4,pi/4000,20);

for i=1:length(dx)
    
   Skewed_right_order_error(i)=skewed_right_fourth_order_approximation(x,dx(i));
   
   Skewed_left_order_error(i)=skewed_left_fourth_order_approximation(x,dx(i));
   
   Central_order_error(i)=central_fourth_order_approximation(x,dx(i));
   
end

figure(1)
loglog(dx,Skewed_right_order_error,'LineWidth',2);

hold on

loglog(dx,Skewed_left_order_error,'r','LineWidth',2)

hold on
loglog(dx,Central_order_error,'g','LineWidth',2)
legend('Skewed right order error','Skewed left order error','central order error')
xlabel('dx values')
ylabel('Absolute Error')
title('Absolute Error for various differencing operation of fourth order derivative')

Function code for Central differencing fourth-order approximation

function out = central_fourth_order_approximation(x,dx)

%given function = exp(x)*cos(x)
analytical_derivaive=-2*exp(x)*sin(x);

central_order_differencing =  (((-1)/12)*(exp(x-2*dx)*cos(x-2*dx))+(4/3)*(exp(x-dx)*cos(x-dx))+((-5)/2)*(exp(x)*cos(x))+(4/3)*(exp(x+dx)*cos(x+dx))+((-1)/12)*(exp(x+2*dx)*cos(x+2*dx)))/((dx)^2);

out=abs(central_order_differencing-analytical_derivaive);

end

Function code for Skewed right differencing fourth-order approximation

function out = skewed_right_fourth_order_approximation(x,dx)

%given function = exp(x)*cos(x)
analytical_derivaive=-2*exp(x)*sin(x);

Skewed_right_order_differencing =  ((15/4)*(exp(x)*cos(x))-(77/6)*(exp(x+dx)*cos(x+dx))+(107/6)*(exp(x+2*dx)*cos(x+2*dx))-13*(exp(x+3*dx)*cos(x+3*dx))+(61/12)*(exp(x+4*dx)*cos(x+4*dx))-(5/6)*(exp(x+5*dx)*cos(x+5*dx)))/((dx)^2);

out=abs(Skewed_right_order_differencing-analytical_derivaive);

end

Function code for Skewed left differencing fourth-order approximation

function out = skewed_left_fourth_order_approximation(x,dx)

%given function = exp(x)*cos(x)
analytical_derivaive=-2*exp(x)*sin(x);

Skewed_left_order_differencing =  ((15/4)*(exp(x)*cos(x))-(77/6)*(exp(x-dx)*cos(x-dx))+(107/6)*(exp(x-2*dx)*cos(x-2*dx))-13*(exp(x-3*dx)*cos(x-3*dx))+(61/12)*(exp(x-4*dx)*cos(x-4*dx))-(5/6)*(exp(x-5*dx)*cos(x-5*dx)))/((dx)^2);

out=abs(Skewed_left_order_differencing-analytical_derivaive);

end

Output Graph:

Conclusion:

  • From the figure, we can conclude that error in central differencing is less as compared to both skewed right and skewed left order differencing
  • But not every time central differencing provides a good result, especially at boundaries where we, can't apply central differencing, there we need to apply either right or left skewed differencing.
  • Central differencing provides a more accurate result than the other two. Skewed right and skewed left provide an almost similar result.
Your portfolio, next

Ship work like this
in your own course.

Talk to a career counsellor

Free 20-minute call. We'll map the right course to where you want to work.

By submitting you agree to our Terms and Privacy Policy.

Talk to a career counsellor

Free 20-minute call. We'll map the right course to where you want to work.

By submitting you agree to our Terms and Privacy Policy.