Higher Order Approximation using Taylor Table Method

Problem Statement : To Derive the following 4th order approximations of the second order derivative of the function f(x) = exp(x)*cos(x)

1. Central difference

2. Skewed right sided difference

3. Skewed left sided difference

1. Provide a plot that compares the absolute error between the above mentioned schemes

2. Compare these errors with FDS and BDS approximations.

3. In your own words, describe why a skewed scheme is useful? What can a skewed scheme do that a CD scheme cannot do?

Explanation : To derive the central difference, skewed right side and skewed left side schemes we first have to come up with the expression to calculate the second order derivative.

For the central differnce scheme the terms to be included should be symmetric about the point at which the function has to be approximated. For central difference scheme the terms are like f(i-2), f(i-1), f(i), f(i+1), f(i+2) and so on.

For the skewed right side scheme the terms only to the right of the point at which the function has to be evaluated are considered. The terms are like f(i), f(i+1), f(i+2), f(i+3) and so on.

For the skewed left side scheme the terms only to the left side of the point at which the fuction has to be evaluated are considered. The terms are like f(i), f(i-1), f(i-2), f(i-3) and so on.

The number of terms to be included in the expression depends upon the 1. the order of accuracy and 2. the order of derivative to be evaluated.

For central difference scheme the number if terms to be included is derived by the expression:

p+q-1 where p is the order of accuracy and q is the order of derivative.

For skewed schems (both left and right) the numcer of terms are derived by the equation

p+q where p is the order of accuracy and q is the order of derivative.

In this challenge the order of accuracy is 4th order and the order of derivative is 2nd order. Hence the number of terms for central difference scheme is p+q-1 = 5.

Hence the formula for Central Difference Scheme is :

d^f/dx^2 = (a*f(i-2) + b*f(i-1) + c*f(i) + d*f(i+1) + e*f(i+2))/h^2

where h is the step size and f(i+1) = f(x+h) and so on.

The number of terms for the skewed right scheme is evaluated as p+q which is equal to 6.

Hence the formula for the skewed right side scheme is :

d^f/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))/h^2

Similarly the number of terms for the skewed left side scheme is 6.

Hence the formula for the skewed left side scheme is :

d^f/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))/h^2

After we have formulated the equation we can individually find the expression for f(i), f(i+1), f(i-1) and so on. And then we can find the values of the coefficients a, b, c, d, e and f.

This above method of approximating the derivatives of the function is called as Taylor Table method.

Code: Below are the codes for the funcitons of CDS, RDS and LDS schemes.

Function code for the CDS scheme.

function f_CDS = CDS_function(x,h)
% Taylor Table Equations
syms a b c d e

eq1 = a + b + c + d + e == 0;
eq2 = -2*a - b + d + 2*e == 0;
eq3 = 2*a + b/2 + d/2 + 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;

%% Matrices from equations

% Coefficients Matrix
A = [1 1 1 1 1;-2 -1 0 1 2;2 1/2 0 1/2 2;-4/3 -1/6 0 1/6 4/3;2/3 1/24 0 1/24 2/3];

% Constants Matrix
C = [0;0;1;0;0];

% Solving the matrices 
B = inv(A)*C;

f1 = (exp(x-(2.*h)).*cos(x-(2.*h)));
f2 = (exp(x-h).*cos(x-h));
f3 = exp(x)*cos(x);
f4 = (exp(x+h).*cos(x+h));
f5 = (exp(x+(2.*h)).*cos(x+(2.*h)));

f_CDS = (B(1).*f1 + B(2).*f2 + B(3)*f3 + B(4).*f4 + B(5).*f5)./(h.^2);
end

This function accepts x and h as an input and evaluates the second order derivative with 4th order accuracy as an output. In this code first the equations have been derived from the taylor table and then those equations have been converted into matrices. If A*X = B and if A is the coefficients matrix and B is the constants matrix and if X is the matrix containing unknowns then X = inv(A)*B will evaluate the unknowns. The same approach has been implemented to calculate the unkowns. Then the unkowns have been put in the final expression.

Below is the function code for the RDS scheme.

function f_RDS = RDS_function(x,h)
%% Taylor Table Equations
syms a b c d e f

eq1 = a + b + c + d + e + f ==0;
eq2 = b + 2*c + 3*d + 4*e + 5*f ==0;
eq3 = b/2 + 2*c + (9/2)*d + 8*e + (25/2)*f == 1;
eq4 = b/6 + (4/3)*c + (9/2)*d + (32/3)*e + (125/6)*f == 0;
eq5 = b/24 + (2/3)*c + (81/24)*d + (256/24)*e + (625/24)*f == 0;
eq6 = b/120 + (4/15)*c + (243/120)*d + (1024/120)*e + (3125/120)*f ==0;

%% Matrices from equations
% Coefficients Matrix
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 4/3 9/2 32/3 125/6; 0 1/24 2/3 81/24 256/24 625/24;0 1/120 4/15 243/120 1024/120 3125/120];

% Constants Matrix
C = [0;0;1;0;0;0];

%% Solving for the variables
B = inv(A)*C;

f1 = exp(x)*cos(x);
f2 = exp(x+h).*cos(x+h);
f3 = exp(x+(2.*h)).*cos(x+(2.*h));
f4 = exp(x+(3.*h)).*cos(x+(3.*h));
f5 = exp(x+(4.*h)).*cos(x+(4.*h));
f6 = exp(x+(5.*h)).*cos(x+(5.*h));

f_RDS = (B(1)*f1 + B(2).*f2 + B(3).*f3 + B(4).*f4 + B(5).*f5 + B(6).*f6)./(h.^2);
end

The same approch as in the previous case has been applied here. The only differnce is in the number of terms (here 6 and in the previous case 5). The function will take x and h as an input and will evaluate the second order derivative with 4th order approximation but will only take the terms to the right of the point in consideration.

Below is the code for the LDS scheme.

function f_LDS = LDS_function(x,h)
%% Taylor Table Equations
syms a b c d e f

eq1 = a + b + c + d + e + f ==0;
eq2 = -b - 2*c - 3*d - 4*e - 5*f ==0;
eq3 = b/2 + 2*c + (9/2)*d + 8*e + (25/2)*f == 1;
eq4 = -b/6 - (4/3)*c - (9/2)*d - (32/3)*e - (125/6)*f == 0;
eq5 = b/24 + (2/3)*c + (81/24)*d + (256/24)*e + (625/24)*f == 0;
eq6 = -b/120 - (4/15)*c - (243/120)*d - (1024/120)*e - (3125/120)*f == 0;

%% Matrices from equations
% Coefficients Matrix
A = [1 1 1 1 1 1;0 -1 -2 -3 -4 -5;0 1/2 2 9/2 8 25/2;0 -1 -4/3 -9/2 -32/3 -125/6; 0 1/24 2/3 81/24 256/24 625/24;0 -1/120 -4/15 -243/120 -1024/120 -3125/120];

% Constants Matrix
C = [0;0;1;0;0;0];

%% Solving for the variables
B = inv(A)*C;

f1 = exp(x)*cos(x);
f2 = exp(x-h).*cos(x-h);
f3 = exp(x-(2.*h)).*cos(x-(2.*h));
f4 = exp(x-(3.*h)).*cos(x-(3.*h));
f5 = exp(x-(4.*h)).*cos(x-(4.*h));
f6 = exp(x-(5.*h)).*cos(x-(5.*h));

f_LDS = (B(1)*f1 + B(2).*f2 + B(3).*f3 + B(4).*f4 + B(5).*f5 + B(6).*f6)./(h.^2);
end

The number of terms here are also equal to 6 and only the terms to the left of the point in consideration have been taken.

Now to compare these values with the analytical solution a seperate code (main code) has been written which provides the values of x and h. It also shows the comparison between the errors of the above schemes as well as with the Forward Difference and Backward Difference approximations.

The exact second order derivative for the function is -2*exp(x)*sin(x).

The formula for the forward differnce (FDS) approximation is :

d^2f/dx^2 = (f(i) - 2*f(i+1) + f(i+2))/h^2

The formula for the backward differnce (BDS) approximation is :

d^2f/dx^2 = (f(i) - 2*f(i-1) + f(i-2))/h^2

Here is the code:

% Computing higher order derivatives with higher order approximations using
% the taylor table method

clear all
close all
clc

%% Constants
x = pi/4; % Considered Point
h = linspace(pi/10,pi/700,1000); % Step Size
analytical_solution = (-2*exp(x)*sin(x));

%% Approximating the derivative of the function using Central Difference Scheme
CDS_approximation = CDS_function(x,h);
CDS_error = abs(analytical_solution - CDS_approximation);
    
%% Approximating the derivative of the function using Forward Difference Scheme
RDS_approximation = RDS_function(x,h);
RDS_error = abs(analytical_solution - RDS_approximation);

%% Approximating the derivative of the function using Backward Difference Scheme
LDS_approximation = LDS_function(x,h);
LDS_error = abs(analytical_solution - LDS_approximation);

%% Forward difference approximation
f1 = exp(x)*cos(x);
f2 = exp(x+h).*cos(x+h);
f3 = exp(x+(2.*h)).*cos(x+(2.*h));
FDS_approximation = (f1 - 2*f2 + f3)./(h.^2);
FDS_error = abs(FDS_approximation - analytical_solution);

%% Backward Difference Appeoximation
f1 = exp(x)*cos(x);
f2 = exp(x-h).*cos(x-h);
f3 = exp(x-(2.*h)).*cos(x-(2.*h));
BDS_approximation = (f1 - 2*f2 + f3)./(h.^2);
BDS_error = abs(BDS_approximation - analytical_solution);

%% Plotting of the results 
% Comparison of errors among the three schemes
figure(1)
loglog(h,RDS_error,\'r\',\'LineWidth\',2);
hold on;
loglog(h,CDS_error,\'b\',\'LineWidth\',2);
hold on;
loglog(h,LDS_error,\'k\',\'LineWidth\',2);
legend(\'RDS Error\',\'CDS Error\',\'LDS Error\',\'location\',\'northwest\');
title(\'Comparison of Error Between Different Discretization Schemes\');
xlabel(\'Step Size (h)\');
ylabel(\'Magnitude of Error\');

% Comparison of errors with the FDS and BDS approximations
figure(2)
loglog(h,CDS_error,\'r\',\'LineWidth\',2);
hold on;
loglog(h,RDS_error,\'b\',\'LineWidth\',2);
hold on;
loglog(h,LDS_error,\'g\',\'LineWidth\',2);
hold on;
loglog(h,FDS_error,\'m\',\'LineWidth\',2);
hold on;
loglog(h,BDS_error,\'k\',\'LineWidth\',2);
title(\'Comparison of errors with FDS and BDS approximations\');
legend(\'CDS Error\',\'RDS Error\',\'LDS Error\',\'FDS Error\',\'BDS Error\',\'location\',\'southeast\');
xlabel(\'Step Size(h)\');
ylabel(\'Magnitude of Error\');

Results : Two plots have been plotted.

1. Comparison of errors among the three schemes (CDS, RDS, BDS).

2. Comparison of errors with the FDS and BDS approximations.

1. Plot 1

Explanation : As we can clearly see that as the step size increases the magnitude of error also increases (positive slope). The error is of the order 10^0 when the steo size is of the order 10^0. The overall error is highest for the skewed left side (LDS) scheme and lowest for the central difference (CDS) scheme. The reason behind this is that in the central difference scheme the terms are symmetric which results in better approximation. Whereas for the skewed schemes the terms are not symmetric hence in any case (left or right) we are not considering the effect of the other sided terms on the current node on our stencil.

2. Plot 2

Explanation: In the above plot we can clearly see that the FDS and BDS approximations are not that accurate as the central scheme or the skewed schemes. This is because the FDS and BDS approximations are first order accurate and the skewed schemes and central schemes are fourth order accurate. The truncation error in FDS and BDS approximation is greater than the CDS or the skewed schemes. Only the skewed left side scheme approximation is comparable to the FDS and BDS approximation.

The skewed schemes are better than the central schemes in many ways:

1. The real world problems where we the values at one end of the boundary is known and we are interested to know the value at the other end the skewed schemes are useful.

2. Problems where the discretization nodes are even in number. the CDS scheme is inaccurate because of loss of symmetry. If the number if terms are 6 (say) then at one end a term less than the other sided term will be taken which will result in loss of symmetry.

3. Though central schemes may provide more accurate results but the skewed schemes will always work in any number of discretization. Skewed schemes are more versatile.

4. Also sometimes when the discretization is large or the equation is complex the central scheme can undergo non-convergent oscillaations.