Computer higher order approximations using skewed and symmetric stencils

Objective -

1. To Derive the 4th order approximations of the second order derivative by using

i. Central difference

ii. Skewed right sided difference

iii. Skewed left sided difference

2. To write a program in Matlab to evaluate the second order derivative of the analytical function exp(x)*cos(x) and compare it with the 3 numerical approximations i) Central difference, ii) Skewed right sided difference, and iii) Skewed left sided difference derived results.

3. To plot a graph that compares the absolute error between the i) Central difference, ii) Skewed right sided difference, and iii) Skewed left sided difference schemes

4. To describe why a skewed scheme is useful and What can a skewed scheme do that a CD scheme cannot do.

Solution:-

1. Fourth order approximation of a second order derivative using central differencing scheme

*---------------*---------------*----------------*---------------*

(i-2) (i-1) (i) (i+1) (i+2)

-Numerical stencil for central differencing scheme

`(del^2u)/(del^2x) ≅ (af(i-2)+bf(i-1)+cf(i)+df(i+1)+ef(i+2))/dx^2 + O(dx^4))`------- Equation(1)

`af(i-2) = af(i) - (af\'(i).(2dx))/(1!) + (af\'\'(i).(2dx)^2)/(2!) +(af\'\'\'(i).(2dx)^3)/(3!) + (af\'\'\'\'(i).(2dx)^4)/(4!)-----`

Expanding the above terms till 4th order.

Similarly, expanding all the terms present in equation(1).

Now these terms are represented in the form of table to solve the coefficients of equation(1) and this can be solved using \'Taylor Table method\'.

Stencil / Taylor Series terms	f(i)	dx.f’(i)	dx2.f’’(i)	dx3.f’’’(i)	dx4.f’’’’(i)
af(i-2)	a	-2a	4a/2	-8a/6	16a/24
bf(i-1)	b	-b	b/2	-b/6	b/24
cf(i)	c	0	0	0	0
df(i+1)	d	d	d/2	d/6	d/24
ef(i+2)	e	2e	4e/2	8e/6	16e/24
Sum of the Coefficients	0	0	1	0	0

From the above table when we sum the derivatives terms we get the equations -

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

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

`2a+b/2+0+d/2+2e = 1`

`(-8a)/16-b/6+0+d/6+(8e)/6 = 0`

`(16a)/24+b/24+0+d/24+(16e)/24 =0`

The above equations have to be solved using matrix inversion method

AX=B

`[[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),(b),(c),(d),(e))=((0),(0),(1),(0),(0))`

X = inversion(A)*B

From the above values we get a=0.0833, b=1.3333, c=-2.5000, d=1.3333, e=-0.0833

substituting the above values in equation(1) for approximations

`(del^2u)/(del^2x) ≅ ((-0.0833)(f(i-2))+(1.3333)(f(i-1))+(-2.5000)(f(i))+(-1.3333)(f(i+1))+(-0.0833)(f(i+2)))/dx^2 + O(dx^4))`

2. Fourth order approximation of a second order derivative using skewed right sided differencing scheme

*---------------*---------------*----------------*---------------*

(i) (i+1) (i+2) (i+3) (i+4)

-Numerical stencil for skewed right sided differencing scheme

`(del^2u)/(del^2x) ≅ (af(i)+bf(i+1)+cf(i+2)+df(i+3)+ef(i+4))/dx^2 + O(dx^4))` ----------- Equation(2)

expanding each of the term in equation(2) and putting them in a Taylor Table, we get

Stencil / Taylor Series terms	f(i)	dx.f’(i)	dx2.f’’(i)	dx3.f’’’(i)	dx4.f’’’’(i)
af(i)	a	0	0	0	0
bf(i+1)	b	b	b/2	b/6	b/24
cf(i+2)	c	2c	4c/2	8c/6	16c/24
df(i+3)	d	3d	9d/2	27d/6	81d/24
ef(i+4)	e	4e	16e/2	64e/6	256e/24
Sum of the Coefficients	0	0	1	0	0

From the above table when we sum the derivatives terms we get the equations -

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

`0+b+2c+3d+4e = 0`

`0+b/2+2c+9/2d+8e = 1`

`0+b/6+8/6c+27/6d+64/6e = 0`

`0+b/24+16/24c+81/24d+256/24e = 0`

The above equations have to be solved using matrix inversion method

AX=B

`[[1,1,1,1,1],[0,1,2,3,4],[0,1/2,2,9/2,8],[0,1/6,8/6,27/6,64/6],[0,1/24,16/24,81/24,256/24]].((a),(b),(c),(d),(e)) = ((0),(0),(1),(0),(0))`

X=inversion(A)*B

From the above values we get a=2.9167, b=1.3333, c=-8.6667, d=9.5000, e=0.9167

substituting the above values in equation(2) for approximations

`(del^2u)/(del^2x) ≅ ((2.9167)(f(i))+(-8.6667)(f(i+1))+(9.5000)(f(i+2))+(0.9167)(f(i+3))+(-0.0833)(f(i+4)))/dx^2 + O(dx^4))`

3. Fourth order approximation of a second order derivative using skewed left sided differencing scheme

*---------------*---------------*----------------*---------------*

(i-4) (i-3) (i-2) (i-1) (i)

-Numerical stencil for skewed left sided differencing scheme

`(del^2u)/(del^2x) ≅ (af(i-4)+bf(i-3)+cf(i-2)+df(i-1)+ef(i))/dx^2 + O(dx^4))` ------- Equation(3)

expanding each of the term in equation(2) and putting them in a Taylor Table, we get

Stencil / Taylor Series terms	f(i)	dx.f’(i)	dx2.f’’(i)	dx3.f’’’(i)	dx4.f’’’’(i)
af(i-4)	a	-4a	8a/2	-64a/6	-256e/24
bf(i-3)	b	-3b	9b/2	-27b/6	81b/24
cf(i-2)	c	-2c	2c	-8c/6	16c/24
df(i-1)	d	-d	d/2	-d/6	d/24
ef(i)	e	0	0	0	0
Sum of the Coefficients	0	0	1	0	0

From the above table when we sum the derivatives terms we get the equations -

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

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

`8a+9/2b+2c+d/2+0 = 1`

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

`256/24a+81/24b+16/24c+d/24+0 = 0`

The above equations have to be solved using matrix inversion method

AX=B

`[[1,1,1,1,1],[-4,-3,-2,-1,0],[8,9/2,2,1/2,0],[-64/6,-27/6,-8/6,-1/6,0],[256/24,81/24,16/24,1/24,0]].((a),(b),(c),(d),(e)) = ((0),(0),(1),(0),(0))`

X = inversion(A)*B

From the above values we get a=0.9167, b=-4.6667, c=9.5000, d=-8.6667, e=2.9167

substituting the above values in equation(3) for approximations

`(del^2u)/(del^2x) ≅ ((0.9167)(f(i-4))+(-4.6667)(f(i-3))+(9.5000)(f(i-2))+(-8.6667)(f(i-1))+(2.9167)(f(i)))/dx^2 + O(dx^4))`

Programme

1. Function Code for central differencing scheme

% defining the function for central differencing method
function error_central_diff = central_difference(x,dx)
%f(x) = e^x*cosx is the given function
% second order derivative for thr given function
second_analytical_derivative = -2*exp(x)*sin(x);

% calculating the central differencing value
central_differencing = (((-0.0833)*((exp(x-(2*dx)))*(cos(x-(2*dx)))))+((1.3333)*((exp(x-dx))*(cos(x-dx))))-((2.5)*((exp(x))*(cos(x))))+((1.3333)*((exp(x+dx))*(cos(x+dx))))+((-0.0833)*((exp(x+(2*dx)))*(cos(x+(2*dx))))))/(dx^2);

% calculating the absolute error
error_central_diff = abs(central_differencing-second_analytical_derivative);
end

2. Function Code for skewed right sided differencing scheme

% defining the function for skewed right differencing method
function error_skewed_right = skewed_right_difference(x,dx)
%f(x) = e^x*cosx is the given function
% second order derivative for thr given function
second_analytical_derivative = -2*exp(x)*sin(x);

% calculating the skewed right sided differencing value
skewed_right_diff = (((2.9167)*((exp(x))*(cos(x)))))-((8.6667)*((exp(x+dx))*(cos(x+dx))))+((9.5000)*((exp(x+(2*dx))*(cos(x+(2*dx))))))-(4.6667)*((exp(x+(3*dx)))*(cos(x+(3*dx))))+((0.9167)*((exp(x+(4*dx)))*(cos(x+(4*dx)))))/(dx^2);

% calculating the absolute error
error_skewed_right = abs(skewed_right_diff-second_analytical_derivative);
end

3. Function Code for skewed left sided differencing scheme

% defining the function for skewed left differencing method
function error_skewed_left = skewed_left_difference(x,dx)
%f(x) = e^x*cosx is the given function
% second order derivative for thr given function
second_analytical_derivative = -2*exp(x)*sin(x);

% calculating the skewed left sided differencing value
skewed_left_diff = ((0.9167)*((exp(x-(4*dx))*(cos(x-(4*dx))))))-((4.6667)*((exp(x-(3*dx)))*(cos(x-(3*dx)))))+((9.5000)*((exp(x-(2*dx))*(cos(x-(2*dx))))))-((8.6667)*((exp(x-dx))*(cos(x-dx))))+((2.9167)*((exp(x))*(cos(x))))/(dx^2);



% calculating the absolute error
error_skewed_left = abs(skewed_left_diff-second_analytical_derivative);
end

Main Code

clear all;
close all;
clc

% Declaring the value of \'x\'
x = pi/3;

% Declaring the variable \'dx\'
dx = linspace(pi/4,pi/4000,20);

% To begin a for loop to calculate the values of different methods

% declaring the variable \'i\'
for i = 1:length(dx)
 central_diff_value(i) = central_difference(x,dx(i));
 skewed_right_value(i) = skewed_right_difference(x,dx(i));
 skewed_left_value(i) = skewed_left_difference(x,dx(i));
 end
 
 % plotting the graph using loglog plot
 figure(1)
 loglog(dx,central_diff_value,\'b\',\'linewidth\',3)
 hold on
 loglog(dx,skewed_right_value,\'g\',\'linewidth\',3)
 hold on
 loglog(dx,skewed_left_value,\'r\',\'linewidth\',3)
 xlabel(\'dx\',\'FontSize\',20)
 ylabel(\'error\',\'FontSize\',20)
 legend(\'central differencing error\',\'skewed right error\',\'skewed left error\')

Plot on loglog graph

Conclusion

From the above plot it can be observed that for the value of \'dx\' Central Differencing Scheme(CDS) gives the least error values when compared to skewed differencing scheme because CDS utilises the information available at neighbouring points on both the sides from the point of consideration to provide the approximation value. CDS results in providing higher order approximations and as the order of approximation increases the value of error decreases.

According to above plotted results initially Skewed Right Differencing Scheme(SRDS) resulting in highest value of error but as \'dx\' value decreases it stands next to Skewed Left Differencing Scheme(SLDS). In Skewed Right Differencing Scheme and Skewed Right Differencing Scheme information is collected from only one side of neighbouring points (either right side neighbouring points in case of SRDS or left side neighbouring points in case of SLDS) from the point of consideration therefore error occurs is higher than other schemes. Skewed scheme is used to compute higher order approximations.

why a skewed scheme is useful?

In the system if the sufficient information is not available on both the sides of the considered point i.e the numerical stencil wont be symmetric, In such cases we cannot use central differencing scheme so Skewed(left/right) scheme is useful .

What can a skewed scheme do that a CD scheme cannot do?

Near the boundaries where there is no access to further points over the boundary, therefore only skewed differencing scheme is used and we cannot use the Central differencing scheme.