Week 3.5 - Deriving 4th order approximation of a 2nd order derivative using Taylor Table method
Objective of this project is to derive fourth order approximation for second order derivative using symmetric and un-symmetric schemes with Taylor's table method and compare the analytical solution with obtained approximation values.
Derivations
1) Central Difference scheme
No. of nodes = p+q-1 = 4+2-1 = 5
p = order of derivative
q = order of approximation
Now, by Taylor's series epansion for all the terms of the equation f(x-2dx), f(x-dx) and all others, we get below tabel which is called Taylor's table
From the above, we get linear equation which will be solved by matrix inversion method
We get coefficints as
a = -0.0833, b= 1.333, c= -2.5, d = 1.333, e = -0.0833
2) Skewed Right hand difference scheme
No. of nodes = p+q = 4+2 = 6
p = order of derivative
q = order of approximation
Derivation: To derive fourth order approximation of second order derivative, we start with below approximation,
Expanding all terms by Taylor series expansionwe get below Taylor's table
From the above, we get linear equation which will be solved by matrix inversion method
We get coefficints as
a = 3.7500, b= -12.8333, c= 17.833, d = -13, e = 5.0833, g= -0.8333
2) Skewed Left hand difference scheme
No. of nodes = p+q = 4+2 = 6
p = order of derivative
q = order of approximation
Derivation: To derive fourth order approximation of second order derivative, we start with below approximation,
Expanding all terms by Taylor series expansionwe get below Taylor's table
From the above, we get linear equation which will be solved by matrix inversion method
We get coefficints as
a = -0.8333, b= 5.0833, c= -13.000, d = 17.8333, e = -12.833, g= 3.75
Matlab code for the same:
clear all
close all
clc
x = pi/3;
dx = linspace(pi/4, pi/400, 200);
% Central Difference scheme
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= inv(A)*B;
% Right handed skew scheme
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] ;
Y = inv(A2)* B2
% Left handed skew scheme
A3 = [1 1 1 1 1 1; -5 -4 -3 -2 -1 0; (25/2) 8 (9/2) 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];
B3 = [0;0;1;0;0;0] ;
Z = inv(A3)* B3
% Analytical derivative
% f''(x) = -2*exp(x)*sin(x)
analytical_derivative = (-2)*exp(x)*sin(x);
for i = 1:length(dx)
central_differencing(i)= fourth_order_central_differencing(x,dx(i),X,analytical_derivative);
right_sided_skew(i) = fourth_order_right_sided_skew(x,dx(i),Y,analytical_derivative);
left_sided_skew(i) = fourth_order_left_sided_skew(x,dx(i),Z,analytical_derivative);
end
% Graphs
figure(1)
loglog(dx,central_differencing, 'r')
hold on
loglog(dx, right_sided_skew, 'g')
hold on
plot(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');
% Code for functions
function error_central_differencing = fourth_order_central_differencing(x,dx,X,analytical_derivative)
% analytical function = exp(x)* cos(x)
f = @(x) exp(x)*cos(x);
% Derivative value numerically
central_differencing = (X(1)*f(x-2*dx) + X(2)*f(x-dx) + X(3)*f(x) + X(4)*f(x+dx) + X(5)*f(x+2*dx))/(dx^2);
% Error
error_central_differencing = abs(analytical_derivative -central_differencing );
end
function error_right_sided_skew = fourth_order_right_sided_skew (x,dx,Y,analytical_derivative)
% analytical function = exp(x)* cos(x)
f = @(x) exp(x)*cos(x);
% Derivative value numerically
% 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.000)*((exp(x+(3*dx)))*(cos(x+(3*dx)))))+((5.0833)*((exp(x+(4*dx)))*(cos(x+(4*dx)))))-((0.8333)*((exp(x+(5*dx))))))/(dx.^2);
right_sided_skew = (Y(1)* f(x) + Y(2)*f(x+dx)+ Y(3)*f(x+2*dx)+ Y(4)*f(x+3*dx)+ Y(5)*f(x+4*dx)+ Y(6)*f(x+5*dx))/(dx^2);
% Error
error_right_sided_skew = abs(analytical_derivative -right_sided_skew);
end
function error_left_sided_skew = fourth_order_left_sided_skew (x,dx,Z,analytical_derivative)
% analytical function = exp(x)* cos(x)
f = @(x) exp(x)*cos(x);
% Derivative value numerically
% left_sided_skew = (((3.75000)*((exp(x-(5*dx)))*(cos(x-(5*dx)))))-((12.8333)*((exp(x-(4*dx)))*(cos(x-(4*dx)))))+((17.8333)*((exp(x-(3*dx)))*(cos(x-(3*dx)))))-((13.000)*((exp(x-(2*dx)))*(cos(x-(2*dx)))))+((5.0833)*((exp(x-dx))*(cos(x-dx))))-((0.8333)*((exp(x))*(cos(x)))))/(dx^2);
left_sided_skew = (Z(1)*f(x-5*dx) + Z(2)* f(x-4*dx) + Z(3)* f(x-3*dx) + Z(4)* f(x-2*dx) + Z(5)* f(x- dx) + Z(6)* f(x))/(dx^2);
% Error
error_left_sided_skew = abs(analytical_derivative - left_sided_skew );
end Result :
Conclusion : Central differencing scheme has very less error compared to other schemes. Hence, it is better to use central differencing scheme because it uses information from both the sides of the point of consideration to give the approximate and accuratye results.
As left and right sided skew methods uses the data only from one side of point of interest it gives less accurate results.
Skewed schemes are used when we are unknown about the data from both the sides of a point. Because we can not use central scheme at that situation.
Difference between central and skewed schemes :
CDS can not be used at the boundary nodes due to unavailability of right, left nodes at the same time. In these cases, skewed schemes will be employed .