Taylor table

We need to derive a fourth order approximation for the second order derivative `f_i\'\'`

1. Central difference: Let us consider a five-point stencil consisting of the points `f_(i-2), f_(i-1), f_i, f_(i+1) and f_(i+2)`.

Let us assume `f_i\'\'approx(af_(i-2)+bf_(i-1)+cf_i+df_(i+1)+ef_(i+2))/h^2`. Taylor table for the approximation is given below:

Hence the system of equations become:

`a+b+c+d+e=0`

`-2a-b+d+2e=0`

`4a+b+0+d+e=2`

`-8a-b+0+d+8e=0`

`16a+b+0+d+16e=0`

Using matlab to solve the above system of equations we get:

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

Therefore:

`f_i\'\'approx(-f_(i-2)+16f_(i-1)-30f_i+16f_(i+1)-f_(i+2))/(12h^2)`

2. Skewed right sided difference: Let us consider a six-point stencil consisting of the points `f_(i-1), f_i, f_(i+1), f_(i+2), f_(i+3) and f_(i+4)`.

Let us assume `f_i\'\'approx(af_(i-1)+bf_i+cf_(i+1)+df_(i+2)+ef_(i+3)+gf_(i+4))/h^2`. Taylor table for the approximation is given below:

Hence the system of equations become:

`a+b+c+d+e+g=0`

`-a+0+c+2d+3e+4g=0`

`a+0+c+4d+9e+16g=2`

`-a+0+c+8d+27e+64g=0`

`a+0+c+16d+81e+256g=0`

`-a+0+c+32d+243e+1024g=0`

Using matlab to solve the above system of equations we get:

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

Therefore:

`f_i\'\'approx(10f_(i-1)-15f_i-4f_(i+1)+14f_(i+2)-6f_(i+3)+f_(i+4))/(12h^2)`

3. Skewed left sided difference: Let us consider a six-point stencil consisting of the points `f_(i-4), f_(i-3), f_(i-2), f_(i-1), f_i and f_(i+1)`.

Let us assume `f_i\'\'approx(af_(i-4)+bf_(i-3)+cf_(i-2)+df_(i-1)+ef_i+gf_(i+1))/h^2`. Taylor table for the approximation is given below:

Hence the system of equations become:

`a+b+c+d+e+g=0`

`-4a-3b-2c-d+0+g=0`

`16a+9b+4c+d+0+g=2`

`-64a-27b-8c-d+0+g=0`

`256a+81b+16c+d+0+g=0`

`-1024a-243b-32c-d+0+g=0`

Using matlab to solve the above system of equations we get:

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

Therefore:

`f_i\'\'approx(f_(i-4)-6f_(i-3)+14f_(i-2)-4f_(i-1)-15f_i+10f_(i+1))/(12h^2)`

(a) Bar graph: Errors incurred due to different approximations

In the above plot errors due to different approximations are plotted at x=pi/3.the grid spacing used is h=0.01. It can be seen that the least error is incurred due to central difference scheme followed by skewed left sided difference scheme and skewed right sided difference scheme respectively.

(b) Output window:

As evident from the output, the errors incurred due to both forward difference scheme and backward difference scheme are several orders of magnitude higher than that incurred due to any of the fourth order approximations, for the same value of grid spacing. This is expected, as both forward difference scheme and backward difference scheme are first order approximations and truncation errors increase with decrease in the order of approximation.

(c) Within a computational domain, at the boundary nodes, derivatives can\'t be approximated using a central difference scheme since no nodes are available on one side of the boundary. Hence the derivatives at the boundary nodes are approximated using skewed difference schemes. Additionally, based on the type of PDE, a skewed difference scheme may offer higher numerical stability as compared to a central difference scheme.

function [y1]=exp_x_cos_x(x1)
y1=exp(x1)*cos(x1);
end

% Main program

clc;
clear all;
close all;
x=pi/3;
% Function: f(x)=exp(x)*cos(x)
% Analytical second order derivative: f\"(x)=-2*exp(x)*sin(x)
Analytical_derivative=-2*exp(x)*sin(x);
% 4th order numerical derivatives
h=0.01;
% 1. Central difference:
% f\"(x)~(-f(x-2h)+16*f(x-h)-30*f(x)+16*f(x+h)-f(x+2h))/(12*h^2)
Fourth_order_central_difference=(-exp_x_cos_x(x-2*h)+16*exp_x_cos_x(x-h)-30*exp_x_cos_x(x)+16*exp_x_cos_x(x+h)-exp_x_cos_x(x+2*h))/(12*h^2);
% 2. Skewed right sided difference:
% f\"(x)~(10*f(x-h)-15*f(x)-4*f(x+h)+14*f(x+2h)-6*f(x+3h)+f(x+4h))/(12*h^2)
Fourth_order_skewed_right_sided_difference=(10*exp_x_cos_x(x-h)-15*exp_x_cos_x(x)-4*exp_x_cos_x(x+h)+14*exp_x_cos_x(x+2*h)-6*exp_x_cos_x(x+3*h)+exp_x_cos_x(x+4*h))/(12*h^2);
% 3. Skewed left sided difference:
% f\"(x)~(f(x-4h)-6*f(x-3h)+14*f(x-2h)-4*f(x-h)-15*f(x)+10*f(x+h))/(12*h^2)
Fourth_order_skewed_left_sided_difference=(exp_x_cos_x(x-4*h)-6*exp_x_cos_x(x-3*h)+14*exp_x_cos_x(x-2*h)-4*exp_x_cos_x(x-h)-15*exp_x_cos_x(x)+10*exp_x_cos_x(x+h))/(12*h^2);
% 4. Forward difference:
% f\"(x)~(f(x)-2*f(x+h)+f(x+2h))/h^2
Forward_difference=(exp_x_cos_x(x)-2*exp_x_cos_x(x+h)+exp_x_cos_x(x+2*h))/h^2;
% 5. Backward difference
% f\"(x)~(f(x-2h)-2f(x-h)+f(x))/h^2
Backward_difference=(exp_x_cos_x(x-2*h)-2*exp_x_cos_x(x-h)+exp_x_cos_x(x))/h^2;
% Calcualtion of absolute error
E(1)=abs(Fourth_order_central_difference-Analytical_derivative);
E(2)=abs(Fourth_order_skewed_right_sided_difference-Analytical_derivative);
E(3)=abs(Fourth_order_skewed_left_sided_difference-Analytical_derivative);
E_forward=abs(Forward_difference-Analytical_derivative);
E_backward=abs(Backward_difference-Analytical_derivative);
% Plotting
e=bar(E);
l=cell(1,3);
l{1}=\'Central\';
l{2}=\'Skewed right\';
l{3}=\'Skewed left\';
set(gca,\'xticklabel\',l);
% Output window
fprintf(\'At x=pi/3 and h=%d \\n \\n\',h);
fprintf(\'Error incurred due to fourth order, central difference scheme =%d \\n \\n\',E(1));
fprintf(\'Error incurred due to fourth order, skewed right sided difference scheme =%d \\n \\n\',E(2));
fprintf(\'Error incurred due to fourth order, skewed left sided difference scheme =%d \\n \\n\',E(3));
fprintf(\'Error incurred due to forward difference scheme =%d \\n \\n\',E_forward);
fprintf(\'Error incurred due to backward difference scheme =%d \\n \\n\',E_backward);