!====================================================================
! 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
!====================================================================
program ch06_diff_richardson
implicit none

double precision x, h
double precision deriv_h, deriv_h2, deriv_r
double precision error_est, exact, actual_error
double precision test_function
integer m, status
external test_function

!
! 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.0d0
h = 0.2d0

!
! Step 2: calculate the first- and second-order derivatives
!
write (*,200)
write (*,201) x, h
write (*,202)

do m=1,2
  call richardson_derivative(test_function, x, h, m, &
       deriv_h, deriv_h2, deriv_r, error_est, status)

  if (status /= 0) then
    write (*,203) m, status
    cycle
  end if

  if (m == 1) then
    exact = 1.0d0/x
  else
    exact = -1.0d0/x**2
  end if

  actual_error = abs(deriv_r-exact)

  write (*,204) m, exact, deriv_h, deriv_h2, deriv_r, &
                error_est, actual_error
end do

200 format ('Centered differences with Richardson extrapolation')
201 format ('f(x) = ln(x),  x = ',f7.3,',  h = ',f7.3,/)
202 format (' order       exact            D(h)          D(h/2)', &
            '       Richardson     estimated error', &
            '        actual error')
203 format ('Derivative of order ',i2, &
            ' could not be evaluated; status = ',i2)
204 format (i6,6es17.8)

stop
end program ch06_diff_richardson


subroutine richardson_derivative(func, x, h, m, deriv_h, deriv_h2, &
                                  deriv_r, error_est, status)
!====================================================================
! 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).
!====================================================================
implicit none

double precision func, x, h
double precision deriv_h, deriv_h2, deriv_r, error_est
double precision correction
double precision centered_difference
integer m, status
external func

deriv_h   = 0.0d0
deriv_h2  = 0.0d0
deriv_r   = 0.0d0
error_est = 0.0d0
status    = 0

! Check the order of the requested derivative
if (m < 1 .or. m > 2) then
  status = 1
  return
end if

! A positive finite-difference step is required
if (h <= 0.0d0) then
  status = 2
  return
end if

! Calculate the centered differences using h and h/2
deriv_h  = centered_difference(func, x, h, m)
deriv_h2 = centered_difference(func, x, h/2.0d0, m)

! Remove the leading O(h**2) error
correction = (deriv_h2-deriv_h)/3.0d0
deriv_r = deriv_h2 + correction

! The Richardson correction provides an error estimate
error_est = abs(correction)

return
end subroutine richardson_derivative


function centered_difference(func, x, h, m)
!====================================================================
! Calculate a first- or second-order derivative using a second-order
! centered-difference formula.
!
! Input ...
! func  - function to be differentiated
! x     - point where the derivative is required
! h     - finite-difference step size
! m     - order of the derivative (1 or 2)
!
! Output ...
! centered_difference - numerical derivative
!====================================================================
implicit none

double precision centered_difference
double precision func, x, h
integer m
external func

if (m == 1) then
  centered_difference = (func(x+h)-func(x-h))/(2.0d0*h)
else
  centered_difference = &
       (func(x+h)-2.0d0*func(x)+func(x-h))/h**2
end if

return
end function centered_difference


function test_function(x)
!====================================================================
! Function used in the driver program. Replace this expression with
! the function required in a practical calculation.
!====================================================================
implicit none

double precision test_function, x

test_function = log(x)

return
end function test_function