Week 3 - Taylor table method and Matlab code

Taylor’s Table method for higher order derivatives

Taylor’s table method is used here to calculate the 4^th order approximations of the second order derivative. The following types of 4^th order approximations are considered.

• Central difference
• Skewed right hand difference
• Skewed left handed difference

Derivation of Central difference Scheme using Taylors Table:

The first task is to determine the number of nodes (n) required for this scheme. Since it’s a central difference scheme,

`n=p+q-1`

Where, p = order of derivative

q = order of approximation

`∴n=2+4-1=5`

The numerical stencil:-

dx dx dx dx

-----o--------o-------o-------o-------o-----

i-2 i-1 i i+1 i+2

`Δx^2 f^'' (x)=a.f_(i-2)+b.f_(i-1)+c.f_i+d.f_(i+1)+e.f_(i+2)`

Where a, b, c, d, e are the coefficients we have to solve using the Taylor’s table.

In the last row of above table the summation of coefficients is taken to zero except the second order derivative term, as we want to find out this particular terms numerical approximation. We dont know the sum of coefficients for fifth and sixth order derivative therfore they are left blank.

Re-arranging the above table in matrix form we get an equation in the form of AX = B;

After solving for `X=A^(-1) B`

We get,

`a= -1/12 ,b=4/3 ,c= -5/2 ,d=4/3 ,e=-1/12`

Determining the order of accuracy:-

Whenever we use numerical derivatives (from Taylor series) to approximate the analytical derivatives, we neglect some terms since Taylor series is a sum of infinite terms. Thus a truncation error is introduced.

For central difference scheme of `Δx^2 f^'' (x)`

rewriting the equation from Taylors table,

`Δx^2 f^'' (x)= a.f_(i-2)+b.f_(i-1)+c.f_i+d.f_(i+1)+e.f_(i+2)+( (-32a)/5!+(-b)/5!+d/5!+32e/5!)Δx^5 f^v (i)+(64a/6!+b/6!+d/6!+64a/6!)Δx^6 f^vi (i)+⋯`

By substituting the values of a,b,c,d and e in above equation the coefficient of gets cancelled and we are left with,

`Δx^2 f^'' (x)= a〖.f〗_(i-2)+b.f_(i-1)+c.f_i+d.f_(i+1)+OΔx^6+⋯`

Dividing by `Δx^2` we get,

`f^'' (x)= a.f_(i-2)+b.f_(i-1)+c.f_i+d.f_(i+1)+OΔx^4`

Where the term `OΔx^4` indicates the order of the approximation is proportional to `x^4` . Thus being 4^th order accurate.

Derivation for Right hand skewed Scheme using Taylors Table:

The first task is to determine the number of nodes (n) required for this scheme. Since it’s a skewed difference scheme,

`n=p+q`

Where, p = order of derivative

q = order of approximation

`∴n=2+4=6`

The numerical stencil:-

dx dx dx dx dx dx

-----o--------o---------o--------o--------o--------o------

i i+1 i+2 i+3 i+4 i+5

`Δx^2 f^'' (x)=a.f_i+b.f_(i+1)+c.f_(i+2)+d.f_(i+3)+e.f_(i+4)+gf_(i+5)`

Where a, b, c, d, e, g are the coefficients we have to solve using the Taylor’s table.

In the last row of above table the summation of coefficients is taken to zero except the second order derivative term, as we want to find out this particular terms numerical approximation. We dont know the sum of coefficients for sixth order derivative therfore it is left blank.

Re-arranging the above table in matrix form we get an equation in the form of AX = B;

After solving for `X=A^(-1) B`

We get,

`a= 15/4 ,b=(-77)/6 ,c= 107/6 ,d=-13,e=61/12 ,g=-5/6`

Determining the order of accuracy:-

For Right hand skewed difference scheme of `Δx^2 f^'' (x)`

rewriting the equation from Taylors table,

`Δx^2 f^'' (x)= a.f_i+b.f_(i+1)+c.f_(i+2)+d.f_(i+3)+e.f_(i+4)+gf_(i+5)+(b/6!+64c/6!+729d/6!+6144e/6!+15625g/6!)Δx^6 f^vi (i)+⋯`

Simplifying and rewritting,

`Δx^2 f^'' (x)= a〖.f〗_i+b.f_(i+1)+c.f_(i+2)+d.f_(i+3)+e.f_(i+4)+gf_(i+5)+OΔx^6+⋯`

Dividing by `Δx^2` we get,

`f^'' (x)= a〖.f〗_i+b.f_(i+1)+c.f_(i+2)+d.f_(i+3)+e.f_(i+4)+gf_(i+5)+OΔx^4`

Where the term `OΔx^4` indicates the order of the accuracy is proportional to `〖 x〗^4`. Thus being 4^th order accurate.

Derivation for Left hand skewed Scheme using Taylors Table:

The first task is to determine the number of nodes (n) required for this scheme. Since it’s a skewed difference scheme,

`n=p+q`

Where, p = order of derivative

q = order of approximation

`∴n=2+4=6`

The numerical stencil:-

dx dx dx dx dx

-----o--------o---------o--------o--------o--------o------

i-5 i-4 i-3 i-2 i-1 i

`Δx^2 f^'' (x)=a.f_i+b.f_(i-1)+c.f_(i-2)+d.f_(i-3)+e.f_(i-4)+gf_(i-5)`

Where a, b, c, d, e, g are the coefficients we have to solve using the Taylor’s table.

In the last row of above table the summation of coefficients is taken to zero except the second order derivative term, as we want to find out this particular terms numerical approximation.We dont know the sum of coefficients for sixth order derivative therfore it is left blank.

Re-arranging the above table in matrix form we get an equation in the form of `AX = B`;

After solving for `X=A^(-1) B`

We get,

`a= 15/4 ,b=(-77)/6 ,c= 107/6 ,d=-13,e=61/12 ,g=-5/6`

Determining the order of accuracy:-

For Left hand skewed difference scheme of `Δx^2 f^'' (x)`

rewriting the equation from Taylors table,

`Δx^2 f^'' (x)= a.f_i+b.f_(i-1)+c.f_(i-2)+d.f_(i-3)+e.f_(i-4)+gf_(i-5)+(b/6!+64c/6!+729d/6!+6144e/6!+15625g/6!)Δx^6 f^vi (i)+⋯`

Simplifying and rewritting,

`Δx^2 f^'' (x)= a.f_i+b.f_(i-1)+c.f_(i-2)+d.f_(i-3)+e.f_(i-4)+gf_(i-5)+OΔx^6+⋯`

Dividing by `Δx^2` we get,

`f^'' (x)= a.f_i+b.f_(i-1)+c.f_(i-2)+d.f_(i-3)+e.f_(i-4)+gf_(i-5)+OΔx^4`

Where the term `OΔx^4` indicates the order of the accuracy is proportional to `x^4` . Thus being 4^th order accurate.

MATLAB CODE FOR COMPARING ERROR OF THE ABOVE MENTIONED SCHEMES

A second order derivative of analytical function `e^x.cos⁡(x)` was compared with right hand skewed, left hand skewed and central difference schemes. Absolute error, between a given scheme and the analytical derivative, was used to compare the errors between the different schemes. The derivatives were found at x = 0

The MATLAB code consists of 1 main program and 3 functions (for each scheme).

Functions for Central diff scheme, Right hand skewed scheme and Left hand skewed scheme.

These codes generate an error value for a particular ‘dx’ term and ‘x’ term. The error value is computed as an absolute difference between the analytical derivative and the respective approximation scheme.

FUNCTION FOR CENTRAL DIFFERENCE SCHEME:-

% Function to calculate the error in 4th order approximation of second order derivate of
% the function exp(x)*cos(x) with a central difference scheme
function out = central_diff_scheme(x,dx)

%Coefficients of central diff scheme
a=-1/12;
b=4/3;
c=-5/2;
d=4/3;
e=-1/12;
% analytical second order derivative
analytical_derivative = -2*exp(x)*(sin(x));
% function
f = @(x)exp(x)*cos(x);
%central difference scheme
Central_difference = (  (a*f(x-2*dx)) + (b*f(x-dx)) + (c*f(x)) + (d*f(x+dx)) + (e*f(x+2*dx)) )/(dx^2);
%error in central differnece scheme
out= abs(analytical_derivative-Central_difference);
end

FUNCTION FOR RIGHT HAND SKEWED SCHEME:-

% Function to calculate the error in 4th order approximation of second order derivate of
% the function exp(x)*cos(x) with a right hand skewed scheme
function out = right_skewed_scheme(x,dx)
%Coefficients of central diff scheme
a=15/4;
b=-77/6;
c=107/6;
d=-13;
e=61/12;
g=-5/6;
% analytical second order derivative
analytical_derivative = -2*exp(x)*(sin(x));
% function
f = @(x)exp(x)*cos(x);
%central difference scheme
Right_skewed_scheme = (  (a*f(x)) + (b*f(x+dx)) + (c*f(x+2*dx)) + (d*f(x+3*dx)) + (e*f(x+4*dx)) + (g*f(x+5*dx)) )/(dx^2);
%error in central differnece scheme
out = abs(analytical_derivative - Right_skewed_scheme  );
end

FUNCTION FOR LEFT HAND SKEWED SCHEME:-

function out = left_skewed_scheme(x,dx)
%Coefficients of central diff scheme
a=15/4;
b=-77/6;
c=107/6;
d=-13;
e=61/12;
g=-5/6;
% analytical second order derivative
analytical_derivative = -2*exp(x)*(sin(x));
% function
f = @(x)exp(x)*cos(x);
%central difference scheme
Left_skewed_scheme = (  (a*f(x)) + (b*f(x-dx)) + (c*f(x-2*dx)) + (d*f(x-3*dx)) + (e*f(x-4*dx)) + (g*f(x-5*dx)) )/(dx^2);
%error in central differnece scheme
out = abs(analytical_derivative - Left_skewed_scheme  );
end

Main MATLAB CODE

This code does the following things:-

• Assigns the value of x and a range of values for dx term
• Calls the functions for different schemes
• Generates an array comprising of error values for different values of ‘dx’ for the different approximation schemes mentioned above
• Plots a log-log graph to compare the error between the different schemes for different values of ‘dx’

%code to compare the errors in the central difference , right and left
%skewed schemes
close all;clear all;clc
% Generating an array for grid spacing points 'dx'
dx = linspace(pi/4,pi/400,20);
% The derivative will be computed at point x = pi/3
x=pi/3;
% Creating a loop to compute the error at each grid point for x = pi/3
for i =1:length(dx)
    %calling the function for central difference,RH skewed and LH skewed scheme scheme to compute error
    %for currnet grid point value
    error_CD(i) = central_diff_scheme(x,dx(i));
    error_RS(i) = right_skewed_scheme(x,dx(i));
    error_LS(i) = left_skewed_scheme(x,dx(i));
end

%plotting the error vs dx
A = figure(1);
%plotting on a log graph to account for difference in the order magnitude
%for different numerical schemes
loglog(dx,error_CD,'marker','o');
hold on
loglog(dx,error_RS,'marker','*');
hold on 
loglog(dx,error_LS,'marker','s','color','k');
legend('Central difference scheme','Right hand skewed scheme','Left hand skewed scheme','location','northwest');
xlabel('dx');
ylabel('Absolute error')
title('Consistency of different schemes');
grid on;
saveas(A,'graph.png');

OUTPUT FROM THE CODE:-

The abscissa shows the values of dx and ordinate depicts the absolute error between the numerical approximation schemes and the analytical derivative. The following things can be concluded from the above figure:-

• The absolute error magnitude decreases as the grid spacing ‘dx’ decreases. Thus it can be inferred that as dx tends to zero , error would vanish. The value of dx for all practical algebraic calculations cannot reach zero. Therfore there will always be some error associated with numerical schemes.
• Error magnitude for the central difference scheme is consistently lower than the other two schemes even though all the schemes are 4^th order accurate. This is beacuse central difference scheme takes value from both sides of the point thus its approximations are better than the other two schemes.
• Left-handed skewed scheme has less error compared to the right-hand skewed scheme.
• As the value of dx reaches zero the maginute of error becomes significantly negligible, thus for practical means the value of any parameter we are measuring will give near-about the same result no matter what scheme we choose.

APPLICATIONS OF SKEWED SCHEMES

The applications of skewed schemes is for boundary nodes where we don’t have nodes either preceding it or succeeding it. Like shown in the below figure where we can use Right-hand skewed scheme for the first node and Left hand skewed scheme for the last node.

1^st Last node

---0----0----0----0----0----0----0----