Week 3 - Taylor table method and Matlab code

FORWARD DIFFERENCE:

Consider a linear interpolation between the current data value `(t_0,x_0)`and the future data value `(t_1,x_1)`. The slope of the secant line between these two points approximates the derivative by the forward (two-point) difference:

For first derivative: f'(`t_0`)=`((x_1-x_0)/(t_1-t_0))`

For instance, the forward difference above predicts the value of `x_1` from the derivative f'(`t_0`) and from the value `x_0`. If the data values are equally spaced with the step size h, the truncation error of the forward difference approximation has the order of O(h).

BACKWARD DIFFERENCE:

Consider a linear interpolation between the current data value `(t0,x0)` and the past data value `(t_-1,x_-1)`. The slope of the secant line between these two points approximates the derivative by the backward (two-point) difference:

For first derivative: f'(`t_0`)=`((x_0-x_-1)/(t_0-t_-1))`

Backward differences are useful for approximating the derivatives if data in the future are not yet available. Moreover, data in the future may depend on the derivatives approximated from the data in the past . If the data values are equally spaced with the step size h, the truncation error of the backward difference approximation has the order of O(h)

CENTRAL DIFFERENCE:

Consider a linear interpolation between the past data value`(t_-1,x_-1)` and the future data value`(t_1,x_1)` . The slope of the secant line between these two points approximates the derivative by the central (three-point) difference:

For first derivative: f'(`t_0`)=`((x_1-x_-1)/(t_1-t_-1))`

If the data values are equally spaced, the central difference is an average of the forward and backward differences. The truncation error of the central difference approximation is order of O(h²), where h is the step size.

In order to increase the error approximation,number of grid points need to increases.

Order of error approximation = (Number of grid points)-(order of derivative)

Here,

Order of error approximation =4

Order of derivative = 2

4 = Number of grid points - 2

Number of grid points = 6

FORWARD DIFFERENCE:

Let grid points are

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

i

`overset(dx)(hArr)`

i+1

`overset(dx)(hArr)`

i+2

`overset(dx)(hArr)`

i+3

`overset(dx)(hArr)`

i+4

`overset(dx)(hArr)`

i+5

`(d^2f)/(dx^2) * dx^2` = `a*f_i + b*f_(i+1)+ c*f_(i+2)+ d*f_(i+3)+ e*f_(i+4)+ g*f_(i+5)`

To calculate a,b,c,d,e,g values, use TAYLOR TABLE method

	`f(i)`	`f'(i)*(dx/(1!))`	`f^2(i)*(dx^2/(2!))`	`f^3(i)*(dx^3/(3!))`	`f^4(i)*(dx^4/(4!))`	`f^5(i)*(dx^5/(5!))`
`f_i`	a	0	0	0	0	0
`f_(i+1)`	b	b	b	b	b	b
`f_(i+2)`	c	2c	4c	8c	16c	32c
`f_(i+3)`	d	3d	9d	27d	81d	243d
`f_(i+4)`	e	4e	16e	64e	256e	1024e
`f_(i+5)`	g	5g	25g	125g	625g	3125g
	0	0	1	0	0	0

a+b+c+d+e+g=0`rarr(1)`

0+b+2c+3d+4e+5g=0`rarr(2)`

0+b+4c+9d+16e+25g=1`rarr(3)`

0+b+8c+27d+64e+125g=0`rarr(4)`

0+b+16c+81d+256e+625g=0`rarr(5)`

0+b+32c+243d+1024e+3125g=0`rarr(6)`

we have 6 equations and 6 unknowns.After solving equations we get

a=1.87500

b=-6.41667

c=8.91667

d=-6.50000

e=2.54167

g=-0.41667

`(d^2f)/(dx^2) ` = `(1.87500*f_i -6.41667*f_(i+1)+ 8.91667*f_(i+2)-6.50000*f_(i+3)+ 2.54167*f_(i+4)-0.41667*f_(i+5))/(dx^2)`

BACKWARD DIFFERENCE:

Let grid points are

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

i-5

`overset(dx)(hArr)`

i-4

`overset(dx)(hArr)`

i-3

`overset(dx)(hArr)`

i-2

`overset(dx)(hArr)`

i-1

`overset(dx)(hArr)`

i

`(d^2f)/(dx^2) * dx^2` = `a*f_i + b*f_(i-1)+ c*f_(i-2)+ d*f_(i-3)+ e*f_(i-4)+ g*f_(i-5)`

To calculate a,b,c,d,e,g values, use TAYLOR TABLE method

	`f(i)`	`f'(i)*(dx/(1!))`	`f^2(i)*(dx^2/(2!))`	`f^3(i)*(dx^3/(3!))`	`f^4(i)*(dx^4/(4!))`	`f^5(i)*(dx^5/(5!))`
`f_i`	a	0	0	0	0	0
`f_(i-1)`	b	-b	b	-b	b	-b
`f_(i-2)`	c	-2c	4c	-8c	16c	-32c
`f_(i-3)`	d	-3d	9d	-27d	81d	-243d
`f_(i-4)`	e	-4e	16e	-64e	256e	-1024e
`f_(i-5)`	g	-5g	25g	-125g	625g	-3125g
	0	0	1	0	0	0

a+b+c+d+e+g=0`rarr(1)`

0-b-2c-3d-4e-5g=0`rarr(2)`

0+b+4c+9d+16e+25g=1`rarr(3)`

0-b-8c-27d-64e-125g=0`rarr(4)`

0+b+16c+81d+256e+625g=0`rarr(5)`

0-b-32c-243d-1024e-3125g=0`rarr(6)`

we have 6 equations and 6 unknowns.After solving equations we get

a=1.87500

b=-6.41667

c=8.91667

d=-6.50000

e=2.54167

g=-0.41667

`(d^2f)/(dx^2) ` = `(1.87500*f_i -6.41667*f_(i-1)+ 8.91667*f_(i-2)-6.50000*f_(i-3)+ 2.54167*f_(i-4)-0.41667*f_(i-5))/(dx^2)`

CENTRAL DIFFERENCE:

Let grid points are

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

i-2

`overset(dx)(hArr)`

i-1

`overset(dx)(hArr)`

i

`overset(dx)(hArr)`

i+1

`overset(dx)(hArr)`

i+2

`overset(dx)(hArr)`

i+3

`(d^2f)/(dx^2) * dx^2` = `a*f_i + b*f_(i+1)+ c*f_(i-1)+ d*f_(i+2)+ e*f_(i-2)+ g*f_(i+3)`

To calculate a,b,c,d,e,g values, use TAYLOR TABLE method

	`f(i)`	`f'(i)*(dx/(1!))`	`f^2(i)*(dx^2/(2!))`	`f^3(i)*(dx^3/(3!))`	`f^4(i)*(dx^4/(4!))`	`f^5(i)*(dx^5/(5!))`
`f_i`	a	0	0	0	0	0
`f_(i+1)`	b	b	b	b	b	b
`f_(i-1)`	c	-c	c	-c	c	-c
`f_(i+2)`	d	2d	4d	8d	16d	32d
`f_(i-2)`	e	-2e	4e	-8e	16e	-32e
`f_(i+3)`	g	3g	9g	27g	81g	243g
	0	0	1	0	0	0

a+b+c+d+e+g=0`rarr(1)`

0+b-c+2d-2e+3g=0`rarr(2)`

0+b+c+4d+4e+9g=1`rarr(3)`

0+b-c+8d-8e+27g=0`rarr(4)`

0+b+c+16d+16e+81g=0`rarr(5)`

0+b-c+32d-32e+243g=0`rarr(6)`

we have 6 equations and 6 unknowns.After solving equations we get

a=-1.25000

b=0.66667

c=0.66667

d=-0.04167

e=-0.04167

g=0.00000

`(d^2f)/(dx^2) ` =`(-1.25000*f_i + 0.66667*f_(i+1)+ 0.66667*f_(i-1)-0.04167*f_(i+2)-0.04167*f_(i-2)+ 0.00000*f_(i+3))/(dx^2)`

MAIN PROGRAM:

clear all;
close all;
clc;

x=2*pi
dx=linspace(pi/4,pi/4000,20)
for i = 1:length(dx)
[skewed_rightside_difference_error(i),skewed_leftside_difference_error(i),central_difference_error(i)]=taylortablefun(x,dx(i));
end
skewed_rightside_difference_error
skewed_leftside_difference_error
central_difference_error

figure(1)
loglog(dx,skewed_rightside_difference_error,'b','linewidth',2)
hold on
loglog(dx,skewed_leftside_difference_error,'g','linewidth',2)
hold on
loglog(dx,central_difference_error,'r','linewidth',2)
xlabel('dx')
ylabel('error_difference')
legend('skewed_rightside_difference_error','skewed_leftside_difference_error','central_difference_error')
title('dx vs error_difference')

SUPPORTING PROGRAM:

function [skewed_rightside_difference_error,skewed_leftside_difference_error,central_difference_error]=taylortablefun(x,dx)
%%Analytical function
%f(x)=e^x*cosx
%%first derivative
%f'(x)=e^x*cosx-e^x*sinx
%% Second derivative
%f''(x)=-2*e^x*sinx

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


%%4th order approximations of the second-order derivative 

%Skewed right-sided difference(forward_difference)
%(d^2f)/(dx^2)=(1.87500*f_i -6.41667*f_(i+1)+ 8.91667*f_(i+2)-6.50000*f_(i+3)+ 2.54167*f_(i+4)-0.41667*f_(i+5))/(dx^2)

skewed_rightside_difference = (1.87500*exp(x)*cos(x) -6.41667*exp(x+dx)*cos(x+dx)+ 8.91667*exp(x+(2*dx))*cos(x+(2*dx))-6.50000*exp(x+(3*dx))*cos(x+(3*dx))+ 2.54167*exp(x+(4*dx))*cos(x+(4*dx))-0.41667*exp(x+(5*dx))*cos(x+(5*dx)))/(dx^2);

%Skewed left-sided difference(backward_difference)
%(d^2f)/(dx^2)=(1.87500*f_i -6.41667*f_(i-1)+ 8.91667*f_(i-2)-6.50000*f_(i-3)+ 2.54167*f_(i-4)-0.41667*f_(i-5))/(dx^2)

skewed_leftside_difference = (1.87500*exp(x)*cos(x) -6.41667*exp(x-dx)*cos(x-dx)+ 8.91667*exp(x-(2*dx))*cos(x-(2*dx))-6.50000*exp(x-(3*dx))*cos(x-(3*dx))+ 2.54167*exp(x-(4*dx))*cos(x-(4*dx))-0.41667*exp(x-(5*dx))*cos(x-(5*dx)))/(dx^2);

%Central difference
%(d^2f)/(dx^2) =(-1.25000*f_i + 0.66667*f_(i+1)+ 0.66667*f_(i-1)-0.04167*f_(i+2)-0.04167*f_(i-2)+ 0.00000*f_(i+3))/(dx^2)

central_difference = (-1.25000*exp(x)*cos(x)+0.66667*exp(x+dx)*cos(x+dx)+ 0.66667*exp(x-(1*dx))*cos(x-(1*dx))-0.04167*exp(x+(2*dx))*cos(x+(2*dx))-0.04167*exp(x-(2*dx))*cos(x-(2*dx))+0.00000*exp(x+(3*dx))*cos(x+(3*dx)))/(dx^2);

%%error
skewed_rightside_difference_error = abs(skewed_rightside_difference-analytical_derivative);
skewed_leftside_difference_error = abs(skewed_leftside_difference-analytical_derivative);
central_difference_error = abs(central_difference-analytical_derivative);


end

RESULTS:

OBSERVATION:

1. It is clear that the central difference gives a much more accurate approximation of the derivative compared to the forward and backward differences.
2. As dx value decreases, error value also decreases but rapidly.
3. The central difference is an average of the forward and backward differences.Hence if we have both forward skewed difference and backward skewed difference values, we can calculate central difference values.But viveversa is not possible