Project

Derivation and analysis of higher-order Finite Difference schemes using Taylor table method

AAngel Gonzalez Llacer29 May 2019

In this project, centered and skewed Finite Differences schemes of 4th order in accuracy will be derived for the 2nd derivative, using the Taylor Table method, and an error comparison will be proposed against lower-order schemes.

The Taylor table method is used to find schemes with a desired accuracy and for a stencil that we choose. It is based on Taylor series.

For example, in the centered stencil that uses `f_(i-2), f_(i-1), f_i, f_(i+1), f_(i+2)`, the Taylor series for each of the terms is:

`f_(i+1) approx f_i + h*f\'_i + h^2/2*f\'\'_i + h^3/6*f_i^((3))+h^4/24*f_i^((4))+h^5/120*f_i^((5)) + O(h^6)`

Substituting 2*h instead of h in the previous formula, we will get the Taylor series of `f_(i+2)`and so on:

`f_(i+2) approx f_i + 2h*f\'_i + (4h^2)/2*f\'\'_i + (8h^3)/6*f_i^((3))+(16h^4)/24*f_i^((4))+(32h^5)/120*f_i^((5)) + O(h^6)`

The centered scheme for 2nd derivative will be of the form: `(a*f_(i-2) + b*f_(i-1) + c*f_i + d*f_(i+1) + e*f_(i+2))/h^2`

So we can use the Taylor expressions for all terms `f_(i+k)`, having for unknowns a,b,c,d,e.

The equations will be:

`f_i = 0`
`f_i\' = 0`
`f_i\'\' = 1`
`f_i^((3)) = 0`
`f_i^((4))=0`

So we have 5 unkwons and 5 equations, the system is determined. As we cancel up to the 4th derivative, and we devide by h^2 to obtain the 2nd derivative, the minimum order of accuracy of the scheme is 3, but we will show that another equation cancelling the 5th derivative will be verified, making the scheme to be of 4th order.

We use the following MATLAB code to solve the coefficients a,b,c,d,e

%% Taylor table for 2nd order derivative schemes
% Derivation and analysis of higher-order Finite Difference schemes using Taylor table method
% @Angel Gonzalez 2019

%% 1. Centered scheme with 4th order error

% fxx = (a*f_i-2 + b*f_i-1 + c*f_i + d*f_i+1 + e*f_i+2)/h^2;

% Equations:
% 1) fi = 0:     a + b + c + d + e = 0
% 2) f'i = 0:    -2*a -1*b + 0*c + 1*d + 2*e = 0
% 3) f''i = 1:   2*a + 1/2*b + 0*c + 1/2*d + 2*e = 1
% 4) f'''i = 0:  -8*a - 1*b + 0*c + 1*d + 8*e = 0
% 5) f''''i = 0: 16*a + 1*b + 0*c + 1*d + 16*e = 0

A1 = [1 1 1 1 1; ...
    -2 -1 0 1 2; ...
    2 1/2 0 1/2 2; ...
    -8 -1 0 1 8; ...
    16 1 0 1 16];

b1 = [0; 0; 1; 0; 0];

coeffsCentered = mldivide(A1,b1); % and coeffsCentered = [a,b,c,d,e]

eq6cent = sum([-32; -1; 0; 1; 32].*coeffsCentered);

% We observe eq5 approx equal to 0 -> it is 4th order (not 3rd order)

We verify at the end that the 6th equation that ensures 4th order, which is `f_i^((5))=0`is also achieved naturally by the centered scheme.

For the skewed schemes, this last equation with 5-point stencils will not be verified, and we will need to add 1 more point and 1 more equation (6th equation), so the stencil will be of 6 points.

The centered scheme for 2nd derivative is hence:

`f_(x x) = (-1f_(i-2) + 16 f_(i-1) -30 f_i + 16 f_(i+1) -1 f_(i+2))/(12 h^2)`

For the skewed right scheme we use the stencil [-1, 0, 1, 2, 3, 4] and the MATLAB code is:

%% 2. Skewed right scheme with 4th order error 

% Needs 6 points for 4th order accuracy
% fxx = (a*f_i-1 + b*f_i + c*f_i+1 + d*f_i+2 + e*f_i+3 + f*f_i+4)/h^2;

% Equations:
% 1) fi = 0:      a + b + c + d + e + f = 0
% 2) f'i = 0:     -1*a 0*b + 1*c + 2*d + 3*e + 4*f = 0
% 3) f''i = 1:    1/2*a + 0*b + 1/2*c + 2*d + 9/2*e + 16/2*f = 1
% 4) f'''i = 0:   -1*a + 0*b + 1*c + 8*d + 27*e + 64*f = 0
% 5) f''''i = 0:  1*a + 0*b + 1*c + 16*d + 81*e + 256*f = 0
% 6) f'''''i = 0: -1*a + 0*b + 1*c + 32*d + 243*e + 1024*f = 0

A2 = [1 1 1 1 1 1; ...
    -1 0 1 2 3 4; ...
    1/2 0 1/2 2 9/2 16/2; ...
    -1 0 1 8 27 64; ...
    1 0 1 16 81 256; ...
    -1 0 1 32 243 1024];

b2 = [0; 0; 1; 0; 0; 0];

coeffsRight = mldivide(A2,b2); % and coeffsRight = [a,b,c,d,e,f]

eq6right = sum([-1; 0; 1; 32; 243; 1024].*coeffsRight);

Equation 6 is forced to be verified by adding it to the system of equations.

The skewed right scheme for 2nd derivative is hence:

`f_(x x) = (10f_(i-1) - 15 f_(i) -4 f_(i+1) + 14 f_(i+2) -6 f_(i+3) + 1 f_(i+4))/(12 h^2)`

For the skewed left scheme we use the stencil [-4, -3, -2, -1, 0, 1] and the MATLAB code is:

%% 3. Skewed left scheme with 4th order error

% Needs 6 points for 4th order accuracy
% fxx = (a*f_i-4 + b*f_i-3 + c*f_i-2 + d*f_i-1 + e*f_i + f*f_i+1)/h^2;


A3 = [1 1 1 1 1 1; ...
    -4 -3 -2 -1 0 1; ...
    16/2 9/2 2 1/2 0 1/2; ...
    -64 -27 -8 -1 0 1; ...
    256 81 16 1 0 1; ...
    -1024 -243 -32 -1 0 1];

b3 = [0; 0; 1; 0; 0; 0];

coeffsLeft = mldivide(A3,b3); % and coeffsLeft = [a,b,c,d,e,f]

eq6left = sum([-1024; -243; -32; -1; 0; 1].*coeffsLeft);

Equation 6 is also forced to be verified by adding it to the system of equations.

The skewed left scheme for 2nd derivative is hence:

`f_(x x) = (1f_(i-4) - 6 f_(i-3) +14 f_(i-2) - 4 f_(i-1) - 15 f_(i) + 10 f_(i+1))/(12 h^2)`

Now that we have derived the schemes, let\'s put them into practice on a particular function, which will be `exp(x)*cos(x)`

We will evaluate the 2nd derivative on the range x = [-2, 2], with dx fixed to 0.1

The exact 2nd derivative is `-2 exp(x) sin(x)`

The following MATLAB code is used to plot the results of different schemes and the error respect to the exact solution:

%% 4. Evaluate on exp(x)*cos(x) and compare

% Parameter definition
fun = @(x) exp(x).*cos(x);
dx = 0.1;
xvec = linspace(-2, 2, 20);

% Compute exact and numerical derivatives
fundotdotexact = @(x) -2*exp(x).*sin(x);
fdotdotexact = fundotdotexact(xvec);

fdotdot1 = centered(fun, xvec, dx, coeffsCentered);
fdotdot2 = skewright(fun, xvec, dx, coeffsRight);
fdotdot3 = skewleft(fun, xvec, dx, coeffsLeft);

% Compute errors
err1 = abs(fdotdot1-fdotdotexact);
err2 = abs(fdotdot2-fdotdotexact);
err3 = abs(fdotdot3-fdotdotexact);

% Plots
figure();
plot(xvec, fdotdotexact, '-');
hold on
plot(xvec, fdotdot1, 'x');
plot(xvec, fdotdot2, '*');
plot(xvec, fdotdot3, 'o');
title('2nd derivative of exp(x)*cos(x)')
xlabel('x');
ylabel('2nd derivative of exp(x)*cos(x)');
legend({'exact', 'centered', 'skewed right', 'skewed left'});

figure()
plot(xvec, err1, 'Marker','x');
hold on
plot(xvec, err2, 'Marker', '*');
plot(xvec, err3, 'Marker', 'o');
title('Error on 2nd derivative of exp(x)*cos(x)')
xlabel('x');
ylabel('Absolute error');
legend({'centered', 'skewed right', 'skewed left'});

%% Scheme definitions
function fdotdot = centered(fun, x, dx, c)
    fdotdot = (fun(x-2*dx).*c(1) + fun(x-dx).*c(2) + fun(x).*c(3) ...
        + fun(x+dx).*c(4) + fun(x+2*dx).*c(5))./dx.^2;
end

function fdotdot = skewright(fun, x, dx, c)
    fdotdot = (fun(x-dx).*c(1) + fun(x).*c(2) + fun(x+dx).*c(3) ...
        + fun(x+2*dx).*c(4) + fun(x+3*dx).*c(5) ...
        + fun(x+4*dx).*c(6))./dx.^2;
end

function fdotdot = skewleft(fun, x, dx, c)
    fdotdot = (fun(x-4*dx).*c(1) + fun(x-3*dx).*c(2) + fun(x-2*dx).*c(3) ...
        + fun(x-dx).*c(4) + fun(x).*c(5) ...
        + fun(x+dx).*c(6))./dx.^2;
end

The results are the following plots:

We can see that, even if all schemes are of the same order of accuracy (4th order), the centered scheme has an advantage for almost all values of x that we have computed. The skewed left and skewed right have similar behavior, with fast increasing error for positive values of x.

Let\'s now compare this schemes to the 1st order FDS and BDS, which are:

FDS: `f_(x x) = (f_i - 2 f_(i+1) + f_(i+2))/h^2`
BDS: `f_(x x) = (f_(i-2) - 2 f_(i-1) + f_(i))/h^2`

We will compare them by varying the stepsize dx for a fixed value x = 2, so that we can verify the order of accuracy of each scheme:

%% 5. Compare to FDS and BDS

% Parameter definition
xval = 2;
dxvec = logspace(-2, 0, 10);

fxxexact = fundotdotexact(xval);
fxx1 = centered(fun, xval, dxvec, coeffsCentered);
fxx2 = skewright(fun, xval, dxvec, coeffsRight);
fxx3 = skewleft(fun, xval, dxvec, coeffsLeft);
fxxFDS = FDS(fun, xval, dxvec);
fxxBDS = BDS(fun, xval, dxvec);

%Compute errors
err1vec = abs(fxx1-fxxexact);
err2vec = abs(fxx2-fxxexact);
err3vec = abs(fxx3-fxxexact);
errFDSvec = abs(fxxFDS-fxxexact);
errBDSvec = abs(fxxBDS-fxxexact);

%Plot
figure()
loglog(dxvec, err1vec, 'Marker', 'x');
hold on
loglog(dxvec, err2vec, 'Marker', '*');
loglog(dxvec, err3vec, 'Marker', 'o');
loglog(dxvec, errFDSvec, 'r--');
loglog(dxvec, errBDSvec, 'g--');
hold off
title('Numerical error in function of step size')
xlabel('dx');
ylabel('Absolute error');
legend({'centered', 'right', 'left', 'FDS', 'BDS'}, ...
    'Location', 'southeast');



%% Scheme definitions
function fdotdot = centered(fun, x, dx, c)
    fdotdot = (fun(x-2*dx).*c(1) + fun(x-dx).*c(2) + fun(x).*c(3) ...
        + fun(x+dx).*c(4) + fun(x+2*dx).*c(5))./dx.^2;
end

function fdotdot = skewright(fun, x, dx, c)
    fdotdot = (fun(x-dx).*c(1) + fun(x).*c(2) + fun(x+dx).*c(3) ...
        + fun(x+2*dx).*c(4) + fun(x+3*dx).*c(5) ...
        + fun(x+4*dx).*c(6))./dx.^2;
end

function fdotdot = skewleft(fun, x, dx, c)
    fdotdot = (fun(x-4*dx).*c(1) + fun(x-3*dx).*c(2) + fun(x-2*dx).*c(3) ...
        + fun(x-dx).*c(4) + fun(x).*c(5) ...
        + fun(x+dx).*c(6))./dx.^2;
end

function fdotdot = FDS(fun, x, dx)
    fdotdot = (fun(x) - 2*fun(x+dx) + fun(x+2*dx))./dx.^2;
end

function fdotdot = BDS(fun, x, dx)
    fdotdot = (fun(x-2*dx) - 2*fun(x-dx) + fun(x))./dx.^2;
end

The resulting plot is:

We verify that the slope of the centered, skewed right and skewed left is 4, and for FDS and BDS is 1, verifying that the orders are according to the theoretical.

We highlight that the centererd scheme has an advantage in accuracy over the skewed schemes, even if they are of same order.

So we rise the question, what could be the advantage of using skewed schemes?

The advantage is that they can be used near boundaries or on the boundaries, where there is no access to further points over the boundary, so that we can still use a high order scheme in those regions, improving overall accuracy. However, for internal points of the mesh, we will tend to use centered schemes as they have an advantage on accuracy.

Full code of the project:

%% Taylor table for 2nd order derivative schemes
% Derivation and analysis of higher-order Finite Difference schemes using Taylor table method
% @Angel Gonzalez 2019

%% 1. Centered scheme with 4th order error

% fxx = (a*f_i-2 + b*f_i-1 + c*f_i + d*f_i+1 + e*f_i+2)/h^2;

% Equations:
% 1) fi = 0:     a + b + c + d + e = 0
% 2) f'i = 0:    -2*a -1*b + 0*c + 1*d + 2*e = 0
% 3) f''i = 1:   2*a + 1/2*b + 0*c + 1/2*d + 2*e = 1
% 4) f'''i = 0:  -8*a - 1*b + 0*c + 1*d + 8*e = 0
% 5) f''''i = 0: 16*a + 1*b + 0*c + 1*d + 16*e = 0

A1 = [1 1 1 1 1; ...
    -2 -1 0 1 2; ...
    2 1/2 0 1/2 2; ...
    -8 -1 0 1 8; ...
    16 1 0 1 16];

b1 = [0; 0; 1; 0; 0];

coeffsCentered = mldivide(A1,b1); % and coeffsCentered = [a,b,c,d,e]

eq6cent = sum([-32; -1; 0; 1; 32].*coeffsCentered);

% We observe eq5 approx equal to 0 -> it is 4th order (not 3rd order)

%% 2. Skewed right scheme with 4th order error 

% Needs 6 points for 4th order accuracy
% fxx = (a*f_i-1 + b*f_i + c*f_i+1 + d*f_i+2 + e*f_i+3 + f*f_i+4)/h^2;

% Equations:
% 1) fi = 0:      a + b + c + d + e + f = 0
% 2) f'i = 0:     -1*a 0*b + 1*c + 2*d + 3*e + 4*f = 0
% 3) f''i = 1:    1/2*a + 0*b + 1/2*c + 2*d + 9/2*e + 16/2*f = 1
% 4) f'''i = 0:   -1*a + 0*b + 1*c + 8*d + 27*e + 64*f = 0
% 5) f''''i = 0:  1*a + 0*b + 1*c + 16*d + 81*e + 256*f = 0
% 6) f'''''i = 0: -1*a + 0*b + 1*c + 32*d + 243*e + 1024*f = 0

A2 = [1 1 1 1 1 1; ...
    -1 0 1 2 3 4; ...
    1/2 0 1/2 2 9/2 16/2; ...
    -1 0 1 8 27 64; ...
    1 0 1 16 81 256; ...
    -1 0 1 32 243 1024];

b2 = [0; 0; 1; 0; 0; 0];

coeffsRight = mldivide(A2,b2); % and coeffsRight = [a,b,c,d,e,f]

eq6right = sum([-1; 0; 1; 32; 243; 1024].*coeffsRight);

%% 3. Skewed left scheme with 4th order error

% Needs 6 points for 4th order accuracy
% fxx = (a*f_i-4 + b*f_i-3 + c*f_i-2 + d*f_i-1 + e*f_i + f*f_i+1)/h^2;


A3 = [1 1 1 1 1 1; ...
    -4 -3 -2 -1 0 1; ...
    16/2 9/2 2 1/2 0 1/2; ...
    -64 -27 -8 -1 0 1; ...
    256 81 16 1 0 1; ...
    -1024 -243 -32 -1 0 1];

b3 = [0; 0; 1; 0; 0; 0];

coeffsLeft = mldivide(A3,b3); % and coeffsLeft = [a,b,c,d,e,f]

eq6left = sum([-1024; -243; -32; -1; 0; 1].*coeffsLeft);


%% 4. Evaluate on exp(x)*cos(x) and compare

% Parameter definition
fun = @(x) exp(x).*cos(x);
dx = 0.1;
xvec = linspace(-2, 2, 20);

% Compute exact and numerical derivatives
fundotdotexact = @(x) -2*exp(x).*sin(x);
fdotdotexact = fundotdotexact(xvec);

fdotdot1 = centered(fun, xvec, dx, coeffsCentered);
fdotdot2 = skewright(fun, xvec, dx, coeffsRight);
fdotdot3 = skewleft(fun, xvec, dx, coeffsLeft);

% Compute errors
err1 = abs(fdotdot1-fdotdotexact);
err2 = abs(fdotdot2-fdotdotexact);
err3 = abs(fdotdot3-fdotdotexact);

% Plots
figure();
plot(xvec, fdotdotexact, '-');
hold on
plot(xvec, fdotdot1, 'x');
plot(xvec, fdotdot2, '*');
plot(xvec, fdotdot3, 'o');
title('2nd derivative of exp(x)*cos(x)')
xlabel('x');
ylabel('2nd derivative of exp(x)*cos(x)');
legend({'exact', 'centered', 'skewed right', 'skewed left'});

figure()
plot(xvec, err1, 'Marker','x');
hold on
plot(xvec, err2, 'Marker', '*');
plot(xvec, err3, 'Marker', 'o');
title('Error on 2nd derivative of exp(x)*cos(x)')
xlabel('x');
ylabel('Absolute error');
legend({'centered', 'skewed right', 'skewed left'});

%% 5. Compare to FDS and BDS

% Parameter definition
xval = 2;
dxvec = logspace(-2, 0, 10);

fxxexact = fundotdotexact(xval);
fxx1 = centered(fun, xval, dxvec, coeffsCentered);
fxx2 = skewright(fun, xval, dxvec, coeffsRight);
fxx3 = skewleft(fun, xval, dxvec, coeffsLeft);
fxxFDS = FDS(fun, xval, dxvec);
fxxBDS = BDS(fun, xval, dxvec);

%Compute errors
err1vec = abs(fxx1-fxxexact);
err2vec = abs(fxx2-fxxexact);
err3vec = abs(fxx3-fxxexact);
errFDSvec = abs(fxxFDS-fxxexact);
errBDSvec = abs(fxxBDS-fxxexact);

%Plot
figure()
loglog(dxvec, err1vec, 'Marker', 'x');
hold on
loglog(dxvec, err2vec, 'Marker', '*');
loglog(dxvec, err3vec, 'Marker', 'o');
loglog(dxvec, errFDSvec, 'r--');
loglog(dxvec, errBDSvec, 'g--');
hold off
title('Numerical error in function of step size')
xlabel('dx');
ylabel('Absolute error');
legend({'centered', 'right', 'left', 'FDS', 'BDS'}, ...
    'Location', 'southeast');



%% Scheme definitions
function fdotdot = centered(fun, x, dx, c)
    fdotdot = (fun(x-2*dx).*c(1) + fun(x-dx).*c(2) + fun(x).*c(3) ...
        + fun(x+dx).*c(4) + fun(x+2*dx).*c(5))./dx.^2;
end

function fdotdot = skewright(fun, x, dx, c)
    fdotdot = (fun(x-dx).*c(1) + fun(x).*c(2) + fun(x+dx).*c(3) ...
        + fun(x+2*dx).*c(4) + fun(x+3*dx).*c(5) ...
        + fun(x+4*dx).*c(6))./dx.^2;
end

function fdotdot = skewleft(fun, x, dx, c)
    fdotdot = (fun(x-4*dx).*c(1) + fun(x-3*dx).*c(2) + fun(x-2*dx).*c(3) ...
        + fun(x-dx).*c(4) + fun(x).*c(5) ...
        + fun(x+dx).*c(6))./dx.^2;
end

function fdotdot = FDS(fun, x, dx)
    fdotdot = (fun(x) - 2*fun(x+dx) + fun(x+2*dx))./dx.^2;
end

function fdotdot = BDS(fun, x, dx)
    fdotdot = (fun(x-2*dx) - 2*fun(x-dx) + fun(x))./dx.^2;
end

Derivation and analysis of higher-order Finite Difference schemes using Taylor table method

Ship work like this
in your own course.

Talk to a career counsellor

Ship work like thisin your own course.

Talk to a career counsellor

Ship work like this
in your own course.