Deriving fourth order scheme for central and skewed difference approximations for a second order derivative using Taylor s table and graphically comparing the errors resulting from these approaches

Objective

To derive the fourth order approximation schemes for the second order derivative of a function using central differencing, skewed right side differencing and skewed left side differencing approach and to draw a comparison between the errors arising from these approaches along with forward and backward differencing approaches. 

Derivation

1. Central Difference

To compute fourth order central difference approximation at point x of a function f, we must have the function defined at 4 other points surrounding x, 2 on the left and 2 on the right. Let those points be (x-2dx), (x-dx), (x+dx) and (x+2dx), dx being the node step size. Then the fourth order central difference approximation for the second derivative of f is given by

`(d^2f)/(dx^2) = a f(x-2dx) + b f(x-dx) + c f(x) + d f(x+dx) + e f(x+2dx)`

The values for the coefficients a, b, c, d, and e can be computed with the help of Taylor's table.

In the above Taylor's table, we're ignoring the higher order derivatives as 5 equations (columns) are enough to calculate the values for 5 unknowns.

a + b + c + d + e = 0          -----(1)

-2a - b + 0 + d + 2e = 0      -----(2)

dx^2(4a + b + 0 + d + 4e)/2 = 1      -----(3)

-8a - b + 0 + d + 8e = 0      -----(4)

16a + b + 0 + d + 16e = 0   -----(5)

Let us omit the term dx^2 in the 3rd equation for convenience as we're using matrix computation. This term will be included in latter steps.

Now, writing these equations in matrix form, we get,

AX = B

Solving this matrix equation will give the values for all the coefficients.

It can be noted from the above computation that a = e and b = d. Also, the above values calculated must be divided by dx^2 to obtain the actual values of a,b,c,d, and e since we had omitted it earlier for our convenience.

Substituting these values, the fourth order central difference approximation for the second order derivative can be obtained,

`(d^2f)/(dx^2) = (-0.0833 f(x-2dx) + 1.3333 f(x-dx) - 2.5 f(x) + 1.3333 f(x+dx) - 0.0833 f(x+2dx))/(dx^2)`

How do we know that the scheme is of fourth order error?

Let us again expand the Taylor table, only this time, including the higher order terms

`(d^2 f)/(dx^2) = (4a + b + 0 + d + 4e) f''(x) (dx^2)/(2!) + (-32a - b + 0 + d + 32e) f^v(x) (dx^5)/(5!) + (64a + b + 0 + d + 64e) f^(vi) dx^6/(6!) + ...`

We know that the coefficient of f"(x) is 1 and since a = e and b = d, the fifth derivative term will be 0. Now we get,

`(d^2 f)/(dx^2) = (4a + b + 0 + d + 4e) f''(x) (dx^2)/(2!) + (64a + b + 0 + d + 64e) f^(vi) dx^6/(6!) + ...`

Now, substituting the initial definition for second derivative

`a f(x-2dx) + b f(x-dx) + c f(x) + d f(x+dx) + e f(x+2dx) = (4a + b + 0 + d + 4e) f''(x) (dx^2)/(2!) + (64a + b + 0 + d + 64e) f^(vi) dx^6/(6!) + ...`

Dividing the equation by dx^2,

`(a f(x-2dx) + b f(x-dx) + c f(x) + d f(x+dx) + e f(x+2dx))/(dx^2) = (4a + b + 0 + d + 4e) (f''(x)) /(2!) + (64a + b + 0 + d + 64e) f^(vi) dx^4/(6!) + ...`

In the RHS of the above equation, we can see that the first dx term has a power of 4, So if we ignore this term with all the higher order terms we get

`(a f(x-2dx) + b f(x-dx) + c f(x) + d f(x+dx) + e f(x+2dx))/(dx^2) = (4a + b + 0 + d + 4e) (f''(x)) /(2!) + O(Delta x^4)`

Therefore, it is assured that our sheme is of fourth order error.

2. Skewed right side difference

 To compute fourth order right skewed difference approximation at point x for the second order derivative of a function f, we must have the function defined at 5 other nodes to the right of x. Let those points be (x+dx), (x+2dx), (x+3dx), (x+4dx) and (x+5dx), dx being the node step size. Then the fourth order right skewed difference approximation for the second derivative of f is given by

`(d^2f)/(dx^2) = a f(x) + b f(x+dx) + c f(x+2dx) + d f(x+3dx) + e f(x+4dx) + h f(x+5dx)`

The values for the coefficients a, b, c, d, e and h can be computed with the help of Taylor's table.

We get 6 coupled algebraic equations for 6 unknowns ignoring the higher order terms from Taylor's table.

a + b + c + d + e + h = 0                        -----(1)

0 + b + 2c + 3d + 4e + 5h = 0                 -----(2)

dx^2(0 + b + 4c + 9d + 16e + 25h)/2 = 1 -----(3)

0 + b + 8c + 27d + 64e + 125h = 0          -----(4)

0 + b + 16c + 81d + 256e + 625h = 0      -----(5)

0 + b + 32c + 243d + 1024e + 3125h = 0 -----(6)

As done earlier, we're omitting the dx^2 term for convenience and will be included in the final step.

Writing these equations in matrix form, we get,

AX = B

Solving this matrix equation will give the values for all the coefficients.

Therefore, the fourth order right skewed difference approximation for the second order derivative is given by

`(d^2f)/(dx^2) = (2 f(x+dx) + 8 f(x+2dx) + 18 f(x+3dx) + 32 f(x+4dx) + 50 f(x+5dx))/(dx^2)`

Proof for the order of error

Expanding the Taylor table 

`(d^2 f)/(dx^2) = (0 + b + 4c + 9d + 15e) f''(x) (dx^2)/(2!) + (0 + b + 64c + 729d + 4096e + 15625 h) f^(vi) dx^6/(6!) + ...`

Substituting the computed values of the constants in the second derivative term, we get,

`(d^2 f)/(dx^2) = f''(x) dx^2 + (0 + b + 64c + 729d + 4096e + 15625 h) f^(vi) dx^6/(6!) + ...`

`a f(x) + b f(x+dx) + c f(x+2dx) + d f(x+3dx) + e f(x+4dx) + h f(x+5dx) = f''(x) dx^2 + (0 + b + 64c + 729d + 4096e + 15625 h) f^(vi) dx^6/(6!) + ...`

Dividing by dx^2, we get

`(a f(x) + b f(x+dx) + c f(x+2dx) + d f(x+3dx) + e f(x+4dx) + h f(x+5dx))/(dx^2) = f''(x) + (0 + b + 64c + 729d + 4096e + 15625 h) f^(vi) dx^4/(6!) + ...`

Ignoring the sixth order and above derivative terms, we get

`(a f(x) + b f(x+dx) + c f(x+2dx) + d f(x+3dx) + e f(x+4dx) + h f(x+5dx))/(dx^2) = f''(x) + O(Delta x^4)`

Therefore, it is assured that the scheme is fourth order accurate.

3. Skewed left side difference

To compute fourth order left skewed difference approximation at point x of a function f, we must have the function defined at 5 other points to the left of x. Let those points be (x-dx), (x-2dx), (x-3dx), (x-4dx) and (x-5dx), dx being the node step size. Then the fourth order left skewed difference approximation for the second derivative of f is given by

`(d^2f)/(dx^2) = a f(x-5dx) + b f(x-4dx) + c f(x-3dx) + d f(x-2dx) + e f(x-dx) + h f(x)`

The values for the coefficients a, b, c, d, e and h can be computed with the help of Taylor\'s table.

We get 5 coupled algebraic equations from Taylor\'s table.

a + b + c + d + e + h = 0                          -----(1)

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

dx^2(25 a + 16b + 9c + 4d + e + 0)/2 = 1 -----(3)

-125a - 64b - 27c - 8d - e + 0 = 0               -----(4)

625a + 256b + 81c + 16d + e + 0 = 0         -----(5)

-3125a - 1024b - 243c - 32d - e + 0 = 0      -----(6)

As done earlier, we're omitting the dx^2 term for convenience and will be included in the final step.

Writing these equations in matrix form, we get,

AX = B

Solving this matrix equation will give the values for all the coefficients.

Therefore, the fourth-order left-skewed difference approximation for the second-order derivative is given by

`(d^2f)/(dx^2) = (-0.83333 f(x-5dx) + 5.08333 f(x-4dx) - 13 f(x-3dx) + 17.83333 f(x-2dx) - 12.83333 f(x-dx) + 3.75 f(x))/dx^2`

From both skewed left and skewed right schemes, we can see that the absolute values of constants are the same just in reverse order. That is expected because both these schemes are of the same order of accuracy since both take information from 5 points to compute the value of the function at a point. Therefore, we do not need to prove that this scheme is of fourth-order error separately because we have already proved for the skewed right scheme.

Matlab Code

% Computation of higher order approximations using symmetric and skewed
% stencils for the function, f(x) = exp(x)*cos(x)

clear all
close all
clc

% Assuming suitable values for x and dx
x = 2*pi;
dx = linspace(pi/4, pi/400000, 60);

% Declaring the size of variables for later use
cd_error = zeros(1, length(dx));
srd_error = zeros(1, length(dx));
sld_error = zeros(1, length(dx));
fwd_error = zeros(1, length(dx));
bwd_error = zeros(1, length(dx));

% Actual or analytical solution for f''(x) = -2*exp(x)*sin(x)
act_solution = -2*exp(x)*sin(x);

% loop to compute various numerical solutions for each value of dx 
for i = 1: length(dx)
    f = exp(x)*cos(x);                      % f(x)- given function
    f1 = exp(x+dx(i))*cos(x+dx(i));         % f(x+1)
    f2 = exp(x+2*dx(i))*cos(x+2*dx(i));     % f(x+2)
    f3 = exp(x+3*dx(i))*cos(x+3*dx(i));     % f(x+3)
    f4 = exp(x+4*dx(i))*cos(x+4*dx(i));     % f(x+4)
    f5 = exp(x+5*dx(i))*cos(x+5*dx(i));     % f(x+5)
    f1m = exp(x-dx(i))*cos(x-dx(i));        % f(x-1)
    f2m = exp(x-2*dx(i))*cos(x-2*dx(i));    % f(x-2)
    f3m = exp(x-3*dx(i))*cos(x-3*dx(i));    % f(x-3)
    f4m = exp(x-4*dx(i))*cos(x-4*dx(i));    % f(x-4)
    f5m = exp(x-5*dx(i))*cos(x-5*dx(i));    % f(x-4)

    % computing solutions of f(x) using various differencing schemes by
    % calling their respective functions
    cd_solution = central_difference(f2m, f1m, f, f1, f2, dx(i));
    srd_solution = skewed_right_difference(f, f1, f2, f3, f4, f5, dx(i));
    sld_solution = skewed_left_difference(f5m, f4m, f3m, f2m, f1m, f, dx(i));
    fwd_solution = forward_difference(f, f1, f2, dx(i));
    bwd_solution = backward_difference(f2, f1, f, dx(i));

    % calculation of absolute error of each differencing scheme 
    cd_error(i) = abs(act_solution - cd_solution);
    srd_error(i) = abs(act_solution - srd_solution);
    sld_error(i) = abs(act_solution - sld_solution);
    fwd_error(i) = abs(act_solution - fwd_solution);
    bwd_error(i) = abs(act_solution - fwd_solution);
end

% plotting graph between dx and various differencing schemes
figure(1)
loglog(dx, cd_error, dx, srd_error, dx, sld_error, dx, fwd_error,...
    dx, bwd_error,'*');
xlabel('Value of dx');
ylabel('Value of error for x = 2*pi');
legend('Central differencing error','Skewed right side differencing error'...
    ,'Skewed left side differencing error','Forward differencing error',...
    'Backward differencing error');

% Function definition for Central Difference Approximation

function cd_solution = central_difference(f2m, f1m, f, f1, f2, dx)

% values derived through Taylor's table
a = -0.0833;
b = 1.3333;
c = -2.5000;
d = 1.3333;
e = -0.0833;

% central difference scheme for 4th order solution to 2nd order derivative
cd_solution = (a*f2m + b*f1m + c*f + d*f1 + e*f2)/dx^2;
end

% Function definition for Skewed Left Side Difference Approximation

function sld_solution = skewed_left_difference(f5m, f4m, f3m, f2m, f1m, f, dx)

% values derived through Taylor's table
a = -0.83333;
b = 5.08333;
c = -13.00000;
d = 17.83333;
e = -12.83333;
h = 3.75000;

% skewed left side difference scheme for 4th order solution to 2nd
% order derivative
sld_solution = (a*f5m + b*f4m + c*f3m + d*f2m + e*f1m + h*f)/dx^2;
end

% Function definition for Skewed Right Side Difference Approximation

function srd_solution = skewed_right_difference(f, f1, f2, f3, f4, f5, dx)

% values derived through Taylor's table
a = 3.75000;
b = -12.83333;
c = 17.83333;
d = -13.00000;
e = 5.08333;
h = -0.83333;

% skewed right side difference scheme for 4th order solution to 2nd
% order derivative
srd_solution = (a*f + b*f1 + c*f2 + d*f3 + e*f4 + h*f5)/dx^2;
end

% Function definition for Backward Difference Approximation

function bwd_solution = backward_difference(f2, f1, f, dx)

% backward difference scheme for 2nd order solution for 2nd order derivative
bwd_solution = (f2 - 2*f1 + f)/dx^2;
end

% Function definition for Forward Difference Approximation

function fwd_solution = forward_difference(f, f1, f2, dx)

% forward difference scheme for 2nd order solution for 2nd order derivative
fwd_solution = (f - 2*f1 + f2)/dx^2;
end

Output

Inference

From the above comparison of different schemes for fourth order approximation of a second-order derivative, we can reassure that the central difference approximation yields more accurate solutions than the right-skewed and left-skewed differencing schemes. Since forward and backward approximations are of second-order error, they give out larger errors than the other fourth-order schemes. 

It may also be noted that if we consider a wider range of dx, at some values of dx, the error produced by higher order schemes is bigger than that of the lower order schemes. This is due to higher roundoff errors at those points. It is important to understand that the value of error at any point is the net sum of numerical error and roundoff error and the round off error increases with increase in the order of error. Therefore, a higher-order scheme may not be always desirable.

Significance of a skewed scheme

A function is always characterized by a domain. A domain is a space or area beyond which a function ceases to exist. Therefore a domain has boundaries and the function exists on only one side of the boundary. Now, if we want to compute a numerical approximation for the function at the boundary, we cannot use a central difference scheme as it requires data from both sides of the boundary. That is where a skewed difference scheme plays a vital role since data from any side can be used to perform numerical approximation at a point.