Week 3 - Taylor table method and Matlab code

INTRODUCTION:

Numerical Solution can only give answers at discrete points in the domain which are called as grid points. The finite difference methods is for replacing the partial derivatives appearing in governing equations of fluid dynamics, with algebraic difference quotients yielding a system of algebraic equations which can be solved for the flow field variables at the specific, discrete grid points in the flow.

SCHEMES FOR FOURTH ORDER APPROXIMATION FOR SECOND ORDER DERIVATIVE:-

We have a set of grid points in our numerical stencil along the X-plane. Each grid point has its own taylor series functions. We first create a taylor table from these sets of grid points as shown below:

1. CENTRAL DIFFERENCE SCHEME:

For Central Difference Scheme all the unknowns are solved using the code given below:

E=[ 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 ];

F = [ 0;0;1;0;0 ];

Z = rats(linsolve(E,F)); 

The output we get for the given lines of code,

2. SKEWED RIGHT DIFFERENCE SCHEME:

For Skewed Right Sided Difference Scheme all the unknowns are solved using the code given below:

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 ];

X = rats(linsolve(A,B)); 

The output we get for the given lines of code,

3.SKEWED LEFT DIFFERENCE SCHEME:

For Skewed Left Sided Difference Scheme all the unknowns are solved using the code given below:

C = [ 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 ];

D = [ 0;0;1;0;0;0 ];

Y = rats(linsolve(C,D)); 

The output we get for the given lines of code,

The partial derivatives are written for schemes which are of 2nd order and 4th order approximations. For approximating a derivative to an arbitrary order of accuracy we use finite difference. For the given problem we say that it is of 4th order accuracy because terms of the order 5 or more are completely neglected. Each column in a tabel has 5 or 6 unknowns. An algebraic equation is formed with these unknowns. If equations are formed for each and every column in a table we end up with a system of equations which is called as a matrix. With the help of a matrix the values of unknowns are calculated.

OBJECTIVE:

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.

MATLAB CODE:

clear all
close all
clc

%Initial Parameters for the given system
x=pi/3;
dx=linspace(pi/4,pi/400,100);

% Analytical Function = exp^(x)*cos(x)

analytical_derivative= -2*exp(x)*sin(x);

for i=1:length(dx)
    
    % Taylor Series Function at various node points
    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))));
    
    %Expression Formulas for various Taylor Series Functions
    %Forward Difference Scheme
    FDS(i)=(g(i)-f(i))/(dx(i));
    
    %Backward Difference Scheme
    BDS(i)=(f(i)-e(i))/(dx(i));
    
    %Central Difference Scheme
    CDS(i)=(-1*(d(i))+16*(e(i))-30*(f(i))+16*(g(i))-1*(h(i)))/(12*(dx(i)^2));
    
    %Skewed Left Sided Difference
    skew_left(i)=(45*(f(i))-154*(e(i))+214*(d(i))-156*(c(i))+61*(b(i))-10*(a(i)))/(12*(dx(i)^2));
    
    %Skewed Right Sided Difference
    skew_right(i)=(45*(f(i))-154*(g(i))+214*(h(i))-156*(j(i))+61*(k(i))-10*(l(i)))/(12*(dx(i)^2));
    
    error_FDS(i)=abs(FDS(i)-analytical_derivative);
    error_BDS(i)=abs(BDS(i)-analytical_derivative);
    error_CDS(i)=abs(CDS(i)-analytical_derivative);
    error_LSD(i)=abs(skew_left(i)-analytical_derivative);
    error_RSD(i)=abs(skew_right(i)-analytical_derivative);
    
end

figure(1)
subplot(1,2,1)
loglog(dx,error_FDS,'linewidth',1,'marker','hexagram')
xlabel('dx')
ylabel('Absolute Error')
title('Absolute Error in Difference Schemes Vs dx')
grid on
hold on
loglog(dx,error_BDS,'linewidth',1,'marker','pentagram')
xlabel('dx')
ylabel('Absolute Error')
legend('Forward Difference Scheme','Backward Difference Scheme')
subplot(1,2,2)
loglog(dx,error_RSD,'linewidth',1,'marker','hexagram')
xlabel('dx')
ylabel('Absolute Error')
title('Absolute Error in various Finite Difference schemes')
grid on
hold on
loglog(dx,error_LSD,'linewidth',1,'marker','pentagram')
xlabel('dx')
ylabel('Absolute Error')
grid on
hold on
loglog(dx,error_CDS,'linewidth',1,'marker','d')
xlabel('dx')
ylabel('Absolute Error')
grid on
legend('Skewed Right Sided Difference','Skewed Left Sided Difference','Central Difference Scheme')

%To find the exact slope of the curve
L_1=log10(dx);
L_2=log10(error_FDS);
L_3=log10(error_BDS);
L_4=log10(error_RSD);
L_5=log10(error_LSD);
L_6=log10(error_CDS);
%g_1=gradient(L_2,L_1);
%g_2=gradient(L_3,L_1);
g_3=gradient(L_4,L_1);
g_4=gradient(L_5,L_1);
g_5=gradient(L_6,L_1);

figure(2)
%plot(g_1)
%hold on 
%plot(g_2)
%hold on
plot(g_3,'linewidth',1,'marker','pentagram')
hold on
plot(g_4,'linewidth',1,'marker','hexagram')
hold on
plot(g_5,'linewidth',1,'marker','d')
legend('Skewed Right Sided Difference','Skewed Left Sided Difference','Central Difference Scheme')
title('Slope of the curve')

OUTPUT:

• In the first plot for schemes that take only one grid point in the forward or backward direction in a numerical stencil along the X-plane, we can tell that as the value of dx keeps on increasing the absolute value of the error will vary at a rapid rate by either increasing or decreasing until the value of dx reaches and after that the value of absolute error becomes unstable for different schemes
• In the second plot the value of absolute error is less for Central Difference Scheme and it increases for the Skewed Schemes. For a range of dx values the plot is linear in Central Difference Schemes and starts at 10^-9 but for the skewed schemes the value of absolute error is a bit higher and starts at 10^-7.
• The linear plot will not be the same if the value of dx becomes greater than 1 and because of which the absolute error values will become unstable.
• The improvement of approximations quality is due to the node(i) being at a closer distance which is observed in case of Central difference scheme seperated by a distance of 2 grid points apart, but for the skewed schemes since the distance to the node (i) is 5 grid points apart which is the reason behind the absolute value for errors being higher.

• The slope of the curve will remain constant at 4 when the log values of dx keep on increasing for the various finite difference schemes.

Why a skewed scheme is useful? What can a skewed scheme do that a CD scheme cannot do?

• When data from the grid points is only available from either the left or right end i.e. at nodes closer to the boundaries, for computations that require a high order approximationsn we make use of skewed schemes. Having a symmetric stencil for computation is not possible in this scenarios like we do for central differencing scheme or having a low order accuracy. For consistency of scheme to be maintained the values of dx are kept as low as possible for providing proper computational plots.

CONCLUSION:

Thus a fourth order approximations are derived for a second order derivative making use of Taylor Table method by providing a suitable MATLAB code.