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

AIM:- Derive the 4th order approximations of the second-order derivative using

1. Central difference (CDS)
2. Skewed right-sided difference (SRSD)
3. Skewed left-sided difference (SLSD)

And write a program in MATLAB to evaluate the second-order derivative of the analytical function and compare it with the derived 3 numerical approximations .

OBJECTIVE:- 1. Create a plot that compares the absolute error between the CDS,SRSD,SLSD.

2. Describe why a skewed scheme is useful? What can a skewed scheme do that a CD scheme cannot?

1. CENTRAL DIFFERENCING SCHEME (CDS):

No. of Nodes for CDS:- n= p+q-1

Where, p is order of derivative.

q is order of approximations.

The no of nodes = 2+4-1 = 5.

For our 5 stencil points are:- (i-2), (i-1),i, (i+1), (i+2)

From Taylor Table we know that fourth order approximation for the second order derivative can be find out using:-

`(∂^2u)/(∂x^2)=(a1*f(i−2)+b1*f(i−1)+c1*f(i)+d1*f(i+1)+e1*f(i+2))/dx^2` ….1

taylor series expansion for , a1*f(i−2) term

`a1f(i−2)=a1*f(i)−(a1f'(i)(2*dx))/(1!)+(a1.f''(i)(2*dx)^2)/(2!)-(a1*f'''(i)(2*dx)^3)/(3!)+(a1*f''''(i)(2*dx^4))/(4!)-……`

taylor series expansion for , b1*f(i−1) term

`b1*f(i−1)=b1*f(i)−(b1*f'(i)(dx))/(1!)+(b1*f''(i)(dx)^2)/(2!)-(b1*f'''(i)(dx)^3)/(3!)+(b1*f''''(i)(dx^4))/(4!)-……`

taylor series expansion for , c1*f(i) term

`c1*f(i)=c1*f(i)`

taylor series expansion for , d1*f(i+1) term

`d1*f(i+1)=d1*f(i)+(d1*f'(i)(dx))/(1!)+(d1*f''(i)(dx)^2)/(2!)+(d1*f'''(i)(dx)^3)/(3!)+(d1*f''''(i)(dx^4))/(4!)+……`

`e1*f(i+2)=e1*f(i)+(e1*f'(i)(2*dx))/(1!)+(e1*f''(i)(2*dx)^2)/(2!)+(e1*f'''(i)(2*dx)^3)/(3!)+(e1*f''''(i)(2*dx^4))/(4!)+……`

Solving the linear equation obtained from the above table using matrix inversion method,

A1.X1 =B1

`[[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],[16/24,1/24,0,1/24,16/24]]*[[a1],[b1],[c1],[d1],[e1]]=[[0],[0],[1],[0],[0]]`

X1 = inv(A1)*B1

a1=-0.0833 ,b1=1.3333 ,c1=-2.5000 ,d1=1.3333 ,e1=-0.0833

putting values of a1,b1,c1,d1,e1 in main eqn of CDS

`(∂^2u)/(∂x^2)=((-0.0833)*f(i−2)+(1.3333)*f(i−1)+(-2.5000)*f(i)+(1.3333)*f(i+1)+(-0.0833)*f(i+2))/dx^2`

2. Skewed right-sided difference (SRSD)

No. of Nodes for SRSD :- n= p+q

Where, p is order of derivative.

q is order of approximations.

The no of nodes = 2+4=6

our 6 stencil points are:- (i), (i+1), (i+2),(i+3), (i+4),(i+5)

From Taylor Table we know that fourth order approximation for the second order derivative can be find out using:-

`(∂^2u)/(∂x^2)=(a2*f(i)+b2*f(i+1)+c2*f(i+2)+d2*f(i+3)+e2*f(i+4)+g2*f(i+5))/dx^2` ….1

taylor series expansion for , a2*f(i) term

`a2f(i)=a2*f(i)`

taylor series expansion for , b2*f(i+1) term

`b2*f(i+1)=b2*f(i)+(b2*f'(i)(dx))/(1!)+(b2*f''(i)(dx)^2)/(2!)+(b2*f'''(i)(dx)^3)/(3!)+(b2*f''''(i)(dx^4))/(4!)+b2*f'''''(i)(dx^5)+……`

taylor series expansion for , c2*f(i+2) term

`c2*f(i+2)=c2*f(i)+(c2*f'(i)(2*dx))/(1!)+(c2*f''(i)(2*dx)^2)/(2!)+(c2*f'''(i)(2*dx)^3)/(3!)+(c2*f''''(i)(2*dx^4))/(4!)+(c2*f'''''(i)(2*dx^5))/(5!)+......`

taylor series expansion for , d2*f(i+3) term

`d2*f(i+3)=d2*f(i)+(d2*f'(i)(3*dx))/(1!)+(d2*f''(i)(3*dx)^2)/(2!)+(d2*f'''(i)(3*dx)^3)/(3!)+(d2*f''''(i)(3*dx^4))/(4!)+(d2*f'''''(i)(3*dx^5))/(5!)……`

`e2*f(i+4)=e2*f(i)+(e2*f'(i)(4*dx))/(1!)+(e2*f''(i)(4*dx)^2)/(2!)+(e2*f'''(i)(4*dx)^3)/(3!)+(e2*f''''(i)(4*dx^4))/(4!)+(e2*f'''''(i)(4*dx^5))/(5!)……`

`g2*f(i+5)=g2*f(i)+(g2*f'(i)(5*dx))/(1!)+(g2*f''(i)(5*dx)^2)/(2!)+(g2*f'''(i)(5*dx)^3)/(3!)+(g2*f''''(i)(5*dx^4))/(4!)+(g2*f'''''(i)(5*dx^5))/(5!)……`

Solving the linear equation obtained from the above table using matrix inversion method,

A2.X2 =B2

`[[1,1,1,1,1,1],[0,1,2,3,4,5],[0,1/2,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]]*[[a2],[b2],[c2],[d2],[e2],[g2]]=[[0],[0],[1],[0],[0],[0]]`

X2 = inv(A2)*B2

a2= 3.7500 ,b2= -12.8333 ,c2=17.8333 ,d2= -13.0000 ,e2= 5.0833 ,g2= -0.8333

putting values of a2,b2,c2,d2,e2 in main eqn of SRSD

`(∂^2u)/(∂x^2)=((3.7500)*f(i)+(-12.8333)*f(i+1)+(17.8333)*f(i+2)+(-13.0000)*f(i+3)+(5.0833)*f(i+4)+(-0.8333)*f(i+5))/dx^2`

3. Skewed left-sided difference (SLSD)

No. of Nodes for SLSD:- n= p+q

Where, p is order of derivative.

q is order of approximations.

The no of nodes = 2+4=6

our 6 stencil points are:- (i), (i-1), (i-2),(i-3), (i-4),(i-5)

From Taylor Table we know that fourth order approximation for the second order derivative can be find out using:-

`(∂^2u)/(∂x^2)=(a3*f(i)+b3*f(i-1)+c3*f(i-2)+d3*f(i-3)+e3*f(i-4)+g3*f(i-5))/dx^2` ….1

taylor series expansion for , a3*f(i) term

`a3f(i)=a3*f(i)`

taylor series expansion for , b3*f(i-1) term

`b3*f(i-1)=b3*f(i)-(b3*f'(i)(dx))/(1!)+(b3*f''(i)(dx)^2)/(2!)-(b3*f'''(i)(dx)^3)/(3!)+(b3*f''''(i)(dx^4))/(4!)-b3*f'''''(i)(dx^5)+……`

taylor series expansion for , c3*f(i-2) term

`c3*f(i-2)=c3*f(i)-(c3*f'(i)(2*dx))/(1!)+(c3*f''(i)(2*dx)^2)/(2!)-(c3*f'''(i)(2*dx)^3)/(3!)+(c3*f''''(i)(2*dx^4))/(4!)-(c3*f'''''(i)(2*dx^5))/(5!)+......`

taylor series expansion for , d3*f(i-3) term

`d3*f(i-3)=d3*f(i)-(d2*f'(i)(3*dx))/(1!)+(d3*f''(i)(3*dx)^2)/(2!)-(d3*f'''(i)(3*dx)^3)/(3!)+(d3*f''''(i)(3*dx^4))/(4!)-(d3*f'''''(i)(3*dx^5))/(5!)+……`

`e3*f(i-4)=e3*f(i)-(e3*f'(i)(4*dx))/(1!)+(e3*f''(i)(4*dx)^2)/(2!)-(e3*f'''(i)(4*dx)^3)/(3!)+(e3*f''''(i)(4*dx^4))/(4!)-(e3*f'''''(i)(4*dx^5))/(5!)+……`

`g3*f(i-5)=g3*f(i)-(g3*f'(i)(5*dx))/(1!)+(g3*f''(i)(5*dx)^2)/(2!)-(g3*f'''(i)(5*dx)^3)/(3!)+(g3*f''''(i)(5*dx^4))/(4!)-(g3*f'''''(i)(5*dx^5))/(5!)+……`

Solving the linear equation obtained from the above table using matrix inversion method,

A3.X3 =B3

`[[1,1,1,1,1,1],[0,-1,-2,-3,-4,-5],[0,1/2,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]]*[[a2],[b2],[c2],[d2],[e2],[g2]]=[[0],[0],[1],[0],[0],[0]]`

X3 = inv(A3)*B3

a3=3.7500 ,b3= -12.8333,c3= 17.8333 ,d3= -13.0000,e3= 5.0833,g3= -0.8333

putting values of a3,b3,c3,d3,e3 in main eqn of SLSD

`(∂^2u)/(∂x^2)=((3.7500)*f(i)+(-12.8333)*f(i-1)+( 17.8333)*f(i-2)+(-13.0000)*f(i-3)+(5.0833)*f(i-4)+(-0.8333)*f(i-5))/dx^2`

Main Code:-

clear all
close all
clc
% defining x 
x=pi/3;
% defining dx
dx=linspace(pi/40000,pi/4,200);
%Finding X1 for central differencing Scheme (CDS)
A1=[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]
B1=[0;0;1;0;0]; 
X1=inv(A1)*B1;
%  Equation for Central_difference_approximation=(a1*(f(i-2))+b1*(f(i-1))+c1*(f(i))+d1*(f(i+1))+e1*(f(i+2)))/(dx^2);
%Finding X2 Skewed Right sided Scheme (SRSD)
A2=[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];
B2=[0;0;1;0;0;0];
X2=(inv(A2))*B2;
% Equation of Skewed_right_sided_approximation=(a2*(f(i))+b2*(f(i+1))+c2*(f(i+2))+d2*(f(i+3))+e2*(f(i+4))+g2*(f(i+5)))/(dx^2);
% Skewed Left sided Scheme (SLSD)
A3=[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];
B3=[0;0;1;0;0;0];
X3=inv(A3)*B3;
%Skewed_left_sided_approximation=(a3*(f(i))+b3*(f(i-1))+c3*(f(i-2))+d3*(f(i-3))+e3*(f(i-4))+g3*(f(i-5)))/(dx^2);

analytical_derivative=(-2*exp(x)*sin(x)); % analytical_function= exp(x)*cos(x);
% Calculating 4th order approximation of second order derivative using, CDS,SRSD,SLSD
for i=1:length(dx)

Central_difference_error(i)=Central_difference_approximation_of_fourth_order(x,dx(i),X1,analytical_derivative);
Skewed_right_sided_error(i)=Skewed_right_sided_approximation_of_fourth_order(x,dx(i),X2,analytical_derivative);
Skewed_left_sided_error(i)=skewed_left_sided_approximation_of_fourth_order(x,dx(i),X3,analytical_derivative);

end

% plotting

loglog(dx,Central_difference_error,'r-+');
hold on 
loglog(dx,Skewed_right_sided_error,'k-o');
hold on 
loglog(dx,Skewed_left_sided_error,'m','marker','s','LineStyle',':')
hold off
xlabel('dx')
ylabel('trunction error')
title('Range of dx Vs error')
legend('Central difference error of 4th order','Skewed right sided error 4th order','Skewed left sided error 4th order')

Function Programs:-

Function for Central difference Error:-

%Function to evaluate  central difference error
function Central_difference_error = Central_difference_approximation_of_fourth_order(x,dx,X1,analytical_derivative)
exp_cos_x=@(a)(exp(a))*(cos(a));  %creating a function to minimize length of Code.

N1=(X1(1)*(exp_cos_x(x-(2*dx)))+X1(2)*(exp_cos_x(x-dx))+X1(3)*(exp_cos_x(x))+X1(4)*(exp_cos_x(x+dx))+X1(5)*(exp_cos_x(x+(2*dx))));

Central_difference_approximation_of_fourth_order=N1/(dx^2);

% analytical_derivative called from main program
Central_difference_error=abs(max(analytical_derivative-Central_difference_approximation_of_fourth_order));
end

Function for Skewed Right sided difference Error:-

function Skewed_right_sided_error = Skewed_right_sided_approximation_of_fourth_order(x,dx,X2,analytical_derivative)

exp_cos_x=@(a)(exp(a))*(cos(a));      %creating a function to minimize length of Code.

N2=(X2(1)*(exp_cos_x(x))+X2(2)*(exp_cos_x(x+dx))+X2(3)*(exp_cos_x(x+(2*dx)))+X2(4)*(exp_cos_x(x+(3*dx)))+X2(5)*(exp_cos_x(x+(4*dx)))+X2(6)*(exp_cos_x(x+(5*dx))));

Skewed_right_sided_approximation_of_fourth_order=N2/(dx^2);

% analytical_derivative called from main program
Skewed_right_sided_error=abs(analytical_derivative-Skewed_right_sided_approximation_of_fourth_order);

end

%Function to evaluate skewed Left sided error
function Skewed_left_sided_error = skewed_left_sided_approximation_of_fourth_order(x,dx,X3,analytical_derivative)

exp_cos_x=@(a)(exp(a))*(cos(a));  %creating a function to minimize length of Code.
N3=(X3(1)*(exp_cos_x(x))+X3(2)*(exp_cos_x(x-dx))+X3(3)*(exp_cos_x(x-(2*dx)))+X3(4)*(exp_cos_x(x-(3*dx)))+X3(5)*(exp_cos_x(x-(4*dx)))+X3(6)*(exp_cos_x(x-(5*dx))));

Skewed_left_sided_approximation_of_fourth_order=N3/(dx^2);

% analytical_derivative called from main program
Skewed_left_sided_error=abs(analytical_derivative-Skewed_left_sided_approximation_of_fourth_order);

end