Degenerate Conic

Algorithms • Modern Fortran Programming • Orbital Mechanics

Jun 14, 2026

Minimization Without Derivatives

The FMIN function at Netlib is based on Richard Brent's classic localmin algorithm from [1]. The method uses a combination of golden section search and successive parabolic interpolation to compute the minimum of a 1D function without the necessity of evaluating any derivatives. There are Fortran 66 and ALGOL 60 versions of the method in the reference. I also have a modernized version here. Let's show an example of how to use it for a simple orbital mechanics application.

Problem setup

Let's solve a trivial problem in orbital mechanics: computing the time of closest approach to a planet for a conic orbit. This is known as periapsis (or, for the Earth, perigee). We don't really need to use a minimizer to compute this, but it's a problem with a known solution that we can test the minimizer on.

We need various algorithms to do this:

Given an initial elliptical orbit state, the Gooding routine pv3els will compute the time since last periapsis passage, so we have the "true" value to compare the results against. The function we pass to fmin to minimize will simply propagate the state (the independent variable is \(\Delta t\)) and return the radius magnitude (\(r\)). The point where the radius magnitude is minimized is the periapsis point we are looking for. The orbital period is used to set the upper time bound for the minimizer (we know there is only one periapsis passage per orbital period).

Example code

In our Fortran Package Manager manifest file, we define the two external dependencies we need:

[dependencies]
fmin = { git="https://github.com/jacobwilliams/fmin.git" }
fortran-astrodynamics-toolkit = { git="https://github.com/jacobwilliams/fortran-astrodynamics-toolkit.git" }

And the main program looks like this:

program main
  use fmin_module, wp => fmin_rk
  use fortran_astrodynamics_toolkit

  implicit none

  real(wp) :: p !! semiparameter [km]
  real(wp) :: period !! orbital period [sec]
  real(wp) :: dt !! time from initial state to periapsis [sec]
  real(wp),dimension(3) :: r, v !! position and velocity vectors [km] and [km/s]
  real(wp),dimension(6) :: rv0 !! initial state vector [km, km/s]
  real(wp),dimension(6) :: e0 !! Gooding orbital elements of the initial state
  integer :: n_evals = 0 !! number of function evaluations

  ! initial conditions (orbital elements):
  real(wp),parameter :: mu   = 398600.4418_wp  !! Earth grav. parameter [km^3/s^2]
  real(wp),parameter :: a    = 7000.0_wp       !! semi-major axis [km]
  real(wp),parameter :: ecc  = 0.2_wp          !! eccentricity
  real(wp),parameter :: inc  = 30.0_wp         !! inclination [deg]
  real(wp),parameter :: raan = 40.0_wp         !! right ascension of ascending node [deg]
  real(wp),parameter :: aop  = 60.0_wp         !! argument of perigee [deg]
  real(wp),parameter :: tru  = -45.0_wp        !! true anomaly [deg]
  real(wp),parameter :: tol  = 1.0e-6_wp       !! tolerance for fmin

  ! get the initial state vector and time from periapsis for the initial state:
  p = a * (1.0_wp - ecc**2)
  period = orbit_period(mu,a)
  call orbital_elements_to_rv(mu, p, ecc, inc, raan, aop, tru, r, v)
  rv0 = [r, v]
  call pv3els (mu, rv0, e0) ! e0(6) is the time from periapsis of the initial state

  ! call the minimizer:
  dt = fmin(func,0.0_wp,period,tol)

  ! print results
  write(*,'(A,F12.6,A)') 'DT to periapsis from fmin        = ', dt, ' sec'
  write(*,'(A,F12.6,A)') 'True initial time from periapsis = ', e0(6), ' sec'
  write(*,'(A,E12.3,A)') 'Error                            = ', dt + e0(6), ' sec'
  write(*,'(A,I0,A)')    'number of function evals         = ', n_evals, ' '

  contains

  function func(x) result(f)
    !! the function is to propagate from the initial state
    !! by the dt and compute the radius value
    real(wp),intent(in) :: x  !! indep. variable (dt)
    real(wp)            :: f  !! function value `f(x)` - radius magnitude
    real(wp),dimension(6) :: rvf
    call propagate(mu, rv0, x, rvf)
    f = norm2(rvf(1:3))
    n_evals = n_evals + 1
  end function func

end program main

Easy peasy!

Results

The result of running this program is:

DT to periapsis from fmin        =   652.983248 sec
True initial time from periapsis =  -652.983240 sec
Error                            =    0.875E-05 sec
number of function evals         = 12

So, it works! Using 12 function evaluations (i.e., kepler propagations of the orbit), the method has located the periapsis time to within about 8 \(\mu s\). This sort of approach can also be used for less-trivial problems where we don't have a good analytical solution, and/or derivatives are unavailable or hard to obtain.

References

  1. R. Brent, "Algorithms for Minimization Without Derivatives", Prentice-Hall, 1973.
posted at 12:00  ·   ·  Fortran  Mathematics  Optimization  ALGOL