Computing higher order approximation by taylor table method for a second order derivative of fourth order approximate using MATLAB
Objective:
- Derive the 4th order approximation of second order derivative using symmetric stencil (central differencing method) and skewed stencil(skewed right sided and skewed left sided differencing method) and prove that the skewed schemes are fourth order approximate.
- Write a matlab code to evaluate the second order derivative for the analytical function f(x)=exp(x)*cos(x) and compare it with central difference scheme, skewed right sided scheme, skewed left sided scheme.
- Plot an error curve between the analytical and the above three numerical schemes.
Theory
The fourth order approximate scheme is calculated by the formula n=p+q, where n is the no of nodes in stencils, p is order of derivative and q is order of approximation
Central differencing scheme
p=2, q=4, n=6
the hypothesis equation can be written as
`(δ^2f)/(δx^2)`= a f(i-2) +b f(i-1)+c f(i)+df(i+1)+e f(i+2)
Writing it in taylor table format with each term expanded in terms of taylor series.
Algebraic equation
a+b+c+d+e=0
-2a-b+0+d+2e=0
2a+b/2+0+d/2+2e=1
-8a/6-b/6+0+d/6+8e/6=0
16a/24+b/24+0+d/24+16e/24=0
Matrix formulation of form AX=B
Evaluate the value of constant
a=-0.0833
b=1.3333
c=-2.5000
d=1.3333
e=-0.0833
`(δ^2f)/(δx^2)`= -0.0833 f(i-2) +1.333 f(i-1)-2.500 f(i)+1.333f(i+1)-0.0833 f(i+2)
Skewed right sided
p=2, q=4, n=6
the hypothesis equation can be written as
`(δ^2f)/(δx^2)`= a f(i)+b f(i+1)+c f(i+2)+df(i+3)+e f(i+4)+g f(i+5)
Writing it in taylor table format with each term expanded in terms of taylor series.
Algebraic equation
a+b+c+d+e+g=0
0+b+2c+3d+4e+5g=0
0+b/2+2c+9d/2+8e+25g/2=1
0+b/6+8c/6+27d/6+64e/6+125g/6=0
0+b/24+16c/24+81d/24+256e/24+625g/24=0
0+b/120+32c/120+243d/120+1024e/120+3125g/120=0
Matrix formulation of form AX=B
Evaluate the value of constant
a=3.7500
b=-12.8333
c=17.8333
d= -13.0000
e=5.0833
g=-0.8333
(δ^2f)/(δx^2)=3.750f(i)-12.833f(i+1)+17.833f(i+2)-13.00f(i+3)+5.0833f(i+4)-0.833f(i+5)
Now let us prove that it is fourth order approximate by expanding each term in taylor series form
Now we know that coefficients of term f(i),f’(i),f’’’(i),fiv(i), fv(i) are all zero except for the term f’’(i) which is equal to 1
Now dividing both sides by x2 we get
(δ^2f)/(δx^2)=f’’(i)+O( x4)
So the truncation error is of order 4 and the scheme is fourth order approximate.
Skewed left sided
p=2, q=4, n=6
the hypothesis equation can be written as
`(δ^2f)/(δx^2)` = a f(i)+b f(i-1)+c f(i-2)+df(i-3)+e f(i-4)+g f(i-5)
Writing it in taylor table format with each term expanded in terms of taylor series.
Algebraic equation
a+b+c+d+e+g=0
0-b-2c-3d-4e-5g=0
0+b/2+2c+9d/2+8e+25g/2=1
0-b/6-8c/6-27d/6-64e/6-125g/6=0
0+b/24+16c/24+81d/24+256e/24+625g/24=0
0-b/120-32c/120-243d/120-1024e/120-3125g/120=0
Matrix formulation of form AX=B
Evaluate the value of constant
a= 3.7500
b=-12.8333
c=17.8333
d=-13.0000
e=5.0833
g=-0.8333
(δ^2f)/(δx^2`=3.750 f(i) -12.833f(i-1)+17.833 f(i-2)-13.00 f(i-3)+5.0833 f(i-4)-0.833f(i-5)
Now let us prove that it is fourth order approximate by expanding each term in taylor series form
Now we know that coefficients of term f(i),f’(i),f’’’(i),fiv(i), fv(i) are all zero except for the term f’’(i) which is equal to 1
Now dividing both sides by x2 we get
(δ^2f)/(δx^2)=f’’(i)+O( x4)
So the truncation error is of order 4 and the scheme is fourth order approximate.
Procedure
- Define the location of x where the analytical and numerical second order derivative is to be calculated.
- Define the no of grid points for a given range of dx(pi/4-pi/400).
- Compute the values of coefficients for central differencing, skewed right and skewed left sided difference scheme by performing inverse matrix operation.
- Then compute the numerical derivative for each of the scheme separately and calculate the error between the analytical derivative and three numerical derivative schemes.
- Perform the entire calculation stated above in three different function block named ‘central_differencing_approximation’, ’skewed_rightsided_differencing_approximation’, ’skewed_leftsided_differencing_approximation’
- In the next step the main code stores the value of errors generated from each of the function block separately.
- A plot between the error and the order of approximation is plotted to observe the rate of accuracy of the given 3 numerical schemes.
The program code for the above procedure is given below.
MATLAB CODE
Main program code
% Matlab main code for taylor series approximation of second order derivative by
% central differencing and skewed differencing
%%%%%%%%%%%%%%%%%%%%
clc
clear all
close all
% location of the point
x=pi/3;
% defining the no of grid points
dx=linspace(pi/10,pi/400000,100);
for i=1:length(dx)
% calling the central differencing approximation function
central_differencing_error(i)=central_difference_approximation(x,dx(i));
% calling the skewed rightsided differencing approximation function
skewed_rightsided_error(i)=skewed_rightsided_difference_approximation(x,dx(i));
% calling the skewed leftsided differencing approximation function
skewed_leftsided_error(i)=skewed_leftsided_difference_approximation(x,dx(i));
end
figure(1)
loglog(dx,central_differencing_error,\'-O\')
hold on
loglog(dx,skewed_rightsided_error,\'-<\')
hold on
loglog(dx,skewed_leftsided_error,\'-*\',\'color\',\'m\')
xlabel(\'grid size dx\')
ylabel(\'error with respect to order of approximation\')
title(\'second order derivative fourth order approximation\')
legend(\'central differencing error\',\'skewed rightsided error\',\'skewed leftsided error\')
function code for central difference scheme
% Matlab function code for central difference scheme for a second order
% derivative of fourth order finite difference approximation
%%%%%%%%%%%%%%%%%%%
function out= central_difference_approximation(x,dx)
%analytical function f(x)=exp(x)*cos(x)
%second order analytical derivative
analytical_derivative=-2*sin(x)*exp(x);
%numerical second order derivative by taylor series method
%central differencing scheme
%f\'\'(x)=af(i-2)+bf(i-1)+cf(i)+df(i+1)+ef(i+2)/dx^2
%computing the coefficients a,b,c,d,e from the linear equations given below
%a+b+c+d+e=0
%-2a-b+0+d+2e=0
%2a+b/2+0+d/2+2e=1
%-8a/6-b/6+0+d/6+8e/6=0
%16a/24+b/24+0+d/24+16e/24=0
%conversion of linear system to matrix system
%defining the coefficient matrix
A=[1 1 1 1 1;-2 -1 0 1 2;2 1/2 0 1/2 2;-8/6 -1/6 0 1/6 8/6;16/24 1/24 0 1/24 16/24];
%defining the constant matrix
B=[0;0;1;0;0];
z=inv(A)*B;
%storing the values of the constants
a=z(1);
b=z(2);
c=z(3);
d=z(4);
e=z(5);
% defining the numerical functions of central differencing scheme for fourth order approximate
function_i=exp(x)*cos(x);
function_i_minus_1=exp(x-dx)*cos(x-dx);
function_i_minus_2=exp(x-(2*dx))*cos(x-(2*dx));
function_i_plus_1=exp(x+dx)*cos(x+dx);
function_i_plus_2=exp(x+(2*dx))*cos(x+(2*dx));
%numerical central differencing scheme
central_difference=[(a*function_i_minus_2)+(b*function_i_minus_1)+(c*function_i)+(d*function_i_plus_1)+(e*function_i_plus_2)]/(dx^2);
%absolute error between analytical and central differencing scheme
out=abs(analytical_derivative-central_difference);
end
functiom code for right sided difference scheme
% Matlab function code for skewed rightsided difference scheme for a second order
% derivative of fourth order finite difference approximation
%%%%%%%%%%%%%%%%%%%
function out= skewed_rightsided_difference_approximation(x,dx)
%analytical function f(x)=exp(x)*cos(x)
%second order analytical derivative
analytical_derivative=-2*sin(x)*exp(x);
%numerical second order derivative by taylor series method
%skewed rightsided differencing scheme
%f\'\'(x)=af(i)+bf(i+1)+cf(i+2)+df(i+3)+ef(i+4)+gf(i+5)
%computing the coefficients a,b,c,d,e from the linear equations given below
%a+b+c+d+e+g=0
%0+b+2c+3d+4e+5g=0
%0+b/2+2c+9d/2+8e+25g/2=1
%0+b/6+8c/6+27d/6+64e/6+125g/6=0
%0+b/24+16c/24+81d/24+256e/24+625g/24=0
%0+b/120+32c/120+243d/120+1024e/120+3125g/120=0
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%conversion of linear system to matrix system
%defining the coefficient matrix
A=[1 1 1 1 1 1;0 1 2 3 4 5;0 1/2 2 9/2 8 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];
%defining the constant matrix
B=[0;0;1;0;0;0];
z=inv(A)*B;
%storing the values of the constants
a=z(1);
b=z(2);
c=z(3);
d=z(4);
e=z(5);
g=z(6);
% defining the numerical functions of skewed rightsided differencing scheme for fourth order approximate
function_i=exp(x)*cos(x);
function_i_plus_1=exp(x+dx)*cos(x+dx);
function_i_plus_2=exp(x+(2*dx))*cos(x+(2*dx));
function_i_plus_3=exp(x+(3*dx))*cos(x+(3*dx));
function_i_plus_4=exp(x+(4*dx))*cos(x+(4*dx));
function_i_plus_5=exp(x+(5*dx))*cos(x+(5*dx));
%numerical skewed rightsided differencing scheme
skewed_rightsided=[(a*function_i)+(b*function_i_plus_1)+(c*function_i_plus_2)+(d*function_i_plus_3)+(e*function_i_plus_4)+(g*function_i_plus_5)]/(dx^2);
%absolute error between analytical and central differencing scheme
out=abs(analytical_derivative-skewed_rightsided);
end
function code for left sided difference scheme
% Matlab function code for skewed leftsided difference scheme for a second order
% derivative of fourth order finite difference approximation
%%%%%%%%%%%%%%%%%%%
function out= skewed_leftsided_difference_approximation(x,dx)
%analytical function f(x)=exp(x)*cos(x)
%second order analytical derivative
analytical_derivative=-2*sin(x)*exp(x);
%numerical second order derivative by taylor series method
%skewed leftsided differencing scheme
%f\'\'(x)=af(i)+bf(i-1)+cf(i-2)+df(i-3)+ef(i-4)+gf(i-5)
%computing the coefficients a,b,c,d,e from the linear equations given below
%a+b+c+d+e+g=0
%0-b-2c-3d-4e-5g=0
%0+b/2+2c+9d/2+8e+25g/2=1
%0-b/6-8c/6-27d/6-64e/6-125g/6=0
%0+b/24+16c/24+81d/24+256e/24+625g/24=0
%0-b/120-32c/120-243d/120-1024e/120-3125g/120=0
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%conversion of linear system to matrix system
%defining the coefficient matrix
A=[1 1 1 1 1 1;0 -1 -2 -3 -4 -5;0 1/2 2 9/2 8 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];
%defining the constant matrix
B=[0;0;1;0;0;0];
z=inv(A)*B;
%storing the values of the constants
a=z(1);
b=z(2);
c=z(3);
d=z(4);
e=z(5);
g=z(6);
% defining the numerical functions of skewed leftsided differencing scheme for fourth order approximate
function_i=exp(x)*cos(x);
function_i_minus_1=exp(x-dx)*cos(x-dx);
function_i_minus_2=exp(x-(2*dx))*cos(x-(2*dx));
function_i_minus_3=exp(x-(3*dx))*cos(x-(3*dx));
function_i_minus_4=exp(x-(4*dx))*cos(x-(4*dx));
function_i_minus_5=exp(x-(5*dx))*cos(x-(5*dx));
%numerical skewed leftsided differencing scheme
skewed_rightsided=[(a*function_i)+(b*function_i_minus_1)+(c*function_i_minus_2)+(d*function_i_minus_3)+(e*function_i_minus_4)+(g*function_i_minus_5)]/(dx^2);
%absolute error between analytical and central differencing scheme
out=abs(analytical_derivative-skewed_rightsided);
end
Results
From the plot it can be inferred that error in cds is much lower compared to right and left skewed scheme. For a better approximate result cds is best. Also the error decreases with decrease in grid size upto a certain value and then it again increases.The error value is lowest for the cds scheme. The maximum error is seen for a right skewed scheme. For a faster convergence cds is better than the other two schemes.
Importance of skewed schemes over CDS
The skewed schemes are useful for computing derivatives at the boundaries. These schemes utilize values either from left or right of the present node points.For e.g. in a fluid flow domain or a heat flow in pipe, in both of these cases the fluid or heat flows in one direction. In these cases, skewed schemes are useful in computing the flow rate derivatives in one particular direction. In these problems cds is not useful as it takes information from both sides. Also in these domain skewed schemes are preferred for faster and accurate solutions.
Conclusion
From the three different scheme central differencing scheme gives least error and is preffered other than two skewed schemes. According to the flow of domain these three schemes should be chosen with proper conditions and obtaining accurate results.