Fourth order approximations using Taylor Table in MATLAB

Taylor Table:-

By creating the table we can find the taylor coefficients for \"forward, backward, central\" differencing methods with equal spacing.

01: Central differencing method:-

function f(x)

Let f\'\'(x)=a*f(i-2)+b*f(i-1)+c*f(i)+d*f(i+1)+e*f(i+2)

`a*f(i-2)=a*f(i)-2*a*Deltax*f\'(i)+(4*a*(Deltax)^2)/(2!)*f\'\'(i)-(8*a*(Deltax)^3)/(3!)*f\'\'\'(i)+(16*a*(Deltax)^4)/(4!)*f\'\'\'\'(i)`

`b*f(i-1)=b*f(i)-b*Deltax*f\'(i)+(b*(Deltax)^2)/(2!)*f\'\'(i)-(b*(Deltax)^3)/(3!)*f\'\'\'(i)+(b*(Deltax)^4)/(4!)*f\'\'\'\'(i)`

c*f(i)=c*f(i)

`d*f(i+1)=d*f(i)+d*Deltax*f\'(i)+(d*(Deltax)^2)/(2!)*f\'\'(i)+(d*(Deltax)^3)/(3!)*f\'\'\'(i)+(d*(Deltax)^4)/(4!)*f\'\'\'\'(i)`

`e*f(i+2)=e*f(i)+2*e*Deltax*f\'(i)+(4*e*(Deltax)^2)/(2!)*f\'\'(i)+(8*e*(Deltax)^3)/(3!)*f\'\'\'(i)+(16*e*(Deltax)^4)/(4!)*f\'\'\'\'(i)`

Transforming these equations into table

	f(i-2)	f(i-1)	f(i)	f(i+1)	f(i+2)
f(i)	a	b	c	d	e	=0
Δx*f\'(i)	-2*a	-b	0	d	2*e	=0
(Δx)2*f\'\'(i)	(4*a)/2	b/2	0	d/2	(4*e)/2	=1
(Δx)3*f\'\'\'(i)	-(8*a)/3!	-b/3!	0	d/3!	(8*e)/3!	=0
(Δx)4*f\'\'\'\'(i)	(16*a)/4!	b/4!	0	d/4!	(16*e)/4!	=0

As this is a second order equation ==> f\'\'(i)=1

Transforming the table into matrix to find the coefficients of taylor table \"a,b,c,d,e\"

A X B

1	1	1	1	1		a		0
-2	-1	0	1	2		b		0
2	1/2	0	1/2	2	*	c	=	1
-8/6	-1/6	0	1/6	8/6		d		0
16/24	1/24	0	1/24	16/24		e		0

By solving A*X=B values of X matrix can be found

`X=A^-1**B`

the MATLAB CODE for this equation is: `X=\"inv\"(A)**B` or`X=\"B\\A\"`

the equation f\'\'(x)=a*f(i-2)+b*f(i-1)+c*f(i)+d*f(i+1)+e*f(i+2) is solved by using the above found values of \"a,b,c,d,e\"

In MATLAB the above step is repeated for a range of `Deltax` and the results are plotted.

`(Deltax)^2*f\'\'(x)+0*(Deltax)^6=a*f(i-2)+b*f(i-1)+c*f(i)+d*f(i+1)+e*f(i+2)`

dividing the above equation by `(Deltax)^2`

`f\'\'(x)+0*(Deltax)^4=(a*f(i-2)+b*f(i-1)+c*f(i)+d*f(i+1)+e*f(i+2))/(Deltax)^2`

the truncation error is of order 4, Hence it is called 4th order approximation of the second order derivative.

02: Right side(Forward) Differencing method:

function f(x)

Let f\'\'(x)=a*f(i)+b*f(i+1)+c*f(i+2)+d*f(i+3)+e*f(i+4)+f*f(i+5)

matrix to find the coefficients of taylor table for forward differencing is:

A X B

1	1	1	1	1	1		a		0
0	1	2	3	4	5		b		0
0	1/2	4/2	9/2	16/2	25/2	*	c	=	1
0	1/6	8/6	64/6	64/6	125/6		d		0
0	1/24	16/24	81/24	256/24	625/24		e		0
0	1/120	32/120	243/120	1024/120	3125/120		f		0

03: Left side(Backward) Differencing method:

function f(x)

Let f\'\'(x)=a*f(i)+b*f(i-1)+c*f(i-2)+d*f(i-3)+e*f(i-4)+f*f(i-5)

matrix to find the coefficients of taylor table for backward differencing is:

A X B

1	1	1	1	1	1		a		0
0	-1	-2	-3	-4	-5		b		0
0	1/2	4/2	9/2	16/2	25/2	*	c	=	1
0	-1/6	-8/6	-64/6	-64/6	-125/6		d		0
0	1/24	16/24	81/24	256/24	625/24		e		0
0	-1/120	-32/120	-243/120	-1024/120	-3125/120		f		0

For skewed methods we have to take the stencil for (i+5) and (i-5)

The MATLAB code is:

clear all
close all
clc
%AX=B
C=[0 0 1 0 0];%X matrix for CDS
C1=[0 0 1 0 0 0];%X matrix for FDS and BDS
A1=[1     1      1         1        1          1; 
    0     1      2         3        4          5; 
    0     1/2    4/2       9/2     16/2       25/2; 
    0     1/6    8/6      27/6     64/6      125/6; 
    0     1/24  16/24     81/24   256/24     625/24; 
    0     1/120 32/120   243/120 1024/120   3125/120];%A matrix for FDS
A2=[1     1       1        1        1;
   -2    -1       0        1        2;
    2     1/2     0        1/2      2; 
   -8/6  -1/6     0        1/6      8/6; 
    16/24 1/24    0        1/24    16/24];%A matrix for CDS
A3=[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     27/8     32/3      625/24;
    0    -1/120 -32/120 -243/120 -1024/120 -3125/120];% A matrix for BDS
B=A1\\C1\';%B=inv(A1)*C1\'    finding the coefficients of taylor table for FDS
D=A2\\C\';%D=inv(A2)*C\'      finding the coefficients of taylor table for CDS
E=A3\\C1\';%E=inv(A3)*C1\'    finding the coefficients of taylor table for BDS
x=pi/4; %solution at pi/4
dx=linspace(pi/4,pi/4000,50);% range of dx
analytical=-2*exp(x)*sin(x);% analytical solution
for i=1:length(dx)%loop for the range of dx
    skewed_right(i)=(B(1)*exp(x)*cos(x)+B(2)*exp(x+dx(i))*cos(x+dx(i))+B(3)*exp(x+2*dx(i))*cos(x+2*dx(i))+B(4)*exp(x+3*dx(i))*cos(x+3*dx(i))+B(5)*exp(x+4*dx(i))*cos(x+4*dx(i))+B(6)*exp(x+5*dx(i))*cos(x+5*dx(i)))/(dx(i)^2);%skewed_right_side
    cds(i)=(D(1)*exp(x-2*dx(i))*cos(x-2*dx(i))+D(2)*exp(x-1*dx(i))*cos(x-1*dx(i))+D(3)*exp(x)*cos(x)+D(4)*exp(x+1*dx(i))*cos(x+1*dx(i))+D(5)*exp(x+2*dx(i))*cos(x+2*dx(i)))/(dx(i)^2);%CDS
    skewed_left(i)=(E(1)*exp(x)*cos(x)+E(2)*exp(x-1*dx(i))*cos(x-1*dx(i))+E(3)*exp(x-2*dx(i))*cos(x-2*dx(i))+E(4)*exp(x-3*dx(i))*cos(x-3*dx(i))+E(5)*exp(x-4*dx(i))*cos(x-4*dx(i))+E(6)*exp(x-5*dx(i))*cos(x-5*dx(i)))/(dx(i)^2);%skewed_left_side
    fbs(i)=(exp(x)*cos(x)-2*exp(x+dx(i))*cos(x+dx(i))+exp(x+2*dx(i))*cos(x+2*dx(i)))/(dx(i)^2);%second_order_forward_Differencing_Scheme
    bds(i)=(exp(x)*cos(x)-2*exp(x-dx(i))*cos(x-dx(i))+exp(x-2*dx(i))*cos(x-2*dx(i)))/(dx(i)^2);%second_order_backward_Differencing_Scheme
    error_cds(i)=abs(cds(i)-analytical);%error_calculation
    error_skewed_right(i)=abs(skewed_right(i)-analytical);%error_calculation
    error_skewed_left(i)=abs(skewed_left(i)-analytical);%error_calculation
    error_fbs(i)=abs(fbs(i)-analytical);%error_calculation
    error_bds(i)=abs(bds(i)-analytical);%error_calculation
end
figure(1)%logarithmic plot comparing errors of CDS FDS BDS
loglog(dx,error_cds,\'r\')
hold on
loglog(dx,error_skewed_right,\'b\')
hold on
loglog(dx,error_skewed_left,\'m\')
legend(\'error cds\',\'error skewed right\',\'error skewed left\',\'location\',\'northwest\')
title(\'variation of error of CDS FDS BDS w.r.to dx\')
xlabel(\'dx\')
ylabel(\'error\')
figure(2)
loglog(dx,error_cds,\'r\')
hold on
loglog(dx,error_skewed_right,\'b\')
hold on
loglog(dx,error_skewed_left,\'m\')
hold on
loglog(dx,error_fbs,\'k\')
hold on
loglog(dx,error_bds,\'g\')
legend(\'error cds\',\'error skewed right\',\'error skewed left\',\'error fbs\',\'error bds\',\'location\',\'northwest\')
title(\'variation of error w.r.to dx\')
xlabel(\'dx\')
ylabel(\'error\')
figure(3)
plot(dx,error_cds,\'r\')
hold on
plot(dx,error_skewed_right,\'b\')
hold on
plot(dx,error_skewed_left,\'m\')
hold on
plot(dx,error_fbs,\'k\')
hold on
plot(dx,error_bds,\'g\')
legend(\'error cds\',\'error skewed right\',\'error skewed left\',\'error fbs\',\'error bds\',\'location\',\'northwest\')
title(\'variation of error w.r.to dx\')
xlabel(\'dx\')
ylabel(\'error\')

The outputs for the code are:

We can see that as the dx increases the error also increases.

The error change is much faster for fourth order approximations compared to forward differencing (fds) and backward diffencing (bds)

After x=0.3 the error for skewed right side increases cosiderably. So skewed right side method can be used from x=7.85*10^-4 to 0.23

When the conditions are flowing from one end of domain to the other end using the central differencing scheme may be difficult(using arbitary nodes). In these conditions the skewed methods are very useful. Such as boundary layer problems where the transfer of conditions is from end to other end.