From 9e8fda2ec08516ef39eff937c7d1f6f4256a27df Mon Sep 17 00:00:00 2001 From: soypat Date: Sat, 3 Jun 2023 20:16:33 -0300 Subject: [PATCH] ode: refactor integrator initialization and documentation --- ode/config.go | 15 --------------- ode/dopri5.go | 36 ++++++++++++++---------------------- ode/ivp.go | 34 +++++++++++++++++++++++++++------- ode/ivp_test.go | 13 +++++++------ ode/ode.go | 10 +++++++--- ode/parameters.go | 38 ++++++++++++++++++++++++++++++++++++++ ode/rkn1210.go | 13 ++++++------- 7 files changed, 99 insertions(+), 60 deletions(-) delete mode 100644 ode/config.go create mode 100644 ode/parameters.go diff --git a/ode/config.go b/ode/config.go deleted file mode 100644 index acce7a1..0000000 --- a/ode/config.go +++ /dev/null @@ -1,15 +0,0 @@ -// 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 - -// Parameters represents the most common configuration parameters -type Parameters struct { - // Permissible tolerance given an adaptive method - AbsTolerance float64 - // Minimum/Maximum step allowed for a single iteration - MinStep, MaxStep float64 -} - -var DefaultParam = Parameters{} diff --git a/ode/dopri5.go b/ode/dopri5.go index c57ece7..11cddd4 100644 --- a/ode/dopri5.go +++ b/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 } diff --git a/ode/ivp.go b/ode/ivp.go index d06cee1..31df57a 100644 --- a/ode/ivp.go +++ b/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 diff --git a/ode/ivp_test.go b/ode/ivp_test.go index 2a7aba5..21a0567 100644 --- a/ode/ivp_test.go +++ b/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 diff --git a/ode/ode.go b/ode/ode.go index 8a27f4f..6320b59 100644 --- a/ode/ode.go +++ b/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) diff --git a/ode/parameters.go b/ode/parameters.go new file mode 100644 index 0000000..e0e15cc --- /dev/null +++ b/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") + } +} diff --git a/ode/rkn1210.go b/ode/rkn1210.go index bb425ef..2c20da4 100644 --- a/ode/rkn1210.go +++ b/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,