SOLVING ODEs USING TAYLOR TABLE METHOD IN MATLAB

AIM:

To derive fourth order approximation of a 2nd order derivative of a function using:

i.Central differencing

ii.Skewed right sided differencing

iii.Skewed left sided differencing

Fourth order approximations:

function: u(x)=exp(x)*cos(x)

1. SKEWED RIGHT SIDED DIFFERENCING SCHEME:

In this scheme, the derivative is calculated by taking into consideration points from the right side of the boundary.This type of differencing scheme is particularly used to obtain higher order accuracy of derivatives when the point considered is at the left boundary.

∂2f/∂x2=af(i)+bf(i+1)+cf(i+2)+df(i+3)+ef(i+4)+f*f(i+5)

where a,b,c,d,e and f are the coefficients we have to determine by solving the taylor table method.

The terms in the eqn above are expanded with the help of Taylor's table:

f(i+1)=f(i)+dx*f'(i)+(dx^2/2!)*f''(i)+(dx^3/3!)*f'''(i)+(dx^4/4!)*f''''(i)+.........

f(i+2)=f(i)+2*dx*f'(i)+(2*dx)^2/2!*f''(i)+(2*dx)^3/3!*f'''(i)+(2*dx)^4/4!*f''''(i)+.........

.

.

f(i+5)=f(i)+5*dx*f'(i)+(5*dx)^2/2!*f''(i)+(5*dx)^3/3!*f'''(i)+(5*dx)^4/4!*f''''(i)+.........

taylor table:

equations:

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

0+b+2c+3d+4e+5f=0

0+b/2+2c+9/2d+8e+25/2f=1 (since we have to find second order dervative)

0+b/6+4/3c+9/2d+32/3e+125/6f=0

0+b/24+2/3c+27/8d+32/3e+625/24f=0

0+b/120+4/15c+81/40d+128/15e+625/24f=0

Solving the equations in matlab we get:

a=3.75

b=-12.833

c=17.8333

d=-13

e=5.0833

f=-0.8333

The equation finally comes to:

∂2u/∂x2=[3.75f(i)-12.833f(i+1)+17.8333f(i+2)-13f(i+3)+5.0833f(i+4)-0.8333*f(i+5)]/(dx)^2

The above equation is for the fourth order approximation of second order derivative.

The snippet for the right sided scheme is:

%code for right sided skewed differencing scheme

function error1=right_error(x,dx)
analytical_derivative=-2*exp(x)*sin(x);
a=3.75;
b=-12.8333;
c=17.8333;
d=-13;
e=5.0833;
f=-0.8333;
B=[a b c d e f];
right_skewed=(B(1)*exp(x)*cos(x)+B(2)*exp(x+dx)*cos(x+dx)+B(3)*exp(x+(2*dx))*cos(x+(2*dx))+B(4)*exp(x+(3*dx))*cos(x+(3*dx))+B(5)*exp(x+(4*dx))*cos(x+(4*dx))+B(6)*exp(x+(5*dx))*cos(x+(5*dx)))/(dx)^2;
error1=abs(right_skewed-analytical_derivative);
end

2. SKEWED LEFT SIDED DIFFERENCING SCHEME:

In this scheme, the derivative is calculated by taking into consideration points from the left side of the boundary.This type of differencing scheme is particularly used to obtain higher order accuracy of derivatives when the point considered is at the right boundary.

∂2f/∂x2=f*f(i-5)+ef(i-4)+df(i-3)+cf(i-2)+bf(i-1)+a*f(i)

where a,b,c,d,e and f are the coefficients we have to determine by solving the taylor table method.

The terms in the eqn above are expanded with the help of Taylor's table.

Similarly after solving the taylor's series:

taylor table:

equations:

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

0-b-2c-3d-4e-5f=0

0+b/2+2c+9/2d+8e+25/2f=1 (since we have to find second order dervative)

0-b/6-4/3c-9/2d-32/3e-125/6f=0

0+b/24+2/3c+27/8d+32/3e+625/24f=0

0-b/120-4/15c-81/40d-128/15e-625/24f=0

Solving the equations in matlab we get:

a=3.75

b=-12.833

c=17.8333

d=-13

e=5.0833

f=-0.8333

The equation finally comes to:

∂2u/∂x2=[3.75f(i)-12.833f(i-1)+17.8333f(i-2)-13f(i-3)+5.0833f(i-4)-0.8333*f(i-5)]/(dx)^2

The above equation is for the fourth order approximation of second order derivative.

The snippet for the left sided skewed scheme:

%code for the skewed left sided differencing scheme

function error1=left_error(x,dx)
analytical_derivative=-2*exp(x)*sin(x);
a=3.75;
b=-12.8333;
c=17.8333;
d=-13;
e=5.0833;
f=-0.8333;
left_skewed=(f*exp(x-(5*dx))*cos(x-(5*dx))+e*exp(x-(4*dx))*cos(x-(4*dx))+d*exp(x-(3*dx))*cos(x-(3*dx))+c*exp(x-(2*dx))*cos(x-(2*dx))+b*exp(x-dx)*cos(x-dx)+a*exp(x)*cos(x))/(dx)^2;
error1=abs(left_skewed-analytical_derivative);
end

3.CENTRAL DIFFERENCING SCHEME:

In this scheme, the derivative is calculated by taking into consideration points from both the sides of the point concerned.This type of differencing scheme is particularly used to obtain higher order accuracy of derivatives when the point considered inside the domain.

∂2f/∂x2=a*f(i-2)+bf(i-1)+cf(i)+df(i+1)+ef(i+2)

where a,b,c,d and e are the coefficients we have to determine by solving the taylor table method.

The terms in the eqn above are expanded with the help of Taylor's table.

Similarly after solving the taylor's series:

taylor table:

equations:

a+b+c+d+e=0

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

2a+b/2+0+d/2+2e=1 (to find the second order derivative)

-8a/6-b/6+0+d/6+8e/6=0

16a/24+b/24+0+d/24+16e/24=0

Solving the equations in matlab we get:

a=-0.0833

b=1.3333

c=-2.5

d=1.3333

e=-0.0833

The equation finally comes to:

∂2u/∂x2=[-0.0833f(i-2)+1.3333f(i-1)-2.5f(i)+1.3333f(i+1)-0.0833f(i+2)]/(dx)^2

The above equation is for the fourth order approximation of second order derivative.

The snippet for the central differencing scheme:

%code for central differencing scheme

function error6=central_error(x,dx)
analytical_derivative=-2*exp(x)*sin(x);
a=-0.0833;
b=1.3333;
c=-2.5;
d=1.3333;
e=-0.0833;
central_difference=(a*exp(x-(2*dx))*cos(x-(2*dx))+b*exp(x-dx)*cos(x-dx)+c*exp(x)*cos(x)+d*exp(x+dx)*cos(x+dx)+e*exp(x+(2*dx))*cos(x+(2*dx)))/(dx)^2;
error6=abs(central_difference-analytical_derivative);
end

In the previous schemes the program has been written to calculate the error between the numerically derived schemes and the exact analytical derivative to find out which derivative scheme has the least error in computing the derivative of the function.

The error differences were also compared with the forward and backward differencing schemes to show that the central and skewed right sided and left sided schemes had much lesser error in the final computed result than the forward and backward differencing scheme.

The program snippet for FDS scheme:

%code for forward differencing scheme

function error_forward=fds_error(x,dx)
 
analytical_dervative=-2*exp(x)*sin(x);
fds_approximation=((exp(x+(2*dx))*cos(x+(2*dx)))-(2*exp(x+dx)*cos(x+dx))+(exp(x)*cos(x)))/(dx)^2;
error_forward=abs(fds_approximation-analytical_dervative);
end

The program snippet for BDS scheme:

%code for backward differencing scheme

function error_backward=bds_error(x,dx)
 
analytical_derivative=-2*exp(x)*sin(x);
bds_approximation=(exp(x-(2*dx))*cos(x-(2*dx))-2*exp(x-dx)*cos(x-dx)+exp(x)*cos(x))/(dx)^2;
error_backward=abs(bds_approximation-analytical_derivative);
end

The main program code which compares the fourth order error values of different numerical schemes:

clear all 
close all
clc
  
x=pi/3;
dx=linspace(pi/4,pi/4000,20);

 for i=1:length(dx)
 central_diff_error(i)=central_error(x,dx(i));
 right_skewed_error(i)=right_error(x,dx(i));
left_skewed_error(i)=left_error(x,dx(i));

forward_error(i)=fds_error(x,dx(i));
backward_error(i)=bds_error(x,dx(i));
 end
 
 loglog(dx,central_diff_error,'color','b','linewidth',3)
 hold on
loglog(dx,right_skewed_error,'color','r','linewidth',3)
 hold on
loglog(dx,left_skewed_error,'color','g','linewidth',3)
hold on
loglog(dx,forward_error,'color','black','linewidth',3)
hold on
loglog(dx,backward_error,'color','magenta','linewidth',3)
xlabel('dx')
ylabel('errors')
title('comparison of errors in central,right sided and left sided skewed with fds and bds scheme')
legend('central diffrence error','right sided skewed error','left sided skewed error','forward error','backward error')

 

The plot comparing the error values for a range of 'dx' of different numerical schemes:

In the plot it can be seen that for the same dx central error at the end of range had the least error while the forward and backward differencing scheme had the highest error.

CONCLUSION:

It can be hence concluded that the skewed scheme as it involves more number of points, it can have much more information than the fds and bds scheme, and hence the error is also less.

In the case of central scheme, it uses information from all the sides and since second order derivative is elliptic in nature; its region of influence and region of dependence are spreaded through the domain.

Hence it uses information from all sides and has got the least error in computing the derivative.