Week 3 - Taylor table method and Matlab code
Fourth order approximations of the second order derivative using the following schemes are derived with the help of programming in MATLAB
- Central difference scheme
- Skewed right sided difference scheme
- Skewed left sided difference scheme
A Matlab code is written to evaluate the second order derivative of the analytical function exp(x)⋅cos(x)exp(x)⋅cos(x) and compared it with the 3 numerical approximations derived.
- Plot which compares the absolute error between the above mentioned schemes is generated.
- Compared with FDS and BDS approximations
- Advantages of the skewed scheme over CD scheme is explained.
Formulation of the Taylor table for CDS.
Taylor Table Method Derivation:
Skewed right sided difference scheme:
`(∂^2 f)/(∂f^2 )=(af_i+bf_(i+1)+cf_(i+2)+df_(i+3)+ef_(i+4)+ff_(i+5))/(∆x^2 )` …..main equation
`f_(i+1)=f_i+(∆xf_i^')/1+(∆x^2 f_i^'')/2+(〖∆x〗^3 f_i^''')/6+(〖∆x〗^4 f_i^'''')/24+(〖∆x〗^5 f_i^''''')/120+higher order terms`....eqn1
`f_(i+2)=f_i+(2∆xf_i^')/1+(4∆x^2 f_i^'')/2+(8〖∆x〗^3 f_i^''')/6+(16〖∆x〗^4 f_i^'''')/24+(32∆x^5 f_i^''''')/120+higher order terms`......eqn2
`f_(i+3)=f_i+(3∆xf_i^')/1+(9∆x^2 f_i^'')/2+(27〖∆x〗^3 f_i^''')/6+(81〖∆x〗^4 f_i^'''')/24+(243∆x^5 f_i^''''')/120+higher order terms`........eqn3
`f_(i+4)=f_i+(4∆xf_i^i)/1+(16∆x^2 f_i^'')/2+(64〖∆x〗^3 f_i^''')/6+(256〖∆x〗^4 f_i^'''')/24+(1024∆x^5 f_i^''''')/120+higher order terms`.......eqn4
`f_(i+5)=f_i+(5∆xf_i^')/1+(25∆x^2 f_i^'')/2+(〖125∆x^3 f_i^''')/6+(625〖∆x〗^4 f_i^'''')/24+(〖3125∆x^5 f_i^''''')/120+higher order terms` ......eqn5`
Substituting equations 1 to 5 in main equation and comparing coefficients of
`f_i,f_i^' ∆x,f_i^'' 〖∆x〗^2,f_i^''' 〖∆x〗^3,f_i^'''' 〖∆x〗^4 and f_i^''''' 〖∆x〗^5`
We get,
`a+b+c+d+e+f=0`
`b+2c+3d+4e+5f=0`
`b/2+2c+9/2 d+8e+25/2 f=1`
`b/6+4/3 c+27/6 d+64/6 e+125/6 f=0`
`b/24+2/3 c+81/24 d+256/24 e+625/24 f=0`
`b/120+32/120 c+243/120 d+1024/120 e+3125/120 f=0`
Solving the above equations we get
a = 3.75, b = -12.833, c = 17.833, d = -13, e = 5.0833, f= -0.833
coding:
% skewed right-sided difference scheme requires information on the right side
[C]= [ 1 1 1 1 1 1;
0 1 2 3 4 5;
0 1/2 2 9/2 8 25/2;
0 1/6 4/3 9/2 32/3 125/6;
0 1/24 2/3 27/8 32/3 625/24;
0 1/120 32/120 243/120 1024/120 3125/120];
[D] = [0;0;1;0;0;0];
%solving matrix C&D to find constants
C*y=D;
y = inv(C)*D
Skewed left sided difference scheme:
`(∂^2 f)/(∂f^2 )=(af_i+bf_(i-1)+cf_(i-2)+df_(i-3)+ef_(i-4)+ff_(i-5))/(∆x^2 )` …..main equation
`f_(i-1)=f_i-(∆xf_i^')/1+(∆x^2 f_i^'')/2-〖∆x〗^3/6 f_i^'''+(〖∆x〗^4 f_i^'''')/24-(〖∆x〗^5 f_i^''''')/120+higher order terms`......eqn1`
`f_(i-2)=f_i-(2∆xf_i^')/1+(4∆x^2 f_i^'')/2-(8〖∆x〗^3 f_i^''')/6+(16〖∆x〗^4 f_i^'''')/24-(〖32∆x^5 f_i^''''')/120+higher order terms`....eqn2`
`f_(i-3)=f_i-(3∆xf_i^')/1+(9∆x^2 f_i^'')/2-(27〖∆x〗^3 f_i^''')/6+(81〖∆x〗^4 f_i^'''')/24-(〖243∆x^5 f_i^''''')/120+higher order terms`.....eqn3`
`f_(i-4)=f_i-(4∆xf_i^i)/1+(16∆x^2 f_i^'')/2-(64〖∆x〗^3 f_i^''')/6+(256〖∆x〗^4 f_i^'''')/24-(〖1024∆x^5 f_i^''''')/120+higher order terms`......eqn4`
`f_(i-5)=f_i-(5∆xf_i^')/1+(25∆x^2 f_i^'')/2-(〖125∆x^3 f_i^''')/6+(625〖∆x〗^4 f_i^'''')/24-(〖3125∆x^5 f_i^''''')/120+higher order terms`....eqn5`
Substituting equations 1 to 5 in main equation and comparing coefficients of
`f_i,f_i^' ∆x,f_i^'' 〖∆x〗^2,f_i^''' 〖∆x〗^3,f_i^'''' 〖∆x〗^4 and f_i^''''' 〖∆x〗^5`
We get,
`a+b+c+d+e+f=0`
`-b-2c-3d-4e-5f=0`
`b/2+2c+9/2 d+8e+25/2 f=1`
`-b/6-4/3 c-27/6 d-64/6 e-125/6 f=0`
`b/24+2/3 c+81/24 d+256/24 e+625/24 f=0`
`-b/120-32/120 c-243/120 d-1024/120 e-3125/120 f=0`
Solving the above equations in octave we get
a = 3.75, b = -12.833, c = 17.833, d = -13, e = 5.0833, f= -0.833
% skewed left-sided difference scheme requires information on the right side
[E]= [ 1 1 1 1 1 1;
0 -1 -2 -3 -4 -5;
0 1/2 2 9/2 8 25/2;
0 -1/6 -4/3 -9/2 -32/3 -125/6;
0 1/24 2/3 27/8 32/3 625/24;0 -1/120 -32/120 -243/120 -1024/120 -3125/120];
[F] = [0;0;1;0;0;0];
%solving matrix C&D to find constants
E*z=F;
z=inv(E)*FCentral difference scheme:
`(∂^2 f)/(∂f^2 )=(a〖f〗_i+bf_(i-1)+c〖f〗_(i-2)+d〖f〗_(i+1)+e〖f〗_(i+2))/(∆x^2 )` …..main equation
`f_(i+1)=f_i+(∆xf_i^')/1+(∆x^2 f_i^'')/2+(〖∆x〗^3 f_i^''')/6+(〖∆x〗^4 f_i^'''')/24+higher order terms`......eqn1
`f_(i+2)=f_i+(2∆xf_i^')/1+(4∆x^2 f_i^'')/2+(8〖∆x〗^3 f_i^''')/6+(16〖∆x〗^4 f_i^'''')/24+higher order terms`.......eqn2
`f_(i-1)=f_i-(∆xf_i^')/1+(∆x^2 f_i^'')/2-〖∆x〗^3/6 f_i^'''+(〖∆x〗^4 f_i^'''')/24-higher order terms`......eqn3
`f_(i-2)=f_i-(2∆xf_i^')/1+(4∆x^2 f_i^'')/2-(8〖∆x〗^3 f_i^''')/6+(16〖∆x〗^4 f_i^'''')/24-higher order terms`.....eqn4
Substituting equations 1 to 4 in main equation and comparing coefficients of
`f_i,f_i^' ∆x,f_i^'' 〖∆x〗^2,f_i^''' 〖∆x〗^3 and f_i^'''' 〖∆x〗^4`
We get,
`a+b+c+d+e=0`
`-b-2c+d+2e=0`
`b/2+2c+d/2+2e=1`
`-b/6-4c/3+d/6+4e/3=0`
`b/24+2/3 c+d/24+2e/3=0`
Solving the above equations in octave we get
a = 2.5, b = 1.334, c = -0.0833, d = 1.334, e = -0.0833
coding:
clear all
close all
clc
% Taylor_table_coeffecients
% central difference scheme requires information from both sides ( symmetric)
% To obtain the 4th order approximated second-order derivative, we need information from 5 equally distributed neighboring points
[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];
[B] = [0;0;1;0;0:0];
A*x=B
x = inv(A)*Bprogram for comparing fds and bds for the given function
clear all
close all
clc
% analytical function f' = exp(x)*cos(x)
%f'' = -2*exp(x)*sin(x)
x = pi/3;
dx = linspace(pi/10,pi/400,100);
analytical_function = exp(x)*cos(x);
analytical_derivative = -2*exp(x)*sin(x);
for i = 1:length(dx);
backward_difference(i) = (exp(x)*cos(x)- exp(x-dx(i))*cos(x-dx(i)))/dx(i);
forward_difference(i) = (exp(x+dx(i))*cos(x+dx(i))-exp(x)*cos(x))/dx(i);
backward_error(i) = abs(backward_difference(i)- analytical_derivative);
forward_error(i) = abs(forward_difference(i)- analytical_derivative);
end
%Error comparison between the above two schemes
figure(1)
loglog(dx,backward_error,'-ks', 'markerfacecolor','b','markersize',10,'linewidth',2)
hold on
loglog(dx,forward_error,'-kd', 'markerfacecolor','g','markersize',10,'linewidth',2)
hold on
legend ('backward error', 'forward error','location','northwest');
xlabel('dx')
ylabel('error')
title ( 'FDS AND BDS ERROR COMPARSION')
program for comparing skewed and central difference scheme
clear all
close all
clc
%Taylor_table_test
% analytical function f' = exp(x)*cos(x)
%f'' = -2*exp(x)*sin(x)
x = pi/3;
dx = linspace(pi/10,pi/400,100);
analytical_function = exp(x)*cos(x);
analytical_derivative = -2*exp(x)*sin(x);
for i = 1:length(dx);
a(i) = (exp(x-(5*dx(i)))*cos(x-(5*dx(i))));
b(i) = (exp(x-(4*dx(i)))*cos(x-(4*dx(i))));
c(i) = (exp(x-(3*dx(i)))*cos(x-(3*dx(i))));
d(i) = (exp(x-(2*dx(i)))*cos(x-(2*dx(i))));
e(i) = (exp(x-dx(i))*cos(x-dx(i)));
f(i) = exp(x)*cos(x)
g(i) = (exp(x+dx(i))*cos(x+dx(i)));
h(i) = (exp(x+(2*dx(i)))*cos(x+(2*dx(i))));
j(i) = (exp(x+(3*dx(i)))*cos(x+(3*dx(i))));
k(i) = (exp(x+(4*dx(i)))*cos(x+(4*dx(i))));
l(i) = (exp(x+(5*dx(i)))*cos(x+(5*dx(i))));
fourth_order_second_derivative_skewed_right(i) = (45*(f(i))-154*(g(i))+214*(h(i))-156*(j(i))+61*(k(i))-10*(l(i)))/(12*(dx(i)^2));
fourth_order_second_derivative_skewed_left(i) = (45*(f(i))-154*(e(i))+214*(d(i))-156*(c(i))+61*(b(i))-10*(a(i)))/(12*(dx(i)^2));
fourth_order_second_derivative_central_difference(i) = (-1*(d(i))+16*(e(i))-30*(f(i))+16*(g(i))-(h(i)))/(12*(dx(i)^2));
fourth_order_skewed_right_error(i) = abs(fourth_order_second_derivative_skewed_right(i)- analytical_derivative)
fourth_order_skewed_left_error(i) = abs(fourth_order_second_derivative_skewed_left(i)- analytical_derivative)
fourth_order_central_difference_error(i) = abs(fourth_order_second_derivative_central_difference(i)- analytical_derivative)
end
%Error comparsion between the above three schemes
figure(1)
loglog(dx,fourth_order_skewed_right_error, 'r')
hold on
loglog(dx,fourth_order_skewed_left_error,'b')
hold on
loglog(dx,fourth_order_central_difference_error,'g')
hold on
legend('fourth order skewed right error','fourth order skewed left error','fourth order central difference error','location','northwest');
xlabel('dx')
ylabel('error')
title( ' Comparision between different schemes of second order derivative of 4th order approximation')
How do these errors compare with FDS and BDS scheme of approximations?
FDS and BDS both acquires information from the current point and the right and left neighbouring points respectively resulting im 1st order approx for 1st order derivative with two points information whereas central diffeerence results im 2nd order approx for 1st order derivative, with two points information.With the same amountb of information Central schemes gives more accuracy due to the reduced truncation error.
FDS and BDS has some constrain, so that the skewed schemes can result in higher order approximation. Due to the perspective of symmetry, CDS schemes are more better than skewed schemes as it gains both the sides information which enhances it's ability to find more accurate derivatives at a certain point.
Advantages of skewed scheme over CD scheme
Skewed difference scheme works well whatever be the number of points available on the grid. For example, if there is only node available at the left and 5 nodes on the right, and we need 4th order approx for second derivative, then at that time only skewed schemes work out well. So, if we need to compute a single or higher order derivative at a corner node, then central scheme is no longer useful and there applies the skewed scheme.