Computer higher order approximations using skewed and symmetric stencils

Problem Statement: Derivation of Central difference, Skewed right-sided difference and Skewed left-sided difference for 4th order approximations of the second-order derivative.

Also, 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.

a. Providing a plot that compares the absolute error between the above-mentioned schemes

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

Theory: The main point is to understand the derivation of skewed and central difference scheme. For skewed difference, the points selected are the sum of order of derivative (p) and the order of error approximation i.e. 6 ( 2(p) +4(q)) but for the central difference scheme, it is one less that it, i.e. 5(2(p)+4(q)-1. Therefore, 6 terms for skewed and 5 terms for central difference scheme is to be taken.

1) Skewed Right-Sided Difference:

Let us assume;

`∂2f/∂x2= a*f(i) + b*f(i + 1) + c*f(i + 2) + d* f(i + 3) + e* f(i + 4) + f* f(i + 5)`

Where;

f(i + 1) =f(i) + ∆x*f’(x) + ∆x²/2! *f’’(x) + ∆x³/3! *f’’’(x) + ∆x⁴/4! *f’’’’(x) + ………

f(i + 2) =f(i) + 2*∆x*f’(x) + 4*∆x²/2! *f’’(x) +8* ∆x³/3! *f’’’(x) +16* ∆x⁴/4! *f’’’’(x) + ………

f(i + 3) =f(i) + 3*∆x*f’(x) +9* ∆x²/2! *f’’(x) + 27*∆x³/3! *f’’’(x) + 81*∆x⁴/4! *f’’’’(x) + ………

f(i + 4) =f(i) + 4*∆x*f’(x) +16* ∆x²/2! *f’’(x) + 64*∆x³/3! *f’’’(x) + 256*∆x⁴/4! *f’’’’(x) + ………

f(i + 5) =f(i) + 5*∆x*f’(x) +25* ∆x²/2! *f’’(x) + 125*∆x³/3! *f’’’(x) + 625*∆x⁴/4! *f’’’’(x) + ………

Therefore;

`∂2f/∂x2`= a*f(i) + b*( f(i) + ∆x*f’(x) + ∆x²/2! *f’’(x) + ∆x³/3! *f’’’(x) + ∆x⁴/4! *f’’’’(x) + ..) + c*( f(i) + 2*∆x*f’(x) + 4*∆x²/2! *f’’(x) +8* ∆x³/3! *f’’’(x) +16* ∆x⁴/4! *f’’’’(x) + …) + d*( f(i) + 3*∆x*f’(x) +9* ∆x²/2! *f’’(x) + 27*∆x³/3! *f’’’(x) + 81*∆x⁴/4! *f’’’’(x) + ..) + e*( f(i) + 4*∆x*f’(x) +16* ∆x²/2! *f’’(x) + 64*∆x³/3! *f’’’(x) + 256*∆x⁴/4! *f’’’’(x) + …) + f*( f(i) + 5*∆x*f’(x) +25* ∆x²/2! *f’’(x) + 125*∆x³/3! *f’’’(x) + 625*∆x⁴/4! *f’’’’(x) + ..)

`∂2f/∂x2`= (a+b+c+d+e+f)*f(i) + (b+2*c+3*d+4*e+5*f) *∆x*f’(x) + (b+4*c+9*d+16*e+25*f)* ∆x²/2! *f’’(x) + (b+8*c+27*d+64*e+125*f)* ∆x³/3! *f’’’(x) + (b+16*c+81*d+256*e+625*f) *∆x⁴/4! *f’’’’(x)………..

We also know that;

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

(b+2*c+3*d+4*e+5*f) *∆x =0

(b+4*c+9*d+16*e+25*f)* ∆x²/2! = 1

(b+8*c+27*d+64*e+125*f)* ∆x³/3! = 0

(b+16*c+81*d+256*e+625*f) *∆x⁴/4! = 0

Solving this we get;

a = 3.75/∆x²

b = -12.83/∆x²

c = 17.83/∆x²

d =-13/∆x²

e= 5.08/∆x²

f = -0.83/∆x²

2) Skewed Left-Sided Difference

Let us assume;

`∂2f/∂x2`= a*f(i) + b*f(i - 1) + c*f(i - 2) + d* f(i - 3) + e* f(i - 4) + f* f(i - 5)

Where;

f(i -1) =f(i) - ∆x*f’(x) + ∆x²/2! *f’’(x) - ∆x³/3! *f’’’(x) + ∆x⁴/4! *f’’’’(x) - ………

f(i - 2) =f(i) - 2*∆x*f’(x) + 4*∆x²/2! *f’’(x) -8* ∆x³/3! *f’’’(x) +16* ∆x⁴/4! *f’’’’(x) -………

f(i - 3) =f(i) - 3*∆x*f’(x) +9* ∆x²/2! *f’’(x) - 27*∆x³/3! *f’’’(x) + 81*∆x⁴/4! *f’’’’(x) - ………

f(i - 4) =f(i) - 4*∆x*f’(x) +16* ∆x²/2! *f’’(x) - 64*∆x³/3! *f’’’(x) + 256*∆x⁴/4! *f’’’’(x) - ………

f(i - 5) =f(i) - 5*∆x*f’(x) +25* ∆x²/2! *f’’(x) - 125*∆x³/3! *f’’’(x) + 625*∆x⁴/4! *f’’’’(x) - ………

and so

∂2f/∂x2 = a*f(i) + b*( f(i) - ∆x*f’(x) + ∆x²/2! *f’’(x) - ∆x³/3! *f’’’(x) + ∆x⁴/4! *f’’’’(x) - ..) + c*( f(i) - 2*∆x*f’(x) + 4*∆x²/2! *f’’(x) -8* ∆x³/3! *f’’’(x) +16* ∆x⁴/4! *f’’’’(x) - …) + d*( f(i) - 3*∆x*f’(x) +9* ∆x²/2! *f’’(x) - 27*∆x³/3! *f’’’(x) + 81*∆x⁴/4! *f’’’’(x) -..) + e*( f(i) - 4*∆x*f’(x) +16* ∆x²/2! *f’’(x) - 64*∆x³/3! *f’’’(x) + 256*∆x⁴/4! *f’’’’(x) - …) + f*( f(i) - 5*∆x*f’(x) +25* ∆x²/2! *f’’(x) - 125*∆x³/3! *f’’’(x) + 625*∆x⁴/4! *f’’’’(x) - ..)

∂2f/∂x2 = (a+b+c+d+e+f)*f(i) - (b+2*c+3*d+4*e+5*f) *∆x*f’(x) + (b+4*c+9*d+16*e+25*f)* ∆x²/2! *f’’(x) -(b+8*c+27*d+64*e+125*f)* ∆x³/3! *f’’’(x) - (b+16*c+81*d+256*e+625*f) *∆x⁴/4! *f’’’’(x)………..

We also know that;

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

(b+2*c+3*d+4*e+5*f) *∆x =0

(b+4*c+9*d+16*e+25*f)* ∆x²/2! = 1

(b+8*c+27*d+64*e+125*f)* ∆x³/3! = 0

(b+16*c+81*d+256*e+625*f) *∆x⁴/4! = 0

Solving this we get;

a = 3.75/∆x²

b = -12.83/∆x²

c = 17.83/∆x²

d =-13/∆x²

e= 5.08/∆x²

f = -0.83/∆x²

^{3) Central Difference :}

Let us assume;

∂2f/∂x2 = a*f(i) + b*f(i - 2) + c*f(i – 1) + d* f(i +1) + e* f(i +2)

Where;

f(i - 2) =f(i) - 2*∆x*f’(x) + 4*∆x²/2! *f’’(x) -8* ∆x³/3! *f’’’(x) +16* ∆x⁴/4! *f’’’’(x) -………

f(i -1) =f(i) - ∆x*f’(x) + ∆x²/2! *f’’(x) - ∆x³/3! *f’’’(x) + ∆x⁴/4! *f’’’’(x) - ………

f(i + 1) =f(i) + ∆x*f’(x) + ∆x²/2! *f’’(x) + ∆x³/3! *f’’’(x) + ∆x⁴/4! *f’’’’(x) + ………

f(i + 2) =f(i) + 2*∆x*f’(x) + 4*∆x²/2! *f’’(x) +8* ∆x³/3! *f’’’(x) +16* ∆x⁴/4! *f’’’’(x) + ………

and so

∂2f/∂x2 = a*f(i) + b*( f(i) - 2*∆x*f’(x) + 4*∆x²/2! *f’’(x) -8* ∆x³/3! *f’’’(x) +16* ∆x⁴/4! *f’’’’(x) -.) + c*( f(i) - ∆x*f’(x) + ∆x²/2! *f’’(x) - ∆x³/3! *f’’’(x) + ∆x⁴/4! *f’’’’(x) - …) + d*( f(i) + ∆x*f’(x) + ∆x²/2! *f’’(x) + ∆x³/3! *f’’’(x) + ∆x⁴/4! *f’’’’(x) +..) + e*( f(i) + 2*∆x*f’(x) + 4*∆x²/2! *f’’(x) +8* ∆x³/3! *f’’’(x) +16* ∆x⁴/4! *f’’’’(x) + …)

∂2f/∂x2 = (a+b+c+d+e)*f(i) - (-2*b-c+d+2*e) *∆x*f’(x) + (4*b+c+d+4*e)*∆x²/2! *f’’(x) -(-8*b-c+d+8*e)* ∆x³/3! *f’’’(x) - (16*b+c+d+16*e) *∆x⁴/4! *f’’’’(x)………..

We also know that;

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

(-2*b-c+d+2*e) *∆x = 0;

(4*b+c+d+4*e)*∆x²/2! = 2;

(-8*b-c+d+8*e)* ∆x³/3! = 0;

(16*b+c+d+16*e) *∆x⁴/4! = 0

Solving this we get;

a = -2.5/∆x²

b = -0.083/∆x²

c = 1.33/∆x²

d = 1.33/∆x²

e = -0.083/∆x²

^Program:

^{1. The first step will be to define the above equation for Skewed Right-Sided Difference and Skewed Left-Sided Difference and calculate the value of constants a,b,c,d,e and f.}

^{2. The first step will be to define the above equation for Central Difference and calculate the value of constants a,b,c,d, and e.}

3. The next step will be to define the function of FDS. The value of constants will be passed as an argument and used in the equation and then calculate the truncation error.

3. Also, the function of BDS will be defined. The value of constants will be passed as an argument and used in the equation and then calculate the truncation error.

4. And the function of BDS will be defined. The value of constants will be passed as an argument and used in the equation and then calculate the truncation error.

6. Also, define a function which will give the value of exp(x)*cos(x).

5. The next step would be to define the values of x and dx and define a loop that will give values of error of FDS, BDS, and CDS.

6. The last step will be to plot all there error values of FDS, BDS and CDS for different values of dx with appropriate labels and legends.

Graph:

The usefulness of Skewed Scheme:

From the graph, it can be noticed that CDS is more accurate than FDS/ BDS, that may be because it used the same no. of iteration on both sides os central points. The central difference scheme is second order in h">h, whereas the other two are first order in h">h. Thus at fixed sufficiently small h">h, the central difference will have a better error than the other two.

But Central Difference lacks in providing Boundedness in several systems (when peclet number > 2, flow is not unidirectional) whereas Skewed System is unconditionally bounded and highly stable in these systems.

Program:

Main:

clear all
close all
clc

%To find a,b,c,d,e for FDS and BDS series.
syms a b c d e f
eqn1=a+b+c+d+e+f == 0;
eqn2=0+b+2*c+3*d+4*e+5*f == 0;
eqn3=0+b+4*c+9*d+16*e+25*f == 2;
eqn4=0+b+8*c+27*d+64*e+125*f == 0;
eqn5=0+b+16*c+81*d+256*e+625*f == 0;
eqn6=0+b+32*c+243*d+256*4*e+625*5*f ==0;    

var=[a,b,c,d,e,f];
[A,B] = equationsToMatrix([eqn1, eqn2, eqn3, eqn4, eqn5, eqn6],var);

Sol= single(linsolve(A,B));


%To fine a,b,c,d,e,f for CDS series
syms a1 b1 c1 d1 e1 f1

eqn1=a1+b1+c1+d1+e1 == 0;
eqn2=0-2*b1-c1+d1+2*e1 == 0;
eqn3=0+4*b1+c1+d1+4*e1== 2;
eqn4=0-8*b1-c1+d1+8*e1 == 0;
eqn5=0+16*b1+c1+d1+16*e1 == 0;

var=[a1,b1,c1,d1,e1];
[A1,B1] = equationsToMatrix([eqn1, eqn2, eqn3, eqn4, eqn5],var);


Sol1= single(linsolve(A1,B1));
x= pi/3;
dx= linspace(pi/50,pi/500,15);
for i = 1:length(dx)
    fds(i)= FDS(x,dx(i),Sol);
    bds(i)= BDS(x,dx(i),Sol);
    cds(i)= CDS(x,dx(i),Sol1);
   
    
end

%Plotting the Graph
figure(1)
loglog(dx,fds)
grid on
hold on
loglog(dx,bds,'r')
hold on

loglog(dx,cds,'g')
hold on
xlabel('dx')
ylabel('Error')
legend('FDS','BDS','CDS')

FDS:

function out = FDS(x, dx, S)

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

fds = (S(1)*F(x)+S(2)*F(x+dx)+S(3)*F(x+2*dx)+S(4)*F(x+3*dx)+S(5)*F(x+4*dx)+S(6)*F(x+5*dx))/(dx^2);
out = abs(Analytical-fds);
end

BDS:

function out = BDS(x, dx, S)

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

bds = (S(1)*F(x)+S(2)*F(x-dx)+S(3)*F(x-2*dx)+S(4)*F(x-3*dx)+S(5)*F(x-4*dx)+S(6)*F(x-5*dx))/(dx^2);
out = abs(Analytical-bds);
end

CDS:

function out = CDS(x,dx,S)

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

cds = (S(1)*F(x)+S(2)*F(x-2*dx)+S(3)*F(x-dx)+S(4)*F(x+dx)+S(5)*F(x+2*dx))/(dx^2);

out = abs(Analytical-cds);

end

F:

function out = F(x)

out=exp(x)*cos(x)
end