Accuracy of fourth order finite difference schemes
Using the Taylor Table method, the fourth order approximate finite difference scheme for the second order derivative has been derived for 3 different stencils.
1. CENTRAL DIFFERENCE SCHEME
To derive the fourth order approximate central difference scheme for the second order derivative, we use the following formula:
`n = p + q - 1`
where, n is the number of nodes in the stencil
p is the order of the derivative
q is the order of approximation
For our case: p = 2, q = 4 which gives n = 5
So, we can write our hypothesis equation for the second order derivative as follows:
`/_\\x^2 (d^2f)/(dx^2) ~~ a f(i-2) + b f(i-1) + c f(i) + d f(i+1) + e f(i+2)`
The Taylor Table to find the co-efficients is shown below:
Now, using the hypothesis equation and the Taylor Table, a set of linear equations involving the co-efficients can be formed which can be represented in matrix form as:
`[[1,1,1,1,1],[-2,-1,0,1,2],[2,1/2,0,1/2,2],[-4/3,-1/6,0,1/6,4/3],[2/3,1/24,0,1/24,2/3]] [[a],[b],[c],[d],[e]] = [[0],[0],[1],[0],[0]]`
By solving the above linear system, using matrix inversion we get the values of the co-efficients as:
`a = -1/12, b = 4/3,c = -5/2,d = 4/3,e = -1/12`
Thus, the the fourth order approximate central difference scheme for the second order derivative can be expressed as:
`f\'\'(i) = (-1/12 f(i-2) + 4/3 f(i-1) -5/2 f(i) + 4/3 f(i+1) - 1/12 f(i+2))/(/_\\x^2)`
2. SKEWED RIGHT SIDED DIFFERENCE SCHEME
To derive the fourth order approximate skewed right sided difference scheme for the second order derivative, we use the following formula:
`n = p + q `
where, n is the number of nodes in the stencil
p is the order of the derivative
q is the order of approximation
For our case: p = 2, q = 4 which gives n = 6
So, we can write our hypothesis equation for the second order derivative as follows:
`/_\\x^2 (d^2f)/(dx^2) ~~ a f(i) + b f(i+1) + c f(i+2) + d f(i+3) + e f(i+4) + gf(i+5)`
The Taylor Table to find the co-efficients is shown below:
Now, using the hypothesis equation and the Taylor Table, a set of linear equations involving the co-efficients can be formed which can be represented in matrix form as:
`[[1,1,1,1,1,1],[0,1,2,3,4,5],[0,1/2,2,9/2,8,25/2],[0,1/6,4/3,9/2,32/3,125/6],[0,1/24,2/3,27/8,32/3,625/24],[0,1/120,4/15,81/40,128/15,625/24]] [[a],[b],[c],[d],[e],[g]] = [[0],[0],[1],[0],[0],[0]]`
By solving the above linear system, using matrix inversion we get the values of the co-efficients as:
`a = 15/4, b = -77/6,c = 107/6,d = -13,e = 61/12,f = -5/6`
Thus, the the fourth order approximate skewed right sided difference scheme for the second order derivative can be expressed as:
`f\'\'(i) = (15/4f(i) -77/6 f(i+1) +107/6 f(i+2) -13 f(i+3) +61/12 f(i+4)-5/6f(i+5))/(/_\\x^2)`
Now, let\'s prove the order of accuracy of this scheme
To do so, we need to expand each term by its Taylor series. On doing so we get:
∆x2 f’’(i) = (a+b+c+d+e+g) f(i) + (b+2c+3d+4e+5g) ∆x f’(i)
+ (b/2 + 4c/2 + 9d/2 + 16e/2 + 25g/2) ∆x2 f’’(i)
+ (b/6 + 8c/6 + 27d/6 + 64e/6 + 125g/6) ∆x3 f’’’(i)
+ (b/24 + 16c/24 +81d/24 + 256e/24 + 625g/24) ∆x4 fIV(i)
+ (b/120 + 32c/120 + 243d/120 + 1024e/120 + 3125g/120) ∆x5 fV(i)
+ (b/720 + 64c/720 + 729c/720 + 4096e/720 + 15625g/720) ∆x6 fVI(i) + …
We know that the linear system used to solve for the values of (a,b,c,d,e,g) was formed by setting the values of the co-efficients of f(i), ∆x f’(i), ∆x3 f’’’(i), ∆x4 fIV(i), ∆x5 fV(i) to zero and the co-efficient of ∆x2 f’’(i) to 1.
Thus, all other terms except ∆x2 f’’(i) will vanish, leaving ∆x6 fVI(i) and other higher order terms.
On dividing both sides by ∆x2 we can write:
f’’(i) ~ f’’(i) + ∆x4 fVI(i) + . . .
which implies the scheme is O(∆x4).
2. SKEWED LEFT SIDED DIFFERENCE SCHEME
To derive the fourth order approximate skewed left sided difference scheme for the second order derivative, we use the following formula:
`n = p + q `
where, n is the number of nodes in the stencil
p is the order of the derivative
q is the order of approximation
For our case: p = 2, q = 4 which gives n = 6
So, we can write our hypothesis equation for the second order derivative as follows:
`/_\\x^2 (d^2f)/(dx^2) ~~ a f(i-5) + b f(i-4) + c f(i-3) + d f(i-2) + e f(i-1) + gf(i)`
The Taylor Table to find the co-efficients is shown below:
Now, using the hypothesis equation and the Taylor Table, a set of linear equations involving the co-efficients can be formed which can be represented in matrix form as:
`[[1,1,1,1,1,1],[5,4,3,2,1,0],[25/2,8,9/2,2,1/2,0],[125/6,32/3,9/2,4/3,1/6,0],[625/24,32/3,27/8,2/3,1/24,0],[625/24,128/15,81/40,4/15,1/120,0]] [[a],[b],[c],[d],[e],[g]] = [[0],[0],[1],[0],[0],[0]]`
By solving the above linear system, using matrix inversion we get the values of the co-efficients as:
`a = -5/6, b = 61/12,c = -13,d = 107/6,e = -77/6,f = 15/4`
Thus, the the fourth order approximate skewed left sided difference scheme for the second order derivative can be expressed as:
`f\'\'(i) = (-5/6f(i-5) +61/12 f(i-4) -13 f(i-3) +107/6 f(i-2) -77/6 f(i-1)+15/4f(i))/(/_\\x^2)`
Now, let\'s prove the order of accuracy of this scheme
To do so, we need to expand each term by its Taylor series. On doing so we get:
∆x2 f’’(i) = (a+b+c+d+e+g) f(i) - (5a+4b+3c+2d+e) ∆x f’(i)
+ (25a/2 + 16b/2 + 9c/2 + 4d/2 + e/2) ∆x2 f’’(i)
- (125a/6 + 64b/6 + 27c/6 + 8d/6 + e/6) ∆x3 f’’’(i)
+ (625a/4 + 256b/24 + 81c/24 + 16d/24 + e/24) ∆x4 fIV(i)
- (3125a/120 + 1024b/120 + 243c/120 + 32d/120 + e/120) ∆x5 fV(i)
+ (15625a/720 + 4096b/720 + 729c/720 + 64d/720 + e/720) ∆x6 fVI(i) + …
We know that the linear system used to solve for the values of (a,b,c,d,e,g) was formed by setting the values of the co-efficients of f(i), ∆x f’(i), ∆x3 f’’’(i), ∆x4 fIV(i), ∆x5 fV(i) to zero and the co-efficient of ∆x2 f’’(i) to 1.
Thus, all other terms except ∆x2 f’’(i) will vanish, leaving ∆x6 fVI(i) and other higher order terms.
On dividing both sides by ∆x2 we can write:
f’’(i) ~ f’’(i) + ∆x4 fVI(i) + . . .
which implies the scheme is O(∆x4).
The MATLAB Code to solve for the co-efficients is shown below:
% Co-efficients for the hypothesis equation
clear all
close all
clc
format rational
i = input(\'1. Central Difference, 2. Right Skewed, 3. Left Skewed\');
%4th order Central difference scheme for 2nd order derivative
if i == 1
A = [ 1 1 1 1 1
-2 -1 0 1 2
2 1/2 0 1/2 2
-4/3 -1/6 0 1/6 4/3
2/3 1/24 0 1/24 2/3]
B = [0
0
1
0
0]
X = inv(A)*B
end
%4th order Skewed right sided difference scheme for 2nd order derivative
if i == 2
A = [1 1 1 1 1 1
0 1 2 3 4 5
0 1/2 2 9/2 8 25/2
0 1/6 4/3 9/2 32/3 125/6
0 1/24 2/3 27/8 32/3 625/24
0 1/120 4/15 81/40 128/15 625/24]
B = [0
0
1
0
0
0]
X = inv(A)*B
end
%4th order Skewed left sided difference scheme for 2nd order derivative
if i == 3
A = [1 1 1 1 1 1
5 4 3 2 1 0
25/2 8 9/2 2 1/2 0
125/6 32/3 9/2 4/3 1/6 0
625/24 32/3 27/8 2/3 1/24 0
625/24 128/15 81/40 4/15 1/120 0]
B = [0
0
1
0
0
0]
X = inv(A)*B
endNow,
In order to analyse the different schemes derived above a test function is chosen for which the analytical derivative can be computed. Then, the schemes are applied for the function and the error is being computed.
The MATLAB Code to compare the accuracy of these schemes is given below:
% Comparison of central difference, right and left skewed schemes
clear all
close all
clc
x = pi/3;
dx = linspace(5e-3,1,50);
% Analytical derivative at x = pi/3
analytical_derivative = derivative(x);
% Numerical derivative at x = pi/3 with varying dx
for i = 1:length(dx)
error_c(i) = central_difference(x,dx(i),analytical_derivative);
error_r(i) = skewed_right(x,dx(i),analytical_derivative);
error_l(i) = skewed_left(x,dx(i),analytical_derivative);
end
figure(1)
loglog(dx,error_c,\'-s\');
hold on
loglog(dx,error_r,\'-*\');
hold on
loglog(dx,error_l,\'-o\');
xlabel(\'dx\');
ylabel(\'Numerical Error\');
xlim([5e-3,1]);
title(\'Comparison of accuracy of 4th order schemes for 2nd order derivative\');
legend(\'Central Difference\',\'Skewed Right\',\'Skewed Left\',\'location\',\'northwest\');
% f(x) = exp(x)*cos(x)
function f = func(x)
f = exp(x)*cos(x);
end
%Analytical 2nd derivative
function f = derivative(x)
f = -2*exp(x)*sin(x);
end
% 4th order central difference scheme for the 2nd order derivative
function e = central_difference(x,dx,d)
point1 = -(1/12)*func(x-(2*dx));
point2 = (4/3)*func(x-dx);
point3 = -(5/2)*func(x);
point4 = (4/3)*func(x+dx);
point5 = -(1/12)*func(x+(2*dx));
d_c = (point1 + point2 + point3 + point4 + point5)/(dx*dx);
e = abs(d-d_c);
end
% 4th order skewed right sided difference scheme for the 2nd order derivative
function e = skewed_right(x,dx,d)
point1 = (15/4)*func(x);
point2 = -(77/6)*func(x+dx);
point3 = (107/6)*func(x+(2*dx));
point4 = -13*func(x+(3*dx));
point5 = (61/12)*func(x+(4*dx));
point6 = - (5/6)*func(x+(5*dx));
d_r = (point1 + point2 + point3 + point4 + point5 + point6)/(dx*dx);
e = abs(d_r-d);
end
% 4th order skewed left sided difference scheme for the 2nd order derivative
function e = skewed_left(x,dx,d)
point1 = 15*func(x)/4;
point2 = -(77/6)*func(x-dx);
point3 = (107/6)*func(x-(2*dx));
point4 = -13*func(x-(3*dx));
point5 = (61/12)*func(x-(4*dx));
point6 = -(5/6)*func(x-(5*dx));
d_l = (point1 + point2 + point3 + point4 + point5 + point6)/(dx*dx);
e = abs(d_l-d);
endUsing the above code, the error in the different 4th order schemes can be calculated. The results have been plotted as shown:
In the above plot, the error in different schemes has been plotted for a range of values of dx.
We can see that the central difference scheme is more accurate compared to the skewed schemes. Moreover, as the value of dx increases, the error in all the schemes increases(with some exceptions in the skewed schemes). It can be seen that the skewed right scheme is less accurate compared to the skewed left scheme.
Why do we need SKEWED SCHEMES?
Skewed schemes utilise the value of the function at either the left or the right nodes of the concerned node points, however, central schemes need both the left and right values for the computation of the derivative.
When we need to compute the derivative at the boundary node points, central schemes won\'t be valid there as we would be either having only the left node value or the right node value. So, we would need skewed schemes in this case.
Moreover, when the dominant flow direction is known, skewed schemes are preffered to get a faster and more accurate solution. On using a skewed scheme which utilizes the nodes in the dominant flow direction, the solution converges quickly.