Computer higher order approximations using skewed and symmetric stencils
Deriving the 4th order approximations for the second-order derivative
1. Central Differencing:
In central differencing, the values are to be taken from both the left and right-hand side of the point of interest.
Hence, The equation can be written as:
`f\'\'(x) = af(i-2) + bf(i-1) + cf(i) + d(fi+1) +ef(i+2)` __________ 1
And from the taylor table we have:
| f(x) | f\'(x)Δx | f\'\'(x)Δx^2 | f\'\'\'(x)Δx^3 | f\'\'\'\'(x)Δx^4 | |
| a*f(i-2) | a | -2a | 2a | -4a/3 | 2a/3 |
| b*f(i-1) | b | -b | b/2 | -b/6 | b/24 |
| c*f(i) | c | 0 | 0 | 0 | 0 |
| d*f(i+1) | d | d | d/2 | d/6 | d/24 |
| e*f(i+2) | e | 2e | 2e | 4e/3 | 2e/3 |
| 0 | 0 | 1 | 0 | 0 |
Adding the rows of each column we get 5 equations from this table
Calculating values of a,b,c,d,e by putting the equations in matrix form and using Matlab we get,
a = -0.0833
b = 1.3333
c = -2.500
d = 1.3333
e = -0.0833
Therefore rewriting equation 1 we get,
`f\'\'(x) = -0.0833*f(x-2dx) + 1.3333*f(x+dx) + -2.5*f(x) + 1.3333*f(x+dx) + -0.0833*f(x+2dx)`
2. Skewed right-sided difference:
In skewed right-sided difference, the values are to be taken only from the right-hand side of the point of interest.
Hence, The equation can be written as:
`f\'\'(x) = af(i) +bf(i+1) +cf(i+2) +df(i+3) + ef(i+4)+gf(i+5)` _______ 2
And from the taylor table we have:
| f(x) | f\'(x)Δx | f\'\'(x)Δx^2 | f\'\'\'(x)Δx^3 | f\'\'\'\'(x)Δx^4 | f\'\'\'\'\'(x)Δx^5 | |
| a*f(i) | a | 0 | 0 | 0 | 0 | 0 |
| b*f(i+1) | b | b | b/2 | b/6 | b/24 | b/120 |
| c*f(i+2) | c | 2c | 2c | 4c/3 | 2c/3 | 4c/15 |
| d*f(i+3) | d | 3d | 9d/2 | 9d/2 | 27d/8 | 81d/140 |
| e*f(i+4) | e | 4e | 8e | 32e/3 | 32e/3 | 128e/15 |
| g*f(i+5) | g | 5g | 25g/2 | 125g/6 | 625g/24 | 625g/24 |
| 0 | 0 | 1 | 0 | 0 | 0 |
Adding the rows of each column we get 6 equations from the table
Calculating the values of a,b,c,d,e,g by putting the equations in matrix form and using Matlab we get,
a = 3.7500
b = -12.8333
c = 17.8333
d = -13
e = 5.0833
g = -0.8333
Therefor rewriting the equation 2 we get,
`f\'\'(x) = 2.3534*f(x) + -5.8506*f(x+dx) + 3.8678*f(x+2dx) + 0.9655*f(x+3dx) + -1.8994*f(x+4dx) + 0.5632*f(x+5dx)`
3. Skewed left sided difference:
In skewed left sided difference the values are to be taken only from the left hand side.
The equation can be written as:
`f\'\'(x) = af(i) +bf(i-1) +cf(i-12) +df(i-3) + ef(i-4)+gf(i-5)` ________ 3
And from the taylor table we have:
| f(x) | f\'(x)Δx | f\'\'(x)Δx^2 | f\'\'\'(x)Δx^3 | f\'\'\'\'(x)Δx^4 | f\'\'\'\'\'(x)Δx^5 | |
| a*f(i) | a | 0 | 0 | 0 | 0 | 0 |
| b*f(i-1) | b | -b | b/2 | -b/6 | b/24 | -b/120 |
| c*f(i-2) | c | -2c | 2c | -4c/3 | 2c/3 | -4c/15 |
| d*f(i-3) | d | -3d | 9d/2 | -9d/2 | 27d/8 | -81d/40 |
| e*f(i-4) | e | -4e | 8e | -32e/3 | 32e/3 | -128e/15 |
| g*f(i-5) | g | -5g | 25g/2 | -125g/6 | 625g/24 | -625g/24 |
| 0 | 0 | 1 | 0 | 0 | 0 |
Adding the rows of each column we get 6 equations from the table
Calculating the values of a,b,c,d,e,g by putting the equations in matrix form and using Matlab we get,
a = 3.7500
b = -12.8333
c = 17.8333
d = -13
e = 5.0833
g = -0.8333
Therefore rewriting equation 3 we get,
`f\'\'(x) = 3.75*f(x) + -12.8333*f(x-dx) + 17.8333*f(x-2dx) + -13*f(x-3dx) + 5.0833*f(x-4dx) + -0.8333*f(x-5dx)`
The program:
clear all
close all
clc
%Given function f(x) = exp(x)*cos(x)
%Second derivative of the function f\'\'(x) = -2*exp(x)*sin(x)
%Defining x
x = pi/3;
%Calculating analytical derivative
analytical_derivative = -2*exp(x)*sin(x);
%Calculating the coefficients for the central difference for fourth order
%approximation
%We have 5 equations whose third equation is equal to 1 and rest are equal
%to zero
A1 = [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]; %Equations of central difference in matrix form
B1 = [0; 0; 1; 0; 0 ];
C1 = A1\\B1; %Values of coefficients of the central difference scheme
%Calculating the coefficients for the skewed right hand sided difference
%for fourth-order approximation
%We have 6 equations whose third equation is equal to 1 and rest are equal
%to zero
A2 = [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]; %Equations of skewed right sided difference in matrix form
B2 = [0; 0; 1; 0; 0 ;0];
C2 = A2\\B2; %Values of coefficients of the skewed right-hand scheme
%Calculating the coefficients for the skewed left hand sided difference
%for fourth-order approximation
%We have 6 equations whose third equation is equal to 1 and rest are equal
%to zero
A3 = [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]; %Equations of skewed left sided difference in matrix form
B3 = [0; 0; 1; 0; 0; 0];
C3 = A3\\B3; %Values of coefficients of the skewed left-hand scheme
%Calculating the numerical approximations for each scheme
%We will use a range of dx to know how the error changes for different
%values of dx
dx = linspace(pi/3,pi/30000,20);
for i = 1:length(dx)
cds(i) = ((C1(1) * ((exp(x-(2*dx(i)))) * (cos(x-(2*dx(i)))))) + (C1(2) * ((exp(x-(dx(i)))) * (cos(x-(dx(i)))))) + (C1(3) * ((exp(x)) * (cos(x)))) + (C1(4) * ((exp(x+(dx(i)))) * (cos(x+(dx(i)))))) + (C1(5) * ((exp(x+(2*dx(i)))) * (cos(x+(2*dx(i)))))))/(dx(i)^2);
srd(i) = ((C2(1) * ((exp(x)) * (cos(x)))) + (C2(2) * ((exp(x+dx(i))) * (cos(x+dx(i))))) + (C2(3) * ((exp(x+(2*dx(i)))) * (cos(x+(2*dx(i)))))) + (C2(4) * ((exp(x+(3*dx(i)))) * (cos(x+(3*dx(i)))))) + (C2(5) * ((exp(x+(4*dx(i)))) * (cos(x+(4*dx(i)))))) + (C2(6) * ((exp(x+(5*dx(i)))) * (cos(x+(5*dx(i)))))))/(dx(i)^2);
sld(i) = ((C3(1) * ((exp(x)) * (cos(x)))) + (C3(2) * ((exp(x-dx(i))) * (cos(x-dx(i))))) + (C3(3) * ((exp(x-(2*dx(i)))) * (cos(x-(2*dx(i)))))) + (C3(4) * ((exp(x-(3*dx(i)))) * (cos(x-(3*dx(i)))))) + (C3(5) * ((exp(x-(4*dx(i)))) * (cos(x-(4*dx(i)))))) + (C3(6) * ((exp(x-(5*dx(i)))) * (cos(x-(5*dx(i)))))))/(dx(i)^2);
fds(i) = (((exp(x+dx(i)))*(cos(x+dx(i)))) - (exp(x)*cos(x)))/(dx(i));
bds(i) = (((exp(x-dx(i)))*(cos(x-dx(i)))) - (exp(x)*cos(x)))/(dx(i));
%Difference between the numerical and actual derivative
cds_error(i) = abs(cds(i) - analytical_derivative);
srd_error(i) = abs(srd(i) - analytical_derivative);
sld_error(i) = abs(sld(i) - analytical_derivative);
fds_error(i) = abs(fds(i) - analytical_derivative);
bds_error(i) = abs(bds(i) - analytical_derivative);
end
loglog(dx,cds_error)
hold on
loglog(dx, srd_error)
hold on
loglog(dx, sld_error)
hold on
loglog(dx, fds_error)
hold on
loglog(dx, bds_error)
legend(\'Central difference\', \'Skewed right sided difference\', \'Skewed left sided difference\', \'Forward difference\', \'Backward difference\')
xlabel(\'Value of dx\')
ylabel(\'Error\')
title(\'Comparison of errors for different scheme\')The plot obtained:
Comparison of fourth-order approximations and FDS and BDS from the plots obtained:
Here we see that these fourth-order approximations plots have significantly lower error than the forward and backward difference scheme. Hence, the fourth-order approximations are much more accurate than the FDS and BDS
Also for higher values of dx, we can see that the error of all fourth-order approximations and first-order approximations(FDS and BDS) are getting closer. Therefore, we can say that for fourth-order approximations the error decreases rapidly compared to first-order approximations, whereas there is no comparable significant change in the first-order approximations.
General Comparison between FDS BDS and fourth-order approximations:
The FDS and BDS are first-order approximations that mean it evaluates by taking information from one neighboring point on the grid.
Whereas the fourth-order approximations take information from multiple points on the grid. From both sided in case of central differencing
This means that the accuracies of the fourth-order approximations will be much higher since more information is been used for evaluation compared to first-order FDS and BDS
Why a skewed system is useful?
In a skewed system, the point of interest is evaluated by taking information from any one side of the POI, whereas in central difference the information is taken from both sides.
Therefore, when there are fewer points available on any one side of the system from the point to be evaluated then CD cannot be used since less information is available
For eg. When calculations are to be made near the boundary, lesser or no information is available on the side of the boundary, therefore, we need to use a skewed system.
From the plots, we can also see that the skewed systems have low difference between their errors although the errors are higher than the central difference system. This point might be helpful wherever applicable.