Week 3.5 - Deriving 4th order approximation of a 2nd order derivative using Taylor Table method
AIM:
To derive the fourth order approximations of a second order derivative using central differencing scheme, skewed right sided difference and skewed leftsided difference with the help of taylor table method and to compare the analytical derivative with obtained values for the given function.
SCHEME DERIVATIONS:
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.
DERIVATION
To derive a fourth order approximation for the second order derivative,we start with these approximation
(NOTE: - To select the number of nodes in consideration, we use following expression:
N=p + q -1 for central differencing
N = p + q for skewed differencing
Here N is no. of nodes, p is order of derivative and q is order of approximation.)
These coefficients will further used in code to get the value of derivative.
1. Skewed right sided difference: - When the data point of the left nodes are not available (at initial boundary conditions) then we have to rely on the available data set.
These coefficients will further used in code to get the value of derivative.
1. Skewed left sided difference: - When the data point of the right nodes are not available (at boundary conditions) then we have to rely on the available data set.
These coefficients will further used in code to get the value of derivative.
MATLAB CODE:
clear all
close all
clc
%analytical function defination and value assignment
x=pi/3;
dx=linspace(pi/10,pi/1000,100);
f=@(x)(exp(x)*cos(x));
g=@(x)(-2*exp(x)*sin(x));
d2f=g(x);
func=f(x);
%for central difference scheme
%coefficients taken from taylor's table
P = [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];
C=[0;0;1;0;0];
%as P*a=C i.e. a=inverse(P)*c(fastest way to getby '')
cf1=PC;
%for right skewwd differnce scheme
Q= [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];
D=[0;0;1;0;0;0];
%as Q*b=s i.e. b=inverse(Q)*s(fastest way to getby '')
cf2=QD;
%for right skewwd differnce scheme
R= [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];
E=[0;0;1;0;0;0];
%as n*b=s i.e. b=inverse(n)*s(fastest way to getby '')
cf3=RE;
%calculating values of derivative using CDS RSDC,LSDC
%defining variable size befor using in loop
central_diff=zeros(1,100);
err_central=zeros(1,100);
right_sdiff=zeros(1,100);
err_right=zeros(1,100);
left_sdiff=zeros(1,100);
err_left=zeros(1,100);
for i=1:length(dx)
%calculating value and error by central difference scheme
central_diff(i)=(cf1(1)*f(x-2*dx(i))+cf1(2)*f(xdx(
i))+cf1(3)*f(x)+cf1(4)*f(x+dx(i))+cf1(5)*f(x+2*dx(i)))/(dx(i))^2;
err_central(i)=abs(central_diff(i)-d2f);
%calculating value and error by right skewed difference scheme
right_sdiff(i)=
(cf2(1)*f(x)+cf2(2)*f(x+dx(i))+cf2(3)*f(x+2*dx(i))+cf2(4)*f(x+3*dx(i))+cf2(5)*f(x+4*dx(i))+cf2(6)*f(x+5*dx(i)))/(dx(i))^2;err_right(i)=abs(right_sdiff(i)-d2f);
%calculating value and error by left skewed difference scheme
left_sdiff(i)=(cf3(1)*f(x-5*dx(i))+cf3(2)*f(x-4*dx(i))+cf3(3)*f(x-3*dx(i))+cf3(4)*f(x-2*dx(i))+cf3(5)*f(xdx(
i))+cf3(6)*f(x))/(dx(i))^2;
err_left(i)=abs(left_sdiff(i)-d2f);
end
%plotting values of derivative
subplot(2,1,1)
hold on
fplot(d2f,'--')
xlim([-0.05 0.35]);
plot(dx,central_diff,'linewidth',2,'color','r');
plot(dx,right_sdiff,'linewidth',2,'color','b');
plot(dx,left_sdiff,'linewidth',2,'color','g');
legend('analytical value at x=pi/3','central diff.','right skewed diff.','left skewed
diff.','location','best')
xlabel('Grid siza(dx)');
ylabel('Value');
title('Comparision bw value of derivative calculated by diff. order of approximation')
hold off
%plotting errors in approximation of derivative
subplot(2,1,2)
loglog(dx,err_central,'linewidth',2,'color','r');
hold on
xlim([.002 0.4]);
loglog(dx,err_right,'linewidth',2,'color','b');
loglog(dx,err_left,'linewidth',2,'color','g');
legend('central diff.','right skewed diff.','left skewed diff.','location','best')
xlabel('Grid siza(dx)');
ylabel('Absolute error');
title('Comparision b/w 4th order approximation method for error in calculating d2f/dx2 ');steps involved: -
1. We assign the value of x and range of dx
2. We define an anonymous function of analytical function and its derivative.
3. Use the matrices from above discussed section and solved it using inversion method i.e.
AX=B ==> X=A B
And stored the coefficients in cf1, cf2 and cf3.
4. Assigned the rows and column to the appropriate matrices.
5. Loop the order of approximation expression in a for loop under the range of dx . Here absolute error is also calculated from the derivative value calculated analytically.
6. The calculated valued are plotted to do a better comparision.
Result:
1. In first plot we can see that the value of all the approximation approached toward the analytically value as the size of dx gets smaller.
2. On the second plot we can see that the error in the central difference is lowest with respect to other two.