Derivation of a Fourth order Accurate schemes for a Second Order Derivative

A skewed scheme is useful in cases where the point that requires computation required a higher order approximation but has data available for data points only on one end of the point that is being computed. The skewed schemes are mostly chosen according to the reliability of the data available at the neighbouring points.If data is available only at one side of the point that is being computed then a skewed scheme would be more useful.

The following project draws a comparison between skewed schemes and central schemes.

The derivative is to be computed for the second order derivative. Writing the taylor series expansion

`(x+dx)=f(x)+f\'(x)dx+f\'\'(x)(dx)^2/(2!)+f\'\'\'(x)(dx)^3/(3!)+f\'\'\'\'(x)(dx)^4/(4!)+...`

Now to construct the Taylors Table the points are considered as shown in the figure below.

For Central Differencing:

`(d^2y)/dx^2=af(i-2)+bf(i-1)+cf(i)+df(i+1)+ef(i+2)+gf(i+3)`

constructing The Taylors table:

The equations to be solved to obtain the coefficients:

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

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

(4/2)a+(1/2)b+(1/2)d+(4/2)e+(9/2)g=1

(-8/6)a+(-1/6)b+(1/6)d+(8/6)e+(27/6)g=0

(16/24)a+(1/24)b+(1/24)d+(16/24)e+(81/24)g=0

(-32/120)a+(-1/120)b+(1/120)d+(32/120)e+(243/120)g=0

Now

af(i-2)+bf(i-1)+cf(i)+df(i+1)+ef(i+2)+gf(i+3)=(4/2)a+(1/2)b+(1/2)d+(4/2)e+(9/2)g f\'\'(x)(dx)^2+Kf\'\'\"\"(x)(dx)^6

Dividing both sides by dx^2

[af(i-2)+bf(i-1)+cf(i)+df(i+1)+ef(i+2)+gf(i+3)]/dx^2=[(4/2)a+(1/2)b+(1/2)d+(4/2)e+(9/2)g]f\'\'(x)+K f\'\'\"\"(x)(dx)^4

As the term corresponding to dx^4 is ignored. It is of fourth order accuracy as the truncation error corresponds to fourth order.

Similarly

Skewed right handed Scheme:

`(d^2y)/dx^2=a f(i)+b f(i+1)+c f(i+2)+df(i+3)+ef(i+4)+gf(i+5)`

equations to be solved are:

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

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

(1/2)b+(4/2)c+(9/2)d+(16/2)e+(25/2)g=1

(1/6)b+(8/6)c+(27/6)d+(64/6)e+(125/6)g=0

(1/24)b+(16/24)c+(81/24)d+(256/24)e+(625/24)g=0

(1/120)b+(32/120)c+(243/120)d+(1024/120)e+(3125/120)g=0

a f(i)+b f(i+1)+c f(i+2)+df(i+3)+ef(i+4)+gf(i+5)=[(1/2)b+(4/2)c+(9/2)d+(16/2)e+(25/2)g]f\'\'(x)(dx)^2+Kf\'\'\"\"(x)(dx)^6

Dividing both sides by dx^2

[a f(i)+b f(i+1)+c f(i+2)+df(i+3)+ef(i+4)+gf(i+5)]/dx^2=[(1/2)b+(4/2)c+(9/2)d+(16/2)e+(25/2)g]f\'\'(x)+K f\'\'\"\"(x)(dx)^4

Skewed Left handed Scheme:

`(d^2y)/dx^2=af(i-5)+bf(i-4)+cf(i-3)+df(i-2)+ef(i-1)+gf(i)`

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

-5a-4b-3c-2d-e=0

(25/2)a+(16/2)b+(9/2)c+(4/2)d+(1/2)e=1

(-125/6)a+(-64/6)b+(-27/6)c+(-8/6)d+(-1/6)e=0

(625/24)a+(256/24)b+(81/24)c+(16/24)d+(1/24)e=0

(-3125/120)a+(-1024/120)b+(-243/120)c+(-32/120)d+-(-1/120)e=0

af(i-5)+bf(i-4)+cf(i-3)+df(i-2)+ef(i-1)+gf(i)=[(25/2)a+(16/2)b+(9/2)c+(4/2)d+(1/2)e]f\'\'(x)(dx)^2+Kf\'\'\"\"(x)(dx)^6

Dividing by dx^2

af(i-5)+bf(i-4)+cf(i-3)+df(i-2)+ef(i-1)+gf(i)=[(25/2)a+(16/2)b+(9/2)c+(4/2)d+(1/2)e]f\'\'(x)+K f\'\'\"\"(x)(dx)^4.

A simple MatLab script is written to find the coefficients using reduced row echlon form command in matlab

clear all
close all
clc

CDS=[1 1 1 1 1 1 0;-2 -1 0 1 2 3 0; 4/2 1/2 0 1/2 4/2 9/2 1; -8/6 -1/6 0 1/6 8/6 27/6 0;16/24 1/24 0 1/24  16/24 32/24 0];
C=rref(CDS);
c=C(:,7)

SLHS=[1 1 1 1 1 1 0;-5 -4 -3 -2 -1 0 0;25/2 16/2 9/2 4/2 1/2 0 1;-125/6 -64/6 -27/6 -8/6 -1/6 0 0 ;625/24 256/24 81/24 16/24 1/24 0 0; -3125/120 -1024/120 -243/120 -32/120 -1/120 0 0];
L=rref(SLHS);
l=L(:,7)
SRHS=[1 1 1 1 1 1 0;0 1 2 3 4 5 0;0 1/2 4/2 9/2 16/2 25/2 1;0 1/6 8/6 27/6 64/6 125/6 0;0 1/24 16/24 81/24 256/24 625/24 0; 0 1/120 32/120 243/120 1024/120 3125/120 0];
R=rref(SRHS);
r=R(:,7)

The equation to be evaluated is exp(x)cos(x)

A comparion of the analytical and numerical values is plotted

clear all
close all
clc

x=pi/4;
dx= linspace(0.5,1,15);

for i=1:length(dx)
    
    central_difference(i)=central_scheme(x,dx(i));
    right_difference(i)=right_scheme_new(x,dx(i));
%     right_diff(i)=right_scheme_new(x,dx(i));
    left_difference(i)=left_scheme(x,dx(i));
    
    
end

loglog(dx,central_difference,dx,right_difference,dx,left_difference)
% loglog(dx,right_difference,dx,right_diff) %,dx,left_difference)
xlabel(\'dx\')
ylabel(\'error\')
legend(\'central difference\',\'skewed right hand\',\'skewed left hand\')

This is evaluated using a function in Matlab instead of having to type it again.

function out=orig(z)
out=exp(z)*cos(z);
end

this function can be called each time the function is to be evaluated.

The central scheme is evaluated in the following function.

function error= central_scheme(x,h)
analytical=-2*exp(x)*sin(x);
a=-0.0833;
b=1.3333;
c=-2.5;
d=1.3333;
e=-0.0833;
m=2*h;
numerical=(a*orig(x-m)+b*orig(x-h)+c*orig(x)+d*orig(x+h)+e*orig(x+m))/(h^2);
error=abs(analytical-numerical);

end

The right Handed skewed scheme:

function error= right_scheme_new(x,h)
analytical=-2*exp(x)*sin(x);
sum=0;
dx=h;
a=[3.75 -12.833 17.8333 -13.0 5.0833 -0.8333];
for i=1:length(a)
    
    sum= sum + a(i)*orig(x+(i-1)*dx);
    
end
numerical=sum/(h^2);

error=abs(analytical-numerical);
end

The Left handed skewed scheme is evaluated in the function.

function error= left_scheme(x,h)
analytical=-2*exp(x)*sin(x);
a=-0.83333;
b=5.0833;
c=-13.0;
d=17.83333;
e=-12.8333;
g=3.75;
m=2*h;
n=3*h;
o=4*h;
p=5*h;
numerical=(a*orig(x-p)+b*orig(x-o)+c*orig(x-n)+d*orig(x-m)+e*orig(x-h)+g*orig(x))/(h^2);
error=abs(analytical-numerical);

The forward difference and backward differences are calculated using the following two functions:

function out=forward_difference(x,h)
analytical=-2*exp(x)*sin(x);
numerical=(orig(x+2*h)-2*orig(x+h)+orig(x))/(h^2);
out=abs(analytical-numerical);
end

function out=backward_difference(x,h)
analytical=-2*exp(x)*sin(x);
numerical=(orig(x)-2*orig(x-h)+orig(x-2*h))/(h^2);
out=abs(analytical-numerical);
end

All the five schemes are plotted using log log plots

The plot shows that the errors reduce as the value of dx is reduced. The central difference has the least error for a given range of dx. It can also be seen that the accuracy of the fourth order schemes are better than the forward difference and backward difference which are second order accurate schemes.