Week 3.5 - Deriving 4th order approximation of a 2nd order derivative using Taylor Table method
AIM and OBJECTIVE:
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.
Write a program in Matlab to evaluate the second-order derivative of the analytical function exp(x)*cos(x) and compare it with the 3 numerical approximations that you have derived.
Provide a plot that compares the absolute error between the above-mentioned schemes
In your own words, describe why a skewed scheme is useful? What can a skewed scheme do that a CD scheme cannot?
CENTRAL DIFFERENCING SCHEME:
NO OF NODES:
The number of nodes = p+q-1
Where, p is order of derivative.
q is order of approximations.
The no of nodes = 2+4-1 =5
These 5 nodes is represented by the below.
Here, we’re computing 4th order approximation at point (i), where we take information from the four-point surrounding (i); i.e. we’re going to take two points on the left of (i) and two points on the right of (i). Let those points be (i-2), (i-1), (i+1), and (i+2) and the distance between each of them is (∆x). Our numerical stencil will look like as shown in the above figure.
Based on the above stencil our hypothesis equation will be:-
Where a, b, c, d, e are coefficients of approximations.
Let ‘dx’ be the discretization size , then from Taylor series expansion of a,b,c,d and e terms.
All the above expansions can be represented in tabular form called Taylor’s table as shown below:
Since only the second derivative of the function is to be found out, the sum of the coefficients of f’’(i)(∆x)2" alone equal to 1 and rest all are zero. We require only 5 equations to find the 5 coefficients.
Therefore considering the coefficients of the five columns f(i)-fIV(i)(∆x)4we have
Converting the above equations, into matrix form:
Solving the matrix by in MATLAB leads to:
a=-0.0833
b=1.3333
c=-2.5000
d=1.3333
e=-0.0833
After putting the values of all the coefficients (a,b,c,d,e) in the hypothesis equation, we will get the fourth-order central difference approximation for the second-order derivative.
SKEWED RIGHT HAND DIFFERENCE SCHEME:
NO OF NODES:
The number of nodes = p+q
Where, p is order of derivative.
q is order of approximations.
The no of nodes = 2+4 =6.
These 6 nodes is represented by the below.
To derive 4th order skewed right-sided difference at point (i) of function. Since we’re using the skewed right-sided method, here we have to define five points on the right of (i). Let those points be (i+1), (i+2), (i+3), (i+4), and (i+5); with a spacing of ∆x each. Our numerical stencil will look like as shown in the above figure.
Based on the above stencil our hypothesis equation will be:-
Where f, g, h, i, j, k are coefficients of approximations.
Let ‘dx’ be the discretization size , then from Taylor series expansion of f, g, h, i, j and k terms.
All the above expansions can be represented in tabular form called Taylor’s table as shown below:
Since only the second derivative of the function is to be found out, the sum of the coefficients of f’’(i)(∆x)2" alone equal to 1 and rest all are zero. We require only 6 equations to find the 6 coefficients.
Therefore considering the coefficients of the five columns f(i)-fV(i)(∆x)5we have
Converting the above equations, into matrix form:
Solving the matrix by in MATLAB leads to:
f=3.7500
g=-12.8333
h=17.8333
i=-13.0000
j=5.0833
k=-0.8333
After putting the values of all the coefficients (f,g,h,i,j,k) in the hypothesis equation, we will get the fourth-order central difference approximation for the second-order derivative.
SKEWED LEFT HAND DIFFERENCE SCHEME:
NO OF NODES:
The number of nodes = p+q
Where, p is order of derivative.
q is order of approximations.
The no of nodes = 2+4 =6.
These 6 nodes is represented by the below.
To derive 4th order skewed right-sided difference at point (i) of function. Since we’re using the skewed right-sided method, here we have to define five points on the left of (i). Let those points be (i-1), (i-2), (i-3), (i-4), and (i-5); with a spacing of ∆x each. Our numerical stencil will look like as shown in the above figure.
Based on the above stencil our hypothesis equation will be:-
Where l, m, n, o, p, q are coefficients of approximations.
Let ‘dx’ be the discretization size , then from Taylor series expansion of l, m, n, o, p and q terms.
All the above expansions can be represented in tabular form called Taylor’s table as shown below:
Since only the second derivative of the function is to be found out, the sum of the coefficients of f’’(i)(∆x)2" alone equal to 1 and rest all are zero. We require only 6 equations to find the 6 coefficients.
Therefore considering the coefficients of the five columns f(i)-fV(i)(∆x)5we have
Converting the above equations, into matrix form:
Solving the matrix by in MATLAB leads to:
l=-0.8333
m=5.0833
n=-13.000
o=17.8333
p=-12.8333
q=3.7500
After putting the values of all the coefficients (l,m,n,o,p,q) in the hypothesis equation, we will get the fourth-order central difference approximation for the second-order derivative.
Analytical derivation of the given function equation.
MATLAB CODE:
After analyzing the graph between ‘Truncation error’ vs. ‘range of (dx)' for central, skewed right-sided, and skewed left-sided difference 4th order error, we can easily notice that central difference approximation produces more accurate solutions, hence produces lesser errors as compared to right-skewed and left-skewed difference schemes.
Since we can see that the central differencing scheme gives us more accurate solutions or lesser errors as compared to the skewed difference scheme. But it is not always viable to use the central differencing to evaluate solutions. Suppose, if we want to compute numerical approximations for a function at the boundary or corner node, we cannot use the central difference scheme at this particular situation because it requires data from both sides of the node or boundary. So, in this situation skewed difference method plays a major role. In the skewed scheme data from one side can be used to compute numerical approximation at a point.
close all
clear all
clc
% Week 3.5 - Deriving 4th order approximation of a 2nd order derivative using Taylor Table method
%Finding the coefficients of central difference scheme taken from taylor's table.
% Matrix A definng
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];
% MAtrix B defining
B=[0;0;1;0;0];
% Finding the coefficients
cds=A\B;
%Finding the coefficients of right skewed differnce scheme taken from taylor's table.
% Matrix C definng
C= [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];
% Matrix D definng
D=[0;0;1;0;0;0];
% Finding the coefficients
skrs=C\D;
%Finding the coefficients of left skewed differnce scheme taken from taylor's table.
% Matrix E definng
E= [1 1 1 1 1 1
-5 -4 -3 -2 -1 0
25/2 16/2 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];
% Matrix F definng
F=[0;0;1;0;0;0];
% Finding the coefficients
skls=E\F;
% Function is given by f(x)=exp(x)*cos(x);
x=pi/3;
dx=linspace(pi/4,pi/4000,20);
% second order Analytical Derivative
Analytical_derivative= -2*exp(x)*sin(x);
for i = 1:length(dx)
% Numerical Derivative is as below
% Central difference scheme equation Δx^2 (ⅆ^2 u)/(ⅆx^2 )=af(ⅈ-2)+bf(ⅈ-1)+cf(ⅈ)+df(ⅈ+1)+ⅇf(ⅈ+2)
% function coefficent values are a=-0.0833 , b=1.3333 ,c=-2.5000, d=1.3333, e=-0.0833
CDS_approximation(i)= (cds(1)*exp(x-(2*dx(i)))*cos(x-(2*dx(i)))+cds(2)*exp(x-(dx(i)))*cos(x-(dx(i)))+cds(3)*exp(x)*cos(x)+cds(4)*exp(x+(dx(i)))*cos(x+(dx(i)))+cds(5)*exp(x+(2*dx(i)))*cos(x+(2*dx(i))))/(dx(i)^2);
CDS_approximation_error(i)=abs(CDS_approximation(i)-Analytical_derivative);
% Skewed right-sided scheme equation Δx^2 (ⅆ^2 u)/(ⅆx^2 )=ff(ⅈ)+gf(ⅈ+1)+hf(ⅈ+2)+if(ⅈ+3)+jf(ⅈ+4)+kf(ⅈ+5)
% function coefficent values are f=3.7500 , g=-12.8333 ,h=17.8333,i=-13.0000 , j=5.0833 , k=-0.8333
Skewed_right_approximation(i)= (skrs(1)*exp(x)*cos(x)+skrs(2)*exp(x+(dx(i)))*cos(x+(dx(i)))+skrs(3)*exp(x+(2*dx(i)))*cos(x+(2*dx(i)))+skrs(4)*exp(x+(3*dx(i)))*cos(x+(3*dx(i)))+skrs(5)*exp(x+(4*dx(i)))*cos(x+(4*dx(i)))+skrs(6)*exp(x+(5*dx(i)))*cos(x+(5*dx(i))))/(dx(i)^2);
Skewed_right_approximation_error(i)=abs(Skewed_right_approximation(i)-Analytical_derivative);
% Skewed left-sided scheme equation Δx^2 (ⅆ^2 u)/(ⅆx^2 )=lf(ⅈ-5)+mf(ⅈ-4)+nf(ⅈ-3)+of(ⅈ-2)+pf(ⅈ-1)+qf(ⅈ)
% function coefficent values are l=-0.8333 , m=5.0833 ,n=-13.0000,o=17.8333 , p=-12.8333 , q=3.75000
Skewed_left_approximation(i)= (skls(1)*exp(x-(5*dx(i)))*cos(x-(5*dx(i)))+skls(2)*exp(x-(4*dx(i)))*cos(x-(4*dx(i)))+skls(3)*exp(x-(3*dx(i)))*cos(x-(3*dx(i)))+skls(4)*exp(x-(2*dx(i)))*cos(x-(2*dx(i)))+skls(5)*exp(x-dx(i))*cos(x-dx(i))+skls(6)*exp(x)*cos(x))/(dx(i)^2);
Skewed_left_approximation_error(i)=abs(Skewed_left_approximation(i)-Analytical_derivative);
end
% plotting the error
figure(1)
loglog(dx,CDS_approximation_error,'b-*')
hold on
loglog(dx,Skewed_right_approximation_error,'r-*')
hold on
loglog(dx,Skewed_left_approximation_error,'g-*')
grid on
legend('CDS_approximation_error','Skewed_right_approximation_error','Skewed_left_approximation_error')
xlabel('range of (dx)')
ylabel('Truncation error')
title('Error for different approximation schemes')