Week 3 - Taylor table method and Matlab code
Aim
- To derive the 4th order approximations for the second order derivatives of a given function using the following methods.
- Central Differencing
- Skewed Right Sided Differencing
- Skewed Left Sided Differencing
- To write a code in Matlab/Octave to study the variations of Skewed left and right approximations and central difference approximation versus dx for a given function.
- Also to compare the above with forward and backward differencing methods
Introduction
Taylor’s series is an important part of CFD. This method is the basic for many problems that we may come across. The Taylor’s table is simply an extension to the Taylor series. Using this method, we can achieve the solution of a problem to the highest accuracy.
The above mentioned methods are methods of approximation and thus give us an approximation of the solution. The actual solution can be obtained by calculating it analytically.
In order to determine the number of nodes to be considered for the approximation, the following equations are used.
For skewed scheme,
N=P+Q
Where N= number of nodes
P= order of derivative
Q=order of approximation
Similarly, for central differencing scheme,
N=P+Q-1
Where N= number of nodes
P= order of derivative
Q=order of approximation
CENTRAL DIFFERENCING SCHEME
In this scheme, we initially calculate the number of nodes using the equation illustrated below
N=P+Q-1
N=2+4-1
N=5
Below is the formula for 4th order central differencing.
`(d^2 y)/(dx^2 )=(af(i-2)+bf(i-1)+cf(i)+df(i+1)+ef(i+2))/dx^2`
|
Taylor series terms |
f(x) |
dx f’(x) |
dx2 f’’(x) |
dx3 f’’’(x) |
dx4 f’’’’(x) |
|
Numerical stencil |
|||||
|
f(i-2) |
a |
-2a |
2a |
-4a/3 |
2a/3 |
|
f(i-1) |
b |
-b |
b/2 |
-b/6 |
b/24 |
|
f(i) |
c |
0 |
0 |
0 |
0 |
|
f(i+1) |
d |
1d |
d/2 |
d/6 |
d/24 |
|
f(i+2) |
e |
2e |
2e |
4e/3 |
2e/3 |
Upon solving the group of equations, the above are the values of the coefficients a, b, c, d and e respectively.
SKEWED DIFFERENCING SCHEME
In this scheme we calculate the number of nodes using the following equation
N=P+Q
N=2+4
N=6
Right scheme
Below is the formula for 4th order Skewed right sided differencing.
`(d^2 y)/(dx^2 )=(af(i)+bf(i+1)+cf(i+2)+df(i+3)+ef(i+4)+gf(i+5))/dx^2`
|
Taylor series terms |
f(x) |
dx f’(x) |
dx2 f’’(x) |
dx3 f’’’(x) |
dx4 f’’’’(x) |
dx5 f’’’’’(x) |
|
Numerical stencil |
||||||
|
f(i) |
a |
0 |
0 |
0 |
0 |
0 |
|
f(i+1) |
b |
b |
b/2 |
b/6 |
b/24 |
b/120 |
|
f(i+2) |
c |
2c |
4c/2 |
8c/6 |
16c/24 |
32c/120 |
|
f(i+3) |
d |
3d |
9d/2 |
27d/6 |
81d/24 |
243d/120 |
|
f(i+4) |
e |
4e |
16e/2 |
64e/6 |
256e/24 |
1024e/120 |
|
f(i+5) |
g |
5g |
25g/2 |
125g/6 |
625g/24 |
3125g/120 |
Upon solving the group of equations, the above are the values of the coefficients a, b, c, d e and g respectively.
Left scheme
Below is the formula for the 4th order approximation for the skewed left scheme
`(d^2 y)/(dx^2 )=(af(i)+bf(i-1)+cf(i-2)+df(i-3)+ef(i-4)+gf(i-5))/dx^2`
|
Taylor series terms |
f(x) |
dx f’(x) |
dx2 f’’(x) |
dx3 f’’’(x) |
dx4 f’’’’(x) |
dx5 f’’’’’(x) |
|
Numerical stencil |
||||||
|
f(i) |
a |
0 |
0 |
0 |
0 |
0 |
|
f(i+1) |
b |
-b |
b/2 |
-b/6 |
b/24 |
-b/120 |
|
f(i+2) |
c |
-2c |
4c/2 |
-8c/6 |
16c/24 |
-32c/120 |
|
f(i+3) |
d |
-3d |
9d/2 |
-27d/6 |
81d/24 |
-243d/120 |
|
f(i+4) |
e |
-4e |
16e/2 |
-64e/6 |
256e/24 |
-1024e/120 |
|
f(i+5) |
g |
-5g |
25g/2 |
-125g/6 |
625g/24 |
-3125g/120 |
Upon solving the group of equations, the above are the values of the coefficients a,b,c,d,e and g respectively.
PROGRAM LISTING
CODE TO SOLVE FOR CENTRAL DIFFERENCING COEFFICIENTS
% program to calculate the values of the coefficients in the central methods
clear all
close all
clc
%values of the coefficients of each derivative
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];
% solving the system of equations to obtain the values of a,b,c,d,e,and g
X=inv(A)*BCODE TO SOLVE FOR SKEWED RIGHT DIFFERENCING COEFFICIENTS
% program to calculate the values of the coefficients in the skewed right method
clear all
close all
clc
%values of the coefficients of each derivative
A=[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];
B=[0;0;1;0;0;0];
% solving the system of equations to obtain the values of a,b,c,d,e,and g
X=inv(A)*BCODE TO SOLVE FOR SKEWED LEFT DIFFERENCING COEFFICIENTS
% program to calculate the values of the coefficients in the skewed left method
clear all
close all
clc
%values of the coefficients of each derivative
A=[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];
B=[0;0;1;0;0;0];
% solving the system of equations to obtain the values of a,b,c,d,e,and g
X=inv(A)*BMAIN PROGRAM
clear all
close all
clc
% f(x)= exp(x)cos(x)
dx=linspace(pi/7,pi/7000,50);
x=pi/3;
analytical=-2*exp(x)*sin(x);
for i = 1:length(dx)
skewed_right_error(i) = fourth_order_skewed_right(x,dx(i));
skewed_left_error(i) = fourth_order_skewed_left(x,dx(i));
central_error(i) = fourth_order_central(x,dx(i));
forward_diff_error(i)=forward_differencing(x,dx(i));
backward_diff_error(i)=backward_differencing(x,dx(i));
end
figure(1)
loglog(dx,central_error,'b','linewidth',2,'marker','o','markersize',5)
hold on
loglog(dx,skewed_right_error,'r','linewidth',2,'marker','d','markersize',5)
hold on
loglog(dx,skewed_left_error,'g','linewidth',2,'marker','p','markersize',5)
hold on
h=legend('central_error','skewed_right_error','skewed_left_error');
legend(h,'location','northwest')
set(h,'fontsize',10);
xlabel("dx",'fontsize',25)
ylabel("error",'fontsize',25)
figure(2)
loglog(dx,forward_diff_error,'y','linewidth',2,'marker','*','markersize',5)
hold on
loglog(dx,backward_diff_error,'c','linewidth',2,'marker','s','markersize',5)
hold on
loglog(dx,central_error,'b','linewidth',2,'marker','o','markersize',5)
hold on
loglog(dx,skewed_right_error,'r','linewidth',2,'marker','d','markersize',5)
hold on
loglog(dx,skewed_left_error,'g','linewidth',2,'marker','p','markersize',5)
hold on
title('first order errors and fourth order errors')
h=legend('forward_diff_error','backward_diff_error','central_error','skewed_right_error','skewed_left_error');
legend(h,'location','northwest')
set(h,'fontsize',10);
xlabel("dx",'fontsize',25)
ylabel("error",'fontsize',25)CODE FOR SKEWED RIGHT SIDED METHOD
function out = fourth_order_skewed_right(x,dx)
%analytical function : exp(x)*cos(x)
% analytical derivative
%f''(x)= -2*exp(x)*sin(x)
analytical=-2*exp(x)*sin(x);
fourth_order_skewed_right=(3.75*(exp(x)*cos(x))+(-12.83333*(exp(x+dx)*cos(x+dx)))+(17.83333*(exp(x+2*dx)*cos(x+2*dx)))+(-13.00000*(exp(x+3*dx)*cos(x+3*dx)))+5.08333*((exp(x+4*dx)*cos(x+4*dx)))+(-0.8333*((exp(x+5*dx)*cos(x+5*dx)))))/dx^2;
out= abs(fourth_order_skewed_right-analytical);
endCODE FOR SKEWED LEFT SIDED METHOD
function out = fourth_order_skewed_left(x,dx)
%analytical function : exp(x)*cos(x)
% analytical derivative
%f''(x)= -2*exp(x)*sin(x)
analytical=-2*exp(x)*sin(x);
fourth_order_skewed_left=(3.75*(exp(x)*cos(x))+(-12.83333*(exp(x-dx)*cos(x-dx)))+(17.83333*(exp(x-2*dx)*cos(x-2*dx)))+(-13.00000*(exp(x-3*dx)*cos(x-3*dx)))+5.08333*((exp(x-4*dx)*cos(x-4*dx)))+(-0.8333*((exp(x-5*dx)*cos(x-5*dx)))))/dx^2;
out= abs(fourth_order_skewed_left-analytical);
endCODE FOR CENTRAL DIFFERENCING METHOD
function out = fourth_order_central(x,dx)
%analytical function : exp(x)*cos(x)
% analytical derivative
%f''(x)= -2*exp(x)*sin(x)
analytical=-2*exp(x)*sin(x);
fourth_order_central=((-0.083333*(exp(x-2*dx)*cos(x-2*dx)))+(1.333333*(exp(x-dx)*cos(x-dx)))+(-2.500000*(exp(x)*cos(x)))+(1.333333*(exp(x+dx)*cos(x+dx)))+(-0.083333*((exp(x+2*dx)*cos(x+2*dx)))))/dx^2;
out= abs(fourth_order_central-analytical);
endCODE FOR FORWARD DIFFERENCING (FIRST ORDER)
function out = forward_differencing(x,dx)
%analytical function : exp(x)*cos(x)
% analytical derivative
%f''(x)= -2*exp(x)*sin(x)
analytical = -2*exp(x)*sin(x);
%numerical derivative
forward_differencing=((exp(x+dx)*sin(x+dx)) - exp(x)*sin(x))/dx;
out= abs(forward_differencing - analytical);
end
CODE FOR BACKWARD DIFFERENCING (FIRST ORDER)
function out = backward_differencing(x,dx)
%analytical function : exp(x)*cos(x)
% analytical derivative
%f''(x)= -2*exp(x)*sin(x)
analytical = -2*exp(x)*sin(x);
%numerical derivative
forward_differencing=((exp(x)*sin(x)) - exp(x-dx)*sin(x-dx))/dx;
out= abs(forward_differencing - analytical);
endRESULTS AND INFERENCES
Plot 1
This is a loglog plot of the fourth orders errors versus dx
Plot 2
THis is a comparison loglog plot of the first and fourth order errors versus dx
Inferences:
- It can be observed from the above graphs that, the central difference scheme has the least error among the rest of the methods.
- This implies the central difference scheme provides a much more accurate result comparing to the other methods of approximation.
- The orders of accuracy to result is central error>skewed left error>skewed right error>backward error>forward error.
- In comparison to the forward and backward schemes, the second order schemes have higher accuracy.
- Also, if dx becomes smaller, there is a probability of round off error due to which the solution may vary.
- Also if dx becomes large, the truncation error will become larger.
Importance of the Skewed scheme and its advantage over the central difference scheme
From the above results, we can determine that the central difference scheme is the most accurate method to obtain the solution. However, this method fails when there is a lack of symmetric nodes. This implies, we cannot apply this method in the extremes of a system. The skewed scheme is useful when symetric nodes are not present. Thus this method can be applied to the extremities of a system of at any location of a system for the desired solution to be obtained.