ODE Initial Value Problems

ODE, Differential Equations

Initial-value problems for ordinary differential equations can be attacked by a great many specialized methods. Numerical analysts agree that there is no best method. Each has situations where it works best and other situations where it fails or is not very good. Also, each technique has numerous parameters, options and defaults.

The default integration method is Bulirsch-Stoer. Usually, the Bulirsch-Stoer algorithm will give better and faster results than others, but there are applications where a quality-controlled trapezoidal method or a quality-controlled 4th order Runge-Kutta method is appropriate. The algorithm used can be set by the user:

(set-ode-integration-method! 'qcrk4)
(set-ode-integration-method! 'bulirsch-stoer)
(set-ode-integration-method! 'qcctrap2)
(set-ode-integration-method! 'explicit-gear)

In SICMUtils we currently only support the Gragg-Bulirsch-Stoer algorithm. If you need any of the others, please file a PR and we can discuss how to get these in for you.

The integration methods all automatically select the step sizes to maintain the error tolerances. But if we have an exceptionally stiff system, or a bad discontinuity, for most integrators the step size will go down to zero and the integrator will make no progress. If you encounter such a disaster try explicit-gear.

We have programs that implement other methods of integration, such as an implicit version of Gear’s stiff solver, and we have a whole language for describing error control, but these features are not available through this interface.

The two main interfaces are evolve and state-advancer.

The procedure state-advancer is used to advance the state of a system according to a system of first order ordinary differential equations for a specified interval of the independent variable. The state may have arbitrary structure, however we require that the first component of the state is the independent variable.

The procedure evolve uses state-advancer to repeatedly advance the state of the system by a specified interval, examining aspects of the state as the evolution proceeds.

In the following descriptions we assume that sysder is a user provided procedure that gives the parametric system derivative. The parametric system derivative takes parameters, such as a mass or length, and produces a procedure that takes a state and returns the derivative of the state. Thus, the system derivative takes arguments in the following way:

((sysder parameter-1 ... parameter-n) state)

There may be no parameters, but then the system derivative procedure must still be called with no arguments to produce the procedure that takes states to the derivative of the state.

For example, if we have the differential equations for an ellipse centered on the origin and aligned with the coordinate axes:

Dx(t) = -a y(t)
Dy(t) = +b x(t)

We can make a parametric system derivative for this system as follows:

(defn ellipse-sysder [a b]
  (fn [[t x y]]
    (up 1				   ;; dt/dt
        (* -1 a y) ;; dx/dt
        (* b x)))) ;; dy/dt

The procedure evolve is invoked as follows:

((evolve sysder . parameters)
 initial-state dt final-t optional-opts)

The user can pass a procedure via the keyword :observe in optional-opts that takes the state as an argument.

The :observe fn is passed successive states of the system as the evolution proceeds. For example it might be used to print the state or to plot some interesting function of the state.

The interval between calls to the monitor is the argument dt. The evolution stops when the independent variable is larger than final-t. The optional keyword argument parameter :epsilon specifies the allowable error.

For example, we can evolve our state forward for 10 seconds:

((evolve ellipse-sysder 0.5 2.0)
 (up 0. .5 .5)	;; initial state
 0.01           ;; step size
 10.0)          ;; final value of t
;;=> (up 9.99999999999992 -0.2835304866702712 -0.9635568769766077)

To take more control of the integration one may use the state advancer directly.

The procedure state-advancer is invoked as follows:

((state-advancer sysder . parameters) start-state dt optional-args)

The state advancer will give a new state resulting from evolving the start state by the increment dt of the independent variable. The allowed local truncation error is specified by the optional keyword argument :epsilon:

For example,

((state-advancer ellipse-sysder 0.5 2.0)
  (up 0 0.5 0.5) 3.0 {:epsilon 1e-10})
;;=> (up 2.999999999999995 -0.530276250315008 -0.35387624023977055)

For a more complex example that shows the use of substructure in the state, consider the two-dimensional harmonic oscillator:

(defn harmonic-sysder [m k]
  (fn [state]
    (let [[x y]   (coordinate state)
          [px py] (momentum state)]
      (up 1                                ;; dt/dt
          (up (/ px m) (/ py m))           ;; dq/dt
          (down (* -1 k x) (* -1 k y)))))) ;; dp/dt

We could monitor the energy (the Hamiltonian):

(defn H [m k]
  (fn [state]
    (+ (/ (square (momentum state))
          (* 2 m))
       (* (/ 1 2) k
          (square (coordinate state))))))

(let [initial-state (up 0
                        (up 0.5 0.5)
                        (down 0.1 0.0))
      monitor (fn [_ state]
                (println
                 (state->t state)
                 "\t"
                 ((H m k) state)))
      step-size 1.0
      final-time 10
      m 0.5
      k 2.0]
  ((evolve harmonic-sysder m k)
   initial-state
   step-size
   final-time
   {:observe monitor}))

;; 0.0                  0.51
;; 0.9999999999999986 	 0.5100000000359725
;; 1.9999999999999933 	 0.5100000001529045
;; 2.9999999999999893 	 0.5100000002789965
;; 3.999999999999993 	 0.5100000004050876
;; 5.000000000000025 	 0.5100000005311794
;; 6.0000000000000515 	 0.5100000006572705
;; 7.000000000000078 	 0.510000000783362
;; 8.000000000000105 	 0.5100000009094526
;; 9.0000000000001 	   0.5100000010355437
;; 10.000000000000103 	 0.510000001191079

Constants

There are a few constants that we find useful, and are thus provided in SICMUtils.

(def pi (* 4 (atan 1 1)))
(def -pi (- pi))

For numerical analysis, we provide the smallest number that when added to 1.0 makes a difference:

(def machine-epsilon
  (loop [e 1.0]
    (if (= 1.0 (+ e 1.0))
      (* e 2.0)
      (recur (/ e 2.0)))))

(def sqrt-machine-epsilon
  (Math/sqrt machine-epsilon))