Taylor Table - 4th Order Approximation of the 2nd Order Derivative

AIM: To derive the 4th approximations of the 2nd derivative using central differencing, skewed right difference and skewed left difference and comparing the truncation errors obtained from the same.

1. Central Differencing : No. of points required = (4+2)-1 = 5.

The following points in the neighbourhood of \'i\' are considered for the required approximation. A CDS uses symetric points on either side.

`(delta^2f)/(deltax^2) = A*f(i-2)+B*f(i-1)+C*f(i)+D*f(i+1)+E*f(i+2)`

The taylor table method involves expanding the terms above in terms of their taylor series expansion & tabulating the coefficients of terms by grouping them on the basis of their order. A set of simultaneous equations needs to be solved with the constaint that we require only the 2nd derivative, hence the sum of coefficients of `f\'(x) , f^3(x),f^4(x) and f^5(x) =0`and sum of coefficients of `f\'\'(x) = 1`.

The set of linear equations is solved to find the 5 unknowns and tabulated as below:

Skewed Right Side Differencing : No. of points required = (4+2) = 6.

The numerical stencil is as shown above.

`(delta^2f)/(deltax^2) = A*f(i)+B*f(i+1)+C*f(i+2)+D*f(i+3)+E*f(i+4)+F*f(i+5)`

Taylor table :

Skewed Left Side Differencing : No. of points required = 4+2 = 6.

`(delta^2f)/(deltax^2) = A*f(i-5)+B*f(i-4)+C*f(i-3)+D*f(i-2)+E*f(i-1)+F*f(i)`

The set of equations are solved in Matlab so as to obtain the values of the unknown variables for each of the schemes.

function [X,Y,Z] = Solving_equations()

% This function returns the solution matrix for CDS, Skewed Right and
% Skewed Left schemes after solving the equations to obtain Taylor
% coefficients. 

% Central Difference scheme.
syms a b c d e % Defining variables of the system of equations. 
% Defining the equations. 
eq1 = a+b+c+d+e==0;
eq2 = -2*a-b+d+2*e==0;
eq3 = 2*a+0.5*b+0.5*d+2*e==1;
eq4 = -(4/3)*a-(b/6)+(d/6)+(4/3)*e==0;
eq5 = (2/3)*a+(b/24)+(d/24)+(2/3)*e==0;
% Converting the equations to matrix form.
[A,B] = equationsToMatrix([eq1,eq2,eq3,eq4,eq5], [a,b,c,d,e]);
% Solving the equations in matrix form simultaneously. 
X=linsolve(A,B);

% Skewed Right sided difference scheme. 
syms a b c d e f % Defining variables of the system of equations. 
% Defining the equations. 
eq1 = a+b+c+d+e+f==0;
eq2 = b+(2*c)+(3*d)+(4*e)+(5*f)==0;
eq3 = 0.5*b+2*c+4.5*d+8*e+12.5*f==1;
eq4 = (b/6)+(8*c/6)+(27/6)*d+(64/6)*e+(125/6)*f==0;
eq5 = (b/24)+(2*c/3)+(81/24)*d+(256/24)*e+(625*f/24)==0;
eq6 = (b/120)+(32*c/120)+(243*d/120)+(1024*e/120)+(3125*f/120)==0;
% Converting the equations to matrix form.
[C,D] = equationsToMatrix([eq1,eq2,eq3,eq4,eq5,eq6], [a,b,c,d,e,f]);
% Solving the equations in matrix form simultaneously. 
Y=linsolve(C,D);

% Skewed Left sided difference scheme. 
syms a b c d e f% Defining variables of the system of equations. 
% Defining the equations. 
eq1 = a+b+c+d+e+f==0;
eq2 = -5*a-(4*b)-(3*c)-(2*d)-e==0;
eq3 = 12.5*a+(8*b)+(4.5*c)+(2*d)+(0.5*e)==1;
eq4 = -(125*a/6)-(64*b/6)-(27/6)*c-(8/6)*d -(e/6)==0;
eq5 = (625*a/24)+(256*b/24)+(81*c/24)+(16*d/24)+(e/24)==0;
eq6 = (-3125*a/120)-(1024*b/120)-(243*c/120)-(32*d/120)-(e/120)==0;
% Converting the equations to matrix form.
[E,F] = equationsToMatrix([eq1,eq2,eq3,eq4,eq5,eq6], [a,b,c,d,e,f]);
% Solving the equations in matrix form simultaneously. 
Z=linsolve(E,F);

end

The out matrices X,Y,Z are as below:

2. MATLAB CODE TO COMPARE TRUNCATION ERRORS

% TAYLOR TABLE & MATLAB CODE CHALLENGE : 
% 4th Order Approximation of 2nd Order Derivative.
% Given function : exp(x)*cos(x)
clear all;
close all; 
clc
x = pi/3; 
dx = linspace(pi/5,pi/200,30); % Range of values for delta x.

[X,Y,Z] = Solving_equations(); % Calling the function that gives coefficients of Taylor table for CDS, Skewed Right & Skewed Left approximations. 

analytical_derivative = -2*exp(x)*sin(x); % Using product rule & differentiating twice to get 2nd derivative.

% Running a loop to compute truncation error at every value of delta x for
% various numerical approximation schemes. 
for i=1:length(dx)
    % Using taylor table and formulae for the FDM schemes; Computing for every \'dx\' and storing in an array. 
    numerical_derivative_central_difference(i)=[(X(1)*exp(x-2*dx(i))*cos(x-2*dx(i)))+ (X(2)*exp(x-dx(i))*cos(x-dx(i))) + (X(3)*exp(x)*cos(x)) + (X(4)*exp(x+dx(i))*cos(x+dx(i))) + (X(5)*exp(x+2*dx(i))*cos(x+dx(i))) ]/(dx(i)^2);
    numerical_derivative_skewedRight_difference(i)= [(Y(1)*exp(x)*cos(x))+ (Y(2)*exp(x+dx(i))*cos(x+dx(i))) + (Y(3)*exp(x+2*dx(i))*cos(x+2*dx(i))) + (Y(4)*exp(x+3*dx(i))*cos(x+3*dx(i))) + (Y(5)*exp(x+4*dx(i))*cos(x+4*dx(i))) + Y(6)*(exp(x+5*dx(i))*cos(x+5*dx(i))) ] /(dx(i)^2);
    numerical_derivative_skewedLeft_difference(i)=[ (Z(1)*exp(x-5*dx(i))*cos(x-5*dx(i))) + (Z(2)*exp(x-4*dx(i))*cos(x-4*dx(i))) + (Z(3)*exp(x-3*dx(i))*cos(x-3*dx(i))) + (Z(4)*exp(x-2*dx(i))*cos(x-2*dx(i))) + (Z(5)*exp(x-dx(i))*cos(x-dx(i))) + Z(6)*exp(x)*cos(x) ] /(dx(i)^2);
    numerical_derivative_forward_difference(i) = [(exp(x+dx(i))*cos(x+dx(i)))-(exp(x)*cos(x))]/dx(i);
    numerical_derivative_backward_difference(i) = [exp(x)*cos(x) - (exp(x-dx(i))*cos(x-dx(i)))]/ dx(i);
    % Computing Error terms for every \'dx\' and storing them in an array.
    error_central_difference(i) = abs(analytical_derivative - numerical_derivative_central_difference(i));
    error_skewedRight_difference(i)= abs(analytical_derivative - numerical_derivative_skewedRight_difference(i));
    error_skewedLeft_difference(i) = abs(analytical_derivative - numerical_derivative_skewedLeft_difference(i));
    error_forward_difference(i)= abs(analytical_derivative - numerical_derivative_forward_difference(i));
    error_backward_difference(i) = abs(analytical_derivative - numerical_derivative_backward_difference(i));
end

% Log log plots for comparing truncation error obtained from various schemes. 
figure(1)
loglog(dx,error_central_difference,\'linewidth\',1.5,\'color\',\'y\');
hold on;
loglog(dx,error_skewedRight_difference,\'linewidth\',1.5,\'color\',\'m\');
hold on;
loglog(dx,error_skewedLeft_difference,\'linewidth\',1.5,\'color\',\'c\');
% Adding legend, x & y axis labels. 
legend(\'Central Difference\',\'Skewed Right Side Difference\',\'Skewed Left Side Difference\');
xlabel(\'\\deltax\');
ylabel(\'Truncation Error\');

% Bar graph to compare values of truncation error at a particular value of \'dx\'.
figure(2)
bardata=[error_central_difference(10),error_skewedRight_difference(10),error_skewedLeft_difference(10),error_forward_difference(10),error_backward_difference(10)];
barlabel={\'CDS\';\'Skew Right\';\'Skew Left\';\'FDS\';\'BDS\'};
bar(bardata,\'facecolor\',\'k\');
set(gca,\'xticklabel\',barlabel)
legend(\'Errors from different approximation schemes\');
ylabel(\'Truncation Error\');
xlabel(\'Numerical Approximation Schemes\');
title(\'Bar chart comparing truncation error values at a particular \\deltax\');

OUTPUT :

The loglog plot comparing the truncation error of the 4th order central and skewed right and left approximation schemes shows us that the central difference error is bigger than that of the other two schemes.

As the value of \'dx\' increases, the truncation errors also increase.

(b) As seen from the bar graph above, CDS, Skewed right and skewed left approximation (4th order ) methods have a truncation error less than that of the first - order FDS & BDS schemes.

(c) A skewed scheme is particularly useful when we are evaluating / approximating corner nodes in a mesh. A CDS scheme requires the mesh to be symmetrical and we might not always have grid points that satisfy this requirement.