Fourth Order Approximation of the Second Derivative

Results and Discussion

FIGURE 1: Error versus grid size for fourth-order approximation for second derivative using Central, Right Skewed, Left Skewed, Forward Difference, and Back Difference Scheme.

In the above plot, Figure 1, we can see that the lower value of delta x, the error decreases. However, upon closer look, we can see that by making delta x too small increases the error, Figure 2.

FIGURE 2: Error vs grid size for different second derivative approximation.

Comparing Forward and Backward scheme with other approximations, the Skewed Left-side, and the Central Difference are more accurate than the other three. Intuitively it makes sense as we increase the size of delta x, the error increases as well. But why is an error for one scheme different than others? That is because of how the way the schemes work. Central Difference takes the approach of using two points on the left and two points on the right as shown in stencil figure 1. Whereas, a skewed scheme uses points from either the front (right) or back (left) of the point of evaluation. Using data from left and right would yield more accurate result than just using the right or left.

However, Central Difference Scheme (CDS) cannot be used in all situation. If there is no data available to the left or right of the "i", CDS would give an inaccurate result. A common example would be near a boundary since outside data is not available, only the interior. In such cases, RS or LS should be used. Therefore, which scheme to use depends on the problem.

Matlab Code

clc; clear; close all;

%% Central Difference
% Nodes:  i-2, i-1, i, i+1, i+2
% Function: f''(i) = af(i-2) + bf(i-1) + cf(i) + df(i+1) + ef(i+2)
Ac = [1     1       1       1       1;
    -2      -1      0       1       2;
    2^2/2   1^2/2   0       1^2/2   2^2/2;
  -2^3/6   -1^3/6   0       1^3/6   2^3/6;
    2^4/24  1^4/24  0       1^4/24  2^4/24;];
    
bc = [0 0 1 0 0]';
xc = Ac\bc; %Solving the linear system

% Coefficients a,b,c,d,e 
c1 = xc(1); c2 = xc(2); c3 = xc(3); c4 = xc(4); c5 = xc(5); 


%% Skewed Right Side
% Nodes: i, i+1, i+2, i+3, i+4, i+5
% Function: f''(i) = af(i) + bf(i+1) + cf(i+2) + df(i+3) + ef(i+4) + gf(i+5)
Ar = [1      1       1       1       1       1;
    0        1       2       3       4       5;
    0       1^2/2   2^2/2   3^2/2   4^2/2   5^2/2;
    0       1^3/6   2^3/6   3^3/6   4^3/6   5^3/6;
    0       1^4/24  2^4/24  3^4/24  4^4/24  5^4/24;
    0       1^5/120 2^5/120 3^5/120 4^5/120 5^5/120;];
br = [0 0 1 0 0 0]';
xr = Ar\br;

% Coefficients a,b,c,d,e,g
r1 = xr(1); r2 = xr(2); r3 = xr(3); r4 = xr(4); r5 = xr(5); r6 = xr(5);


%% Skewed Left Side
% Nodes: i, i-1, i-2, i-3, i-4, i-5
% Function: f''(i) = af(i) + bf(i-1) + cf(i-2) + df(i-3) + ef(i-4) + gf(i-5)
Al = [1      1       1       1       1       1;
    0       -1       -2       -3       -4       -5;
    0       1^2/2   2^2/2   3^2/2   4^2/2   5^2/2;
    0      -1^3/6  -2^3/6  -3^3/6   -4^3/6   -5^3/6;
    0       1^4/24  2^4/24  3^4/24  4^4/24  5^4/24;
    0     -1^5/120 -2^5/120 -3^5/120 -4^5/120 -5^5/120;];
bl = [0 0 1 0 0 0]';
xl = Al\bl;

% Coefficients a,b,c,d,e,g
l1 = xl(1); l2 = xl(2); l3 = xl(3); l4 = xl(4); l5 = xl(5); l6 = xl(6);


%% Forward Difference Scheme
% Nodes: i, i+1, i+2
% Function: f''(i) = af(i) + bf(i+1) + cf(i+2);
Af = [1     1       1;
    0       -1      -2;
    0       1/2     2;];
bf = [0 0 1]';
xf = Af\bf;

% Coefficients a,b,c
f1 = xf(1); f2 = xf(2); f3 = xf(3);


%% Backward Difference Scheme
% Nodes: i, i-1, i-2
% Function: f''(i) = af(i) + bf(i-1) + cf(i-2);
Ab = [1     1       1;
    0       -1      -2;
    0       1/2     2;];
bb = [0 0 1]';
xb = Ab\bb;

% Coefficients a,b,c
b1 = xb(1); b2 = xb(2); b3 = xb(3);


%% Function
f = @(x) exp(x)*cos(x); % Given
x = 0; % From previous challenge
% spacing = linspace(1E-6,10,251); % dx range from 1E-6 to 10
spacing = 1E-5:0.01:100;

% Analytically solution 
d2dx2 = -2*exp(x)*sin(x); % Second derivative of the function

for i = 1:length(spacing)
    dX = spacing(i);
    
    d2dx2_1 = (c1*f(x-2*dX) + c2*f(x-1*dX) + c3*f(x) + c4*f(x+1*dX) + c5*f(x+2*dX))/(dX^2);
    d2dx2_2 = (r1*f(x) + r2*f(x+1*dX) + r3*f(x+2*dX) + r4*f(x+3*dX) + r5*f(x+4*dX) + r6*f(x+5*dX))/(dX^2);
    d2dx2_3 = (l1*f(x) + l2*f(x-1*dX) + l3*f(x-2*dX) + l4*f(x-3*dX) + l5*f(x-4*dX) + l6*f(x-5*dX))/(dX^2);
    d2dx2_4 = (f1*f(x) + f2*f(x+1*dX) + f3*f(x+2*dX))/(dX^2);
    d2dx2_5 = (b1*f(x) + b2*f(x-1*dX) + b3*f(x-2*dX))/(dX^2);
    
    
    Error(1,i) = abs(d2dx2_1 - d2dx2);
    Error(2,i) = abs(d2dx2_2 - d2dx2);
    Error(3,i) = abs(d2dx2_3 - d2dx2);
    Error(4,i) = abs(d2dx2_4 - d2dx2);
    Error(5,i) = abs(d2dx2_5 - d2dx2);
     
end