Week 3 - Taylor table method and Matlab code

The code for the given challenge is enclosed below:

This part of the code is the primary part of the code where function calling and for loop is used:

 close all 
clear all
clc

dx=linspace(pi/4,pi/4000,20); % creating the mesh
x=pi/3; % the point of evaluation


for i=1:length(dx)
    fourth_error_CDS(i)= fourth_order_approximation_CDS(x,dx(i));
     fourth_error_Right_skewed(i)= fourth_order_approximation_Right_skewed(x,dx(i));
     fourth_error_Left_skewed(i)=fourth_order_approximation_Left_skewed(x,dx(i));
end
figure(1)
loglog(dx,fourth_error_CDS,'r');
hold on
loglog(dx,fourth_error_Right_skewed,'color','b');
hold on
loglog(dx,fourth_error_Left_skewed,'color','g');
legend('fourth_order_error_CDS','fourth_order_error_Right_skewed','fourth_order_error_Left_skewed');
xlabel('mesh spacing');
ylabel('error');

For the CDS scheme:

function  fourth_error_CDS= fourth_order_approximation_CDS(x,dx)
A2=[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 ];
[X2]=A2\B';
% to compare analytical & numerical derivative exp(x)*cos(x)
analytical_second_derivative= -2*exp(x)*sin(x);
% taking the first order in forward diffrencing
 central_diffrence_fouth_order=(X2(1)*exp(x-2*dx)*cos(x-2*dx)+X2(2)*exp(x-dx)*cos(x-dx)+X2(3)*exp(x)*cos(x)+X2(4)*exp(x+dx)*cos(x+dx)+X2(5)*exp(x+2*dx)*cos(x+2*dx))/(dx)^2;
 fourth_error_CDS=abs(central_diffrence_fouth_order-analytical_second_derivative);
end

For the Left scheme:

function fourth_order_error_Left_skewed= fourth_order_approximation_Left_skewed(x,dx)
A1=[1 1 1 1 1 1; 0 -1 -2 -3 -4 -5; 0 1/2 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];
[X1]=A1\B';
% to compare analytical & numerical derivative exp(x)*cos(x)
analytical_second_derivative= -2*exp(x)*sin(x);
% taking the first order in forward diffrencing
Left_skewed_fouth_order= (X1(1)*exp(x)*cos(x)+X1(2)*exp(x-dx)*cos(x-dx)+X1(3)*exp(x-2*dx)*cos(x-2*dx)+X1(4)*exp(x-3*dx)*cos(x-3*dx)+X1(5)*exp(x-4*dx)*cos(x-4*dx)+X1(6)*exp(x-5*dx)*cos(x-5*dx))/(dx)^2;
 fourth_order_error_Left_skewed=abs(Left_skewed_fouth_order-analytical_second_derivative);

For the right skewed :

function fourth_error_Right_skewed= fourth_order_approximation_Right_skewed(x,dx)
A=[1 1 1 1 1 1; 0 1 2 3 4 5; 0 1/2 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]=A\B';

% to compare analytical & numerical derivative exp(x)*cos(x)
analytical_second_derivative= -2*exp(x)*sin(x);
% taking the first order in forward diffrencing
Right_skewed_fouth_order= (X(1)*exp(x)*cos(x)+X(2)*exp(x+dx)*cos(x+dx)+X(3)*exp(x+2*dx)*cos(x+2*dx)+X(4)*exp(x+3*dx)*cos(x+3*dx)+X(5)*exp(x+4*dx)*cos(x+4*dx)+X(6)*exp(x+5*dx)*cos(x+5*dx))/(dx)^2;
 fourth_error_Right_skewed=abs(Right_skewed_fouth_order-analytical_second_derivative);

The plot between the different schemes can be figured as:

On evaluating the above plot we can see that for greater mesh spacing CDS scheme has more error whereas Right and Left schemes do have a lesser error. On lower mesh spacing the CDS scheme has more error, whereas the Right and Left schemes have similar error profiles though lesser, with little variance.

The CDS scheme uses 5 terms, P+Q-1 so has a higher error since there are lesser points in consideration, whereas the skewed schemes have 6 terms P+Q, hence lesser truncation error, though care is taken to avoid roundoff error as much as possible.

Skewed schemes are better than CDS schemes because it becomes particularly useful when the flow direction is known eg. ( left to right) since skewed schemes take information in a particular direction. On the other hand CDS, the scheme takes information from the left and right points. 

The derivation of the schemes can be shown as:

d^2(u)/dx^2= a*f(i-2)+b*f(i-1)+c*f(i)+d*f(i+1)+e*f(i+2) (This is for cds scheme)

TERMS           f(i)     delta*X*f'(x)    delta^2*f''(x)   delta^3*f'''(x)    delta^4*f''''(x)     delta^5*f'''''(x)

a*f(i-2)          a             -2a                  2a                      -8/6*a             16/24*a            -32/120*a

b*f(i-1)          b               -b                   b/2                       -b/6                b/24                  -b/120

c*f(i)              c                 0                    0                             0                  0                         0

d*f(i+1)          d                 d                    d/2                          d/6               d/24                    d/120

e*f(i+2)           e                 2e                   2e                         8/6*e               16/24*e      32/120*e 

Skewed right:

d^2(u)/dx^2= a*f(i)+b*f(i+1)+c*f(i+2)+d*f(i+3)+e*f(i+4)+f*f(i+5)

TERMS           f(i)     delta*X*f'(x)    delta^2*f''(x)   delta^3*f'''(x)    delta^4*f''''(x)     delta^5*f'''''(x)

a*f(i)            a                 0                      0                        0                      0                       0

b*f(i+1)          b               b                  b/2                       b/6                 b/24                   b/120

c*f(i+2)           c                 2c                  2c                         8c/6                16c/24               32c/120

d*f(i+3)          d               3d                  9d/2                        27d/6             81 d/24            243d/120

e*f(i+4)           e               4e                   8e                         64/6*e           256/24*e      1024/120*e 

 f*f(i+5)           f               5f                 25f/2                       125f/6             625f/24        3125f/120

Left skewed:

d^2(u)/dx^2= af(i)+bf(i-1)+cf(i-2)+df(i-3)+ef(i-4)+f*f(i-6) (left skewed scheme)

TERMS           f(i)     delta*X*f'(x)    delta^2*f''(x)   delta^3*f'''(x)    delta^4*f''''(x)     delta^5*f'''''(x)

a*f(i)            a                 0                      0                        0                      0                       0

b*f(i-1)          b               -b                  b/2                       -b/6                 b/24                  -b/120

c*f(i-2)           c               -2c                  2c                        -8c/6                16c/24              -32c/120

d*f(i-3)          d               -3d                  9d/2                      -27d/6             81 d/24           -243d/120

e*f(i-4)           e               -4e                  8e                       -64/6*e           256/24*e      -1024/120*e 

 f*f(i-5)           f               -5f                 25f/2                     -125f/6             625f/24        -3125f/120

 These coefficients are derived through the Taylor series and then matrix inversion in MATLAB is used to compute the required coefficients.