ode: refactor integrator initialization and documentation

ode/dopri5.go

@@ -27,9 +27,8 @@ type DoPri5 struct {
 // Set implements the Integrator interface. Initializes a Dormand Prince method
 // for use as a solver
 func (dp *DoPri5) Init(ivp IVP) {
-
+	ivp.mustValidate()
 	dp.fx = ivp.Func
-
 	t0, x0 := ivp.T0, mat.VecDenseCopyOf(ivp.Y0)
 	dp.dom = t0
 	nx := x0.Len()
@@ -52,8 +51,9 @@ func (dp *DoPri5) State(s *State) {
 }
 
 // Step implements Integrator interface. Advances solution by step h. If algorithm
-// is set to adaptive then h is just a suggestion
-func (dp *DoPri5) Step(h float64) (float64, error) {
+// is set to adaptive then h is just a suggestion. The returned stepTaken is the
+// actual step taken by the algorithm.
+func (dp *DoPri5) Step(h float64) (stepTaken float64, _ error) {
 	const c20, c21 = 1. / 5., 1. / 5.
 	const c30, c31, c32 = 3. / 10., 3. / 40., 9. / 40.
 	const c40, c41, c42, c43 = 4. / 5., 44. / 45., -56. / 15., 32. / 9.
@@ -151,26 +151,18 @@ SOLVE:
 	return h, nil
 }
 
-// NewDormandPrince5 returns a adaptive-ready Dormand Prince solver of order 5
-//
+// NewDormandPrince5 returns a adaptive-ready Dormand Prince solver of order 5.
 // To enable step size adaptivity minimum step size must be set and
-// absolute tolerance must be set. i.e:
-//
-//  NewDormandPrince5(ConfigScalarTolerance(0, 0.1), ConfigStepLimits(1, 1e-3))
+// absolute tolerance must be set. A zero value of Parameters is a valid
+// configuration.
 //
 // If a invalid configuration is passed the function panics.
-func NewDormandPrince5(cfg Parameters) *DoPri5 {
-	dp := new(DoPri5)
-	if cfg.AbsTolerance == 0 {
-		cfg.AbsTolerance = 1e-3
+func NewDormandPrince5(params Parameters) *DoPri5 {
+	params.mustValidate()
+	return &DoPri5{
+		minStep:  params.MinStep,
+		maxStep:  params.MaxStep,
+		atol:     params.AbsTolerance,
+		adaptive: params.AbsTolerance > 0,
 	}
-	if cfg.MinStep == 0 {
-		cfg.MinStep = 1e-4
-	}
-	if cfg.MaxStep == 0 {
-		cfg.MaxStep = 1e3
-	}
-	dp.atol = cfg.AbsTolerance
-
-	return dp
 }

ode/ivp.go

@@ -14,8 +14,9 @@ import (
 // IVP defines a multivariable, initial value problem represented by a system of ordinary differential equations.
 //
 // These problems have the form
-//  y'(t) = f(t, y(t))
-//  y(t_0) = y_0
+//
+//	y'(t) = f(t, y(t))
+//	y(t_0) = y_0
 //
 // Where:
 // t is a scalar representing the integration domain, which is time for most physical problems.
@@ -39,10 +40,19 @@ type IVP struct {
 	Func func(dst *mat.VecDense, t float64, y mat.Vector)
 }
 
+func (ivp IVP) mustValidate() {
+	if ivp.Y0 == nil || ivp.Func == nil {
+		panic("ode: nil value in IVP")
+	}
+	if math.IsNaN(ivp.T0) || math.IsInf(ivp.T0, 0) {
+		panic("ode: infinite or NaN initial domain value")
+	}
+}
+
 // NewModel returns a IVP given initial conditions (x0,u0), differential equations (xeq) and
 // input functions for non-autonomous ODEs (ueq).
 func NewIVP(t0 float64, y0 mat.Vector, f func(y *mat.VecDense, dom float64, x mat.Vector)) (IVP, error) {
-	if y0 == nil || math.IsNaN(t0) || f == nil {
+	if y0 == nil || math.IsNaN(t0) || math.IsInf(t0, 0) || f == nil {
 		return IVP{}, errors.New("bad model value")
 	}
 	return IVP{Func: f, Y0: y0, T0: t0}, nil
@@ -52,9 +62,10 @@ func NewIVP(t0 float64, y0 mat.Vector, f func(y *mat.VecDense, dom float64, x ma
 // second order system of ordinary differential equations.
 //
 // These problems have the form
-//  y''(t)   = f(t, y(t))
-//  y(t_0)   = y_0
-//  y'(t_0)  = y'_0
+//
+//	y''(t)   = f(t, y(t))
+//	y(t_0)   = y_0
+//	y'(t_0)  = y'_0
 type IVP2 struct {
 	Y0  mat.Vector
 	DY0 mat.Vector
@@ -64,9 +75,18 @@ type IVP2 struct {
 	Func func(dst *mat.VecDense, t float64, y mat.Vector)
 }
 
+func (ivp IVP2) mustValidate() {
+	if ivp.Y0 == nil || ivp.DY0 == nil || ivp.Func == nil {
+		panic("ode: nil value in IVP2")
+	}
+	if math.IsNaN(ivp.T0) || math.IsInf(ivp.T0, 0) {
+		panic("ode: infinite or NaN initial domain value")
+	}
+}
+
 // NewIVP2 returns a second order initial value problem.
 func NewIVP2(t0 float64, y0, dy0 mat.Vector, f func(y *mat.VecDense, dom float64, x mat.Vector)) (IVP2, error) {
-	if y0 == nil || dy0 == nil || math.IsNaN(t0) || f == nil || y0.Len() != dy0.Len() {
+	if y0 == nil || dy0 == nil || math.IsNaN(t0) || math.IsInf(t0, 0) || f == nil || y0.Len() != dy0.Len() {
 		return IVP2{}, errors.New("bad model value")
 	}
 	return IVP2{Func: f, Y0: y0, DY0: dy0, T0: t0}, nil

ode/ivp_test.go

@@ -24,7 +24,7 @@ func TestSolve(t *testing.T) {
 	// domain start.
 	quad := quadTestModel(t)
 	ivp := quad.IVP
-	solver := ode.NewDormandPrince5(ode.DefaultParam)
+	solver := ode.NewDormandPrince5(ode.Parameters{})
 	solver.Init(ivp)
 	stepsize := 0.5 / 15.
 	end := 0.5
@@ -78,7 +78,7 @@ func Example_solve() {
 		log.Fatal(err)
 	}
 	// Here we choose our algorithm. Runge-Kutta 4th order is used
-	var solver ode.Integrator = ode.NewDormandPrince5(ode.DefaultParam)
+	var solver ode.Integrator = ode.NewDormandPrince5(ode.Parameters{})
 	// Solve function makes it easy to integrate a problem without having
 	// to implement the `for` loop. This example integrates the IVP with a step size
 	// of 0.1 over a domain of 10. arbitrary units, in this case, 10 seconds.
@@ -113,9 +113,10 @@ func quadTestModel(t *testing.T) *TestModel {
 }
 
 // exponential unidimensional model may be used for future algorithms
-//  y'(t) = -15*y(t)
-//  y(t=0) = 1
-//  solution: y(t) = exp(-15*t)
+//
+//	y'(t) = -15*y(t)
+//	y(t=0) = 1
+//	solution: y(t) = exp(-15*t)
 func exp1DTestModel(t *testing.T) *TestModel {
 	tau := -2.
 	t0 := 0.0
@@ -158,7 +159,7 @@ func TestRK1210_fallingball(t *testing.T) {
 			dst.SetVec(0, gravity)
 		},
 	}
-	integrator := ode.NewRKN1210(ode.DefaultParam)
+	integrator := ode.NewRKN1210(ode.Parameters{})
 	integrator.Init(ivp)
 	state := ode.State2{Y: mat.NewVecDense(1, nil), DY: mat.NewVecDense(1, nil)}
 	var time, position, velocity float64 = 0, initPos, initVel

ode/ode.go

@@ -15,11 +15,15 @@ import (
 // system of ordinary differential equations (ODEs).
 type Integrator interface {
 	// Init initializes the integrator and sets the initial condition.
+	// It should be called before the first call to Step or State.
 	Init(IVP)
 
-	// Step advances the current state by taking at most the given step. It returns a proposed step size
-	// for the next step and an error indicating whether the step was successful.
-	Step(step float64) (stepNext float64, err error)
+	// Step advances the current state by taking at most the argument step. The step taken
+	// will depend on whether the integrator implementation is adaptive or not.
+	// Non-adaptive methods will always take the argument step.
+	// If an error is returned by Step the integrator should be considered invalidated
+	// and not be used again.
+	Step(step float64) (stepTaken float64, err error)
 
 	// State stores the current state of the integrator in-place in dst.
 	State(dst *State)

ode/parameters.go

@@ -0,0 +1,38 @@
+// Copyright ©2022 The Gonum Authors. All rights reserved.
+// Use of this source code is governed by a BSD-style
+// license that can be found in the LICENSE file.
+
+package ode
+
+import "math"
+
+// Parameters is a generic configuration object for adaptive methods.
+type Parameters struct {
+	// Permissible error tolerance given an adaptive method.
+	AbsTolerance float64
+	// Minimum/Maximum step allowed for a single iteration.
+	MinStep, MaxStep float64
+}
+
+func (p Parameters) mustValidate() {
+	if math.IsInf(p.MaxStep, 0) || math.IsNaN(p.MaxStep) {
+		panic("ode: max step must be finite and not NaN")
+	}
+	if math.IsInf(p.MinStep, 0) || math.IsNaN(p.MinStep) {
+		panic("ode: min step must be finite and not NaN")
+	}
+	if math.IsInf(p.AbsTolerance, 0) || math.IsNaN(p.AbsTolerance) {
+		panic("ode: tolerance must be finite and not NaN")
+	}
+	if p.AbsTolerance < 0 || p.MinStep < 0 || p.MaxStep < 0 {
+		panic("ode: negative step parameters or tolerance")
+	}
+
+	isadaptive := p.AbsTolerance > 0
+	if isadaptive && (p.MaxStep == 0 || p.MinStep == 0) {
+		panic("ode: max and min step must be set when using adaptive methods")
+	}
+	if isadaptive && p.MinStep >= p.MaxStep {
+		panic("ode: min step must be greater than max step")
+	}
+}

ode/rkn1210.go

@@ -33,6 +33,7 @@ const rk1210Len = 17
 // Init initializes the Runge-Kutta-Nyström integrator with a second order
 // ODE initial value problem.
 func (rk *RKN1210) Init(ivp IVP2) {
+	ivp.mustValidate()
 	rk.fx = ivp.Func
 	rk.dom, rk.y = ivp.T0, mat.VecDenseCopyOf(ivp.Y0)
 	rk.dy = mat.VecDenseCopyOf(ivp.DY0)
@@ -66,8 +67,8 @@ func (rk *RKN1210) SetState(s State2) {
 
 // Step implements Integrator interface. Advances solution by step h. If algorithm
 // is set to adaptive then h is just a suggestion.
-func (rk *RKN1210) Step(h float64) (step float64, err error) {
-	step = h
+func (rk *RKN1210) Step(h float64) (stepTaken float64, err error) {
+	stepTaken = h
 	adaptive := rk.atol > 0
 
 	for {
@@ -114,7 +115,7 @@ func (rk *RKN1210) Step(h float64) (step float64, err error) {
 				continue
 			}
 			// The error is within tolerance and we may suggest the user use a larger step.
-			step = hnew
+			stepTaken = hnew
 		}
 		break
 	}
@@ -125,14 +126,12 @@ func (rk *RKN1210) Step(h float64) (step float64, err error) {
 	rk.y.AddScaledVec(rk.y, h, rk.aux)
 	rk.dy.AddVec(rk.dy, rk.hFDbhat)
 	rk.dom += h
-	return step, nil
+	return stepTaken, nil
 }
 
 // NewRKN1210 configures a new RKN1210 instance.
 func NewRKN1210(cfg Parameters) *RKN1210 {
-	if (cfg.AbsTolerance != 0 && cfg.MaxStep <= 0) || cfg.MaxStep < cfg.MinStep || cfg.MinStep < 0 {
-		panic("invalid parameter supplied")
-	}
+	cfg.mustValidate()
 	return &RKN1210{
 		atol:    cfg.AbsTolerance,
 		minStep: cfg.MinStep,