Week 3 - Taylor table method and Matlab code

When we want to solve PDE's or ODE's in the 2D domain, we have to discretize the equations using Taylor's table method To find the derivative of PDEs or ODE's using Taylor's table method we have to consider stencil. The value of forward or backward or central scheme points from a given initial condition can be calculated by this method.

There are different schemes that are used to calculate the values viz- Central differencing, Foward differencing, and Backward differencing.

Forward differencing is when the points are considered in the forward direction of the flow.

Backward differencing is when the points are considered in the backward or the opposite direction of the flow.

Central differencing is when the points are considered on either side of the selected location.

For the given problem, we have to find the Taylor series expansion for all the three schemes. This solution needed is a fourth-order approximation, and the order of the equation is second order.

For this condition, we calculated that the no. Terms/ nodes that are needed to be considered in forward and the backward difference is six, and for the central difference, it is 5.

We define the derivative of the equation as a series of terms with unknown variables.

FDS = `a*f(i) + bf(i+1) + c*f(i+2) + d*f(i+3) + e*f(i+4) + g*f(i+5)`

BDS = `a*f(i) + bf(i-1) + c*f(i-2) + d*f(i-3) + e*f(i-4) + g*f(i-5)`

CDS = `a*f(i-2) + bf(i-1) + c*f(i) + d*f(i+1) + e*f(i+2)`

Now using Taylor's table method the terms of the expansions are calculated. Then the final derivative of the given derivative term is expressed.

FDS = `{a*f(x) + bf(x+dx) + c*f(x+2dx) + d*f(x+3dx) + e*f(x+4dx) + g*f(x+5dx)}/dx^2`

BDS = `{a*f(x) + bf(x-dx) + c*f(x-2dx) + d*f(x-3dx) + e*f(x-4dx) + g*f(x-5dx)}/dx^2`

CDS = `{a*f(x-2dx) + bf(x-dx) + c*f(x) + d*f(x+dx) + e*f(x+2dx) }/dx^2`

The terms a,b,c,d,e,g for BDS and FDS are solved by making matrices of the terms calculated from Taylor's table. Then these variables are substituted back to the equations and derivative of the given function is calculated.

Similarly, the terms in CDS a,b,c,d,e are solved by making matrices of the terms calculated from Taylor's table. Then these variables are substituted back to the equations and derivative of the given function is calculated.

clear all
close all
clc

% input conditions 
x= pi/4;
dx = linspace(pi/4,pi/400,10);

% analytical derivative - after solving the given fucntion manually the
% fucntion is expressed below 

analytical_derivative = -2*exp(x)*sin(x);
%% matrix calculations
% Coefficients calculation
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];

C = [0;0;1;0;0;0];
B = inv(A)*C;

% skewed left side difference 

D = [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];
E = [0;0;1;0;0;0];
G = inv(D)*E;
% central differencing 

H = [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];
I=[0;0;1;0;0];
J = inv(H)*I;
%% the main calcuations of the error and differences 
for i=1:length(dx)
% central differencing 
central_difference(i) = ((J(1)*f(x-(2*dx(i))))+ J(2)*f(x-dx(i)) +(J(3)*f(x)) + J(4)*f(x+dx(i))+ (J(5)*f(x+(2*dx(i)))))/(dx(i)^2);

% skewed left side differencing 
left_differencing(i) = ((3.75*f(x)) -(12.8333*f(x-dx(i))) + (17.8333*f(x-(2*dx(i)))) -(13*f(x-(3*dx(i))))+(5.0833*f(x-(4*dx(i)))) -(0.8333*f(x-(5*dx(i)))))/(dx(i)^2);

% skewed right side differencing 
right_differencing(i) = ((3.75*f(x)) -(12.8333*f(x+dx(i))) +(17.8333*f(x+(2*dx(i)))) -(13*f(x+(3*dx(i)))) + (5.0833*f(x+(4*dx(i)))) -(0.833*f(x+(5*dx(i)))) )/(dx(i)^2);

% error 
central_error(i) = abs(analytical_derivative - central_difference(i));
left_error(i) = abs(analytical_derivative - left_differencing(i));
right_error(i) = abs(analytical_derivative - right_differencing(i));
end

% plotting commands 
loglog(dx,central_error)
hold on
loglog(dx,left_error)
loglog(dx,right_error)
grid on
legend('central error','left error','right error')

the function that I have defined in the main code is below

function y=f(x)
y= exp(x)*cos(x);
end

The error graph is depicted below

The central error is lowest than right and left side error. The Right side error gives more error when compared to left side error. The central difference is preferred usually.