Week 3.5 - Deriving 4th order approximation of a 2nd order derivative using Taylor Table method
AIM:- Deriving 4th order approximation of a 2nd order derivative using the Taylor Table method
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 i.e the thumb rule given in
HINT (ii).:-
Once you have done this, 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?
HINT :
i) The order of accuracy P = Number of nodes in the stencil N - order of derivative Q + 1 for symmetric/central schemes
ii) The order of accuracy P = Number of nodes in the stencil N - order of derivative Q for one-sided/skewed schemes
iii) Solve the system of equations using the matrix inversion method in Matlab to obtain the Taylor table coefficients and use them directly in the numerical approximation statements to avoid input errors.
introduction:-
The differential equation is an equation.
which contains one or more terms that are in value derivative of the dependent variable.
An ordinary differential equation is one independent variable equation.
Q = -kAdT/dx.
dT = depend varible
dx = independent varible
ODE = Initial value problem
Partial diffrention Equation
Partial differential equations are boundary value problems.
Discretization:-
in mathematics, discretization concerns the process of transferring continuous models and equations into discrete computer parts
This process is usually carried out as a first step toward making them suitable for numerical evaluation and implementation of a digital computer.
Different scheme:-
∂f/∂x, C. D. S . F. D. S. B. D. S.
Central Differencing Forward Differencing Backward Differencing
Scheme Scheme Scheme
Concept of stencil:-
In mathematics especially in the areas of numerical solutions of the partial differencing equation
Stencil is a generic arrangement of a node group that relates to the point of interest by using a numerical approximation routine.
Central differencing:- Here the nodes selected are symmetrical i.e. the data is taken from both sides of the point of focus.
Taylor's method, we get
∂f/∂x= af(i - 2) + b(i - 1) + cf(i) + df(i +2)
Here there are nodes taken from both sides of the point of focus i.e.(i).
(NOTE):- To select the number of nodes in consideration, we use the following expression:-
N = P+q -1 for central differencing
N = P + q for skewed differencing
Here N is no. of nodes, p is the order of derivative and q order of approximation.)
Taylor table for central difference:-
f(x) dx dx^2 dx^3 dx^4 d^5
f'(x) f''(x) f'''(x) f''''(x) f'''''(x)
af(i-1) a -2a 4a/2 -8a/6 16a/24 -32a/120
b(i-2) 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 d/120
ef(i-2) e 2e 4e/2 8e/6 16e/24 32e/120
af(i-2) + bf(i-1)
+cf(i) +df(i+1)
+ef(i+2) 0 0 1 0 0 ?
it look at hypothetic term 1 to 5 term
This tabel transform into matrix of type Ax = B
a + b + c + d + e = 0 [1 1 1 1 1] [a] [0]
-2a - b + 0 + d + 2e = 0 [-2 -1, 0 1 2] [b] [0]
2a + b/2 + 0 + d/2 + 2e = 1 [ 2 1/2 0 1/2 2] [c] [0]
-8a/6 -b/6 + 0 + d/6 + 8e/6 = 0 [ 8/6 -1/6 0 1/6 8/6] [d] [0]
16/24 + b/24 + 0 + d/24 + 16e/24 =0 [ 16/24 1/24 0 1/24 16/24] [e] [0]
get the value of the unknown
i.e. values = a = -0.0833 b = 1.3333 c = 2.5000
d = 1.3333 e = 0.0833
These coefficients will further in code to get the value of derivative.
Skewed right sided difference:-
when the data point of the lift nodes are not aviable (at intial boundary conditions) .
we have to rely on the aviable data set
∂^2f/∂x^2 = af(i) + bf(i+1) +cf(1+i) + df(i+3) + ef(i+4) + ff(i+5)
Taylor's table for skewed right-sided differences
f(x) dxf'(x) dx^2f"(x) dx^3f'''(x) dx^4f''''(x) dx^5f'''''(x)
af(i) a 0 0 0 0 0
bf(i+1) b b b/2 b/6 b/24 b/120
cf(i+2) c 2c 4c/2 8c/6 16c/24 32c/120
df(i+3) d 3d 9d/2 27d/6 81d/24 243d/120
ef(i+4) e 4e 16e/2 64e/6 256e/24 1024e/120
ff(i+5) f 5f 25f/2 125f/6 625f/24 3125f/120
af(i)+bf(i+1)+
cf(i+2)+ df(i+3)
ef(i+4)+ff(i+5)
o 0 1 0 0 0
Tjis table transfer into matrix of the AX = B .
1 1 1 1 1 1 a 0
0 1 2 3 4 5 b 0
0 1/2 2 9/2 16/2 25/2 * c = 1
0 1/6 8/6 27/6 64/6 125/6 d 0
0 1/24 16/24 81/24 256/24 625/24 e 0
0 1/120 32/120 243/120 1024/120 3125/120 0
on further solving we get the value of the unkown values = a = 3.75000
b = -12.8333 c = 17.8333 d = 13.0000 e = 5.0833 f = 0.8333
These coefficients will futher used in code to set the value of derivative
Sekwed left sided difference.
When the data point of the left nodes is not available (at boundary conditions then we have to rely on the available data set).
∂^2f/∂x^2 = af(i - 5) + bf(i-4) + cf(i -3) + df(i-2) + ef(i-1)+ ff(i)
Taylor table for skewed right-sided difference
f(x) dxf(x) dx^2f''(x) dx^3f'''(x) dx^4f''''(x) dx^5f'''''(x)
af(i - 5) a -5a 25a/2 -125a/6 625a/24 3125a/120
bf(i-4) b -4b 16b/2 -64b/6 256b/24 -1024b/120
cf(i-3) c -3c 9c/2 -27c/6 81c/24 -243c/120
df(i-2) d -2d 4d/2 -8d/6 16d/24 -32d/120
ef(i-1) e -e e/2 e/6 e/24 -e/120
ff(i) f 0 0 0 0 0
af(i-5)+bf(i-4)+
cf(i-3) + df(i-2)+
ef(i-1)+ff(i)
0 0 1 0 0 0
This tabel tansform into matrix of type AX = B.
[ 1 1 1 1 1 1] a 0
[ -5 -4 -3 -2 -1 0] b 0
[ 25/2 16/2 9/2 2 1/2 0] * c = 1
[ -125/6 -64/6 -27/6 -8/6 -1/6 0] d 0
[ 625/24 256/24 81/24 16/24 7/24 0] e 0
[3121/120 -1o24/120 -243/120 32/120 1/120 0] f 0
unknown i.e. values = a= 0.8333 b = 5.08333 c = 13.0000 d = 17.8333 e = -12.8333
f = 3.7500
These coefficients will further use incode to get the value of derivative
step involved:-
1 we design 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 the above-discussed section and solved it using the inversion method
i.e.
Ax = B = = X = A-1 B
And stored the coefficients in cf1 cf2 and cf3
4. Assigned the rows and columns 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 anal6ytically
6. The calculated values are plotted to do a better comparison
Matlab_code:-
clear all
close all
%analytical function defination and value assigment
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.
%coefficient 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];
pc = inv(p)*c;
%as p*a = c i.e. a inverse (p)*c('fastest way to set by')
cf1 = pc;
%for right skewed difference 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/4 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];
qD = inv(q)*D;
%as q*b = s i.e. b = inverse(q)*s (fastest way to get by")
cf2 = qD;
%for left skewed different 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;];
RE = inv(R)*E;
%as n*b = s i.e. b = inverse(b)*s(fastest way to set by ")
cf3 = RE;
%calculated values of derivatives using C.D.S. R.S.D.C. L.S.D.C.
%defining variable size before 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(x-dx(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_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+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);
%calculate value and error by left skewed difference skewed.
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(x-dx(i))+cf3(6)*f(i))/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,central_diff,'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 calcuated by diff. order of approximation')
hold off
%plotting errors in approximatiion of derivative
subplot(2,1,2)
loglog(dx,err_central,'linewidth',2,'color','b');
hold on
xlim([.002 0.4]);
loglog(dx,err_right,'linewidth',2,'color','b');
loglog(dx,err_left,'linewidth',2,'color','g');
legend('cental diff.','right skewed diff.','left skewed diff.','location','best')
xlabel('grid siza(dx)');
ylabel('Absoulte error');
title('comparission b/w 4th order approximation method for error in calculating d2f/dx2');
Results:-
1 in the first plot we can see that the value of all the approximations approaches toward the analytically calculated value as the size of dx gets smaller.
2. on the second plot we can see that the error in the coefficient difference is lowest with respect to the other two