CFD || Derivation and Analysis of Higher Order Finite Difference with various Difference Schemes

PROBLEM STATEMENT

1. Derive the following 4th order approximations of the second order derivative.
2. Write a program in Matlab to evaluate the second order derivative of the analytical function exp(x)*cos(x) and compare it with various difference schemes (Central Difference, Skewed Right Sided Difference, Skewed Left Sided Difference).
3. Compares the absolute error between the above mentioned schemes

OVERVIEW

• Analytical Function = exp(x)*cos(x)
• Analytical Derivative = exp(x)*sin(x)

Since we need to calculate the fourth order approximation we need atleast 4 points for consideration.

• Skewed Right Side Difference

Let `(del^2f)/(delx^2)=af_i + bf_(i+1)+cf_(i+2)+df_(i+3)+ef_(i+4)`......(1)

According to Taylor's Table;

Since we are finding `(del^2f)/(delx^2)` the co-efficient of f''(i) = 1 and others will be zero;

`[a+b+c+d+e]=0`....(i)

`[b+2c+3d+4e]Deltax=0`....(ii)

`[b+4c+9d+16e](Deltax^2)/2=1`....(iii)

`[b+8c+27d+64e](Deltax^3)/6=0`....(iv)

`[b+16c+81d+256e](Deltax^4)/24=0`....(v)

Solving the above equation we get;

`a = 35/12; b=-104/12;c=114/12;d=-56/12;e=11/12`

Substituting the above values in equation ......(1) we get;

`bb((del^2f)/(delx^2)= (35f_i-104f_(i+1)+114f_(i+2)-56f_(i+3)+11f_(i+4))/(12Deltax^2))`

• Skewed Left Side Difference

Let `(del^2f)/(delx^2)=af_i + bf_(i-1)+cf_(i-2)+df_(i-3)+ef_(i-4)`......(2)

According to Taylor's Table;

Since we are finding `(del^2f)/(delx^2)` the co-efficient of f''(i) = 1 and others will be zero;

`[a+b+c+d+e]=0`....(i)

`[-b-2c-3d-4e]Deltax=0`....(ii)

`[b+4c+9d+16e](Deltax^2)/2=1`....(iii)

`[-b-8c-27d-64e](Deltax^3)/6=0`....(iv)

`[b+16c+81d+256e](Deltax^4)/24=0`....(v)

Solving the above equation we get;

`a = 35/12; b=-104/12;c=114/12;d=-56/12;e=11/12`

Substituting the above values in equation ......(2) we get;

`bb((del^2f)/(delx^2)= (35f_i-104f_(i-1)+114f_(i-2)-56f_(i-3)+11f_(i-4))/(12Deltax^2))`

• Central Difference

Let `(del^2f)/(delx^2)=af_i + bf_(i+1)+cf_(i+2)+df_(i-1)+ef_(i-2)`......(3)

According to Taylor's Table;

Since we are finding `(del^2f)/(delx^2)` the co-efficient of f''(i) = 1 and others will be zero;

`[a+b+c+d+e]=0`....(i)

`[b+2c-d-2e]Deltax=0`....(ii)

`[b+4c+d+4e](Deltax^2)/2=1`....(iii)

`[b+8c-d-8e](Deltax^3)/6=0`....(iv)

`[b+16c+d+16e](Deltax^4)/24=0`....(v)

Solving the above equation we get;

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

Substituting the above values in equation ......(3) we get;

`bb((del^2f)/(delx^2)= (-30f_i+16f_(i+1)-f_(i+2)+16f_(i-1)-f_(i-2))/(12Deltax^2))`

PROGRAM

%################  MAIN FUNCTION ####################

clear all
close all
clc

%%_____Analytical function = exp(x)*cosx


x = pi/3;
dx=linspace(1e-8,1,1000);


%%_______Finding Error


for i=1:length(dx)
    
SRD_error(i) = srsdfn(x,dx(i));
SLD_error(i) = slsdfn(x,dx(i));
CD_error(i) = cdfn(x,dx(i));

end


%%______Plotting

 
loglog(dx,SRD_error,'r')
hold on
loglog(dx,SLD_error,'k')
hold on
loglog(dx,CD_error,'b')
legend('Skewed RHD','Skewed LHD','Central Difference')
xlabel('log(dx)')
ylabel('log(error)')






%########################### SKEWED RIGHT SIDE SCHEME  ###############

function error=srsdfn(x,dx)

analytical_second_derivative=(-2*exp(x)*sin(x));
srsd=((35*exp(x)*cos(x))-(104*exp(x+dx)*cos(x+dx))+(114*exp(x+2*dx)*cos(x+2*dx))-(56*exp(x+3*dx)*cos(x+3*dx))+(11*exp(x+4*dx)*cos(x+4*dx)))/(12*(dx^2));
 
error=abs(analytical_second_derivative-srsd);
end


%############################  SKEWED LEFT SIDE SCHEME  ################

function error=slsdfn(x,dx)

analytical_second_derivative=(-2*exp(x)*sin(x));
slsd=((35*exp(x)*cos(x))-(104*exp(x-dx)*cos(x-dx))+(114*exp(x-2*dx)*cos(x-2*dx))-(56*exp(x-3*dx)*cos(x-3*dx))+(11*exp(x-4*dx)*cos(x-4*dx)))/(12*(dx^2));

 
error=abs(analytical_second_derivative-slsd);
end

%############################  CENTRAL DIFFERENCE SCHEME  ##################

function error=cdfn(x,dx)

analytical_second_derivative=(-2*exp(x)*sin(x));
cd=(-(exp(x+2*dx)*cos(x+2*dx))+(16*exp(x+dx)*cos(x+dx))-(30*exp(x)*cos(x))+(16*exp(x-dx)*cos(x-dx))-(exp(x-2*dx)*cos(x-2*dx)))/(12*(dx^2));

error=abs(analytical_second_derivative-cd);
end

RESULTS & CONCLUSION

1. Its clear from the above plot that the central difference scheme is more accurate compared to the skewed schemes;because it gathers information from both sides equally.
2. The cental difference is more stable compared to the skewed difference scheme.
3. Since difference schemes provide only approximate solution there are two errors involved in it;truncation error and round-off error. As grid size reduces the truncation error decreases and round-off error increases proportionally, but the total error decreases. After a certain point, reduction in grid size results in drastic increase of round-off error but the truncation error decreases propotional, hence total error increases.
4. When there are no symmetrc stencil available the central difference cannot be applied hence skewed scheme are used to provide solution.(eg;Skewed schemes are used in boundaries because points are not available on either sides of the corner node)