Deriving 4th order approximations of 2nd order derivatives using taylor table for central skewed right skewed left differences and MATLAB code for the same

clear all
close all
clc

x=pi/4;
dx=linspace(pi/40,pi/60);

% f(x)=exp(x)*cos(x);
% exact_derivative_second_order =  -2*exp(x)*sin(x);

exact_derivative = -2*exp(x)*sin(x);

% calculating the coefficients of the equation in forward differencing
% scheme
% obtaining 6 equations for 6 unknowns from the taylor table, we use matrix
% form

eq_1_FDS = [1,1,1,1,1,1];
eq_2_FDS = [0,1,2,3,4,5];
eq_3_FDS = (1/factorial(2)) * [0,1,2^2,3^2,4^2,5^2];
eq_4_FDS = (1/factorial(3)) * [0,1,2^3,3^3,4^3,5^3];
eq_5_FDS = (1/factorial(4)) * [0,1,2^4,3^4,4^4,5^4];
eq_6_FDS = (1/factorial(5)) * [0,1,2^5,3^5,4^5,5^5];

% only equation 3 shall be equated to 1 as the coefficient of f\'\'(x) equals
% unity. Rest all shall be zero.

A_FDS = [eq_1_FDS;eq_2_FDS;eq_3_FDS;eq_4_FDS;eq_5_FDS;eq_6_FDS];
B_FDS = [0;0;1;0;0;0];
X_FDS = A_FDS\\B_FDS;

% Coefficents for Backward differencing scheme

eq_1_BDS = [1,1,1,1,1,1];
eq_2_BDS = -1 * [0,1,2,3,4,5];
eq_3_BDS = (1/factorial(2)) * [0,1,2^2,3^2,4^2,5^2];
eq_4_BDS = -1 * (1/factorial(3)) * [0,1,2^3,3^3,4^3,5^3];
eq_5_BDS = (1/factorial(4)) * [0,1,2^4,3^4,4^4,5^4];
eq_6_BDS = -1 * (1/factorial(5)) * [0,1,2^5,3^5,4^5,5^5];

A_BDS = [eq_1_BDS;eq_2_BDS;eq_3_BDS;eq_4_BDS;eq_5_BDS;eq_6_BDS];
B_BDS = [0;0;1;0;0;0];
X_BDS = A_BDS\\B_BDS;

% Coefficents for Central differencing scheme

eq_1_CDS = [1,1,1,1,1];
eq_2_CDS = [-2,-1,0,1,2];
eq_3_CDS = (1/factorial(2)) * [2^2,1,0,1,2^2];
eq_4_CDS = (1/factorial(3)) * [-2^3,-1,0,1,2^3];
eq_5_CDS = (1/factorial(4)) * [2^4,1,0,1,2^4];

A_CDS = [eq_1_CDS;eq_2_CDS;eq_3_CDS;eq_4_CDS;eq_5_CDS];
B_CDS = [0;0;1;0;0];
X_CDS = A_CDS\\B_CDS;

% Calculating 2nd order derivative with 4th order approximation by FDS
% (a*f(x)+b*f(x+dx)+c*f(x+2*dx)+d*f(x+3*dx)+e*f(x+4*dx)+h*f(x+5*dx))/dx^2

% Calculating 2nd order derivative with 4th order approximation by BDS
% (a*f(x)+b*f(x-dx)+c*f(x-2*dx)+d*f(x-3*dx)+e*f(x-4*dx)+h*f(x-5*dx))/dx^2

% Calculating 2nd order derivative with 4th order approximation by CDS
% (a*f(x-dx)+b*f(x-2*dx)+c*f(x)+d*f(x+dx)+e*f(x+2*dx))/dx^2


second_order_FDS=zeros(1,length(dx));
second_order_BDS=zeros(1,length(dx));
second_order_CDS=zeros(1,length(dx));


for i=1:length(dx)
    second_order_FDS(i)=(1/(dx(i)^2))*((X_FDS(1)*f(x))+(X_FDS(2)*f(x+dx(i)))+(X_FDS(3)*f(x+2*dx(i)))+(X_FDS(4)*f(x+3*dx(i)))+(X_FDS(5)*f(x+4*dx(i)))+(X_FDS(6)*f(x+5*dx(i))));
    second_order_BDS(i)=(1/(dx(i)^2))*((X_BDS(1)*f(x))+(X_BDS(2)*f(x-dx(i)))+(X_BDS(3)*f(x-2*dx(i)))+(X_BDS(4)*f(x-3*dx(i)))+(X_BDS(5)*f(x-4*dx(i)))+(X_BDS(6)*f(x-5*dx(i))));
    second_order_CDS(i)=(1/(dx(i)^2))*((X_CDS(1)*f(x-2*dx(i)))+(X_CDS(2)*f(x-dx(i)))+(X_CDS(3)*f(x))+(X_CDS(4)*f(x+dx(i)))+(X_CDS(5)*f(x+2*dx(i))));
end

% Error

error_FDS = abs(second_order_FDS - exact_derivative);
error_BDS = abs(second_order_BDS - exact_derivative);
error_CDS = abs(second_order_CDS - exact_derivative);

% Plotting the errors of the three schemes against the step size

figure(1)
loglog(dx,error_FDS)
hold on
loglog(dx,error_BDS)
hold on
loglog(dx,error_CDS,\'color\',\'g\')
hold off
grid on
xlabel(\'dx/step size\');
ylabel(\'Magnitude of Error\');
legend(\'FDS\',\'BDS\',\'CDS\',\'location\',\'northwest\');

% Defining the function

function y = f(x)
y=exp(x)*cos(x);
 end 

Step size vs error

1. Deriving by a Central difference scheme, the 4^th order approximation for a 2^nd order derivative.

Central Difference Method

`f^\'\' (x_i ) ~ (af_(i-2)+bf_(i-1)+cf_i+df_(i+1)+ef_(i+2))/(∆x^2 )`

Where,

a, b, c, d and e are the coefficients obtained after solving the equation by Taylors expansion of each term.

Expanding each term about the point ‘i’ and forming a table,

Taylors table

	f(i)	f I (i)*(∆x)1	f ii (i)*(∆x)2	f iii (i)*(∆x)3	f iv (i)*(∆x)4
	a	-2*a	a*22/2!	-a*23/3!	a*24/4!
	b	-b	b/2!	-b/3!	b/4!
	c	0	0	0	0
	d	d	d/2!	d/3!	d/4!
	e	2*e	e*22/2!	e*23/3!	e*24/4!

Each column of the table forms an equation. We have taken 2 points from each side of the central node ‘i’. We have 5 unknowns and 5 equations; thus, the unknowns are solved using a matrix form.

Equation 1; `a+b+c+d+e=0`

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

Equation 3; `2a+0.5b+0.5d+2e=1`

Equation 4; `-1.33a-0.17b+0.17d+1.33e=0`

Equation 5; `0.67a+0.0417b+0.0417d+0.67e=0`

Note: Equation 3 is equated to 1 as we need the coefficient corresponding to the second derivative of ‘f’ equation be equal to 1 and the remaining coefficients be 0.

The equation in the matrix form is,

AX=B

2. Deriving by a Forward difference scheme, the 4^th order approximation for a 2^nd order derivative (Skewed right difference)

In this derivation, we will be considering 6 nodes in all starting from the i^th node.

Reason for considering 6 nodes in the stencil,

When we take first order approximation in second order derivative, we take 3 nodes.

When we take second order approximation in second order derivative, we take 4 nodes.

® requires 3 nodes for O(h) accuracy

® requires 4 nodes for O(h²) accuracy

® requires 5 nodes for O(h³) accuracy

® requires 6 nodes for O(h⁴) accuracy

The parent equation for the forward difference scheme is as below.

`f^\'\' (x_0 ) ~ (af_x+bf_(x+dx)+cf_(x+2dx)+df_(x+3dx)+ef_(x+4dx)+hf_(x+5dx))/(∆x^2 )`

`f_(x+∆x)= f(x)+f^\' (x)∆x+f^\'\' (x) 〖∆x〗^2⁄2!+f^\'\'\' (x) 〖∆x〗^3⁄3!+ f^iv 〖∆x〗^4⁄4!+f^v 〖∆x〗^5⁄5!`

`f_(x+2∆x)= f(x)+f^\' (x)(2.∆x)+f^\'\' (x) ((2.〖∆x〗^2))⁄2!+f^\'\'\' (x) ((2.〖∆x〗^3))⁄3!+ f^iv (〖(2∆x〗^4))⁄4!+f^v ((2.〖∆x〗^5))⁄5!`

`f_(x+3∆x)= f(x)+f^\' (x)(3.∆x)+f^\'\' (x) ((3〖∆x〗^2))⁄2!+f^\'\'\' (x) ((3〖∆x〗^3))⁄3!+ f^iv ((〖3∆x〗^4))⁄4!+f^v ((3〖∆x〗^5))⁄5!`

Where,

a, b, c, d and e are the coefficients to be determined.

Expanding each term about the point ‘i’ and forming a table,

Taylors table

f(i)	f I (i)*(∆x)1	f ii (i)*(∆x)2	f iii (i)*(∆x)3	f iv (i)*(∆x)4	f v (i)*(∆x)5	f vi (i)*(∆x)6
a	0	0	0	0	0	0
b	b	b/2!	b/3!	b/4!	b/5!	b/6!
c	2*c	c*22/2!	c*23/3!	c*24/4!	c*25/5!	c*26/6!
d	3*d	d*32/2!	d*33/3!	d*34/4!	d*35/5!	d*36/6!
e	4*e	e*42/2!	e*43/3!	e*44/4!	e*45/5!	e*46/6!
h	5*h	h*52/2!	h*53/3!	h*54/4!	h*55/5!	h*56/6!

Forming equations and writing in the matrix form to solve for the unknowns, we get,

`Equation 1; a+b+c+d+e+h=0`

`Equation 2; b+2c+3d+4e+5h=0`

`Equation 3; 0.5b+2c+4.5d+8e+12.5h=1`

`Equation 4; 0.17b+1.33c+4.5d+10.67e+20.83h=0`

`Equation 5; 0.0417b+0.67c+3.38d+10.67e+26.04h=0`

`Equation 6; 0.01b+0.27c+2.025d+8.533e+26.04h=0`

3. Deriving by a Backward difference scheme, the 4^th order approximation for a 2^nd order derivative (Skewed left difference)

The parent equation for the forward difference scheme is as below.

`f^\'\' (x) ~ (af_x+bf_(x-dx)+cf_(x-2dx)+df_(x-3dx)+ef_(x-4dx)+hf_(x-5dx))/(∆x^2 )`

`f_(x-∆x)= f(x)-f^\' (x)∆x+f^\'\' (x) 〖∆x〗^2⁄2!-f^\'\'\' (x) 〖∆x〗^3⁄3!+ f^iv 〖∆x〗^4⁄4!-f^v 〖∆x〗^5⁄5!`

`f_(x-2∆x)= f(x)-f^\' (x)(2.∆x)+f^\'\' (x) ((2.〖∆x〗^2 ))⁄2!-f^\'\'\' (x) ((2.〖∆x〗^3 ))⁄3!+ f^iv (〖(2∆x〗^4))⁄4!-f^v ((2.〖∆x〗^5))⁄5!`

`f_(x-3∆x)= f(x)-f^\' (x)(3.∆x)+f^\'\' (x) ((3〖∆x〗^2 ))⁄2!-f^\'\'\' (x) ((3〖∆x〗^3 ))⁄3!+ f^iv ((〖3∆x〗^4 ))⁄4!-f^v ((3〖∆x〗^5))⁄5!`

And so on…

Where,

a, b, c, d, e and h are the coefficients to be determined.

Expanding each term about the point ‘i’ and forming a table,

Taylors table

f(i)	f I (i)*(∆x)1	f ii (i)*(∆x)2	f iii (i)*(∆x)3	f iv (i)*(∆x)4	f v (i)*(∆x)5
a	0	0	0	0	0
b	-b	b/2!	-b/3!	b/4!	-b/5!
c	-2*c	c*22/2!	-c*23/3!	c*24/4!	-c*25/5!
d	-3*d	d*32/2!	-d*33/3!	d*34/4!	-d*35/5!
e	-4*e	e*42/2!	-e*43/3!	e*44/4!	-e*45/5!
h	-5*h	h*52/2!	-h*53/3!	h*54/4!	-h*55/5!

Forming equations and writing in the matrix form to solve for the unknowns, we get,

Solving the matrix form in MATLAB, using the equation X=A\\B, we get the corresponding values of X. These are the coefficients of the derivatives of the equation.

Our parent equation takes the form,

`f^\'\' (x_0 ) ~ (3.75f_x-12.833f_(x-dx)+17.833f_(x-2dx)-13f_(x-3dx)+5.0833f_(x-4dx)-0.833f_(x-5dx))/(∆x^2 )+O(h^4)`

1. Why a skewed scheme is useful? What can a skewed scheme do that a CD scheme cannot do?

In the FDS method aka skewed right sided difference methods, we take information from all or majority of the nodes to the right of the node we need to analyse.

Similarly, in BDS method aka skewed left sided difference methods, as the name suggests, we acquire the function information from all or majority of the nodes to the left of the current node. Whereas in the central differencing scheme, we take information from equal number of nodes left and right to the current node. However, this method cannot be applied to the nodal points falling on the edges of a meshed component. Therein we use the FDS and BDS methods accordingly.

The errors in forward difference and backward difference are roughly proportional to O(h) for a first order derivative.

For the same order of the derivative, a CDS system has a higher degree of approximation.

The fact that the central difference estimate has a dependency of order of h² rather than O(h), makes it a better estimate than forward and backward differencing.

2. How do these errors compare with FDS and BDS approximations?

From the plot, and the fact that CDS has a dependency of the order of square of the step size, we understand that the magnitude of absolute error, for a range of step sizes, between the analytic and the exact derivative will always be lower for a central system as compared to the Forward and Backward difference scheme.

3. Detailed explanation of results, plots and codes

Code

The code is self-explanatory as each and every step in the script is explained as a commented note.

The script was started with x and dx definitions.

x is taken as pi/4 as the function has cos(x) which is in radians.

dx is taken in the range pi/40 to pi/60, as the value cannot be too high or too low.

Too low a value would not affect the error and the whole reason of understanding the variation of error between an analytical and exact solution would be lost.

On the other hand, too large a value of dx would result in greater magnitudes of error and this is not acceptable during discretization.

In short, extreme values of dx will not be able to relate with the particular value of x.

The exact derivative is defined in the script file.

In order to calculate the coefficients for a forward differencing scheme, following method was employed.

1. From the Taylor’s table, 6 equations are formed corresponding to the 6 unknowns.

2. These 6 equations are then assigned to a Matrix A.

3. The results of the equations are assigned to the Matrix B

4. The unknown coefficients are assigned to the Matrix X and is calculated in MATLAB using the relation X=A\\B which is derived from AX=B.

5. The complete process is repeated for the Backward and Central differencing schemes.

The numerical value corresponding to each value of the step size is stored in an array form in all the three schemes using the equation we derived above.

An array of absolute error is now created by difference between the numerical value array and the exact derivative.

Errors corresponding to the three schemes are now plotted against the step sizes.

Finally, the function is defined, which will eventually be used in the equation to calculate the second order derivatives in all the three schemes.

Temporary information

For forward / rearward differencing

· Order of derivative + order of approximation = number of nodal info required

For central differencing

· Order of derivative + order of approximation -1 = number of nodal info required

A very brief explaination of the script file,

The step size is \'set up\' using the linspace command.

exact 2nd order derivative is defined.

An array of size 1 by x is formed where x is the number of unknowns (different for central and skewed schemes). The constant coefficients are assigned in each array. x number of equations are formed. These arrays are then stacked in an x by x matrix following the equation AX=B

X=A\\B. This is the command to calculate A inverse times B.

B matrix is the answer we get after equating all.

Similar procedure is followed for the rest of the differences.

A null array is formed again to reduce computing time. This was explained in the previous challenges.

In the iteration, each value is defined as explained and derived in the code. The function handler is called for inside the iterations which is defined at the end of the script file as is required by MATLAB.

A figure window is created plotting the absolute error vs the step size for Forward , backward and central differencing schemes.

From the plot we understand that the error magnitude is lower in central differencing as compared to forward and backward differencing.