Numerical Methods

Interpolation (Lagrange, cubic splines), least-squares fitting, FFT (or call FFTW), Monte Carlo integration with random_number, finite-difference heat/wave equations, conjugate gradient, power iteration for eigenvalues, Gauss–Legendre quadrature, adaptive s…

37.1 Root finding — Newton–Raphson & bisection

module roots
  use, intrinsic :: iso_fortran_env, only: dp => real64
  implicit none
  abstract interface
    pure real(dp) function fun(x)
      import :: dp
      real(dp), intent(in) :: x
    end function fun
  end interface
contains
  function newton(f, df, x0, tol, maxit) result(x)
    procedure(fun) :: f, df
    real(dp), intent(in) :: x0, tol
    integer, intent(in) :: maxit
    real(dp) :: x, dx
    integer :: it
    x = x0
    do it = 1, maxit
      dx = f(x) / df(x)
      x = x - dx
      if (abs(dx) < tol) return
    end do
    error stop "newton: no convergence"
  end function newton

  function bisect(f, a0, b0, tol) result(c)
    procedure(fun) :: f
    real(dp), intent(in) :: a0, b0, tol
    real(dp) :: a, b, c
    a = a0; b = b0
    if (f(a) * f(b) > 0) error stop "bisect: root not bracketed"
    do while (b - a > tol)
      c = 0.5_dp * (a + b)
      if (f(a) * f(c) <= 0) then
        b = c
      else
        a = c
      end if
    end do
    c = 0.5_dp * (a + b)
  end function bisect
end module roots

37.2 Numerical integration — Simpson's rule

function simpson(f, a, b, n) result(s)
  procedure(fun) :: f
  real(dp), intent(in) :: a, b
  integer, intent(in) :: n          ! must be even
  real(dp) :: s, h
  integer :: i
  h = (b - a) / n
  s = f(a) + f(b)
  do i = 1, n - 1
    s = s + merge(4.0_dp, 2.0_dp, mod(i, 2) == 1) * f(a + i * h)
  end do
  s = s * h / 3.0_dp
end function simpson

37.3 ODE solver — classical Runge–Kutta (RK4) for systems

module rk4_mod
  use, intrinsic :: iso_fortran_env, only: dp => real64
  implicit none
  abstract interface
    pure function rhs(t, y) result(dydt)
      import :: dp
      real(dp), intent(in) :: t, y(:)
      real(dp) :: dydt(size(y))
    end function rhs
  end interface
contains
  subroutine rk4_step(f, t, y, h)
    procedure(rhs) :: f
    real(dp), intent(in) :: t, h
    real(dp), intent(inout) :: y(:)
    real(dp), dimension(size(y)) :: k1, k2, k3, k4
    k1 = f(t, y)
    k2 = f(t + h / 2, y + h / 2 * k1)
    k3 = f(t + h / 2, y + h / 2 * k2)
    k4 = f(t + h, y + h * k3)
    y = y + h / 6 * (k1 + 2 * k2 + 2 * k3 + k4)
  end subroutine rk4_step

  pure function pendulum(t, y) result(dydt)      ! y = [theta, omega]
    real(dp), intent(in) :: t, y(:)
    real(dp) :: dydt(size(y))
    dydt(1) = y(2)
    dydt(2) = -9.81_dp * sin(y(1))
  end function pendulum
end module rk4_mod

program run_pendulum
  use rk4_mod
  implicit none
  real(dp) :: y(2), t, h
  integer :: i
  y = [1.0_dp, 0.0_dp]; t = 0; h = 0.01_dp
  do i = 1, 1000
    call rk4_step(pendulum, t, y, h)
    t = t + h
    if (mod(i, 100) == 0) print '(f8.3, 2f12.6)', t, y
  end do
end program run_pendulum

37.4 Linear systems — Gaussian elimination with partial pivoting

subroutine gauss_solve(A, b, x)
  real(dp), intent(inout) :: A(:, :), b(:)
  real(dp), intent(out) :: x(:)
  integer :: n, k, i, p
  real(dp) :: factor
  n = size(b)
  do k = 1, n - 1
    p = maxloc(abs(A(k:n, k)), dim=1) + k - 1          ! pivot row
    if (p /= k) then
      A([k, p], :) = A([p, k], :)                       ! swap rows via vector subscripts
      b([k, p]) = b([p, k])
    end if
    do i = k + 1, n
      factor = A(i, k) / A(k, k)
      A(i, k:n) = A(i, k:n) - factor * A(k, k:n)
      b(i) = b(i) - factor * b(k)
    end do
  end do
  do i = n, 1, -1                                       ! back substitution
    x(i) = (b(i) - dot_product(A(i, i+1:n), x(i+1:n))) / A(i, i)
  end do
end subroutine gauss_solve

37.5 Iterative solvers & PDEs — 2-D Jacobi for Laplace's equation

program laplace
  use, intrinsic :: iso_fortran_env, only: dp => real64
  implicit none
  integer, parameter :: n = 200
  real(dp) :: u(0:n+1, 0:n+1), unew(0:n+1, 0:n+1), err
  integer :: it, i, j
  u = 0.0_dp
  u(:, n+1) = 1.0_dp                 ! hot top boundary
  unew = u
  do it = 1, 100000
    do concurrent (j = 1:n, i = 1:n)
      unew(i, j) = 0.25_dp * (u(i-1, j) + u(i+1, j) + u(i, j-1) + u(i, j+1))
    end do
    err = maxval(abs(unew(1:n, 1:n) - u(1:n, 1:n)))
    u(1:n, 1:n) = unew(1:n, 1:n)
    if (err < 1.0e-6_dp) exit
  end do
  print '(a, i0, a, es10.3)', "converged after ", it, " iterations, err = ", err
end program laplace