Computer higher order approximations using skewed and symmetric stencils

The objective is to Derive the 4th order approximations of the second order derivative for

1. Central difference
2. Skewed right sided difference
3. Skewed left sided difference

Also write a MATLAB program calculate 2^nd order derivative of the given analytical function

F(x)= exp(x)*cos(x)

And compare it with the three numerical approximations that were derived. Provide a plot that compares the absolute error between the above mentioned schemes. Describe why a skewed scheme is useful? What can a skewed scheme do that a CD scheme cannot do?

Theory

Taylor Series: Represents a function as an infinite sum of terms that are calculated from the values of function derivatives at a single point

Where n! Is factorial of n and nth derivative of f evaluated at point a. When a= 0 the above series is called as maclaurin series.

Taylor Table: Determines the coefficient of terms which minimizes the taylor series error and the order of accuracy of the method.

Derivation

Central Differencing Scheme:

For fourth order approximation 5 points stencil is selected

_________*_________*_________*_________*_________*_________

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

Information is taken from both the forward and backward position in a space to compute the value of i. Therefore the central differencing of 4^th order approximations of 2^nd order derivative is given below

`(del^2f)/(delx^2)= (af(i-2)+bf(i-1)+cf(i)+df(i+1)+ef(i+2))/dx^2 +O(dx^4)`-----------------------------(1)`

Applying taylor series for coefficient a b c d and e

`ef(i-2)=ef(i)-e(2Deltax)/(1!) f^i(i)+e(2Deltax)^2/(2!) f^(ii)(i)-e(2Deltax)^3/(3!) f^(iii)(i)+e(2Deltax)^4/(4!) f^(iv)(i)-e(2Deltax)^5/(5!) f^(v)(i)+e(2Deltax)^6/(6!) f^(vi)(i)`

`bf(i-a)=bf(i)-b(Deltax)/(1!)f^i(i)+b(Deltax)^2/(2!)f^(ii)(i)-b(Deltax)^3/(3!)f^(iii)(i)+b(Deltax)^4/(4!)f^(iv)(i)-b(Deltax)^5/(5!)f^(v)(i)+b(Deltax)^6/(6!)f^(vi)(i)`

`cf(i)=cf(i)`

`df(i+1)=bf(i)+d(Deltax)/(1!)f^i(i)+d(Deltax)^2/(2!)f^(ii)(i)+d(Deltax)^3/(3!)f^(iii)(i)+d(Deltax)^4/(4!)f^(iv)(i)+d(Deltax)^5/(5!)f^(v)(i)+d(Deltax)^6/(6!)f^(vi)(i)`

`ef(i+2)=ef(i)+e(2Deltax)/(1!) f^i(i)+e(2Deltax)^2/(2!) f^(ii)(i)+e(2Deltax)^3/(3!) f^(iii)(i)+e(2Deltax)^4/(4!) f^(iv)(i)+e(2Deltax)^5/(5!) f^(v)(i)+e(2Deltax)^6/(6!) f^(vi)(i)`

Calculating the values of a b c d and e by putting the equation in matrix form AX=B

`[[1,1,1,1,1],[-2,-1,0,1,2],[2,1/2,0,1/2,2],[-8/6,-1/6,0,8/6,1/6],[16/24,1/24,0,16/24,1/24]]*((a),(b),(c),(d),(e))=((0),(0),(1),(0),(0))`

X = inversion (A)*B

a =e= -0.0833333

b=d= 1.333333

c= -2.5

Substituting all the above values in equation (1)

`(del^2f)/(delx^2)= ((1.333333(f(i-1)+f(i+1))-2.5f(i))-0.0833333e(f(i+2)+f(i-2)))/dx^2 +O(dx^4)`

Skewed Right Side Differencing Scheme

For fourth order approximations again 5 points stencils is considered

_________*_________*_________*_________*_________*_________

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

In the skewed right side differencing information is taken from the right i.e. forward positions in space to compute the value of i

`(del^2f)/(delx^2)= (af(i)+bf(i+1)+cf(i+2)+df(i+3)+ef(i+4))/dx^2 +O(dx^4)`---------------------------------------------(2)`

Similar to the central differencing applying taylor series to coefficient a b c d and e and further simplifying we get the taylor table for skewed right side differencing scheme

the matrix will be AX= B

`[[1,1,1,1,1],[0,1,2,3,4],[2,1/2,0,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

a=2.9167 b=-8.6667 c= 9.5000 d= -4.6667 e= 0.9167

Substituting all the above values in equation (2)

`del^2f)/(delx^2)= (2.9167f(i)-8.667f(i+1)+9.5000f(i+2)-4.6667f(i+3)+0.9167f(i+4))/dx^2 +O(dx^4)`

similarly

Skewed left Differencing Scheme

For fourth order approximations again 5 points stencils is considered

_________*_________*_________*_________*_________*_________

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

In the skewed left side differencing information is taken from the left i.e. backward positions in space to compute the value of i

`(del^2f)/(delx^2)= (af(i-4)+bf(i-3)+cf(i-2)+df(i+1)+ef(i))/dx^2 +O(dx^4)`----------------------------------------------------------------(3)

applying taylor series to coefficient a b c d and e and further simplifying we get the taylor table for skewed left side differencing scheme

the matrix will be AX= B

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

a=0.9167 b=-4.6667 c= 9.5000 d= -8.6667 e= 2.9167

Substituting all the above values in equation (2)

`del^2f)/(delx^2)= (0.9167f(i)-4.667f(i+1)+9.5000f(i+2)-8.6667f(i+3)+2.9167f(i+4))/dx^2 +O(dx^4)`

Approximations for forward differencing scheme

`(d^2f(x))/dx^2=2*(f(x+a)-f(x)-f^\'(x)dx)/dx^2`

Approximations for backrward differencing scheme

`(d^2f(x))/dx^2=2*(f(x+a)-f(x)+f^\'(x)dx)/dx^2`

The main program depends on the other function code which is given below.

Main Program

clear all
close all
clc

%varible decleration
x = pi/3;
dx = linspace(pi/10,pi/40,10);
 %loop to calculate a different scheme
 for  i = 1:length(dx)  % declaring variable i
     %calling a function of various scheme
     central_diff_value(i) = central_diff(x,dx(i)); 
     skewed_right_diff_value(i) = skewed_right_diff(x,dx(i));
     skewed_left_diff_value(i) = skewed_left_diff(x,dx(i));
     forward_diff_value(i) = forward_diff(x,dx(i));
     backward_diff_value(i)= backward_diff(x,dx(i));
 end
 %Plot
 loglog(dx,central_diff_value,\'linewidth\',1.5)
 hold on
 loglog(dx,skewed_right_diff_value,\'linewidth\',1.5)
 hold on
 loglog(dx,skewed_left_diff_value,\'linewidth\',1.5)
 hold on
 loglog(dx,forward_diff_value,\'linewidth\',1.5)
 hold on
 loglog(dx,backward_diff_value,\'linewidth\',1.5)
 
 xlabel(\'dx\');
 ylabel(\'Error\');
 legend(\'CD\',\'SRD\',\'SLD\',\'FD\',\'BD\');

Function code for Central differencing Method

%functional code for central differencing
function error_CD = central_diff(x,dx)
%given function f(x) = e^x*cos(x)
%second order derivative of function
scecond_order_derivative = -2*exp(x)*sin(x);
% difference calculation
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);
                       
% error calculation
error_CD = abs(central_differencing-scecond_order_derivative);
end                   

Function code for Skewed right differencing Method

%functional code for central differencing
function error_SR = skewed_right_diff(x,dx)
%given function f(x) = e^x*cos(x)
%second order derivative of function
scecond_order_derivative = -2*exp(x)*sin(x);
% difference calculation
skewed_right_differencing = (((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);

% error calculation
error_SR = abs(skewed_right_differencing-scecond_order_derivative);
end

Function code for Skewed left differencing Method

%functional code for central differencing
function error_SL = skewed_left_diff(x,dx)
%given function f(x) = e^x*cos(x)
%second order derivative of function
scecond_order_derivative = -2*exp(x)*sin(x);
% difference calculation
skewed_left_differencing = ((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);
                       % error calculation
error_SL = abs(skewed_left_differencing-scecond_order_derivative);
end

Function code for forward differencing Method

%functional code for central differencing
function error_FD = forward_diff(x,dx)
%given function f(x) = e^x*cos(x)
%second order derivative of function
scecond_order_derivative = -2*exp(x)*sin(x);
% difference calculation
forward_differencing = (2*((exp(x+dx)*cos(x+dx))-(exp(x)*cos(x))...
                       -((exp(x)*(cos(x)-sin(x)))*dx)))/(dx^2);
%error calculation
error_FD = abs(forward_differencing-scecond_order_derivative);
end

Function code for backward differencing Method

%functional code for central differencing
function error_BD = backward_diff(x,dx)
 %given function f(x) = e^x*cos(x)
 %second order derivative of function
scecond_order_derivative = -2*exp(x)*sin(x);
% difference calculation
backward_differencing = (2*((exp(x-dx)*cos(x-dx))-(exp(x)*cos(x))+((exp(x)...
                        *(cos(x)-sin(x)))*dx)))/(dx^2);
%error calculation
error_BD = abs(backward_differencing-scecond_order_derivative);
end

Plot of Comparision of Error of different Scheme

The above graph shows that for the value of dx the Central differencing scheme gives the very less error as compare to skewed right differencing and skewed left differencing scheme. also SRDS and SLDS are quit same.

Conclusion

• above graph shows that for the value of dx the Central differencing scheme gives the very less error as compare to skewed right differencing and skewed left differencing scheme. as CDS utilises avilable information from left and right both the sides. CDS results in giving higher order approximation and as the order of approximation increases the error decreases.
• Backward and forward scheme are used to compute the first order approximation. in this nebouring information is utilise to get the value.
• In skewed scheme in both right and left the information is collected from successive nebouring points therefore the error occurs is higher than all others scheme. skewed scheme is used to compute higher order approximations.

why a skewed scheme is useful? What can a skewed scheme do that a CD scheme cannot do?

Skewed scheme utilises one side of neighbouring points to compute point of intrest and central differencing utilises both right and left side of points to compute point of interest thus the error is higher in skewed scheme than central differencing scheme.

When there is fewer points available at the either side of point of intrest the central differencing schem can not be use due to insufficient points thus there we use skewed (right or left) scheme.

generally at or near the boundary when calculations are to be made then in such conditions lesser or no information is available at the either side of boundary there skewed scheme is used.

also from the plot it can be seen that skewed scheme has low differences in there as compaire to central differencing although the error is higher between skewed and central scheme.

At last this is been concluded by this project that although the CD gives the lesser error but skewed scheme has their own significance when ti comes to compute the things at boundary or insufficient availabity of information.