Week 3 - Taylor table method and Matlab code

Problem Statement: To derive the numerical 4th order accurate schemes for the second-order derivative using

1) A left sided skewed stencil

2) A right sided skewed stencil

3) Symmetric Scheme

and plot their respective absolute errors.

Procedure:

For the skewed schemes, the number of points in the stencil so that the scheme is 4th order accurate is given by the formula `(N+P)`, where `N` is the order of the derivative and `P` is the order of accuracy. Thus, `2+4 = 6` points must be chosen.

For the symmetric scheme, the number of points in the stencil so that the scheme is 4th order accurate is given by the formula `(N+P-1)`, where `N` is the order of the derivative and `P` is the order of accuracy. Thus, `2+4-1 = 5` points must be chosen.

Now, consider the example of the right-sided skewed scheme with points `(i, i+1, i+2, i+3, i+4, i+5)`, a hypothesis function can then be obtained: `(trianglex)^2(d^2y)/(dx^2) = (af(i) + bf(i+1) + cf(i+2) + df(i+3) + ef(i+4) + gf(i+5)``

Now, once the coefficients for this hypothesis equation are determined, that would provide the scheme for the second order derivative that is 4th order accurate, using a right sided skewed scheme.

The co-efficients can be obtained using the Taylor Table method.

Firstly, each term in the hypothesis equation is expanded using Taylor Series expansion :

`af(i) = a[f(i) + 0trianglexf'(i) + (0trianglex^2f''(i))/(2!) + (0trianglex^3f'''(i))/(3!) + (0trianglexf''''(i))/(4!) + (0trianglexf'''''(i))/(5!)]`

`bf(i+1) = b[f(i) + trianglexf'(i) + (trianglex^2f''(i))/(2!) + (trianglex^3f'''(i))/(3!) + (trianglexf''''(i))/(4!) + (trianglexf'''''(i))/(5!)]`

`cf(i+2) = c[f(i) + 2trianglexf'(i) + (2trianglex^2f''(i))/(2!) + (2trianglex^3f'''(i))/(3!) + (2trianglexf''''(i))/(4!) + (2trianglexf'''''(i))/(5!)]`

`df(i+3) = d[f(i) + 3trianglexf'(i) + (3trianglex^2f''(i))/(2!) + (3trianglex^3f'''(i))/(3!) + (3trianglexf''''(i))/(4!) + (3trianglexf'''''(i))/(5!)]`

`ef(i+4) = e[f(i) + 4trianglexf'(i) + (4trianglex^2f''(i))/(2!) + (4trianglex^3f'''(i))/(3!) + (4trianglexf''''(i))/(4!) + (4trianglexf'''''(i))/(5!)]`

`gf(i+5) = g[f(i) + 5trianglexf'(i) + (5trianglex^2f''(i))/(2!) + (5trianglex^3f'''(i))/(3!) + (5trianglexf''''(i))/(4!) + (5trianglexf'''''(i))/(5!)]`

Now, once the co-efficients have been calculated, they can be tracked by placing them inside the Taylor Table.

Now, since only the second derivative term must exist, the remaining columns can be added up to sum up to 0.

Hence, the following matrices will be obtained.

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

`x = [[a],[b],[c],[d],[e],[g]]`

`b = [[0],[0],[1],[0],[0],[0]]`

Now, these matrices can be solved using matrix inversion to give the corresponding values of the co-efficients that constitute the hypothesis equation.

`ax = b`

`x = a^-1b`

This process can be repeated for the left sided stencil, and the symmetric stencil.

Why are these schemes 4th order accurate?

Since we have expanded each term in the hypothesis equation to the 5th order derivative using Fourier series expansion, all terms from the 6th order derivative are truncated, we get `(d^2y)/(dx^2) = (...)f''(x)(trianglex)^2 + f(O(trianglex)^6)`. Now, on dividing both sides of the hypothesis equation by `(trianglex)^2`, we get a truncation error of the order `(trianglex)^4`

Once the schemes have been developed, the program can be written to calculate the absolute errors of the different schemes with the analytical solution of the given function.

% Making running the code multiple times eeasy
clc
clear all
close all

% Using Taylor Table Method done by hand calculations, matrices are inputted here - RSS
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];
b = [0;0;1;0;0;0];

% Using Taylor Table Method - CDS
c = [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];
d = [0;0;1;0;0];

% % Using Taylor Table Method - LSS
e = [1 1 1 1 1 1; -5 -4 -3 -2 -1 0; 25/2 8 9/2 2 1/2 0; -125/6 -64/6 -27/6 -8/6 -1/6 0; 625/24 256/24 81/24 16/24 1/24 0; -3125/120 -1024/120 -243/120 -32/120 -1/120 0];
f = [0;0;1;0;0;0];

% Matrix Inversion
z = inv(a)* b;
y = inv(c)* d;
m = inv(e)* f;

% Creating dx array for plotting
%dx=linspace(pi/40,pi/4000,200);
dx=linspace(10^(-2),10^(-0.05),200);
x = pi/3;

% Analytical Soln
Analytical = -2 * exp(x) * sin(x);

% Looping through the values in dx, using the derived schemes
for i=1:length(dx)
%Right sided skewed scheme:
    RSS(i) = (z(1)* u(x) + z(2) * u(x+dx(i)) + z(3)* u(x+(2*dx(i))) + z(4)* u(x+(3*dx(i))) +z(5)* u(x+(4*dx(i))) + z(6)* u(x+(5*dx(i))))/(dx(i))^2;
    
%Central Diff Scheme    
    CDS(i) = (y(1)* u(x-(2*dx(i))) + y(2)* u(x-dx(i)) + y(3)* u(x) + y(4)* u(x+dx(i)) + y(5)* u(x+(2*dx(i))))/(dx(i)^2);
    
% Left sided skewed scheme:
    LSS(i) = (m(1)* u(x-(5*dx(i))) + m(2)* u(x-(4*dx(i))) + m(3)* u(x-(3*dx(i))) + m(4)* u(x-(2*dx(i))) + m(5)* u(x-(dx(i))) + m(6)* u(x))/(dx(i))^2;

% Calculating Absolute errors
    A(i) = abs(Analytical - RSS(i));
    B(i) = abs(Analytical - LSS(i));
    C(i) = abs(Analytical - CDS(i));
   
end

% Plotting using loglog because, the errors can be visualised better with
% loglog plot
loglog(dx,A,'LineWidth',3)
hold on
loglog(dx,B,'LineWidth',3)
loglog(dx,C,'LineWidth',3)

legend('Right Skewed', 'Left Skewed','Central Symmetric')
xlabel('dx')
ylabel('Abs Error')
hold off

Output: dx 10^0 to 10^-2

When the dx range is limited to 10^-2 from 10^0, we see a clear trend. That is, The central schemes are more accurate by a significant factor.

The growth in error for the right skewed scheme is greater that that for the left symmetric scheme, as dx increases. These errors are primarily due to the effect of truncation errors becoming more prominent as the dx size increases.

There is also a very signigicant growth in absolute errors as dx becomes bigger.

Output: When dx is 10^0 to 10^-4

It might seem logical that the absolute errors would go on decreasing as dx continues to decrease. However, there is a value for dx, where the absolute errors are at a minima. Then they start to increase with decrease in dx.

This is due to round-off errors in MATLAB, wherein MATLAB rounds off the small dx values as it is capable of storing the data for dx with 16bit storage only. In case, there are too many decimals involved, MATLAB would round it off.

One advantage of skewed schemes over a symmetric scheme are at nodes that lie on the boundary of a domain. A central scheme requires nodes on either side of the current scheme. This condition is not met at boundaries.