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

Aim:-

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

Theory:-

Taylor Series

The Taylor series of a function is an infinite sum of terms that are expressed in terms of the function's derivatives at a single point.

f(x)=f(a)+f'(a)1!(x−a)+f''(a)2!(x−a)2+f'''(a)3!(x−a)3+...,

1) Central Difference

The number of Nodes required for 4th order approximation of a second-order derivative can be given by the equation

N=P+Q−1 where

N is the number of nodes

Q is the order of derivative

P is the order of approximation

So, here we need 5 nodes in order to find the 2nd order derivative with the 4th order of approximation.

Consider the value of ∂^2y/∂x^2 with the 4th order of approximation is given by the equation

∂^2y/∂x^2=af(i−2)+bf(i−1)+cf(i)+df(i+1)+ef(i+2)

Consider Taylor Table for central difference

Taylor Table
Assumed Function		Taylor Series Expansion
		f(i)	Δxf'(i)	Δx2f''(i)	Δx3f'''(i)	Δx4f''''(i)	Δx5f'''''(i)
Stencil	af(i−2)	a	−2a	4a/2	−8a/6	16a/24	−32a/120
	bf(i−1)	b	−b	b/2	−b/6	b/24	−b/120
	cf(i)	c	0	0	0	0	0
	df(i+1)	d	d	d/2	d/6	d/24	−c/120
	ef(i+2)	e	2e	4e/2	−8e/6	16e/24	32e//120
Sum of Coefficients		0	0	1	0	0	?

Solving using MATLAB to find the values of the coefficients

>> 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=AB

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 =

     0
     0
     1
     0
     0


X =

   -0.0833
    1.3333
   -2.5000
    1.3333
   -0.0833

Now we can write the equation using the Taylors Table as

∂^2y/∂x^2=af(i−2)+bf(i−1)+cf(i)+df(i+1)+ef(i+2)/Δx^2

2) Skewed right-sided difference

The number of Nodes required for 4th order approximation of a second-order derivative can be given by the equation

N=P+Q where

N is the number of nodes

Q is the order of derivative

P is the order of approximation

So, here we need 6 nodes in order to find the 2nd order derivative with the 4th order of approximation.

Consider the value of ∂^2y/∂x^2 with the 4th order of approximation is given by the equation

∂2y∂x2=af(i)+bf(i+1)+cf(i+2)+df(i+3)+ef(i+4)+gf(i+5)

Now, using taylor series to find the values of discrete poits with repect to point i

af(i)=af(i)

bf(i+1)=bf(i)+bΔxf'(i)+(bΔx^2/2!)f''(i)+(bΔx^3/3!)f'''(i)+(bΔx^4/4!)f''''(i)+(bΔx^5/5!)f'''''(i)+(bΔx^6/6!)f''''''(i)+...

cf(i+2)=cf(i)+c(2Δx)f'(i)+(c(2Δx)^2/2!)f''(i)+(c(2Δx)^3/3!)f'''(i)+(c(2Δx)^4/4!)f''''(i)+(c(2Δx)^5/5!)f'''''(i)+(c(2Δx)^6/6!)f''''''(i)+...

df(i+3)=df(i)+d(3Δx)f'(i)+(d(3Δx)^2/2!)f''(i)+(d(3Δx)^3/3!)f'''(i)+(d(3Δx)^4/4!)f''''(i)+(d(3Δx)^5/5!)f'''''(i)+(d(3Δx)^6/6!)f''''''(i)+...

ef(i+4)=ef(i)+e(4Δx)f'(i)+(e(4Δx)^2/2!)f''(i)+(e(4Δx)^3/3!)f'''(i)+(e(4Δx)^4/4!)f''''(i)+(e(4Δx)^5/5!)f'''''(i)+(e(4Δx)^6/6!)f''''''(i)+...

gf(i+5)=gf(i)+g(5Δx)f'(i)+(g(5Δx)^2/2!)f''(i)+(g(5Δx)^3/3!)f'''(i)+(g(5Δx)^4/4!)f''''(i)+(g(5Δx)^5/5!)f'''''(i)+(g(5Δx)^6/6!)f''''''(i)+...

Adding the LHS and RHS of the above equations we get

af(i)+bf(i+1)+cf(i+2)+df(i+3)+ef(i+4)+gf(i+5)=f(i)[a+b+c+d+e+g]+Δxf'(i)[b+2c+3d+4e+5g]+Δx2f''(i)[b2+4c2+9d2+16e2+25g2]+Δx3f'''(i)[b6+8c6+27d6+64e6+125g6]+Δx4f''''(i)[b24+16c24+81d24+256e24+625g24]+Δx5f'''''(i)[b120+32c120+243d120+1024e120+3125g120]

In the above equations, we want only the coefficients of Δx2f''(i) and we assume the coefficient to be one and for the remaining terms we assume the coefficients to be zero.

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

Taylor Table
Assumed Function		Taylor Series Expansion
		f(i)	Δf'(i)	Δx2f''(i)	Δx3f'''(i)	Δx4f''''(i)	Δx5f'''''(i)	Δx6f''''''(i)
Stencil	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	4c/2	8c6	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	64e/6	256e/24	1024e/120	4096e/720
	gf(i+5)	g	5g	25g/2	125g/6	625g/24	3125g/120	15625g/720
Sum of Coefficients		0	0	1	0	0	0	?

Solving using MATLAB to find the values of the coefficients

>> 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]
B=[0;0;1;0;0;0]
X=AB

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


X =

    3.7500
  -12.8333
   17.8333
  -13.0000
    5.0833
   -0.8333

Now we can write the equation using the Taylors Table as ∂^2y/∂x^2=af(i)+bf(i+1)+cf(i+2)+df(i+3)+ef(i+4)+gf(i+5)/Δx^2

3) Skewed left-sided difference

The number of Nodes required for 4th order approximation of a second-order derivative can be given by the equation

N=P+Qwhere

N is the number of nodes

Q is the order of derivative

P is the order of approximation

So, here we need 6 nodes in order to find the 2nd order derivative with the 4th order of approximation.

Consider the value of ∂^2y/∂x^2 with the 4th order of approximation is given by the equation

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

Consider Taylor Table for left-sided difference

Taylor Table
Assumed Function		Taylor Series Expansion
		f(i)	Δxf'(i)	Δx2f''(i)	Δx3f'''(i)	Δx4f''''(i)	Δx5f'''''(i)	Δx6f''''''(i)
Stencil	af(i)	a	00	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	4c/2	−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	−64e/6	256e/24	−1024e/120	4096e/720
	gf(i−5)	g	−5g	25g/2	−125g/6	625g/24	−3125g/120	15625g/720
Sum of Coefficients		0	0	1	0	0	0	?

Solving using MATLAB to find the values of the coefficients

>> 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]
B=[0;0;1;0;0;0]
X=AB

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


X =

    3.7500
  -12.8333
   17.8333
  -13.0000
    5.0833
   -0.8333

Now we can write the equation using the Taylors Table as ∂^2y/∂x^2=af(i)+bf(i−1)+cf(i−2)+df(i−3)+ef(i−4)+gf(i−5)

MATLAB CODE:

MAIN CODE:

clear all
close all
clc

%value of x
x= pi/3;

%range of dx
dx= linspace(pi/4,pi/4000,20);

%for loop
for i=1:length(dx)
    
    central_differencing(i) = fourth_order_central_differencing(x,dx(i));
    
    right_sided_skew(i) = fourth_order_right_sided_skew(x,dx(i));
    
    left_sided_skew(i) = fourth_order_left_sided_skew(x,dx(i));
    
end

%plotting
figure(1)
loglog(dx,central_differencing,'r')
hold on
loglog(dx,right_sided_skew,'g')
hold on
loglog(dx,left_sided_skew,'b')
xlabel('range of dx')
ylabel('truncation error')
title('range of dx vs error')
legend('central differencing','skewed right sided difference','skewed left sided difference','location','northwest')


function Error_central_differencing = fourth_order_central_differencing(x,dx)

% analytic function = exp(x)*cos(x)

% analytic derivative
%f''(x) = -2*exp(x)*sin(x)

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

%numerical derivative

%central differencing
central_differencing = (((-0.0833)*((exp(x-(2*dx)))*(cos(x-(2*dx)))))+((1.3333)*((exp(x-dx))*(cos(x-dx))))-((2.5000)*((exp(x))*(cos(x))))+((1.3333)*((exp(x+dx))*(cos(x+dx))))-((0.0833)*((exp(x+(2*dx)))*(cos(x+(2*dx))))))/(dx^2)

%error
Error_central_differencing = abs(central_differencing- analytical_derivative)
end


function Error_right_sided_skew = fourth_order_right_sided_skew(x,dx)

% analytic function = exp(x)*cos(x)

% analytic derivative
%f''(x) = -2*exp(x)*sin(x)

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

%numerical derivative

%central differencing
right_sided_skew = (((3.7500)*((exp(x))*(cos(x))))-((12.8333)*((exp(x+dx))*(cos(x+dx))))+((17.8333)*((exp(x+(2*dx)))*(cos(x+(2*dx)))))-((13.0000)*((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)

%error
Error_right_sided_skew = abs(right_sided_skew- analytical_derivative)
end


function Error_left_sided_skew = fourth_order_left_sided_skew(x,dx)

% analytic function = exp(x)*cos(x)

% analytic derivative
%f''(x) = -2*exp(x)*sin(x)

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

%numerical derivative

%central differencing
left_sided_skew = (((-0.8333)*((exp(x-(5*dx)))*(cos(x-(5*dx)))))+((5.0833)*((exp(x-(4*dx)))*(cos(x-(4*dx)))))-((13.0000)*((exp(x-(3*dx)))*(cos(x-(3*dx)))))+((17.8333)*((exp(x-(2*dx)))*(cos(x-(2*dx)))))-((12.8333)*((exp(x-dx))*(cos(x-dx))))+((3.7500)*((exp(x))*(cos(x)))))/(dx^2)

%error
Error_left_sided_skew = abs(left_sided_skew- analytical_derivative)
end

RIGHT SIDE SKEWED :

LEFT SIDE SKEWED :

CENTRAL DIFFRENCING :

PLOT :

SNAPSHOT :

Conclusion:-

The comparison of absolute error for the left, right and central differencing scheme can be seen in the below graph. The graph shows that the central scheme is most accurate compared to the skewed differencing. The curves also show that when the value of dx is very small or large the approximating function becomes unstable. This instability is caused due to the combination of various errors such as roundoff error, residual error and truncation error.

From the below graph it is seen that the range of dx where the function is stable ranges from 0.1 to 0.01, beyond this range we see an abrupt rise or fall.