Let `f(x)` be the function under question. Let `dx` be the grid size of the discretization. The first derivative of the function would be denoted as `f\'(x)`, the second derivative as `f\'\'(x)` and so on.
- 4th order approximation of the 2nd derivative in Central Difference Scheme:
For a central difference scheme, the numerical stencil would include values of the function at grid points from either side of the grid point where the derivative needs to be calculated. The numerical stencil looks as follows:
`obrace(underset(i-2)(ast)--underset(i-1)(ast))^(dx)--underset(i)(ast)--underset(i+1)(ast)--underset(i+2)(ast)`
Hence,
`f\'\' (x)=af(x-2dx)+bf(x-dx)+cf(x)+df(x+dx)+ef(x+2dx)`
Writing the taylor series expansion for each of the terms, we get,
`f\'\' (x)=a[f(x)-f\' (x)*2dx+f\'\' (x)*(2dx)^2/(2!)- f\'\'\' (x)*(2dx)^3/(3!)+ f\'\'\'\' (x)*(2dx)^4/(4!)]+b[f(x)-f\' (x)*dx+f\'\' (x)*(dx^2)/(2!)- f\'\'\' (x)*(dx^3)/(3!)+ f\'\'\'\' (x)*(dx^4)/(4!)]+c[f(x)]+d[f(x)+f\' (x)*dx+f\'\' (x)*(dx^2)/(2!)+ f\'\'\' (x)*(dx^3)/(3!)+ f\'\'\'\' (x)*(dx^4)/(4!)]+e[f(x)+f\' (x)*2dx+f\'\' (x)*(2dx)^2/(2!)+ f\'\'\' (x)*(2dx)^3/(3!)+ f\'\'\'\' (x)*(2dx)^4/(4!)]`
Since we are looking for the 4th order approximation, we neglect any terms in the expansion of an order higher than that and thus expand each term only until the 4th order term.
Rearranging this gives us,
`f\'\' (x)=(a+b+c+d+e)*f(x)+(-2a-b+d+2e)*f\' (x)*dx+(2a+b/2+d/2+2e)*f\'\' (x)*dx^2+(-4a/3-b/6+d/6+4e/3)*f\'\'\' (x)*dx^3+(2a/3+b/24+d/24+2e/3)*f\'\'\'\'(x)*dx^4`
Using this, we equate coefficients on LHS with the ones on RHS. This gives us the following 5 equations:
`a+b+c+d+e=0`
`-2a-b+d+2e=0`
`2a+(b/2)+(d/2)+2e=dx^2`
`-(4a/3)-(b/6)+(d/6)+(4e/3)=0`
`(2a/3)+(b/24)+(d/24)+(2e/3)=0`
These equations were solved using Wolfram Equation solver. [They can also be solved using a simple matrix inversion method in Matlab] Solving these equations, we get,
`a =(-1)/(12*dx^2); b = 4/(3*dx^2); c (-5)/(2*dx^2); d = 4/(3*dx^2); e =(-1)/(12*dx^2)`
Hence, using a central difference scheme,
`f\'\' (x)=(1/dx^2)*((-1)/12 f(x-2dx)+4/3 f(x-dx)-5/2 f(x)+4/3 f(x+dx)-1/12 f(x+2dx))`
- Skewed right sided difference:
We use the same process as before, except that the numerical stencil is now right-side skewed. The numerical stencil looks as follows:
`obrace(underset(i)(ast)--underset(i+1)(ast))^(dx)--underset(i+2)(ast)--underset(i+3)(ast)--underset(i+4)(ast)`
Hence,
`f\'\' (x)=af(x)+bf(x+dx)+cf(x+2dx)+df(x+3dx)+ef(x+4dx)`
This gives,
`f\'\' (x)=a[f(x)]+b[f(x)+f\' (x)*dx+f\'\' (x)*(dx^2)/(2!)+ f\'\'\' (x)*(dx^3)/(3!)+ f\'\'\'\' (x)*(dx^4)/(4!)]+c[f(x)+f\' (x)*2dx+f\'\' (x)*(2dx)^2/(2!)+ f\'\'\' (x)*(2dx)^3/(3!)+ f\'\'\'\' (x)*(2dx)^4/(4!)]+d[f(x)+f\' (x)*3dx+f\'\' (x)*(3dx)^2/(2!)+ f\'\'\' (x)*(3dx)^3/(3!)+ f\'\'\'\' (x)*(3dx)^4/(4!)]+e[f(x)+f\' (x)*4dx+f\'\' (x)*(4dx)^2/(2!)+ f\'\'\' (x)*(4dx)^3/(3!)+ f\'\'\'\' (x)*(4dx)^4/(4!)]`
`=>f\'\' (x)=(a+b+c+d+e)*f(x)+(b+2c+3d+4e)*f\' (x)*dx+(b/2+2c+9d/2+8e)*f\'\' (x)*dx^2+(b/6+4c/3+9d/2+32e/3)*f\'\'\' (x)*dx^3+(b/24+2c/3+27d/8+32e/3)*f\'\'\'\'(x)*dx^4`
This gives us the following 5 equations:
`a+b+c+d+e=0`
`b+2c+3d+4e=0`
`b/2+2c+9d/2+8e=dx^2`
`b/6+4c/3+9d/2+32e/3=0`
`b/24+2c/3+27d/8+32e/3=0`
Solving these gives,
`a=35/(12*dx^2);b=-26/(3*dx^2);c=19/(2*dx^2);d=-14/(3*dx^2);e=11/(12*dx^2)`
Hence, using skewed right-side difference,
`f\'\' (x)=(1/dx^2)*((35/12)f(x)-(26/3)f(x+dx)+(19/2)f(x+2dx)-(14/3)f(x+3dx)+(11/12)f(x+4dx))`
- Skewed left side Difference:
Using a left-side skewed difference, the numerical stencil looks as follows:
`obrace(underset(i-4)(ast)--underset(i-3)(ast))^(dx)--underset(i-2)(ast)--underset(i-1)(ast)--underset(i)(ast)`
Hence,
`f\'\' (x)=af(x)+bf(x-dx)+cf(x-2dx)+df(x-3dx)+ef(x-4dx)`
This gives,
`f\'\' (x)=a[f(x)]+b[f(x)-f\' (x)*dx+f\'\' (x)*(dx^2)/(2!)- f\'\'\' (x)*(dx^3)/(3!)+ f\'\'\'\' (x)*(dx^4)/(4!)]+c[f(x)-f\' (x)*2dx+f\'\' (x)*(2dx)^2/(2!)- f\'\'\' (x)*(2dx)^3/(3!)+ f\'\'\'\' (x)*(2dx)^4/(4!)]+d[f(x)-f\' (x)*3dx+f\'\' (x)*(3dx)^2/(2!)- f\'\'\' (x)*(3dx)^3/(3!)+ f\'\'\'\' (x)*(3dx)^4/(4!)]+e[f(x)-f\' (x)*4dx+f\'\' (x)*(4dx)^2/(2!)- f\'\'\' (x)*(4dx)^3/(3!)+ f\'\'\'\' (x)*(4dx)^4/(4!)]`
`=>f\'\' (x)=(a+b+c+d+e)*f(x)-(b+2c+3d+4e)*f\' (x)*dx+(b/2+2c+9d/2+8e)*f\'\' (x)*dx^2-(b/6+4c/3+9d/2+32e/3)*f\'\'\' (x)*dx^3+(b/24+2c/3+27d/8+32e/3)*f\'\'\'\'(x)*dx^4`
This gives us the following 5 equations:
`a+b+c+d+e=0`
`b+2c+3d+4e=0`
`b/2+2c+9d/2+8e=dx^2`
`b/6+4c/3+9d/2+32e/3=0`
`b/24+2c/3+27d/8+32e/3=0`
Solving these gives,
`a=35/(12*dx^2);b=-26/(3*dx^2);c=19/(2*dx^2);d=-14/(3*dx^2);e=11/(12*dx^2)`
Hence, using skewed right-side difference,
`f\'\' (x)=(1/dx^2)*((35/12)f(x)-(26/3)f(x-dx)+(19/2)f(x-2dx)-(14/3)f(x-3dx)+(11/12)f(x-4dx))`
- Plot Comparing Error of Various Schemes:
The plot shows the variation of error with respect to changing grid-size on a logarithmic scale and compares central difference method to skewed right sided and skewed left sided difference methods. With increase in dx, all methods give greater error. This means that with increase in grid-size, accuracy of a computational solution falls, no matter the scheme used. It can be seen that the central difference scheme has a less error or is more accurate than the skewed schemes. Between the skewed schemes, a skewed left-sided difference scheme seems to be slightly more accurate than a skewed right-sided scheme at greater grid-sizes.
- Advantages of a Skewed scheme over Central Difference:
As can be seen from the plot above, a central difference scheme is better than any skewed scheme of the same order with respect to the error. However, this comes with its own set of restrictions. A central difference can only be used if the numerical stencil is available on either side of the location where the differential is to be calculated. This is usually not available on the discretized nodes/elements close to the boundary of the computational mesh. Hence, at locations where numerical stencil is unavailable on one side, due to geometric restrictions of the problem, a skewed scheme needs to be used instead of a central difference scheme.
We assume a value for x, in this case, `pi/3`.
clc;
clear all;
close all;
x = pi/3;
dx = linspace(pi/1000, pi/50, 500);
% Analytical Derivative
analytical_derivative = -2*exp(x)*sin(x);
Analytical Derivative for the function `f(x) = (e^x)*cos(x)` would be `f\'(x) = (e^x)*(cos(x)-sin(x))` and its second derivative would be:
`f\'\'(x) = (e^x)*(cos(x)-sin(x)-sin(x)-cos(x)) = -2*(e^x)*sin(x)`
Using the equations derived in the previous section, we compute the numerical derivatives. Usually, it would be better to define a function for `f(x)` and call the function at different values in the main Matlab code. However, since in this case we use a simple function, directly writing the equations in the main function is computationally less expensive.
% Numerical Derivative
% Central Difference Scheme
central_difference = (1./dx.^2).*((-1./12).*(exp(x-2.*dx).*cos(x-2.*dx)) + (4./3).*(exp(x-dx).*cos(x-dx)) - (5./2).*(exp(x).*cos(x)) + (4./3).*(exp(x+dx).*cos(x+dx)) - (1./12).*(exp(x+2.*dx).*cos(x+2.*dx)));
% Skewed Right-Sided Difference Scheme
skewed_right = (1./dx.^2).*((35./12).*(exp(x).*cos(x)) - (26./3).*(exp(x+dx).*cos(x+dx)) + (19./2).*(exp(x+2.*dx).*cos(x+2.*dx)) - (14./3).*(exp(x+3.*dx).*cos(x+3.*dx)) + (11./12).*(exp(x+4.*dx).*cos(x+4.*dx)));
% Skewed Left-Sided Difference Scheme
skewed_left = (1./dx.^2).*((35./12).*(exp(x).*cos(x)) - (26./3).*(exp(x-dx).*cos(x-dx)) + (19./2).*(exp(x-2.*dx).*cos(x-2.*dx)) - (14./3).*(exp(x-3.*dx).*cos(x-3.*dx)) + (11./12).*(exp(x-4.*dx).*cos(x-4.*dx)));
Element-wise operators are used wherever multiplication(`*`), division(`/`) or power(`^`) operators are used. This is because the grid-size, `dx` is a vector. We then compute the error for each numerical scheme.
% Error Computation
error_central = abs(analytical_derivative - central_difference);
error_right = abs(analytical_derivative - skewed_right);
error_left = abs(analytical_derivative - skewed_left);
To plot the error curves in a logarithmic plot,
plot1 = loglog(dx, error_central, dx, error_left, dx, error_right);
set(plot1, \"linewidth\", 4);
title(\"Error - Central Difference vs Skewed Left-Sided Difference vs Skewed Right-Sided Difference\", \"fontsize\", 28);
xlabel(\"dx\", \"fontsize\", 28);
ylabel(\"Error\", \"fontsize\", 28);
grid on;
legend1 = legend(\"Central Difference Derivative\", \"Skewed Left-Sided Derivative\", \"Skewed Right-Sided Derivative\", \"location\", \'northwest\');
set(legend1, \"fontsize\", 14);
This gives us the previously shown plot.