%=====================================================================
% Driver program: Calculates first- and second-order derivatives
% Method: centered differences with Richardson extrapolation
%
% The derivatives are first calculated with step sizes h and h/2.
% Richardson extrapolation then removes the leading O(h^2) error.
%
% Alexander Godunov, September 2026
%=====================================================================
clear;
clc;

%---------------------------------------------------------------------
% Step 1: specify the evaluation point and initial step size
%
% The example uses
%
%       f(x) = ln(x),
%
% for which f'(x) = 1/x and f''(x) = -1/x^2.
%---------------------------------------------------------------------
x = 2.0;
h = 0.2;

% Function to be differentiated
func = @testFunction;

%---------------------------------------------------------------------
% Step 2: calculate the first- and second-order derivatives
%---------------------------------------------------------------------
fprintf('Centered differences with Richardson extrapolation\n');
fprintf('f(x) = ln(x),  x = %7.3f,  h = %7.3f\n\n',x,h);

fprintf([' order       exact            D(h)          D(h/2)', ...
         '       Richardson     estimated error', ...
         '        actual error\n']);

for m=1:2
    [deriv_h, deriv_h2, deriv_r, error_est, status] = ...
        richardsonDerivative(func,x,h,m);

    if status ~= 0
        fprintf(['Derivative of order %d could not be evaluated; ', ...
                 'status = %d\n'],m,status);
        continue;
    end

    if m == 1
        exact = 1.0/x;
    else
        exact = -1.0/x^2;
    end

    actual_error = abs(deriv_r-exact);

    fprintf('%6d%17.8e%17.8e%17.8e%17.8e%17.8e%17.8e\n', ...
            m,exact,deriv_h,deriv_h2,deriv_r,error_est,actual_error);
end


function [deriv_h,deriv_h2,deriv_r,error_est,status] = ...
         richardsonDerivative(func,x,h,m)
%=====================================================================
% Calculate a first- or second-order derivative using a centered
% difference and one level of Richardson extrapolation.
%
% Input ...
% func       - function to be differentiated
% x          - point where the derivative is required
% h          - initial finite-difference step size
% m          - order of the derivative (1 or 2)
%
% Output ...
% deriv_h    - centered difference calculated with h
% deriv_h2   - centered difference calculated with h/2
% deriv_r    - Richardson-extrapolated derivative
% error_est  - magnitude of the Richardson correction
% status     - 0: success
%              1: m is not 1 or 2
%              2: h is not positive
%
% A centered difference has an error proportional to h^2. For step
% refinement by a factor of 2, Richardson extrapolation gives
%
%       D_R = D(h/2) + [D(h/2)-D(h)]/(2^2-1).
%=====================================================================
deriv_h = 0.0;
deriv_h2 = 0.0;
deriv_r = 0.0;
error_est = 0.0;
status = 0;

% Check the order of the requested derivative
if m < 1 || m > 2
    status = 1;
    return;
end

% A positive finite-difference step is required
if h <= 0.0
    status = 2;
    return;
end

% Calculate the centered differences using h and h/2
deriv_h = centeredDifference(func,x,h,m);
deriv_h2 = centeredDifference(func,x,h/2.0,m);

% Remove the leading O(h^2) error
correction = (deriv_h2-deriv_h)/3.0;
deriv_r = deriv_h2 + correction;

% The Richardson correction provides an error estimate
error_est = abs(correction);
end


function value = centeredDifference(func,x,h,m)
%=====================================================================
% Calculate a first- or second-order derivative using a second-order
% centered-difference formula.
%=====================================================================
if m == 1
    value = (func(x+h)-func(x-h))/(2.0*h);
else
    value = (func(x+h)-2.0*func(x)+func(x-h))/h^2;
end
end


function value = testFunction(x)
%=====================================================================
% Function used in the driver program. Replace this expression with
% the function required in a practical calculation.
%=====================================================================
value = log(x);
end