Program on deriving fourth order approximation for different difference schemes for second order derivative

In this programming, I have derived the four order approximation of second order derivative using finite difference method.

In Finite difference method,we generally convert ordinary differential equation into difference equation. In order to achieve so, we use simple three techniques i.e Forward differencing, Backward differencing and central differencing scheme.

So as a part of challenge, I have to use central difference(symmetrical difference), skewed right difference(forward difference), skewed left difference(backward difference) in order to find out the fourth order approximation of second order derivative.

So as to derive the difference equation for all these schemes for this challenge, I have used taylor table method which is nothing else but the combination of taylor series expansion for each node/grid point where we represent the expression for each grid point/node as follow

TAYLOR SERIES EXPANSION:-

 

 

TAYLOR TABLE METHOD:-

 

 

 

Using the same technique I have derived the equation as follow below:-

1. CENTRAL DIFFERENCE SCHEME:-

    In this scheme, equal number of stencils are taken from both the right and left direction.Hence, in this scheme stencils are in symmetrical way.

    In central differnce scheme, we consider (n-1) stencils for deriving the equation.

    Number of stencils=N-1

    for,4th order approximation:- number of stencils=6-1=5,This is because  here we have to derive the fourth order approximation.

    Taylor series for stencils:-

    

 

     Taylor series table method for this scheme is:-

 

    

 

    Hence set of algeebric equations obtained for fouth order approximation are:-

 

   a+b+c+d+e=0

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

  -4a/2-b/2+0+d/2+4e/2=1

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

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

 

    using matrix multipliclation method, coeffcients  obtained are:-

    a=-1/12            b=16/12         c=-30/12          d=16/12               e=-1/12

 

    hence, difference equation obtained is:-

   

 

2.  SKEWED RIGHT SIDED DIFFERENCE

    In this scheme, the stencils are arranged in  positive direction such that each stencil is next by one.

    In this scheme, we consider 'N' of  stencils

    so, number of stencils=N=6.

    here there are 6 stencils because we are deriving the difference equation for four order approximation through skewed stencil not by symmetrical one.

    stencils derived using taylor series:-

   

 

   Taylor table is obtained below from expansion of taylor series above.

 

   

 

   therefore, set of algebraic equations are:-

   a+b+c+d+e+g=0.

   0+b+2c+3d+4e+5g=0.

   0+b/2+2c+9d/2+16e/2+25g/2=1.

   0+b/6+8c/6+27d/6+64e/6+125g/6=0.

   0+b/24+16c/24+81d/24+256e/24+625g/24=0.

   0+b/120+32c/120+243d/120+1024e/120+3125g/120=0.

 

   so, coeffcients are:-

   a=3.75    b=-77/6    c=107/6    d=-13   e=61/12    g=-5/6.

   

   hence equation is:-

   

     

 

3. SKEWED LEFT SIDED DIFFERENCE:-

    In this scheme, stencils are arranged in negative direction such that each stencil is one less than the previous one.

    similarly like forward differencing, we will consider here 'N' number of stencils.

    so, number of stencils =6.

    again reason is same.

    Taylor series expansion for this scheme is:-

    

 

    Taylor table method for this scheme is:-

    

 

   Therefore, set of algebraic equation are:-

   a+b+c+d+e+g=0.

   0+b+2c+3d+4e+5g=0.

   0+b/2+2c+9d/2+16e/2+25g/2=1.

   0+b/6+8c/6+27d/6+64e/6+125g/6=0.

   0+b/24+16c/24+81d/24+256e/24+625g/24=0.

   0+b/120+32c/120+243d/120+1024e/120+3125g/120=0.

   

   so, coeffcients are:-

   a=3.75    b=-77/6    c=107/6    d=-13   e=61/12    g=-5/6.

 

   

 

 

clear all
close all
clc

x=pi/3;

dx=linspace(pi/30,pi/300,20);

% given:- f(x)=exp(x)*cos(x);

% second order derivative is given below:-

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

for i=1:length(dx)
   
% second_order_numerical_derivative with central differencing:-

% second_order_numerical_derivative=((-30/12)*f(x)+(16/12)*f(x+dx)+(16)/12*f(x-dx)-(1/12)(f(x+2*dx)-(1/12)f(x-2*dx)))/(dx^2);
   
  numerical_derivative_CDS(i)=((-30/12)*(exp(x)*cos(x))+(16/12)*(exp(x+dx(i))*cos(x+dx(i))+exp(x-dx(i))*cos(x-dx(i)))-(1/12)*(exp(x-2*dx(i))*cos(x-2*dx(i))+exp(x+2*dx(i))*cos(x+2*dx(i))))/((dx(i)^2));

% second_order_numerical_derivative with right_skewed differencing scheme:-
% second_order_numerical_derivative=(3.75*f(x)-(77/6)*f(x+dx)+(107/6)*f(x+2*dx)-13*f(x+3*dx)+(61/12)*f(x+4*dx)-(5/6)*f(x+5*dx)/(dx^2);
 
 numerical_derivative_skew_right(i)=((3.75)*(exp(x)*cos(x))-(77/6)*(exp(x+dx(i))*cos(x+dx(i)))+(107/6)*(exp(x+2*dx(i))*cos(x+2*dx(i)))-13*(exp(x+3*dx(i))*cos(x+3*dx(i)))+(61/12)*(exp(x+4*dx(i))*cos(x+4*dx(i)))-(5/6)*exp(x+5*dx(i))*cos(x+5*dx(i)))/((dx(i)^2));
 
 
% second_order_numerical_derivative with left_sided differencing scheme:-
% second_order_numerical_derivative=(3.75*f(x)-(77/6)*f(x-dx)+(107/6)*f(x-2*dx)-13*f(x-3*dx)+(61/12)*f(x-4*dx)-(5/6)*f(x-5*dx)/(dx^2);
 
 numerical_derivative_skew_left(i)=((3.75)*(exp(x)*cos(x))-(77/6)*(exp(x-dx(i))*cos(x-dx(i)))+(107/6)*(exp(x-2*dx(i))*cos(x-2*dx(i)))-13*(exp(x-3*dx(i))*cos(x-3*dx(i)))+(61/12)*(exp(x-4*dx(i))*cos(x-4*dx(i)))-(5/6)*exp(x-5*dx(i))*cos(x-5*dx(i)))/((dx(i)^2));
  
% second_order_numerical_derivative with forward differencing scheme:-
% f''(x)=(f(x+2*dx)+f(x+dx)-2*f(x))/((dx)^2);
   
  numerical_derivative_FDS(i)=(exp(x+dx(i))*cos(x+dx(i))+exp(x+2*dx(i))*cos(x+2*dx(i))-2*exp(x)*cos(x))/(dx(i)^2);

% second_order_numerical_derivative with backward differencing scheme:-
% f''(x)=(f(x-2*dx)+f(x-dx)-2*f(x))/((dx)^2);
   
  numerical_derivative_BDS(i)=(exp(x-dx(i))*cos(x-dx(i))+exp(x-2*dx(i))*cos(x-2*dx(i))-2*exp(x)*cos(x))/(dx(i)^2);
 
 end 

% error produced due to difference between analytical_derivative and second_order_numerical_derivative with CDS :-

 error_CDS=abs(analytical_derivative-numerical_derivative_CDS);
 
% error produced due to difference between analytical_derivative and second_order_numerical_derivative with skew_right :-
 
 error_skew_right=abs(analytical_derivative-numerical_derivative_skew_right);
 
% error produced due to difference between analytical_derivative and second_order_numerical_derivative with skew_left :- 
 
  error_skew_left=abs(analytical_derivative-numerical_derivative_skew_left);
% error produced due to difference between analytical_derivative and second_order_numerical_derivative with FDS :-  
 
  error_FDS=abs(analytical_derivative-numerical_derivative_FDS);

% error produced due to difference between analytical_derivative and second_order_numerical_derivative with BDS :-  
 
  error_BDS=abs(analytical_derivative-numerical_derivative_BDS);
 
figure(1)
loglog(dx,error_CDS,"color",'r',"linewidth",3)

hold on
loglog(dx,error_skew_right,"color",'b',"linewidth",3)

loglog(dx,error_skew_left,"color",'g',"linewidth",3)
title('fourth order approximation with different difference scheme')
ylabel('fourth order error')
xlabel('different value of dx') 
legend('central difference','skewed right sided difference','skewed left sided difference',4)

figure(2)
loglog(dx,error_FDS,"color","black","linewidth",3)
hold on
loglog(dx,error_BDS,"color","magenta","linewidth",3)
hold on
loglog(dx,error_CDS,"color",'r',"linewidth",3)

hold on
loglog(dx,error_skew_right,"color",'b',"linewidth",3)

loglog(dx,error_skew_left,"color",'g',"linewidth",3)
title('fourth order approximation with different difference scheme')
ylabel('fourth order error')
xlabel('different value of dx') 
legend('forward difference','backward difference','central difference','skewed right sided difference','skewed left sided difference',1)

 

 

(a) In the above plot in the form of log, you can easily see that central scheme is having the least value of error as compared to other two schemes i.e BDS and CDS. After central scheme,backward differencing scheme is having the least value followed by forward differencing scheme.so it give the order of scheme in following way,

                central scheme(symmetrical scheme)

  central scheme has less error because it has both sided stencil into consideration compare to other schemes. While in skewed sided, there is only one sided stencil considerations.Hence creates less reduction in error.In this plot, the central scheme has error ranging around 10^-9 to 10^-5 as we increase the small value of 'dx' from 10^-2 to 10^-1 over the plot.where as,  both skewed scheme has error around 10^-6 to 10^-2 as we increase the value of dx.

 

 

(b)

   

  The forward and backward scheme has totally different curve apart from other three schemes.And It has error around 10^3 to 10^1 although it still tries to decrease which simply implies that they are unable to fit in this second order    approximation due to its such a large error. Hence it proves that central and skewed sided schemes has more accuracy in approximation.Therefore these three schemes fits in such types of approximation.

 

 (c) Skewed sided difference i.e left/right difference is sometimes more useful than when compare to          symmetrical sided difference. This is because in skewed differencing, we genrally have 'N' of              stencils which increases more grid point.Hence provides more information about the nearest grid          points.while in symmetrical stencils, there is no such case as there is always one less stencil than        thee skewed stencil due to its symmetric nature.Hence, there is less information about nearest            grid point.

      Also, there is one more advantage in skewed stencils is that it takes corner nodes/stencils into              consideration for more information which is generally not in the case of symmetrical stencils.                Hence creates backlog for symmetric sided difference.


Projects by Shyam Rai

THEORY:- BOUNDARY CONDITION :- IN PHYSICS ,THESE ARE THOSE CONDITIONS WHICH ARE APPLIED AT THE BOUNDARY OF THE DOMAIN . IT SPECIFIES THE CONDITION AT BOUNDARY ON THE BASIS OF WHAT THE GOVERNING EQUATION DEFINES THE CONDITION OF DOMAIN ALONG WITH THE HELP OF INTIAL COND Read more

OBJECTIVE:- STEADY STATE SIMULATION OF FLOW OVER A THROTTLE BODY . SOFTWARE USED :- 1. CONVERGE STUDIO:- TO SETUP THE MODEL.  2. CYGWIN:- FOR RUNNING THE SIMULATION. 3. OPENFOAM:- TO VISUALIZE THE DIFFERENT POST-PROCESSED RESULT.   THEORY:- HERE, WE HAV Read more

OBJECTIVE:- STEADY STATE SIMULATION OF FLOW OVER A THROTTLE BODY . SOFTWARE USED :- 1. CONVERGE STUDIO:- TO SETUP THE MODEL.  2. CYGWIN:- FOR RUNNING THE SIMULATION. 3. OPENFOAM:- TO VISUALIZE THE DIFFERENT POST-PROCESSED RESULT.   THEORY:- HERE, WE HAV Read more

AIM:- IN THIS PROJECT SIMULATION IS DONE ON THE FLUID FLOW THROUGH PIPE USING OPENFOAM SOFTWARE AND MATLAB . THE TOPIC GIVEN BELOW IS COVERED IN THIS                      PROJECT. TO MAKE A PROGRAM IN MATLAB THAT Read more

AIM:- IN THIS PROJECT SIMULATION IS DONE ON THE FLUID FLOW THROUGH PIPE USING OPENFOAM SOFTWARE AND MATLAB . THE TOPIC GIVEN BELOW IS COVERED IN THIS PROJECT. TO MAKE A PROGRAM IN MATLAB THAT CAN GENERATE MESH AUTOMATICALLY FOR ANY WEDGE ANGLE AND GRADING SCHEME. TO Read more

AIM :- IN THIS PROJECT, I HAVE WORKED ON SIMULATING THE FLOW OF FLUID THROUGH  THE DOMAIN (CONTROL VOLUME). DESCRIPTION:- THE DOMAIN IS NOTHING ELSE, IT IS JUST A CONTROL VOLUME THROUGH WHICH I HAVE SIMULATED THE RESULT OBTAINED FROM THE GOVERNING EQUATION. THE Read more

  IN THIS PROJECT, I HAVE WORKED ON THE ONE DIMENSIONAL SUPERSONIC NOZZLE FLOW. IN THIS FLOW, I HAVE ASSUMED THE FLOW TO BE STEADY,ISENTROPIC. IN THIS NOZZLE, I HAVE CONSIDERED THE FLOW AT INLET SECTION COMES FROM RESERVOIR WHERE THE PRESSURE, TEMPERATURE ARE DENO Read more

  In this project, I have worked on the simulation of conduction of temperature heat in two dimensional space (2D domain).  The aim was to determine the rate of flow of  heat due to temperature difference over the space domain under steady condition and Read more

In this report, I have simulated the oscillation of 2D pendulum and also generate a plot of \'angular_displacement\' & \'angular_velocity\' vs \"time\" in octave. Consider a pendulum which is having a string of length of \'1metre\' connected to a ball of mass,m=1kg Read more


Loading...

The End