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

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

AAmol Patel19 Jun 2021
MATLAB

Fourth Order Approximation of Second Order Derivative using Central Differencing

We will derive the fourth order approximation for second order derivative using a taylors table method:

The fourth order approximation uses four point nearby the point where we have to find the derivative , it can be seen in the following figure that to find the derivative at point `i` we use two point to the left and two point from the right of point `i`, this indicated that we have used central differencing.

let the spacing between each nodes be `Deltax`

the derivative of the fourth order approximation for the second order derivative at point `i` can be given using

`Deltax^2(del^2f)/(delx^2) = a f(i-2) + b f(i-1) + cf(i) + df(i+1) + ef(i+2)` ........(1`

now using taylor table method we can find the coefficients a, b, c, d, e.

`f(i)` `Deltax .(delf)/(delx)` `Deltax^2.(del^2f)/(delx^2)` `Deltax^3.(del^3f)/(delx^3)` `Deltax^4.(del^4f)/(delx^4)` `Deltax^5.(del^5f)/(delx^5)` `Deltax^6.(del^6f)/(delx^6)`

`af(i-2)`

 `a` `-2.a` `4/(2!).a` `-8/(3!).a` `16/(4!).a` `-32/(5!).a` `64/(6!).a`
`bf(i-1)` `b` `-b` `b/(2!)` `-b/(3!)` `b/(4!)` `-b/(5!)` `b/(6!)`

`cf(i)`

 `c` `0` `0` `0` `0` `0` `0`

`df(i+1)`

 `d` `d` `d/(2!)` `d/(3!)` `d/(4!)` `d/(5!)` `d/(6!)`

`ef(i+2)`

 `e` `2.e` `4/(2!).e` `8/(3!).e` `16/(4!).e` `32/(5!).e` `64/(6!).e`

sum

 `0` `0` `1` `0` `0` `?` `?`

now we can easily find the coefficients a, b, c, d, e using the 5 equation that we get from the Taylors Table

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

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

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

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

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

Now using the coefficient matrix we can find the values of a,b,c,d,e in MATLAB

>> A

A =

    1.0000    1.0000    1.0000    1.0000    1.0000
   -2.0000   -1.0000         0    1.0000    2.0000
    2.0000    0.5000         0    0.5000    2.0000
   -1.3333   -0.1667         0    0.1667    1.3333
    0.6667    0.0417         0    0.0417    0.6667

>> B

B =

     0
     0
     1
     0
     0

>> inv(A)*B

ans =

   -0.0833
    1.3333
   -2.5000
    1.3333
   -0.0833

form here we get a = -0.0833, b = 1.3333, c = -2.5000, d = 1.3333, e = -0.0833.

If we consider the sum of coefficients of the terms with `Deltax^5` in the central differencing scheme we get

`-32a/(5!)-b/(5!) +0.c+d/(5!)+32e/(5!)`

puting the values of a,b,c,d,e from the central differencing case we get `0`.

and the sum of the coefficient of the terms with `Deltax^6` we get

`(64a+b+0.c +d+64e)/(6!) = -103.9958`

i.e the term with `Deltax^6`is non zero. so the order of approximation = smallest order of term in error - order of derivative

`therefore sf" order of approximation" = 6 - 2 = 4`

so our equation (1) becomes

`(del^2f)/(delx^2) =( -0.0833f(i-2) +1.3333 f(i-1) -2.5000f(i) +1.3333f(i+1) -0.0833f(i+2))/(Deltax^2)` ......(2)`

Thus, we have derived the fourth order approxiamtion of the second order derivative.

Now , Let `N` be the number of nodes used in the derivative . As we have 5 nodes , `N=5`

also the order of derivation is second order, let it be denoted by `Q` ,so we have ` Q=2`

Let the order of accuracy be denoted by `P` so we have `P =4`

we can now give the relation between the order of accuracy, the order of derivative and the nu,ber of nodes used as

`P = N - Q + 1`

i.e for symmetry/central differencing ` sf" order of accuracy" =sf"number of nodes" -sf"order of derivation" +1`

Skewed right-sided diference for Fourth order approximation of second order derivative

Consider the following figure which shows skewed arrangement of nodes on the right side or in forward differencing.

We can write the approximation at nose `i` using all the 5 nodes at the rigth side and it can be represented as

`Deltax^2(del^2f)/(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)` ........(3)

we can expand each term of the above in the form of taylor series and we will write it in the form of taylors table

`f(i)` `Deltax.(delf)/(delx)` `Deltax^2.(del^2f)/(delx^2)` `Deltax^3.(del^3f)/(delx^3)` `Deltax^4.(del^4f)/(delx^4)` `Deltax^5.(del^5f)/(delx^5)` `Deltax^6.(del^6f)/(delx^6)`
`a.f(i)` a 0 0 0 0 0 0
`b.f(i+1)` b b b/(2!) b/(3!) b/(4!) b/(5!) b/(6!)
`c.f(i+2)` c 2c 4c/(2!) 8c/(3!) 16c/(4!) 32c/(5!) 64c/(6!)
`d.f(i+3)` d 3d 9d/(2!) 27d/(3!) 81d/(4!) 243d/(5!) 729d/(6!)
`e.f(i+4)` e 4e 16e/(2!) 64e/(3!) 256e/(4!) 1024e/(5!) 4096e/(6!)
`g.f(i+5)` g 5g 25g/(2!) 125g/(3!) 625g/(5!) 3125g/(5!) 15625g/(6!)
sum 0 0 1 0 0 0 ?

as we have 6 nodes we will have 6 equations

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

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

`b/2+4c/2+9d/2+16e/2+25g/5=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`

solving the coefficient matrix

>> 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]

A =

    1.0000    1.0000    1.0000    1.0000    1.0000    1.0000
         0    1.0000    2.0000    3.0000    4.0000    5.0000
         0    0.5000    2.0000    4.5000    8.0000   12.5000
         0    0.1667    1.3333    4.5000   10.6667   20.8333
         0    0.0417    0.6667    3.3750   10.6667   26.0417
         0    0.0083    0.2667    2.0250    8.5333   26.0417

>> B = [0;0;1;0;0;0]

B =

     0
     0
     1
     0
     0
     0

>> X=inv(A)*B

X =

    3.7500
  -12.8333
   17.8333
  -13.0000
    5.0833
   -0.8333

we have a = 3.75, b = -12.8333, c = 17.8333, d = -13.0000, e = 5.0833, g = -0.8333

the coefficient of term having `Deltax^6` becomes

`(-12.8333+64*17.8333-729*13.0000+4096*5.0833-15625*0.8333)/720 = -0.7606`

the coefficient is non-zero so the largest term of truncating error is of order 4 ( `Deltax^6` will be divided by `Deltax^2` so we will be remained with `Deltax^4`)

So our equation for the foruth order accurate approximation of second order derivative becomes :

`(del^2f)/(delx^2) = (3.7500f(i) -12.8333f(i+1) +17.8333f(i+2)-13.0000f(i+3)+5.0833f(i+4)-0.8333f(i+5))/(Deltax^2)` ........(4

Here the number of nodes, N = 6 ,

the order of dervative, Q = 2

order of accuracy , P = number of nodes, N - order of derivaitve, Q

P = N - Q = 6 - 2 = 4

so we have our four order approximation.

Skewed left-sided diference for Fourth order approximation of second order derivative

Consider the following figure which shows skewed arrangement of nodes on the left side or in backward differencing.

We can write the approximation at nose `i` using all the 5 nodes at the left side and it can be represented as

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

we can expand each term of the above in the form of taylor series and we will write it in the form of taylors table

`f(i)` `Deltax.(delf)/(delx)` `Deltax^2.(del^2f)/(delx^2)` `Deltax^3.(del^3f)/(delx^3)` `Deltax^4.(del^4f)/(delx^4)` `Deltax^5.(del^5f)/(delx^5)` `Deltax^6.(del^6f)/(delx^6)`
`a.f(i-5)` a -5a 25a/(2!) -125a/(3!) 625a/(4!) -3125a/(5!) 15625a/(6!)
`b.f(i-4)` b -4b 16b/(2!) -64b/(3!) 256b/(4!) -1024b/(5!) 4096b/(6!)
`c.f(i-3)` c -3c 9c/(2!) -27c/(3!) 81c/(4!) -243c/(5!) 729c/(6!)
`d.f(i-2)` d -2d 4d/(2!) -8d/(3!) 16d/(4!) -32d/(5!) 64d/(6!)
`e.f(i-1)` e -e e/(2!) -e/(3!) e/(4!) -e/(5!) e/(6!)
`g.f(i)` g 0 0 0 0 0 0
sum 0 0 1 0 0 0 ?

as we have 6 nodes we will have 6 equations

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

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

`25a/2+16b/2+9c/2+4d/2+e/2=1`

`-125a/6-64b/6-27c/6-8d/6-e/6=0`

`625a/24+256b/24+81c/24+16d/24+e/24=0`

`-3125a/120-1024b/120-243c/120-32d/120-e/120=0`

solving the coefficient matrix

>> A=[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]

A =

    1.0000    1.0000    1.0000    1.0000    1.0000    1.0000
   -5.0000   -4.0000   -3.0000   -2.0000   -1.0000         0
   12.5000    8.0000    4.5000    2.0000    0.5000         0
  -20.8333  -10.6667   -4.5000   -1.3333   -0.1667         0
   26.0417   10.6667    3.3750    0.6667    0.0417         0
  -26.0417   -8.5333   -2.0250   -0.2667   -0.0083         0

>> B=[0;0;1;0;0;0]

B =

     0
     0
     1
     0
     0
     0

>> X=inv(A)*B

X =

   -0.8333
    5.0833
  -13.0000
   17.8333
  -12.8333
    3.7500

we have a = -0.8333, b = 5.0833, c = -13.0000, d = 17.8333, e = -12.8333, g = 3.7500

the coefficient of term having `Deltax^6` becomes

`(-15625*0.8333 + 4096*5.0833 -729*13 +64*17.8333 -12.8333+0)/720= -0.7606`

the coefficient is non-zero so the largest term of truncating error is of order 4 ( `Deltax^6` will be divided by `Deltax^2` so we will be remained with `Deltax^4`)

So our equation for the foruth order accurate approximation of second order derivative becomes :

`(del^2f)/(delx^2) = (-0.8333*f(i-5) +5.0833*f(i-4) -13.0000*f(i-3)+17.8333*f(i-2)-12.8333*f(i-1)+3.7500*f(i))/(Deltax^2)` ........(6)

Here the number of nodes, N = 6 ,

the order of dervative, Q = 2

`sf"order of accuracy , P = number of nodes, N - order of derivaitve, Q"`

P = N - Q = 6 - 2 = 4

so we have our four order approximation.

It is clear that we have to use one less node in symetry scheme than in the skewed scheme

so for symmety we have the formula : `sf"order of accuracy,P = number of nodes,N - order of derivative,Q +1"`

Solving the Program in MATLAB to compare the error between the analytical and numerical approxiamtions

Now we will write a program in MATLAB to compare the analytical derivative with the 3 numerical approximations that we have derived

First we need to write a function for each approximation so that we can call that functions in our program to find the error between the analytical and the numerical values.

the first function that we write is for the fourth order approximation of the second order derivative using central difference method.

we have named this function as fourth_order_approximation_second_derivative_cd and it is the function of x and dx.

function [out]= fourth_order_approximation_second_derivative_cd(x,dx)

% analyitcal funtion = exp(x)*cos(x)
% second order derivative of analytical function = -2exp(x)*sin(x)
analytical_derivative = (-2*exp(x)*sin(x))

% coefficient matrix
A= [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]

% solution matrix
B = [0;0;1;0;0]

% result matrix
X = inv(A)*B

% fourth orer approximation ofsecond order numerical derivative using central difference method
central_difference = (X(1)*(exp(x-2*dx)*cos(x-2*dx))+X(2)*(exp(x-1*dx)*cos(x-1*dx))+X(3)*(exp(x)*cos(x))+X(4)*(exp(x+1*dx)*cos(x+1*dx))+X(5)*(exp(x+2*dx)*cos(x+2*dx)))/(dx^2)

% absolute error beteeen analytical and numerical derivaitve
out = abs(analytical_derivative - central_difference)


end

next we write the function for the fourth order approximation of the second order derivative using forward difference method.

we have named this function as fourth_order_approximation_second_derivative_fd and it is the function of x and dx.

function [out]= fourth_order_approximation_second_derivative_fd(x,dx)

% analyitcal funtion = exp(x)*cos(x)
% second order derivative of analytical function = -2exp(x)*sin(x)
analytical_derivative = (-2*exp(x)*sin(x))

% coefficeint matrix
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]

% solution matrix
B = [0;0;1;0;0;0]

% result matrix
X=inv(A)*B

% fourth order approximation of second order numerical derivaitve using forward difference method
forward_difference = (X(1)*(exp(x)*cos(x))+X(2)*(exp(x+1*dx)*cos(x+1*dx))+X(3)*(exp(x+2*dx)*cos(x+2*dx))+X(4)*(exp(x+3*dx)*cos(x+3*dx))+X(5)*(exp(x+4*dx)*cos(x+4*dx))+X(6)*(exp(x+5*dx)*cos(x+5*dx)))/(dx^2)

% absolute error between analytical and numerical derivatives
out = abs(analytical_derivative - forward_difference) 


end

next we write the function for the fourth order approximation of the second order derivative using backward difference method.

we have named this function as fourth_order_approximation_second_derivative_bd and it is the function of x and dx.

function [out]= fourth_order_approximation_second_derivative_bd(x,dx)

% analyitcal funtion = exp(x)*cos(x)
% second order derivative of analytical function = -2exp(x)*sin(x)
analytical_derivative = (-2*exp(x)*sin(x))

% coefficient matrix
A=[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]

% solution matrix
B=[0;0;1;0;0;0]

% result matrix
X=inv(A)*B

% fourth order approximation of second order numerical derivative using backward difference method
backward_difference = (X(1)*(exp(x-5*dx)*cos(x-5*dx))+X(2)*(exp(x-4*dx)*cos(x-4*dx))+X(3)*(exp(x-3*dx)*cos(x-3*dx))+X(4)*(exp(x-2*dx)*cos(x-2*dx))+X(5)*(exp(x-1*dx)*cos(x-1*dx))+X(6)*(exp(x)*cos(x)))/(dx^2)

% absolute error between analytical and numerical derivaitve
out = abs(analytical_derivative - backward_difference)

end

we have calculated to coefficeints for our approximation within the function to reduce the error while rounding off in all the 3 differencing methods.

now we will write our main program 

clear all
close all
clc

% assigning the initial value for x and a range of values for dx
x = pi/3;
dx = linspace(pi/4000,pi/4,40)

% using for loop for each value of dx 
for i = 1:length(dx)
    % using variables to store the value of error 
    error_central_differencing(i) = fourth_order_approximation_second_derivative_cd(x,dx(i))
    error_forward_differencing(i) = fourth_order_approximation_second_derivative_fd(x,dx(i))
    error_backward_differencing(i) = fourth_order_approximation_second_derivative_bd(x,dx(i))
    
end

%plotting the loglog graph for error vs the value of dx in each differencing schemes
figure(1)
loglog(dx,error_central_differencing,'-or')
hold on
loglog(dx,error_forward_differencing,'-db')
hold on
loglog(dx,error_backward_differencing,'-sg')
hold on
% decorating the graph 
legend('fourth order approximation of second order derivative using central differencing','fourth order approxiamtion of second order derivative using forward differencing','fourth order approximation of second order derivative using backward differencing','location','northwest')
xlabel('dx')
ylabel('error')
title('errors in central, forward and backward differencing schemes')

we can see the output figure of the loglog graph as following

for the figure we can see that the central differencing method has the lowest error, also we can obeserve the after a certain value of dx the round off errors dominates the solution.

We now know that the central differencing has the least error, but there are places where we cannot use this methods such as regions like boundary layer where we have nodes only at one side of our initial node so there we cant use central differencing method in such place skewed schemes either left or right handed can be used depending on the situation.

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.