From 6071bd986b0b7c53e771ee866172dd7e2fd679ca Mon Sep 17 00:00:00 2001 From: soypat Date: Fri, 21 May 2021 16:13:39 -0300 Subject: [PATCH 01/10] add ivp package to solve multivariable, non-autonomous ODEs --- ivp/README.md | 3 + ivp/doc.go | 7 +++ ivp/ivp.go | 108 +++++++++++++++++++++++++++++++++ ivp/ivp_test.go | 76 +++++++++++++++++++++++ ivp/model.go | 46 ++++++++++++++ ivp/rungekutta4.go | 74 +++++++++++++++++++++++ ivp/rungekutta4_test.go | 129 ++++++++++++++++++++++++++++++++++++++++ 7 files changed, 443 insertions(+) create mode 100644 ivp/README.md create mode 100644 ivp/doc.go create mode 100644 ivp/ivp.go create mode 100644 ivp/ivp_test.go create mode 100644 ivp/model.go create mode 100644 ivp/rungekutta4.go create mode 100644 ivp/rungekutta4_test.go diff --git a/ivp/README.md b/ivp/README.md new file mode 100644 index 0000000..1c501d0 --- /dev/null +++ b/ivp/README.md @@ -0,0 +1,3 @@ +# Gonum ivp [![GoDoc](https://godoc.org/gonum.org/v1/gonum/ivp?status.svg)](https://godoc.org/gonum.org/v1/gonum/ivp) + +Package ivp provides numerical methods solving multivariable ODE initial value problems for the Go programming language. diff --git a/ivp/doc.go b/ivp/doc.go new file mode 100644 index 0000000..a336e44 --- /dev/null +++ b/ivp/doc.go @@ -0,0 +1,7 @@ +// Copyright ©2021 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 ivp provides interfaces and structures to solve +// multivariable ODE initial value problems. +package ivp // import "gonum.org/v1/gonum/ivp" diff --git a/ivp/ivp.go b/ivp/ivp.go new file mode 100644 index 0000000..3da320c --- /dev/null +++ b/ivp/ivp.go @@ -0,0 +1,108 @@ +// Copyright ©2021 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 ivp + +import ( + "errors" + "fmt" + + "gonum.org/v1/gonum/mat" +) + +// IVP defines a multivariable, non-autonomous initial value problem. +// It is worth mentioning system does not necessarily need +// to be non-autonomous. https://en.wikipedia.org/wiki/Autonomous_system_(mathematics) +// +// +// These problems have the form (capital letters are vectors in this example) +// X'(t) = F(t, X, U) +// X(0) = X_0 +// +// Where X' is the vector of first derivatives of the state vector X. t is +// a scalar representing the integrations domain, which is usually time +// for most physical problems. U is a vector which is a function of the current +// state. Put simply: +// U = F_u(t, X, U) +// +// An initial value problem is characterized by boundary conditions imposed +// on the state vector X at the beginning of the integration domain. These +// boundary conditions are returned by the IV() method for the state vector +// as x0 and for the input vector as u0. +// +// The term "state vector" and "state variables" are used interchangeably +// throughout the code and refer to X vector of independent variables. +type IVP interface { + // Initial values vector for state variables x and inputs u. x0 defines + // the first values the state vector takes when integration begins. + IV() (x0, u0 mat.Vector) + // Equations returns the coupled, non-linear algebraic differential + // equations (xequations) for the state variables (x) and the functions for inputs (u). + // The input functions (ufunc) are not differential equations but rather + // calculated directly from a given x vector and current input vector. + // Results are stored in y which are the length of x and u, respectively. + // The scalar (float64) argument is the domain over which the + // problem is integrated, which is usually time for most physical problems. + // + // If problem has no input functions then u supplied and ufunc returned + // may be nil. x equations my not be nil. + Equations() (xequations func(y []float64, t float64, x, u []float64), ufunc func(u_next []float64, t float64, x, u []float64)) + // Dimensions of x state vector and u inputs. + Dims() (nx, nu int) +} + +// Integrator abstracts algorithm specifics. For anyone looking to +// implement it, Set(ivp) should be called first to initialize the IVP with +// initial values. Step will calculate the next x values and store them in y +// u values should not be stored as they can easily be obtained if one has +// x values. Integrator should store 1 or more (depending on algorithm used) +// of previously calculated x values to be able to integrate. +type Integrator interface { + // Set initializes an initial value problem. First argument + // is the initial domain integration point, is usually zero. + Set(float64, IVP) error + // Dimensions of x state vector and u inputs. Set must be called first. + Dims() (nx, nu int) + // Step integrates IVP and stores result in y. step is a suggested step + // for the algorithm to take. The algorithm may decide that it is not sufficiently + // small or big enough (these are adaptive algorithms) and take a different step. + // The resulting step is returned as the first parameter + Step(y []float64, step float64) (float64, error) +} + +type result = struct { + DomainStartOffset float64 + X []float64 +} + +// Solve solves an already initialized Integrator returning state vector results. +// Returns error upon needing to allocate over 2GB of memory +func Solve(solver Integrator, stepsize, domainLength float64) (results []result, err error) { + const maxAllocGB = 2 + integrated := 0. + expectedLength := int(domainLength/stepsize) + 1 + nx, _ := solver.Dims() + if nx == 0 { + return nil, errors.New("state vector length can not be equal to zero. Has ivp been set?") + } + size := 8 * (nx + 1) * expectedLength + if size > maxAllocGB*1e9 { + return nil, fmt.Errorf("solution exceeds %dGB or not initialized (size is %dMB)", maxAllocGB, size/1e6) + } + results = make([]result, 0, expectedLength) + + for integrated < domainLength { + res := make([]float64, nx) + stepsize, err = solver.Step(res, stepsize) + if err != nil { + return results, err + } + if stepsize <= 0 { + return results, errors.New("got zero or negative step size from Integrator") + } + integrated += stepsize + results = append(results, result{DomainStartOffset: integrated, X: res}) + } + return results, nil +} diff --git a/ivp/ivp_test.go b/ivp/ivp_test.go new file mode 100644 index 0000000..3a7be4f --- /dev/null +++ b/ivp/ivp_test.go @@ -0,0 +1,76 @@ +// Copyright ©2021 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 ivp_test + +import ( + "fmt" + "log" + "math" + "testing" + + "gonum.org/v1/gonum/ivp" + "gonum.org/v1/gonum/mat" +) + +func TestSolve(t *testing.T) { + quad := quadTestModel(t) + solver := new(ivp.RK4) + solver.Set(0.0, quad) + stepsize := 0.1 + results, err := ivp.Solve(solver, stepsize, 10) + quadsol := quad.solution + sol, _ := quadsol.Equations() + if err != nil { + t.Fatal(err) + } + nxs, _ := quadsol.Dims() + solresults := make([]float64, nxs) + for i := range results { + sol(solresults, results[i].DomainStartOffset, results[i].X, nil) + for j := range solresults { + got := math.Abs(solresults[j] - results[i].X[j]) + expect := quad.err(stepsize, float64(i)) + if true || got > expect { + t.Errorf("error %g is greater than permitted tolerance %g", got, expect) + } + } + } +} + +func Example_solve() { + const ( + g = -10. // gravity field [m.s^-2] + ) + // we declare our physical model in the following function + ballModel, err := ivp.NewModel(mat.NewVecDense(2, []float64{100., 0.}), + nil, func(yvec []float64, _ float64, xvec, _ []float64) { + // this anonymous function defines the physics. + // The first variable xvec[0] corresponds to position + // second variable xvec[1] is velocity. + Dx := xvec[1] + // yvec represents change in xvec, or derivative with respect to domain + // Change in position will be equal to velocity, which is the second variable: + // thus yvec[0] = xvec[1], which is the same as saying "change in xvec[0]" is equal to xvec[1] + yvec[0] = Dx + // change in velocity is acceleration. We suppose our ball is on earth accelerating at `g` + yvec[1] = g + }, nil) + if err != nil { + log.Fatal(err) + } + // Here we choose our algorithm. Runge-Kutta 4th order is used + var solver ivp.Integrator = new(ivp.RK4) + // Before integrating the IVP is passed to the integrator (a.k.a solver). Domain (time) starts at 0 + dom := 0. + err = solver.Set(dom, ballModel) + if err != nil { + log.Fatal(err) + } + // 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. + results, err := ivp.Solve(solver, 0.1, 10.) + fmt.Println(results) +} diff --git a/ivp/model.go b/ivp/model.go new file mode 100644 index 0000000..46491ba --- /dev/null +++ b/ivp/model.go @@ -0,0 +1,46 @@ +// Copyright ©2021 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 ivp + +import "gonum.org/v1/gonum/mat" + +// Model implements IVP interface. X vector and equations can not be nil or zero length. +type Model struct { + x0, u0 mat.Vector + xeq, ins func(y []float64, dom float64, x, u []float64) +} + +// IV returns initial values of the IVP. First returned parameter is the +// starting x vector and second parameter are inputs when solving non-autonomous +// ODEs. +func (m *Model) IV() (mat.Vector, mat.Vector) { return m.x0, m.u0 } + +// Equations returns differential equations relating to state vector x and input functions +// for non-autonomous ODEs. +// +// Input functions may be nil (ueq). +func (m *Model) Equations() (xeq, ueq func(y []float64, dom float64, x, u []float64)) { + return m.xeq, m.ins +} + +// Dims returns dimension of state and input (x length and u length, respectively). +// Dims has some dimension checking involved and can serve as a preliminary system verifier. +// +// Dims panics if x vector is nil. +func (m *Model) Dims() (nx, nu int) { + if m.x0 == nil { + panic("x vector can not be nil") + } + if m.u0 == nil { + nu = 0 + } + return m.x0.Len(), nu +} + +// NewModel returns a IVP given initial conditions (x0,u0), differential equations (xeq) and +// input functions for non-autonomous ODEs (ueq). +func NewModel(x0, u0 mat.Vector, xeq, ueq func(y []float64, dom float64, x, u []float64)) (*Model, error) { + return &Model{xeq: xeq, ins: ueq, x0: x0, u0: u0}, nil +} diff --git a/ivp/rungekutta4.go b/ivp/rungekutta4.go new file mode 100644 index 0000000..6af223d --- /dev/null +++ b/ivp/rungekutta4.go @@ -0,0 +1,74 @@ +// Copyright ©2021 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 ivp + +import ( + "errors" + + "gonum.org/v1/gonum/floats" +) + +// RK4 implements Integrator interface for Runke-Kutta +// 4th order multivariable method (non adaptive) +type RK4 struct { + dom float64 + x, u []float64 + a, b, c, d []float64 + fx, fu func(y []float64, t float64, x, u []float64) +} + +// Step implements Integrator interface +func (rk *RK4) Step(y []float64, h float64) (float64, error) { + const overSix = 1. / 6. + t := rk.dom + // set a, b, c, d (in some literatures these are k1,k2,k3,k4) + rk.fx(rk.a, t, rk.x, rk.u) + rk.fx(rk.b, t+0.5*h, floats.AddScaledTo(rk.b, rk.x, 0.5*h, rk.a), rk.u) + rk.fx(rk.c, t+0.5*h, floats.AddScaledTo(rk.c, rk.x, 0.5*h, rk.b), rk.u) + rk.fx(rk.d, t+h, floats.AddScaledTo(rk.d, rk.x, h, rk.c), rk.u) + + floats.Add(rk.a, rk.d) + floats.Add(rk.b, rk.c) + floats.AddScaled(rk.a, 2, rk.b) + floats.AddScaledTo(y, rk.x, h*overSix, rk.a) + // finished integrating. Now we update RK4 structure for future Steps + copy(rk.x, y) // store results + rk.dom += h + if len(rk.u) > 0 { + rk.fu(rk.u, t, rk.x, rk.u) + } + return h, nil +} + +// Return size of state vector and input vector. Solver must be set beforehand. +func (rk *RK4) Dims() (nx, nu int) { return len(rk.x), len(rk.u) } + +// Set implements integrator interface. All RK4 data +// is reset when calling Set. +func (rk *RK4) Set(d0 float64, ivp IVP) error { + if ivp == nil { + return errors.New("IVP is nil") + } + nx, nu := ivp.Dims() + rk.dom = d0 //set start domain + rk.a = make([]float64, nx) + rk.b = make([]float64, nx) + rk.c = make([]float64, nx) + rk.d = make([]float64, nx) + rk.x = make([]float64, nx) + rk.u = make([]float64, nu) + // set initial values + x0, u0 := ivp.IV() + for i := 0; i < x0.Len(); i++ { + rk.x[i] = x0.AtVec(i) + } + if nu > 0 { + for i := 0; i < u0.Len(); i++ { + rk.u[i] = u0.AtVec(i) + } + } + rk.fx, rk.fu = ivp.Equations() + return nil +} diff --git a/ivp/rungekutta4_test.go b/ivp/rungekutta4_test.go new file mode 100644 index 0000000..e4d7192 --- /dev/null +++ b/ivp/rungekutta4_test.go @@ -0,0 +1,129 @@ +// Copyright ©2021 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 ivp_test + +import ( + "fmt" + "log" + "math" + "testing" + + "gonum.org/v1/gonum/ivp" + "gonum.org/v1/gonum/mat" +) + +type TestModel struct { + *ivp.Model + solution *ivp.Model + err func(h, i float64) float64 +} + +func TestQuadratic(t *testing.T) { + Quadratic := quadTestModel(t) + var solver ivp.Integrator = new(ivp.RK4) + + err := solver.Set(0, Quadratic.Model) + if err != nil { + log.Fatal(err) + } + + nx, _ := Quadratic.Model.Dims() + steps := 10 + dt := 0.1 + + results := make([]float64, nx) + + solmodel := Quadratic.solution + soleq, _ := solmodel.Equations() + solDims, _ := solmodel.Dims() + solution := make([]float64, solDims) + for i := 1.; i < float64(steps+1); i++ { + dom := i * dt + solver.Step(results, dt) + soleq(solution, dom, results, nil) + for j := range results { + got := math.Abs(solution[j] - results[j]) + expected := Quadratic.err(dt, i) + if got > expected { + t.Errorf("error %e greater than threshold %e", got, expected) + } + + } + } +} + +func Example_fallingBall() { + const ( + g = -10. // gravity field [m.s^-2] + ) + // we declare our physical model in the following function + ballModel, err := ivp.NewModel(mat.NewVecDense(2, []float64{100., 0.}), + nil, func(yvec []float64, _ float64, xvec, _ []float64) { + // this anonymous function defines the physics. + // The first variable xvec[0] corresponds to position + // second variable xvec[1] is velocity. + Dx := xvec[1] + // yvec represents change in xvec, or derivative with respect to domain + // Change in position will be equal to velocity, which is the second variable: + // thus yvec[0] = xvec[1], which is the same as saying "change in xvec[0]" is equal to xvec[1] + yvec[0] = Dx + // change in velocity is acceleration. We suppose our ball is on earth accelerating at `g` + yvec[1] = g + }, nil) + if err != nil { + log.Fatal(err) + } + // Here we choose our algorithm. Runge-Kutta 4th order is used + var solver ivp.Integrator = new(ivp.RK4) + // Before integrating the IVP is passed to the integrator (a.k.a solver). Domain (time) starts at 0 + dom := 0. + err = solver.Set(dom, ballModel) + if err != nil { + log.Fatal(err) + } + + // we define the domain over which we integrate: 10 steps with step length of 0.1 + steps := 10 // number of steps + dt := 0.1 // step length + // we will store position in xvec and domain (time) in dvec + xvec := mat.NewVecDense(steps, nil) + dvec := mat.NewVecDense(steps, nil) + // results is an auxiliary vector that stores integration results + nx, _ := ballModel.Dims() + results := make([]float64, nx) + + for i := 0; i < steps; i++ { + // Step integrates the IVP. Each step advances the solution + step, _ := solver.Step(results, dt) // for non-adaptive algorithms step == dt + dom += step // calculate domain at current step + xvec.SetVec(i, results[0]) // set x value + dvec.SetVec(i, dom) + } + // print results + fmt.Println(mat.Formatted(xvec), "\n\n", mat.Formatted(dvec)) +} + +// Quadratic model may be used for future algorithms +func quadTestModel(t *testing.T) *TestModel { + Quadratic := new(TestModel) + quad, err := ivp.NewModel(mat.NewVecDense(2, []float64{0, 0}), + nil, func(y []float64, dom float64, x, u []float64) { + y[0], y[1] = x[1], 1. + }, nil) + if err != nil { + t.Fatal(err) + } + Quadratic.Model = quad + quadsol, err := ivp.NewModel(mat.NewVecDense(2, []float64{0, 0}), + nil, func(y []float64, dom float64, x, u []float64) { + y[0], y[1] = dom*dom/2., dom + }, nil) + if err != nil { + t.Fatal(err) + } + Quadratic.solution = quadsol + Quadratic.err = func(h, i float64) float64 { return math.Pow(h*i, 4) } + return Quadratic +} From e47b2cfd3b118ee5bd63f101f6d2ff43da0fc591 Mon Sep 17 00:00:00 2001 From: soypat Date: Fri, 21 May 2021 18:21:36 -0300 Subject: [PATCH 02/10] add draft DormandPrince solver --- ivp/dopri5.go | 144 ++++++++++++++++++++++++++++++++++++++++++++++++ ivp/ivp.go | 9 ++- ivp/ivp_test.go | 2 +- 3 files changed, 151 insertions(+), 4 deletions(-) create mode 100644 ivp/dopri5.go diff --git a/ivp/dopri5.go b/ivp/dopri5.go new file mode 100644 index 0000000..5887987 --- /dev/null +++ b/ivp/dopri5.go @@ -0,0 +1,144 @@ +package ivp + +import ( + "math" + + "gonum.org/v1/gonum/floats" +) + +type DoPri5 struct { + maxError, minStep, maxStep float64 + adaptive bool + x, u []float64 + // solutions + y5, err45 []float64 + // integration coefficients (for Butcher Tableau) + k1, k2, k3, k4, k5, k6, k7 []float64 + dom float64 + fx, fu func(y []float64, t float64, x, u []float64) +} + +// Set implements the Integrator interface. Initializes a Dormand Prince method +// for use as a solver +func (dp *DoPri5) Set(t0 float64, ivp IVP) error { + dp.dom = t0 + xequations, ufunc := ivp.Equations() + dp.fx = xequations + dp.fu = ufunc + nx, nu := ivp.Dims() + dp.k1, dp.k2, dp.k3 = make([]float64, nx), make([]float64, nx), make([]float64, nx) + dp.k4, dp.k5, dp.k6 = make([]float64, nx), make([]float64, nx), make([]float64, nx) + dp.k7 = make([]float64, nx) + dp.x = make([]float64, nx) + dp.y5 = make([]float64, nx) + dp.err45 = make([]float64, nx) + dp.u = make([]float64, nu) + return nil +} + +// 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(y4 []float64, h float64) (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 + const c50, c51, c52, c53, c54 = 8. / 9., 19372. / 6561., -25360. / 2187., 64448. / 6561., -212. / 729. + const c60, c61, c62, c63, c64, c65 = 1., 9017. / 3168., -355. / 33., 46732. / 5247., 49. / 176., -5103. / 18656. + const c70, c71, c72, c73, c74, c75, c76 = 1., 35. / 384., 0., 500. / 1113., 125. / 192., -2187. / 6784., 11. / 84. + // Alternate solution for error calculation + const a1, a3, a4, a5, a6, a7 = 5179. / 57600., 7571. / 16695., 393. / 640., -92097. / 339200., 187. / 2100., 1. / 40. + // Fourth order solution is used to advance + const b1, b3, b4, b5, b6 = 35. / 384., 500. / 1113., 125. / 192., -2187. / 6784., 11. / 84. + + // prettier variable names, also someone once said pointers are equivalent to JMP + x, u := dp.x, dp.u + F := dp.fx + t := dp.dom + // if adaptive then this tag will be used until h==DoPri5.MinStep +SOLVE: + F(dp.k1, t, x, u) + floats.Scale(h, dp.k1) + + floats.AddScaledTo(dp.k2, x, c21, dp.k1) + F(dp.k2, t+c20*h, x, u) + floats.Scale(h, dp.k2) + + floats.AddScaledTo(dp.k3, x, c31, dp.k1) + floats.AddScaled(dp.k3, c32, dp.k2) + F(dp.k3, t+c30*h, x, u) + floats.Scale(h, dp.k3) + + floats.AddScaledTo(dp.k4, x, c41, dp.k1) + floats.AddScaled(dp.k4, c42, dp.k2) + floats.AddScaled(dp.k4, c43, dp.k3) + F(dp.k4, t+c40*h, x, u) + floats.Scale(h, dp.k4) + + floats.AddScaledTo(dp.k5, x, c51, dp.k1) + floats.AddScaled(dp.k5, c52, dp.k2) + floats.AddScaled(dp.k5, c53, dp.k3) + floats.AddScaled(dp.k5, c54, dp.k4) + F(dp.k5, t+c50*h, x, u) + floats.Scale(h, dp.k5) + + floats.AddScaledTo(dp.k6, x, c61, dp.k1) + floats.AddScaled(dp.k6, c62, dp.k2) + floats.AddScaled(dp.k6, c63, dp.k3) + floats.AddScaled(dp.k6, c64, dp.k4) + floats.AddScaled(dp.k6, c65, dp.k5) + F(dp.k6, t+c60*h, x, u) + floats.Scale(h, dp.k6) + + // fourth order approximation used to advance solution + floats.AddScaledTo(y4, x, b1, dp.k1) + floats.AddScaled(y4, b3, dp.k3) + floats.AddScaled(y4, b4, dp.k4) + floats.AddScaled(y4, b5, dp.k5) + floats.AddScaled(y4, b6, dp.k6) + if dp.adaptive { + // calculation of order 5 coefficients + floats.AddScaledTo(dp.k7, x, c71, dp.k1) + floats.AddScaled(dp.k7, c72, dp.k2) + floats.AddScaled(dp.k7, c73, dp.k3) + floats.AddScaled(dp.k7, c74, dp.k4) + floats.AddScaled(dp.k7, c75, dp.k5) + floats.AddScaled(dp.k7, c76, dp.k6) + F(dp.k7, t+h*c70, x, u) + floats.Scale(h, dp.k7) + + floats.AddScaledTo(dp.y5, x, a1, dp.k1) + floats.AddScaled(dp.y5, a3, dp.k3) + floats.AddScaled(dp.y5, a4, dp.k4) + floats.AddScaled(dp.y5, a5, dp.k5) + floats.AddScaled(dp.y5, a6, dp.k6) + floats.AddScaled(dp.y5, a7, dp.k7) + + // compute error between fifth order solution and fourth order solution + errRatio := dp.maxError / floats.Norm(floats.SubTo(dp.err45, y4, dp.y5), math.Inf(1)) + hnew := math.Min(math.Max(0.9*h*math.Pow(errRatio, .2), dp.minStep), dp.maxStep) + // given "bad" error ratio, algorithm will recompute with a smaller step, given it is not the minimum permissible step or less + if errRatio < 1 && h > dp.minStep { + h = hnew + goto SOLVE + } + } + // advance solution with fourth order solution + copy(dp.x, y4) + if dp.fu != nil { + dp.fu(u, t+h, dp.x, u) + } + return h, nil +} + +// NewDormandPrince5 returns a adaptive-ready Dormand Prince solver of order 5 +func NewDormandPrince5(maxError, maxStep, minStep float64) *DoPri5 { + dp := new(DoPri5) + if maxError > 0 && maxStep > minStep && minStep > 0 { + dp.adaptive = true + dp.maxError = maxError + dp.minStep, dp.maxStep = minStep, maxStep + } else { + panic("NewDormandPrince5 generator requires proper arguments: maxError > 0 && maxStep > minStep && minStep > 0") + } + return dp +} diff --git a/ivp/ivp.go b/ivp/ivp.go index 3da320c..4194c22 100644 --- a/ivp/ivp.go +++ b/ivp/ivp.go @@ -20,12 +20,15 @@ import ( // X'(t) = F(t, X, U) // X(0) = X_0 // -// Where X' is the vector of first derivatives of the state vector X. t is +// Where X' is the vector of first derivatives of the state vector X. +// F would be xequations as returned by Equations(). t is // a scalar representing the integrations domain, which is usually time // for most physical problems. U is a vector which is a function of the current -// state. Put simply: -// U = F_u(t, X, U) +// state. Put simply, the next input is a function of all current state variables +// and, possibly, current input as well. +// U_next = F_u(t, X, U) // +// Where F_u is ufunc as returned by Equations() // An initial value problem is characterized by boundary conditions imposed // on the state vector X at the beginning of the integration domain. These // boundary conditions are returned by the IV() method for the state vector diff --git a/ivp/ivp_test.go b/ivp/ivp_test.go index 3a7be4f..e96b1b7 100644 --- a/ivp/ivp_test.go +++ b/ivp/ivp_test.go @@ -32,7 +32,7 @@ func TestSolve(t *testing.T) { for j := range solresults { got := math.Abs(solresults[j] - results[i].X[j]) expect := quad.err(stepsize, float64(i)) - if true || got > expect { + if got > expect { t.Errorf("error %g is greater than permitted tolerance %g", got, expect) } } From 06913c30f71770751155edde626121f17ac1622e Mon Sep 17 00:00:00 2001 From: soypat Date: Fri, 21 May 2021 22:48:35 -0300 Subject: [PATCH 03/10] ivp now generic. input vector no longer part of IVP --- ivp/_ignore.go | 24 ++++++++++++++++++++++++ ivp/dopri5.go | 36 +++++++++++++++++++----------------- ivp/ivp.go | 40 +++++++++++----------------------------- ivp/ivp_test.go | 9 +++++---- ivp/model.go | 14 +++++++------- ivp/rungekutta4.go | 35 ++++++++++++++--------------------- ivp/rungekutta4_test.go | 12 ++++++------ 7 files changed, 86 insertions(+), 84 deletions(-) create mode 100644 ivp/_ignore.go diff --git a/ivp/_ignore.go b/ivp/_ignore.go new file mode 100644 index 0000000..d881490 --- /dev/null +++ b/ivp/_ignore.go @@ -0,0 +1,24 @@ +package ivp + +/* idea drafting for possible future shape of IVPs with an input vector. + +// U is a vector which is a function of the current +// state. Put simply, the next input is a function of all current state variables +// and, possibly, current input as well. +// U_next = F_u(t, X, U) +type NAIVP interface { + // The input functions (ufunc) are not differential equations but rather + // calculated directly from a given x vector and current input vector. + IVP + // Where F_u is ufunc as returned by Equations() + // An initial value problem is characterized by boundary conditions imposed + // on the state vector X at the beginning of the integration domain. These + // boundary conditions are returned by the IV() method for the state vector + // as x0 and for the input vector as u0. + Inputs() (ufunc func(y []float64, t float64, x []float64)) + // Dimensions of x state vector and u inputs. // If problem has no input functions then u supplied and ufunc returned + // may be nil. x equations my not be nil. + Dims() (nx, nu int) +} + +*/ diff --git a/ivp/dopri5.go b/ivp/dopri5.go index 5887987..10cb496 100644 --- a/ivp/dopri5.go +++ b/ivp/dopri5.go @@ -9,30 +9,32 @@ import ( type DoPri5 struct { maxError, minStep, maxStep float64 adaptive bool - x, u []float64 + x []float64 // solutions y5, err45 []float64 // integration coefficients (for Butcher Tableau) k1, k2, k3, k4, k5, k6, k7 []float64 dom float64 - fx, fu func(y []float64, t float64, x, u []float64) + fx func(y []float64, t float64, x []float64) } // Set implements the Integrator interface. Initializes a Dormand Prince method // for use as a solver func (dp *DoPri5) Set(t0 float64, ivp IVP) error { dp.dom = t0 - xequations, ufunc := ivp.Equations() + xequations := ivp.Equations() dp.fx = xequations - dp.fu = ufunc - nx, nu := ivp.Dims() + // dp.fu = ufunc + nx := ivp.IV().Len() dp.k1, dp.k2, dp.k3 = make([]float64, nx), make([]float64, nx), make([]float64, nx) dp.k4, dp.k5, dp.k6 = make([]float64, nx), make([]float64, nx), make([]float64, nx) dp.k7 = make([]float64, nx) dp.x = make([]float64, nx) + for i := range dp.x { + dp.x[i] = ivp.IV().AtVec(i) + } dp.y5 = make([]float64, nx) dp.err45 = make([]float64, nx) - dp.u = make([]float64, nu) return nil } @@ -51,34 +53,34 @@ func (dp *DoPri5) Step(y4 []float64, h float64) (float64, error) { const b1, b3, b4, b5, b6 = 35. / 384., 500. / 1113., 125. / 192., -2187. / 6784., 11. / 84. // prettier variable names, also someone once said pointers are equivalent to JMP - x, u := dp.x, dp.u + x := dp.x F := dp.fx t := dp.dom // if adaptive then this tag will be used until h==DoPri5.MinStep SOLVE: - F(dp.k1, t, x, u) + F(dp.k1, t, x) floats.Scale(h, dp.k1) floats.AddScaledTo(dp.k2, x, c21, dp.k1) - F(dp.k2, t+c20*h, x, u) + F(dp.k2, t+c20*h, x) floats.Scale(h, dp.k2) floats.AddScaledTo(dp.k3, x, c31, dp.k1) floats.AddScaled(dp.k3, c32, dp.k2) - F(dp.k3, t+c30*h, x, u) + F(dp.k3, t+c30*h, x) floats.Scale(h, dp.k3) floats.AddScaledTo(dp.k4, x, c41, dp.k1) floats.AddScaled(dp.k4, c42, dp.k2) floats.AddScaled(dp.k4, c43, dp.k3) - F(dp.k4, t+c40*h, x, u) + F(dp.k4, t+c40*h, x) floats.Scale(h, dp.k4) floats.AddScaledTo(dp.k5, x, c51, dp.k1) floats.AddScaled(dp.k5, c52, dp.k2) floats.AddScaled(dp.k5, c53, dp.k3) floats.AddScaled(dp.k5, c54, dp.k4) - F(dp.k5, t+c50*h, x, u) + F(dp.k5, t+c50*h, x) floats.Scale(h, dp.k5) floats.AddScaledTo(dp.k6, x, c61, dp.k1) @@ -86,7 +88,7 @@ SOLVE: floats.AddScaled(dp.k6, c63, dp.k3) floats.AddScaled(dp.k6, c64, dp.k4) floats.AddScaled(dp.k6, c65, dp.k5) - F(dp.k6, t+c60*h, x, u) + F(dp.k6, t+c60*h, x) floats.Scale(h, dp.k6) // fourth order approximation used to advance solution @@ -103,7 +105,7 @@ SOLVE: floats.AddScaled(dp.k7, c74, dp.k4) floats.AddScaled(dp.k7, c75, dp.k5) floats.AddScaled(dp.k7, c76, dp.k6) - F(dp.k7, t+h*c70, x, u) + F(dp.k7, t+h*c70, x) floats.Scale(h, dp.k7) floats.AddScaledTo(dp.y5, x, a1, dp.k1) @@ -124,12 +126,12 @@ SOLVE: } // advance solution with fourth order solution copy(dp.x, y4) - if dp.fu != nil { - dp.fu(u, t+h, dp.x, u) - } return h, nil } +// XLen returns length of x state vector +func (dp *DoPri5) XLen() (nx int) { return len(dp.x) } + // NewDormandPrince5 returns a adaptive-ready Dormand Prince solver of order 5 func NewDormandPrince5(maxError, maxStep, minStep float64) *DoPri5 { dp := new(DoPri5) diff --git a/ivp/ivp.go b/ivp/ivp.go index 4194c22..572904a 100644 --- a/ivp/ivp.go +++ b/ivp/ivp.go @@ -12,47 +12,29 @@ import ( ) // IVP defines a multivariable, non-autonomous initial value problem. -// It is worth mentioning system does not necessarily need -// to be non-autonomous. https://en.wikipedia.org/wiki/Autonomous_system_(mathematics) -// +// It is worth mentioning the system need not be non-autonomous. https://en.wikipedia.org/wiki/Autonomous_system_(mathematics) // // These problems have the form (capital letters are vectors in this example) -// X'(t) = F(t, X, U) +// X'(t) = F(t, X) // X(0) = X_0 // // Where X' is the vector of first derivatives of the state vector X. // F would be xequations as returned by Equations(). t is // a scalar representing the integrations domain, which is usually time -// for most physical problems. U is a vector which is a function of the current -// state. Put simply, the next input is a function of all current state variables -// and, possibly, current input as well. -// U_next = F_u(t, X, U) -// -// Where F_u is ufunc as returned by Equations() -// An initial value problem is characterized by boundary conditions imposed -// on the state vector X at the beginning of the integration domain. These -// boundary conditions are returned by the IV() method for the state vector -// as x0 and for the input vector as u0. +// for most physical problems. // // The term "state vector" and "state variables" are used interchangeably // throughout the code and refer to X vector of independent variables. type IVP interface { - // Initial values vector for state variables x and inputs u. x0 defines + // Initial values vector for state variables x. x0 defines // the first values the state vector takes when integration begins. - IV() (x0, u0 mat.Vector) + IV() (x0 mat.Vector) // Equations returns the coupled, non-linear algebraic differential - // equations (xequations) for the state variables (x) and the functions for inputs (u). - // The input functions (ufunc) are not differential equations but rather - // calculated directly from a given x vector and current input vector. - // Results are stored in y which are the length of x and u, respectively. + // equations (xequations) for the state variables (x) + // Results are stored in y which are the length of x. // The scalar (float64) argument is the domain over which the // problem is integrated, which is usually time for most physical problems. - // - // If problem has no input functions then u supplied and ufunc returned - // may be nil. x equations my not be nil. - Equations() (xequations func(y []float64, t float64, x, u []float64), ufunc func(u_next []float64, t float64, x, u []float64)) - // Dimensions of x state vector and u inputs. - Dims() (nx, nu int) + Equations() (xequations func(y []float64, t float64, x []float64)) } // Integrator abstracts algorithm specifics. For anyone looking to @@ -65,8 +47,8 @@ type Integrator interface { // Set initializes an initial value problem. First argument // is the initial domain integration point, is usually zero. Set(float64, IVP) error - // Dimensions of x state vector and u inputs. Set must be called first. - Dims() (nx, nu int) + // Length of state vector x + XLen() (nx int) // Step integrates IVP and stores result in y. step is a suggested step // for the algorithm to take. The algorithm may decide that it is not sufficiently // small or big enough (these are adaptive algorithms) and take a different step. @@ -85,7 +67,7 @@ func Solve(solver Integrator, stepsize, domainLength float64) (results []result, const maxAllocGB = 2 integrated := 0. expectedLength := int(domainLength/stepsize) + 1 - nx, _ := solver.Dims() + nx := solver.XLen() if nx == 0 { return nil, errors.New("state vector length can not be equal to zero. Has ivp been set?") } diff --git a/ivp/ivp_test.go b/ivp/ivp_test.go index e96b1b7..974ef63 100644 --- a/ivp/ivp_test.go +++ b/ivp/ivp_test.go @@ -10,7 +10,8 @@ import ( "math" "testing" - "gonum.org/v1/gonum/ivp" + "gonum.org/v1/exp/ivp" + "gonum.org/v1/gonum/mat" ) @@ -21,14 +22,14 @@ func TestSolve(t *testing.T) { stepsize := 0.1 results, err := ivp.Solve(solver, stepsize, 10) quadsol := quad.solution - sol, _ := quadsol.Equations() + sol := quadsol.Equations() if err != nil { t.Fatal(err) } nxs, _ := quadsol.Dims() solresults := make([]float64, nxs) for i := range results { - sol(solresults, results[i].DomainStartOffset, results[i].X, nil) + sol(solresults, results[i].DomainStartOffset, results[i].X) for j := range solresults { got := math.Abs(solresults[j] - results[i].X[j]) expect := quad.err(stepsize, float64(i)) @@ -45,7 +46,7 @@ func Example_solve() { ) // we declare our physical model in the following function ballModel, err := ivp.NewModel(mat.NewVecDense(2, []float64{100., 0.}), - nil, func(yvec []float64, _ float64, xvec, _ []float64) { + nil, func(yvec []float64, _ float64, xvec []float64) { // this anonymous function defines the physics. // The first variable xvec[0] corresponds to position // second variable xvec[1] is velocity. diff --git a/ivp/model.go b/ivp/model.go index 46491ba..08b05db 100644 --- a/ivp/model.go +++ b/ivp/model.go @@ -8,21 +8,21 @@ import "gonum.org/v1/gonum/mat" // Model implements IVP interface. X vector and equations can not be nil or zero length. type Model struct { - x0, u0 mat.Vector - xeq, ins func(y []float64, dom float64, x, u []float64) + x0, u0 mat.Vector + xeq func(y []float64, dom float64, x []float64) } // IV returns initial values of the IVP. First returned parameter is the // starting x vector and second parameter are inputs when solving non-autonomous // ODEs. -func (m *Model) IV() (mat.Vector, mat.Vector) { return m.x0, m.u0 } +func (m *Model) IV() mat.Vector { return m.x0 } // Equations returns differential equations relating to state vector x and input functions // for non-autonomous ODEs. // // Input functions may be nil (ueq). -func (m *Model) Equations() (xeq, ueq func(y []float64, dom float64, x, u []float64)) { - return m.xeq, m.ins +func (m *Model) Equations() (xeq func(y []float64, dom float64, x []float64)) { + return m.xeq } // Dims returns dimension of state and input (x length and u length, respectively). @@ -41,6 +41,6 @@ func (m *Model) Dims() (nx, nu int) { // NewModel returns a IVP given initial conditions (x0,u0), differential equations (xeq) and // input functions for non-autonomous ODEs (ueq). -func NewModel(x0, u0 mat.Vector, xeq, ueq func(y []float64, dom float64, x, u []float64)) (*Model, error) { - return &Model{xeq: xeq, ins: ueq, x0: x0, u0: u0}, nil +func NewModel(x0, u0 mat.Vector, xeq, ueq func(y []float64, dom float64, x []float64)) (*Model, error) { + return &Model{xeq: xeq, x0: x0, u0: u0}, nil } diff --git a/ivp/rungekutta4.go b/ivp/rungekutta4.go index 6af223d..4283980 100644 --- a/ivp/rungekutta4.go +++ b/ivp/rungekutta4.go @@ -14,9 +14,9 @@ import ( // 4th order multivariable method (non adaptive) type RK4 struct { dom float64 - x, u []float64 + x []float64 a, b, c, d []float64 - fx, fu func(y []float64, t float64, x, u []float64) + fx func(y []float64, t float64, x []float64) } // Step implements Integrator interface @@ -24,10 +24,10 @@ func (rk *RK4) Step(y []float64, h float64) (float64, error) { const overSix = 1. / 6. t := rk.dom // set a, b, c, d (in some literatures these are k1,k2,k3,k4) - rk.fx(rk.a, t, rk.x, rk.u) - rk.fx(rk.b, t+0.5*h, floats.AddScaledTo(rk.b, rk.x, 0.5*h, rk.a), rk.u) - rk.fx(rk.c, t+0.5*h, floats.AddScaledTo(rk.c, rk.x, 0.5*h, rk.b), rk.u) - rk.fx(rk.d, t+h, floats.AddScaledTo(rk.d, rk.x, h, rk.c), rk.u) + rk.fx(rk.a, t, rk.x) + rk.fx(rk.b, t+0.5*h, floats.AddScaledTo(rk.b, rk.x, 0.5*h, rk.a)) + rk.fx(rk.c, t+0.5*h, floats.AddScaledTo(rk.c, rk.x, 0.5*h, rk.b)) + rk.fx(rk.d, t+h, floats.AddScaledTo(rk.d, rk.x, h, rk.c)) floats.Add(rk.a, rk.d) floats.Add(rk.b, rk.c) @@ -36,39 +36,32 @@ func (rk *RK4) Step(y []float64, h float64) (float64, error) { // finished integrating. Now we update RK4 structure for future Steps copy(rk.x, y) // store results rk.dom += h - if len(rk.u) > 0 { - rk.fu(rk.u, t, rk.x, rk.u) - } return h, nil } -// Return size of state vector and input vector. Solver must be set beforehand. -func (rk *RK4) Dims() (nx, nu int) { return len(rk.x), len(rk.u) } - // Set implements integrator interface. All RK4 data // is reset when calling Set. func (rk *RK4) Set(d0 float64, ivp IVP) error { if ivp == nil { return errors.New("IVP is nil") } - nx, nu := ivp.Dims() + x0 := ivp.IV() + nx := x0.Len() rk.dom = d0 //set start domain rk.a = make([]float64, nx) rk.b = make([]float64, nx) rk.c = make([]float64, nx) rk.d = make([]float64, nx) rk.x = make([]float64, nx) - rk.u = make([]float64, nu) // set initial values - x0, u0 := ivp.IV() + for i := 0; i < x0.Len(); i++ { rk.x[i] = x0.AtVec(i) } - if nu > 0 { - for i := 0; i < u0.Len(); i++ { - rk.u[i] = u0.AtVec(i) - } - } - rk.fx, rk.fu = ivp.Equations() + + rk.fx = ivp.Equations() return nil } + +// XLen returns length of x state vector +func (rk *RK4) XLen() (nx int) { return len(rk.x) } diff --git a/ivp/rungekutta4_test.go b/ivp/rungekutta4_test.go index e4d7192..bc38940 100644 --- a/ivp/rungekutta4_test.go +++ b/ivp/rungekutta4_test.go @@ -10,7 +10,7 @@ import ( "math" "testing" - "gonum.org/v1/gonum/ivp" + "gonum.org/v1/exp/ivp" "gonum.org/v1/gonum/mat" ) @@ -36,13 +36,13 @@ func TestQuadratic(t *testing.T) { results := make([]float64, nx) solmodel := Quadratic.solution - soleq, _ := solmodel.Equations() + soleq := solmodel.Equations() solDims, _ := solmodel.Dims() solution := make([]float64, solDims) for i := 1.; i < float64(steps+1); i++ { dom := i * dt solver.Step(results, dt) - soleq(solution, dom, results, nil) + soleq(solution, dom, results) for j := range results { got := math.Abs(solution[j] - results[j]) expected := Quadratic.err(dt, i) @@ -60,7 +60,7 @@ func Example_fallingBall() { ) // we declare our physical model in the following function ballModel, err := ivp.NewModel(mat.NewVecDense(2, []float64{100., 0.}), - nil, func(yvec []float64, _ float64, xvec, _ []float64) { + nil, func(yvec []float64, _ float64, xvec []float64) { // this anonymous function defines the physics. // The first variable xvec[0] corresponds to position // second variable xvec[1] is velocity. @@ -109,7 +109,7 @@ func Example_fallingBall() { func quadTestModel(t *testing.T) *TestModel { Quadratic := new(TestModel) quad, err := ivp.NewModel(mat.NewVecDense(2, []float64{0, 0}), - nil, func(y []float64, dom float64, x, u []float64) { + nil, func(y []float64, dom float64, x []float64) { y[0], y[1] = x[1], 1. }, nil) if err != nil { @@ -117,7 +117,7 @@ func quadTestModel(t *testing.T) *TestModel { } Quadratic.Model = quad quadsol, err := ivp.NewModel(mat.NewVecDense(2, []float64{0, 0}), - nil, func(y []float64, dom float64, x, u []float64) { + nil, func(y []float64, dom float64, x []float64) { y[0], y[1] = dom*dom/2., dom }, nil) if err != nil { From 37348c2dcedd2b44fe53c2f0ef99cc78b8057e37 Mon Sep 17 00:00:00 2001 From: soypat Date: Sat, 22 May 2021 09:52:00 -0300 Subject: [PATCH 04/10] add Configuration/parameters to modify solver behaviour --- ivp/config.go | 65 +++++++++++++++++++++++++++++++++++++++++++++++++++ ivp/dopri5.go | 34 +++++++++++++++++---------- 2 files changed, 87 insertions(+), 12 deletions(-) create mode 100644 ivp/config.go diff --git a/ivp/config.go b/ivp/config.go new file mode 100644 index 0000000..383416e --- /dev/null +++ b/ivp/config.go @@ -0,0 +1,65 @@ +package ivp + +import "errors" + +// Configuration represents a change in default parameters in the solving +// of IVPs. Parameters are often needed to enable step size adaptivity and to +// overwrite solver default values. +type Configuration func(*parameters) error + +type parameters struct { + // step length + minStep, maxStep float64 + // Relative and absolute error tolerances + rtol, atol float64 + // parameters for step size selection. + // The new step size is subject to restriction minStepRatio <= hnew/hold <= maxStepRatio + minStepRatio, maxStepRatio float64 +} + +// ConfigStepLimits sets minimum step length and max step length +// an integrator can take when using an adaptive method. +func ConfigStepLimits(minStep, maxStep float64) Configuration { + return func(p *parameters) error { + if minStep <= 0 || minStep > maxStep { + return errors.New("minimum step length too small or greater than max step length") + } + p.maxStep, p.minStep = maxStep, minStep + return nil + } +} + +// ConfigScalarTolerance Sets relative and absolute tolerances for controlling step error. +// If value passed is zero the corresponding tolerance is not modified. +func ConfigScalarTolerance(rtol, atol float64) Configuration { + return func(p *parameters) error { + if rtol < 0 || atol < 0 { + return errors.New("negative error tolerance") + } + if rtol > 0 { + p.rtol = rtol + } + if atol > 0 { + p.atol = atol + } + return nil + } +} + +// ConfigAdaptiveRatioLimits sets parameters for step size selection. +// The new step size is subject to restriction min <= hnew/hold <= max. +// If min or max are zero then value is not changed +func ConfigAdaptiveRatioLimits(min, max float64) Configuration { + return func(p *parameters) error { + if min < 0 || max < min { + return errors.New("bad step size selection ratio") + } + if min > 0 { + p.minStepRatio = min + } + if max > 0 { + p.maxStepRatio = max + } + return nil + } +} diff --git a/ivp/dopri5.go b/ivp/dopri5.go index 10cb496..332de2d 100644 --- a/ivp/dopri5.go +++ b/ivp/dopri5.go @@ -7,9 +7,9 @@ import ( ) type DoPri5 struct { - maxError, minStep, maxStep float64 - adaptive bool - x []float64 + param parameters + adaptive bool + x []float64 // solutions y5, err45 []float64 // integration coefficients (for Butcher Tableau) @@ -116,10 +116,10 @@ SOLVE: floats.AddScaled(dp.y5, a7, dp.k7) // compute error between fifth order solution and fourth order solution - errRatio := dp.maxError / floats.Norm(floats.SubTo(dp.err45, y4, dp.y5), math.Inf(1)) - hnew := math.Min(math.Max(0.9*h*math.Pow(errRatio, .2), dp.minStep), dp.maxStep) + errRatio := dp.param.atol / floats.Norm(floats.SubTo(dp.err45, y4, dp.y5), math.Inf(1)) + hnew := math.Min(math.Max(0.9*h*math.Pow(errRatio, .2), dp.param.minStep), dp.param.maxStep) // given "bad" error ratio, algorithm will recompute with a smaller step, given it is not the minimum permissible step or less - if errRatio < 1 && h > dp.minStep { + if errRatio < 1 && h > dp.param.minStep { h = hnew goto SOLVE } @@ -133,14 +133,24 @@ SOLVE: func (dp *DoPri5) XLen() (nx int) { return len(dp.x) } // NewDormandPrince5 returns a adaptive-ready Dormand Prince solver of order 5 -func NewDormandPrince5(maxError, maxStep, minStep float64) *DoPri5 { +// +// 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)) +// +// If a invalid configuration is passed the function panics. +func NewDormandPrince5(configs ...Configuration) *DoPri5 { dp := new(DoPri5) - if maxError > 0 && maxStep > minStep && minStep > 0 { + for i := range configs { + err := configs[i](&dp.param) + if err != nil { + panic(err) + } + } + + if dp.param.atol > 0 && dp.param.minStep > 0 { dp.adaptive = true - dp.maxError = maxError - dp.minStep, dp.maxStep = minStep, maxStep - } else { - panic("NewDormandPrince5 generator requires proper arguments: maxError > 0 && maxStep > minStep && minStep > 0") } return dp } From f124a23d0e9d4a5a6386c19b477044f128c432b5 Mon Sep 17 00:00:00 2001 From: soypat Date: Sat, 22 May 2021 19:27:10 -0300 Subject: [PATCH 05/10] add draft of implicit Newton-Raphson method --- ivp/config.go | 15 ++++++ ivp/ivp.go | 1 - ivp/newton.go | 115 +++++++++++++++++++++++++++++++++++++++++++++ ivp/newton_test.go | 63 +++++++++++++++++++++++++ 4 files changed, 193 insertions(+), 1 deletion(-) create mode 100644 ivp/newton.go create mode 100644 ivp/newton_test.go diff --git a/ivp/config.go b/ivp/config.go index 383416e..e629dfb 100644 --- a/ivp/config.go +++ b/ivp/config.go @@ -15,6 +15,8 @@ type parameters struct { // parameters for step size selection. // The new step size is subject to restriction minStepRatio <= hnew/hold <= maxStepRatio minStepRatio, maxStepRatio float64 + // relaxation factor offers stability for implicit solvers + relax float64 } // ConfigStepLimits sets minimum step length and max step length @@ -63,3 +65,16 @@ func ConfigAdaptiveRatioLimits(min, max float64) Configuration { return nil } } + +// ConfigRelaxation sets the relaxation factor for methods +// like Newton-Raphson which can be unstable in certain cases +// without it. +func ConfigRelaxation(alpha float64) Configuration { + return func(p *parameters) error { + if alpha < 0 || alpha >= 1. { + return errors.New("relaxation factor must be in [0,1)") + } + p.relax = alpha + return nil + } +} diff --git a/ivp/ivp.go b/ivp/ivp.go index 572904a..43ff7a6 100644 --- a/ivp/ivp.go +++ b/ivp/ivp.go @@ -76,7 +76,6 @@ func Solve(solver Integrator, stepsize, domainLength float64) (results []result, return nil, fmt.Errorf("solution exceeds %dGB or not initialized (size is %dMB)", maxAllocGB, size/1e6) } results = make([]result, 0, expectedLength) - for integrated < domainLength { res := make([]float64, nx) stepsize, err = solver.Step(res, stepsize) diff --git a/ivp/newton.go b/ivp/newton.go new file mode 100644 index 0000000..8f477e5 --- /dev/null +++ b/ivp/newton.go @@ -0,0 +1,115 @@ +package ivp + +import ( + "math" + + "gonum.org/v1/exp/linsolve" + "gonum.org/v1/gonum/diff/fd" + "gonum.org/v1/gonum/floats" + "gonum.org/v1/gonum/mat" +) + +type NewtonRaphson struct { + param parameters + x, b []float64 + //guess contains best guess for next X vector which zeros the residual function + guess, err []float64 + + dom float64 + fx func(y []float64, t float64, x []float64) + residualers func(h, t float64, xnow []float64) func(y, xnext []float64) + // Residual functions jacobian + J *mat.Dense + Jband *mat.BandDense + result *linsolve.Result +} + +func (nr *NewtonRaphson) Step(y []float64, h float64) (float64, error) { + var err error + if nr.param.atol <= 0 { + panic("Newton method requires error") + } + guess := nr.guess + + nx := len(nr.x) + copy(guess, nr.x) + t := nr.dom + h + iter := 0 + // while step length does not change in iter loop, this may remain outside + F := nr.residualers(h, t, nr.x) + aerr := math.Inf(1) + for iter < 10 && aerr > nr.param.atol { + // set b vector with "old" x vector + F(nr.b, guess) + b := mat.NewVecDense(nx, nr.b) + + fd.Jacobian(nr.J, F, guess, nil) + for i := 0; i < nx; i++ { + for j := 0; j < nx; j++ { + nr.Jband.SetBand(i, j, nr.J.At(i, j)) // we clone Jacobian Dense to our Jacobian BandDense for MulVecToer + } + } + nr.result, err = linsolve.Iterative(nr.Jband, b, &linsolve.GMRES{}, &linsolve.Settings{MaxIterations: 2}) + if err != nil { + return 0, err + } + resx := nr.result.X.RawVector().Data + // X_(i+1) = X_(i) - alpha * F(X_(g)) / J(X_(g)) where g are guesses, and alpha is the relaxation factor + // resx now contains the next guess. we proceed to calculate error. + floats.AddScaledTo(resx, guess, nr.param.relax-1., resx) + // error calculation if enabled + if nr.param.atol > 0 { + floats.SubTo(nr.err, guess, resx) + aerr = floats.Norm(nr.err, math.Inf(1)) + } + // advance solution + copy(guess, resx) + iter++ + } + copy(nr.x, guess) + copy(y, guess) + return h, nil +} + +func (nr *NewtonRaphson) Set(t0 float64, ivp IVP) error { + nr.dom = t0 + xequations := ivp.Equations() + nr.fx = xequations + // nr.fu = ufunc + nx := ivp.IV().Len() + + nr.x = make([]float64, nx) + for i := range nr.x { + nr.x[i] = ivp.IV().AtVec(i) + } + nr.err = make([]float64, nx) + nr.guess = make([]float64, nx) + + nr.b = make([]float64, nx) + nr.J = mat.NewDense(nx, nx, nil) + nr.Jband = mat.NewBandDense(nx, nx, nx-1, nx-1, nil) + + // First propose residual functions such that + // F(X_(i+1)) = 0 = X_(i+1) - X_(i) - step * f(X_(i+1)) + // where f is the vector of differential equations. We + // seek the zero of the function F evaluated at xnext a.k.a X_(i+1) + nr.residualers = func(h, t float64, xnow []float64) func(r, xnext []float64) { + return func(r, xnext []float64) { + // store results of function in r + nr.fx(r, t, xnext) + floats.SubTo(r, xnext, floats.AddTo(r, xnow, floats.ScaleTo(r, h, r))) + } + } + return nil +} +func (nr *NewtonRaphson) XLen() int { return len(nr.x) } +func NewNewtonRaphson(configs ...Configuration) *NewtonRaphson { + nr := new(NewtonRaphson) + for i := range configs { + err := configs[i](&nr.param) + if err != nil { + panic(err) + } + } + return nr +} diff --git a/ivp/newton_test.go b/ivp/newton_test.go new file mode 100644 index 0000000..b8ed56e --- /dev/null +++ b/ivp/newton_test.go @@ -0,0 +1,63 @@ +package ivp_test + +import ( + "math" + "testing" + + "gonum.org/v1/exp/ivp" + "gonum.org/v1/gonum/mat" +) + +// This is a stiff equation https://en.wikipedia.org/wiki/Stiff_equation +// y'(t) = -15*y(t) +// solution: y(t) = exp(-15*t) +func exponential1DTestModel(t *testing.T) *TestModel { + tau := -15. + Stiff := new(TestModel) + stiff, err := ivp.NewModel(mat.NewVecDense(1, []float64{1.}), + nil, func(y []float64, dom float64, x []float64) { + y[0] = tau * x[0] + }, nil) + if err != nil { + t.Fatal(err) + } + Stiff.Model = stiff + stiffsol, err := ivp.NewModel(mat.NewVecDense(1, []float64{0}), + nil, func(y []float64, dom float64, x []float64) { + y[0] = math.Exp(tau * dom) + }, nil) + if err != nil { + t.Fatal(err) + } + Stiff.solution = stiffsol + Stiff.err = func(h, i float64) float64 { return 2 * h * i } + return Stiff +} + +func TestNewton1DStiff(t *testing.T) { + stiff1D := exponential1DTestModel(t) + solver := ivp.NewNewtonRaphson( + ivp.ConfigScalarTolerance(0, 1e-5), + ) + solver.Set(0.0, stiff1D) + stepsize := 1. / 30. + results, err := ivp.Solve(solver, stepsize, 0.5) + if err != nil { + t.Fatal(err) + } + + quadsol := stiff1D.solution + sol := quadsol.Equations() + nxs, _ := quadsol.Dims() + solresults := make([]float64, nxs) + for i := range results { + sol(solresults, results[i].DomainStartOffset, results[i].X) + for j := range solresults { + got := math.Abs(solresults[j] - results[i].X[j]) + expect := stiff1D.err(stepsize, float64(i)) + if got > expect { + t.Errorf("error %g is greater than permitted tolerance %g", got, expect) + } + } + } +} From f05c7f89673ebfd6cf1580de9c14cc8fd49d982a Mon Sep 17 00:00:00 2001 From: soypat Date: Sat, 22 May 2021 19:54:48 -0300 Subject: [PATCH 06/10] change IVP description. thanks @acanalis --- ivp/ivp.go | 16 ++++++++++------ 1 file changed, 10 insertions(+), 6 deletions(-) diff --git a/ivp/ivp.go b/ivp/ivp.go index 43ff7a6..11b49e3 100644 --- a/ivp/ivp.go +++ b/ivp/ivp.go @@ -18,19 +18,23 @@ import ( // X'(t) = F(t, X) // X(0) = X_0 // -// Where X' is the vector of first derivatives of the state vector X. -// F would be xequations as returned by Equations(). t is -// a scalar representing the integrations domain, which is usually time -// for most physical problems. +// Where: +// t is a scalar representing the integration domain, which is time for most physical problems. +// X' is the vector of first derivatives of the state vector X. +// F would be xequations as returned by Equations(). +// An initial value problem is characterized by the initial conditions imposed +// on the state vector X at the beginning of the integration domain. These +// initial conditions are returned by the IV() method for the state vector +// as x0 and for the input vector as u0. // // The term "state vector" and "state variables" are used interchangeably // throughout the code and refer to X vector of independent variables. type IVP interface { - // Initial values vector for state variables x. x0 defines + // Initial values vector for state vector x. x0 defines // the first values the state vector takes when integration begins. IV() (x0 mat.Vector) // Equations returns the coupled, non-linear algebraic differential - // equations (xequations) for the state variables (x) + // equations (xequations) for the state vector (x) // Results are stored in y which are the length of x. // The scalar (float64) argument is the domain over which the // problem is integrated, which is usually time for most physical problems. From d2407cbaf9235ff9f3c566e5c5ae94a3f59195ae Mon Sep 17 00:00:00 2001 From: soypat Date: Tue, 25 May 2021 11:51:31 -0300 Subject: [PATCH 07/10] accept change proposed by @vladimir-ch --- ivp/_ignore.go | 24 ------ ivp/ivp_test.go | 77 ------------------ ivp/model.go | 46 ----------- ivp/newton.go | 115 --------------------------- ivp/newton_test.go | 63 --------------- ivp/rungekutta4.go | 67 ---------------- ivp/rungekutta4_test.go | 129 ------------------------------ {ivp => ode}/README.md | 2 +- {ivp => ode}/config.go | 2 +- {ivp => ode}/doc.go | 4 +- {ivp => ode}/dopri5.go | 167 +++++++++++++++++++------------------- {ivp => ode}/ivp.go | 57 ++++++------- ode/ivp_test.go | 172 ++++++++++++++++++++++++++++++++++++++++ ode/model.go | 41 ++++++++++ 14 files changed, 332 insertions(+), 634 deletions(-) delete mode 100644 ivp/_ignore.go delete mode 100644 ivp/ivp_test.go delete mode 100644 ivp/model.go delete mode 100644 ivp/newton.go delete mode 100644 ivp/newton_test.go delete mode 100644 ivp/rungekutta4.go delete mode 100644 ivp/rungekutta4_test.go rename {ivp => ode}/README.md (73%) rename {ivp => ode}/config.go (99%) rename {ivp => ode}/doc.go (66%) rename {ivp => ode}/dopri5.go (50%) rename {ivp => ode}/ivp.go (72%) create mode 100644 ode/ivp_test.go create mode 100644 ode/model.go diff --git a/ivp/_ignore.go b/ivp/_ignore.go deleted file mode 100644 index d881490..0000000 --- a/ivp/_ignore.go +++ /dev/null @@ -1,24 +0,0 @@ -package ivp - -/* idea drafting for possible future shape of IVPs with an input vector. - -// U is a vector which is a function of the current -// state. Put simply, the next input is a function of all current state variables -// and, possibly, current input as well. -// U_next = F_u(t, X, U) -type NAIVP interface { - // The input functions (ufunc) are not differential equations but rather - // calculated directly from a given x vector and current input vector. - IVP - // Where F_u is ufunc as returned by Equations() - // An initial value problem is characterized by boundary conditions imposed - // on the state vector X at the beginning of the integration domain. These - // boundary conditions are returned by the IV() method for the state vector - // as x0 and for the input vector as u0. - Inputs() (ufunc func(y []float64, t float64, x []float64)) - // Dimensions of x state vector and u inputs. // If problem has no input functions then u supplied and ufunc returned - // may be nil. x equations my not be nil. - Dims() (nx, nu int) -} - -*/ diff --git a/ivp/ivp_test.go b/ivp/ivp_test.go deleted file mode 100644 index 974ef63..0000000 --- a/ivp/ivp_test.go +++ /dev/null @@ -1,77 +0,0 @@ -// Copyright ©2021 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 ivp_test - -import ( - "fmt" - "log" - "math" - "testing" - - "gonum.org/v1/exp/ivp" - - "gonum.org/v1/gonum/mat" -) - -func TestSolve(t *testing.T) { - quad := quadTestModel(t) - solver := new(ivp.RK4) - solver.Set(0.0, quad) - stepsize := 0.1 - results, err := ivp.Solve(solver, stepsize, 10) - quadsol := quad.solution - sol := quadsol.Equations() - if err != nil { - t.Fatal(err) - } - nxs, _ := quadsol.Dims() - solresults := make([]float64, nxs) - for i := range results { - sol(solresults, results[i].DomainStartOffset, results[i].X) - for j := range solresults { - got := math.Abs(solresults[j] - results[i].X[j]) - expect := quad.err(stepsize, float64(i)) - if got > expect { - t.Errorf("error %g is greater than permitted tolerance %g", got, expect) - } - } - } -} - -func Example_solve() { - const ( - g = -10. // gravity field [m.s^-2] - ) - // we declare our physical model in the following function - ballModel, err := ivp.NewModel(mat.NewVecDense(2, []float64{100., 0.}), - nil, func(yvec []float64, _ float64, xvec []float64) { - // this anonymous function defines the physics. - // The first variable xvec[0] corresponds to position - // second variable xvec[1] is velocity. - Dx := xvec[1] - // yvec represents change in xvec, or derivative with respect to domain - // Change in position will be equal to velocity, which is the second variable: - // thus yvec[0] = xvec[1], which is the same as saying "change in xvec[0]" is equal to xvec[1] - yvec[0] = Dx - // change in velocity is acceleration. We suppose our ball is on earth accelerating at `g` - yvec[1] = g - }, nil) - if err != nil { - log.Fatal(err) - } - // Here we choose our algorithm. Runge-Kutta 4th order is used - var solver ivp.Integrator = new(ivp.RK4) - // Before integrating the IVP is passed to the integrator (a.k.a solver). Domain (time) starts at 0 - dom := 0. - err = solver.Set(dom, ballModel) - if err != nil { - log.Fatal(err) - } - // 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. - results, err := ivp.Solve(solver, 0.1, 10.) - fmt.Println(results) -} diff --git a/ivp/model.go b/ivp/model.go deleted file mode 100644 index 08b05db..0000000 --- a/ivp/model.go +++ /dev/null @@ -1,46 +0,0 @@ -// Copyright ©2021 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 ivp - -import "gonum.org/v1/gonum/mat" - -// Model implements IVP interface. X vector and equations can not be nil or zero length. -type Model struct { - x0, u0 mat.Vector - xeq func(y []float64, dom float64, x []float64) -} - -// IV returns initial values of the IVP. First returned parameter is the -// starting x vector and second parameter are inputs when solving non-autonomous -// ODEs. -func (m *Model) IV() mat.Vector { return m.x0 } - -// Equations returns differential equations relating to state vector x and input functions -// for non-autonomous ODEs. -// -// Input functions may be nil (ueq). -func (m *Model) Equations() (xeq func(y []float64, dom float64, x []float64)) { - return m.xeq -} - -// Dims returns dimension of state and input (x length and u length, respectively). -// Dims has some dimension checking involved and can serve as a preliminary system verifier. -// -// Dims panics if x vector is nil. -func (m *Model) Dims() (nx, nu int) { - if m.x0 == nil { - panic("x vector can not be nil") - } - if m.u0 == nil { - nu = 0 - } - return m.x0.Len(), nu -} - -// NewModel returns a IVP given initial conditions (x0,u0), differential equations (xeq) and -// input functions for non-autonomous ODEs (ueq). -func NewModel(x0, u0 mat.Vector, xeq, ueq func(y []float64, dom float64, x []float64)) (*Model, error) { - return &Model{xeq: xeq, x0: x0, u0: u0}, nil -} diff --git a/ivp/newton.go b/ivp/newton.go deleted file mode 100644 index 8f477e5..0000000 --- a/ivp/newton.go +++ /dev/null @@ -1,115 +0,0 @@ -package ivp - -import ( - "math" - - "gonum.org/v1/exp/linsolve" - "gonum.org/v1/gonum/diff/fd" - "gonum.org/v1/gonum/floats" - "gonum.org/v1/gonum/mat" -) - -type NewtonRaphson struct { - param parameters - x, b []float64 - //guess contains best guess for next X vector which zeros the residual function - guess, err []float64 - - dom float64 - fx func(y []float64, t float64, x []float64) - residualers func(h, t float64, xnow []float64) func(y, xnext []float64) - // Residual functions jacobian - J *mat.Dense - Jband *mat.BandDense - result *linsolve.Result -} - -func (nr *NewtonRaphson) Step(y []float64, h float64) (float64, error) { - var err error - if nr.param.atol <= 0 { - panic("Newton method requires error") - } - guess := nr.guess - - nx := len(nr.x) - copy(guess, nr.x) - t := nr.dom + h - iter := 0 - // while step length does not change in iter loop, this may remain outside - F := nr.residualers(h, t, nr.x) - aerr := math.Inf(1) - for iter < 10 && aerr > nr.param.atol { - // set b vector with "old" x vector - F(nr.b, guess) - b := mat.NewVecDense(nx, nr.b) - - fd.Jacobian(nr.J, F, guess, nil) - for i := 0; i < nx; i++ { - for j := 0; j < nx; j++ { - nr.Jband.SetBand(i, j, nr.J.At(i, j)) // we clone Jacobian Dense to our Jacobian BandDense for MulVecToer - } - } - nr.result, err = linsolve.Iterative(nr.Jband, b, &linsolve.GMRES{}, &linsolve.Settings{MaxIterations: 2}) - if err != nil { - return 0, err - } - resx := nr.result.X.RawVector().Data - // X_(i+1) = X_(i) - alpha * F(X_(g)) / J(X_(g)) where g are guesses, and alpha is the relaxation factor - // resx now contains the next guess. we proceed to calculate error. - floats.AddScaledTo(resx, guess, nr.param.relax-1., resx) - // error calculation if enabled - if nr.param.atol > 0 { - floats.SubTo(nr.err, guess, resx) - aerr = floats.Norm(nr.err, math.Inf(1)) - } - // advance solution - copy(guess, resx) - iter++ - } - copy(nr.x, guess) - copy(y, guess) - return h, nil -} - -func (nr *NewtonRaphson) Set(t0 float64, ivp IVP) error { - nr.dom = t0 - xequations := ivp.Equations() - nr.fx = xequations - // nr.fu = ufunc - nx := ivp.IV().Len() - - nr.x = make([]float64, nx) - for i := range nr.x { - nr.x[i] = ivp.IV().AtVec(i) - } - nr.err = make([]float64, nx) - nr.guess = make([]float64, nx) - - nr.b = make([]float64, nx) - nr.J = mat.NewDense(nx, nx, nil) - nr.Jband = mat.NewBandDense(nx, nx, nx-1, nx-1, nil) - - // First propose residual functions such that - // F(X_(i+1)) = 0 = X_(i+1) - X_(i) - step * f(X_(i+1)) - // where f is the vector of differential equations. We - // seek the zero of the function F evaluated at xnext a.k.a X_(i+1) - nr.residualers = func(h, t float64, xnow []float64) func(r, xnext []float64) { - return func(r, xnext []float64) { - // store results of function in r - nr.fx(r, t, xnext) - floats.SubTo(r, xnext, floats.AddTo(r, xnow, floats.ScaleTo(r, h, r))) - } - } - return nil -} -func (nr *NewtonRaphson) XLen() int { return len(nr.x) } -func NewNewtonRaphson(configs ...Configuration) *NewtonRaphson { - nr := new(NewtonRaphson) - for i := range configs { - err := configs[i](&nr.param) - if err != nil { - panic(err) - } - } - return nr -} diff --git a/ivp/newton_test.go b/ivp/newton_test.go deleted file mode 100644 index b8ed56e..0000000 --- a/ivp/newton_test.go +++ /dev/null @@ -1,63 +0,0 @@ -package ivp_test - -import ( - "math" - "testing" - - "gonum.org/v1/exp/ivp" - "gonum.org/v1/gonum/mat" -) - -// This is a stiff equation https://en.wikipedia.org/wiki/Stiff_equation -// y'(t) = -15*y(t) -// solution: y(t) = exp(-15*t) -func exponential1DTestModel(t *testing.T) *TestModel { - tau := -15. - Stiff := new(TestModel) - stiff, err := ivp.NewModel(mat.NewVecDense(1, []float64{1.}), - nil, func(y []float64, dom float64, x []float64) { - y[0] = tau * x[0] - }, nil) - if err != nil { - t.Fatal(err) - } - Stiff.Model = stiff - stiffsol, err := ivp.NewModel(mat.NewVecDense(1, []float64{0}), - nil, func(y []float64, dom float64, x []float64) { - y[0] = math.Exp(tau * dom) - }, nil) - if err != nil { - t.Fatal(err) - } - Stiff.solution = stiffsol - Stiff.err = func(h, i float64) float64 { return 2 * h * i } - return Stiff -} - -func TestNewton1DStiff(t *testing.T) { - stiff1D := exponential1DTestModel(t) - solver := ivp.NewNewtonRaphson( - ivp.ConfigScalarTolerance(0, 1e-5), - ) - solver.Set(0.0, stiff1D) - stepsize := 1. / 30. - results, err := ivp.Solve(solver, stepsize, 0.5) - if err != nil { - t.Fatal(err) - } - - quadsol := stiff1D.solution - sol := quadsol.Equations() - nxs, _ := quadsol.Dims() - solresults := make([]float64, nxs) - for i := range results { - sol(solresults, results[i].DomainStartOffset, results[i].X) - for j := range solresults { - got := math.Abs(solresults[j] - results[i].X[j]) - expect := stiff1D.err(stepsize, float64(i)) - if got > expect { - t.Errorf("error %g is greater than permitted tolerance %g", got, expect) - } - } - } -} diff --git a/ivp/rungekutta4.go b/ivp/rungekutta4.go deleted file mode 100644 index 4283980..0000000 --- a/ivp/rungekutta4.go +++ /dev/null @@ -1,67 +0,0 @@ -// Copyright ©2021 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 ivp - -import ( - "errors" - - "gonum.org/v1/gonum/floats" -) - -// RK4 implements Integrator interface for Runke-Kutta -// 4th order multivariable method (non adaptive) -type RK4 struct { - dom float64 - x []float64 - a, b, c, d []float64 - fx func(y []float64, t float64, x []float64) -} - -// Step implements Integrator interface -func (rk *RK4) Step(y []float64, h float64) (float64, error) { - const overSix = 1. / 6. - t := rk.dom - // set a, b, c, d (in some literatures these are k1,k2,k3,k4) - rk.fx(rk.a, t, rk.x) - rk.fx(rk.b, t+0.5*h, floats.AddScaledTo(rk.b, rk.x, 0.5*h, rk.a)) - rk.fx(rk.c, t+0.5*h, floats.AddScaledTo(rk.c, rk.x, 0.5*h, rk.b)) - rk.fx(rk.d, t+h, floats.AddScaledTo(rk.d, rk.x, h, rk.c)) - - floats.Add(rk.a, rk.d) - floats.Add(rk.b, rk.c) - floats.AddScaled(rk.a, 2, rk.b) - floats.AddScaledTo(y, rk.x, h*overSix, rk.a) - // finished integrating. Now we update RK4 structure for future Steps - copy(rk.x, y) // store results - rk.dom += h - return h, nil -} - -// Set implements integrator interface. All RK4 data -// is reset when calling Set. -func (rk *RK4) Set(d0 float64, ivp IVP) error { - if ivp == nil { - return errors.New("IVP is nil") - } - x0 := ivp.IV() - nx := x0.Len() - rk.dom = d0 //set start domain - rk.a = make([]float64, nx) - rk.b = make([]float64, nx) - rk.c = make([]float64, nx) - rk.d = make([]float64, nx) - rk.x = make([]float64, nx) - // set initial values - - for i := 0; i < x0.Len(); i++ { - rk.x[i] = x0.AtVec(i) - } - - rk.fx = ivp.Equations() - return nil -} - -// XLen returns length of x state vector -func (rk *RK4) XLen() (nx int) { return len(rk.x) } diff --git a/ivp/rungekutta4_test.go b/ivp/rungekutta4_test.go deleted file mode 100644 index bc38940..0000000 --- a/ivp/rungekutta4_test.go +++ /dev/null @@ -1,129 +0,0 @@ -// Copyright ©2021 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 ivp_test - -import ( - "fmt" - "log" - "math" - "testing" - - "gonum.org/v1/exp/ivp" - "gonum.org/v1/gonum/mat" -) - -type TestModel struct { - *ivp.Model - solution *ivp.Model - err func(h, i float64) float64 -} - -func TestQuadratic(t *testing.T) { - Quadratic := quadTestModel(t) - var solver ivp.Integrator = new(ivp.RK4) - - err := solver.Set(0, Quadratic.Model) - if err != nil { - log.Fatal(err) - } - - nx, _ := Quadratic.Model.Dims() - steps := 10 - dt := 0.1 - - results := make([]float64, nx) - - solmodel := Quadratic.solution - soleq := solmodel.Equations() - solDims, _ := solmodel.Dims() - solution := make([]float64, solDims) - for i := 1.; i < float64(steps+1); i++ { - dom := i * dt - solver.Step(results, dt) - soleq(solution, dom, results) - for j := range results { - got := math.Abs(solution[j] - results[j]) - expected := Quadratic.err(dt, i) - if got > expected { - t.Errorf("error %e greater than threshold %e", got, expected) - } - - } - } -} - -func Example_fallingBall() { - const ( - g = -10. // gravity field [m.s^-2] - ) - // we declare our physical model in the following function - ballModel, err := ivp.NewModel(mat.NewVecDense(2, []float64{100., 0.}), - nil, func(yvec []float64, _ float64, xvec []float64) { - // this anonymous function defines the physics. - // The first variable xvec[0] corresponds to position - // second variable xvec[1] is velocity. - Dx := xvec[1] - // yvec represents change in xvec, or derivative with respect to domain - // Change in position will be equal to velocity, which is the second variable: - // thus yvec[0] = xvec[1], which is the same as saying "change in xvec[0]" is equal to xvec[1] - yvec[0] = Dx - // change in velocity is acceleration. We suppose our ball is on earth accelerating at `g` - yvec[1] = g - }, nil) - if err != nil { - log.Fatal(err) - } - // Here we choose our algorithm. Runge-Kutta 4th order is used - var solver ivp.Integrator = new(ivp.RK4) - // Before integrating the IVP is passed to the integrator (a.k.a solver). Domain (time) starts at 0 - dom := 0. - err = solver.Set(dom, ballModel) - if err != nil { - log.Fatal(err) - } - - // we define the domain over which we integrate: 10 steps with step length of 0.1 - steps := 10 // number of steps - dt := 0.1 // step length - // we will store position in xvec and domain (time) in dvec - xvec := mat.NewVecDense(steps, nil) - dvec := mat.NewVecDense(steps, nil) - // results is an auxiliary vector that stores integration results - nx, _ := ballModel.Dims() - results := make([]float64, nx) - - for i := 0; i < steps; i++ { - // Step integrates the IVP. Each step advances the solution - step, _ := solver.Step(results, dt) // for non-adaptive algorithms step == dt - dom += step // calculate domain at current step - xvec.SetVec(i, results[0]) // set x value - dvec.SetVec(i, dom) - } - // print results - fmt.Println(mat.Formatted(xvec), "\n\n", mat.Formatted(dvec)) -} - -// Quadratic model may be used for future algorithms -func quadTestModel(t *testing.T) *TestModel { - Quadratic := new(TestModel) - quad, err := ivp.NewModel(mat.NewVecDense(2, []float64{0, 0}), - nil, func(y []float64, dom float64, x []float64) { - y[0], y[1] = x[1], 1. - }, nil) - if err != nil { - t.Fatal(err) - } - Quadratic.Model = quad - quadsol, err := ivp.NewModel(mat.NewVecDense(2, []float64{0, 0}), - nil, func(y []float64, dom float64, x []float64) { - y[0], y[1] = dom*dom/2., dom - }, nil) - if err != nil { - t.Fatal(err) - } - Quadratic.solution = quadsol - Quadratic.err = func(h, i float64) float64 { return math.Pow(h*i, 4) } - return Quadratic -} diff --git a/ivp/README.md b/ode/README.md similarity index 73% rename from ivp/README.md rename to ode/README.md index 1c501d0..68459d5 100644 --- a/ivp/README.md +++ b/ode/README.md @@ -1,3 +1,3 @@ # Gonum ivp [![GoDoc](https://godoc.org/gonum.org/v1/gonum/ivp?status.svg)](https://godoc.org/gonum.org/v1/gonum/ivp) -Package ivp provides numerical methods solving multivariable ODE initial value problems for the Go programming language. +package ode provides numerical methods solving multivariable ODE initial value problems for the Go programming language. diff --git a/ivp/config.go b/ode/config.go similarity index 99% rename from ivp/config.go rename to ode/config.go index e629dfb..1975bd6 100644 --- a/ivp/config.go +++ b/ode/config.go @@ -1,4 +1,4 @@ -package ivp +package ode import "errors" diff --git a/ivp/doc.go b/ode/doc.go similarity index 66% rename from ivp/doc.go rename to ode/doc.go index a336e44..b647ba7 100644 --- a/ivp/doc.go +++ b/ode/doc.go @@ -2,6 +2,6 @@ // Use of this source code is governed by a BSD-style // license that can be found in the LICENSE file. -// Package ivp provides interfaces and structures to solve +// package ode provides interfaces and structures to solve // multivariable ODE initial value problems. -package ivp // import "gonum.org/v1/gonum/ivp" +package ode // import "gonum.org/v1/gonum/ode" diff --git a/ivp/dopri5.go b/ode/dopri5.go similarity index 50% rename from ivp/dopri5.go rename to ode/dopri5.go index 332de2d..16089cd 100644 --- a/ivp/dopri5.go +++ b/ode/dopri5.go @@ -1,49 +1,50 @@ -package ivp +package ode import ( "math" - "gonum.org/v1/gonum/floats" + "gonum.org/v1/gonum/mat" ) type DoPri5 struct { param parameters adaptive bool - x []float64 + x, aux *mat.VecDense // solutions - y5, err45 []float64 + y5, err45 *mat.VecDense // integration coefficients (for Butcher Tableau) - k1, k2, k3, k4, k5, k6, k7 []float64 + k1, k2, k3, k4, k5, k6, k7 *mat.VecDense dom float64 - fx func(y []float64, t float64, x []float64) + fx func(y *mat.VecDense, t float64, x mat.Vector) } // Set implements the Integrator interface. Initializes a Dormand Prince method // for use as a solver -func (dp *DoPri5) Set(t0 float64, ivp IVP) error { +func (dp *DoPri5) Set(ivp IVP) { + + dp.fx = ivp.Func + + t0, x0 := ivp.IV() dp.dom = t0 - xequations := ivp.Equations() - dp.fx = xequations - // dp.fu = ufunc - nx := ivp.IV().Len() - dp.k1, dp.k2, dp.k3 = make([]float64, nx), make([]float64, nx), make([]float64, nx) - dp.k4, dp.k5, dp.k6 = make([]float64, nx), make([]float64, nx), make([]float64, nx) - dp.k7 = make([]float64, nx) - dp.x = make([]float64, nx) - for i := range dp.x { - dp.x[i] = ivp.IV().AtVec(i) - } - dp.y5 = make([]float64, nx) - dp.err45 = make([]float64, nx) - return nil + nx := x0.Len() + dp.k1, dp.k2, dp.k3 = mat.NewVecDense(nx, nil), mat.NewVecDense(nx, nil), mat.NewVecDense(nx, nil) + dp.k4, dp.k5, dp.k6 = mat.NewVecDense(nx, nil), mat.NewVecDense(nx, nil), mat.NewVecDense(nx, nil) + dp.k7 = mat.NewVecDense(nx, nil) + dp.aux = mat.NewVecDense(nx, nil) + dp.x = mat.NewVecDense(nx, nil) + dp.x.CloneFromVec(x0) + + dp.y5 = mat.NewVecDense(nx, nil) + dp.err45 = mat.NewVecDense(nx, nil) + } // 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(y4 []float64, h float64) (float64, error) { +func (dp *DoPri5) Step(y4 *mat.VecDense, h float64) (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 + const c40, c41, c42, c43 = 4. / 5., 44. / 45., -56. / 15., 32. / 9. const c50, c51, c52, c53, c54 = 8. / 9., 19372. / 6561., -25360. / 2187., 64448. / 6561., -212. / 729. const c60, c61, c62, c63, c64, c65 = 1., 9017. / 3168., -355. / 33., 46732. / 5247., 49. / 176., -5103. / 18656. const c70, c71, c72, c73, c74, c75, c76 = 1., 35. / 384., 0., 500. / 1113., 125. / 192., -2187. / 6784., 11. / 84. @@ -54,69 +55,75 @@ func (dp *DoPri5) Step(y4 []float64, h float64) (float64, error) { // prettier variable names, also someone once said pointers are equivalent to JMP x := dp.x + // auxiliary vector for storing results of calling ivp.Func + aux := dp.aux F := dp.fx t := dp.dom // if adaptive then this tag will be used until h==DoPri5.MinStep SOLVE: - F(dp.k1, t, x) - floats.Scale(h, dp.k1) - - floats.AddScaledTo(dp.k2, x, c21, dp.k1) - F(dp.k2, t+c20*h, x) - floats.Scale(h, dp.k2) - - floats.AddScaledTo(dp.k3, x, c31, dp.k1) - floats.AddScaled(dp.k3, c32, dp.k2) - F(dp.k3, t+c30*h, x) - floats.Scale(h, dp.k3) - - floats.AddScaledTo(dp.k4, x, c41, dp.k1) - floats.AddScaled(dp.k4, c42, dp.k2) - floats.AddScaled(dp.k4, c43, dp.k3) - F(dp.k4, t+c40*h, x) - floats.Scale(h, dp.k4) - - floats.AddScaledTo(dp.k5, x, c51, dp.k1) - floats.AddScaled(dp.k5, c52, dp.k2) - floats.AddScaled(dp.k5, c53, dp.k3) - floats.AddScaled(dp.k5, c54, dp.k4) - F(dp.k5, t+c50*h, x) - floats.Scale(h, dp.k5) - - floats.AddScaledTo(dp.k6, x, c61, dp.k1) - floats.AddScaled(dp.k6, c62, dp.k2) - floats.AddScaled(dp.k6, c63, dp.k3) - floats.AddScaled(dp.k6, c64, dp.k4) - floats.AddScaled(dp.k6, c65, dp.k5) - F(dp.k6, t+c60*h, x) - floats.Scale(h, dp.k6) - + F(aux, t, x) + dp.k1.ScaleVec(h, aux) + + dp.k2.AddScaledVec(x, c21, dp.k1) + F(aux, t+c20*h, dp.k2) + dp.k2.ScaleVec(h, aux) + + dp.k3.AddScaledVec(x, c31, dp.k1) + dp.k3.AddScaledVec(dp.k3, c32, dp.k2) + F(aux, t+h*c30, dp.k3) + dp.k3.ScaleVec(h, aux) + + dp.k4.AddScaledVec(x, c41, dp.k1) + dp.k4.AddScaledVec(dp.k4, c42, dp.k2) + dp.k4.AddScaledVec(dp.k4, c43, dp.k3) + F(aux, t+h*c40, dp.k4) + dp.k4.ScaleVec(h, aux) + + dp.k5.AddScaledVec(x, c51, dp.k1) + dp.k5.AddScaledVec(dp.k5, c52, dp.k2) + dp.k5.AddScaledVec(dp.k5, c53, dp.k3) + dp.k5.AddScaledVec(dp.k5, c54, dp.k4) + F(aux, t+h*c50, dp.k5) + dp.k5.ScaleVec(h, aux) + + dp.k6.AddScaledVec(x, c61, dp.k1) + dp.k6.AddScaledVec(dp.k6, c62, dp.k2) + dp.k6.AddScaledVec(dp.k6, c63, dp.k3) + dp.k6.AddScaledVec(dp.k6, c64, dp.k4) + dp.k6.AddScaledVec(dp.k6, c65, dp.k5) + F(aux, t+h*c60, dp.k6) + dp.k6.ScaleVec(h, aux) + + y4.AddScaledVec(x, b1, dp.k1) + y4.AddScaledVec(y4, b3, dp.k3) + y4.AddScaledVec(y4, b4, dp.k4) + y4.AddScaledVec(y4, b5, dp.k5) + y4.AddScaledVec(y4, b6, dp.k6) // fourth order approximation used to advance solution - floats.AddScaledTo(y4, x, b1, dp.k1) - floats.AddScaled(y4, b3, dp.k3) - floats.AddScaled(y4, b4, dp.k4) - floats.AddScaled(y4, b5, dp.k5) - floats.AddScaled(y4, b6, dp.k6) + if dp.adaptive { + y5 := dp.y5 + dp.k7.AddScaledVec(x, c71, dp.k1) + dp.k7.AddScaledVec(dp.k7, c72, dp.k2) + dp.k7.AddScaledVec(dp.k7, c73, dp.k3) + dp.k7.AddScaledVec(dp.k7, c74, dp.k4) + dp.k7.AddScaledVec(dp.k7, c75, dp.k5) + dp.k7.AddScaledVec(dp.k7, c76, dp.k6) // calculation of order 5 coefficients - floats.AddScaledTo(dp.k7, x, c71, dp.k1) - floats.AddScaled(dp.k7, c72, dp.k2) - floats.AddScaled(dp.k7, c73, dp.k3) - floats.AddScaled(dp.k7, c74, dp.k4) - floats.AddScaled(dp.k7, c75, dp.k5) - floats.AddScaled(dp.k7, c76, dp.k6) - F(dp.k7, t+h*c70, x) - floats.Scale(h, dp.k7) - - floats.AddScaledTo(dp.y5, x, a1, dp.k1) - floats.AddScaled(dp.y5, a3, dp.k3) - floats.AddScaled(dp.y5, a4, dp.k4) - floats.AddScaled(dp.y5, a5, dp.k5) - floats.AddScaled(dp.y5, a6, dp.k6) - floats.AddScaled(dp.y5, a7, dp.k7) + + F(aux, t+h*c70, dp.k7) + dp.k7.ScaleVec(h, aux) + + y5.AddScaledVec(x, a1, dp.k1) + y5.AddScaledVec(y5, a3, dp.k3) + y5.AddScaledVec(y5, a4, dp.k4) + y5.AddScaledVec(y5, a5, dp.k5) + y5.AddScaledVec(y5, a6, dp.k6) + y5.AddScaledVec(y5, a7, dp.k7) + dp.err45.SubVec(y4, y5) // compute error between fifth order solution and fourth order solution - errRatio := dp.param.atol / floats.Norm(floats.SubTo(dp.err45, y4, dp.y5), math.Inf(1)) + errRatio := dp.param.atol / mat.Norm(dp.err45, math.Inf(1)) hnew := math.Min(math.Max(0.9*h*math.Pow(errRatio, .2), dp.param.minStep), dp.param.maxStep) // given "bad" error ratio, algorithm will recompute with a smaller step, given it is not the minimum permissible step or less if errRatio < 1 && h > dp.param.minStep { @@ -125,13 +132,11 @@ SOLVE: } } // advance solution with fourth order solution - copy(dp.x, y4) + dp.x.CloneFromVec(y4) + dp.dom += h return h, nil } -// XLen returns length of x state vector -func (dp *DoPri5) XLen() (nx int) { return len(dp.x) } - // NewDormandPrince5 returns a adaptive-ready Dormand Prince solver of order 5 // // To enable step size adaptivity minimum step size must be set and diff --git a/ivp/ivp.go b/ode/ivp.go similarity index 72% rename from ivp/ivp.go rename to ode/ivp.go index 11b49e3..4ac174b 100644 --- a/ivp/ivp.go +++ b/ode/ivp.go @@ -2,11 +2,11 @@ // Use of this source code is governed by a BSD-style // license that can be found in the LICENSE file. -package ivp +package ode import ( "errors" - "fmt" + "math" "gonum.org/v1/gonum/mat" ) @@ -32,13 +32,13 @@ import ( type IVP interface { // Initial values vector for state vector x. x0 defines // the first values the state vector takes when integration begins. - IV() (x0 mat.Vector) - // Equations returns the coupled, non-linear algebraic differential - // equations (xequations) for the state vector (x) - // Results are stored in y which are the length of x. + IV() (t0 float64, x0 mat.Vector) + // Func is the coupled, non-linear algebraic differential + // equations for the state vector (x) + // Results are stored in dst which are the length of x. // The scalar (float64) argument is the domain over which the // problem is integrated, which is usually time for most physical problems. - Equations() (xequations func(y []float64, t float64, x []float64)) + Func(dst *mat.VecDense, t float64, x mat.Vector) } // Integrator abstracts algorithm specifics. For anyone looking to @@ -50,38 +50,39 @@ type IVP interface { type Integrator interface { // Set initializes an initial value problem. First argument // is the initial domain integration point, is usually zero. - Set(float64, IVP) error - // Length of state vector x - XLen() (nx int) - // Step integrates IVP and stores result in y. step is a suggested step + Set(IVP) + // Step integrates IVP and stores result in dst. step is a suggested step // for the algorithm to take. The algorithm may decide that it is not sufficiently // small or big enough (these are adaptive algorithms) and take a different step. // The resulting step is returned as the first parameter - Step(y []float64, step float64) (float64, error) + Step(dst *mat.VecDense, step float64) (float64, error) } type result = struct { - DomainStartOffset float64 - X []float64 + T float64 + X *mat.VecDense } // Solve solves an already initialized Integrator returning state vector results. -// Returns error upon needing to allocate over 2GB of memory -func Solve(solver Integrator, stepsize, domainLength float64) (results []result, err error) { - const maxAllocGB = 2 - integrated := 0. - expectedLength := int(domainLength/stepsize) + 1 - nx := solver.XLen() +func Solve(p IVP, solver Integrator, stepsize, tend float64) (results []result, err error) { + t0, x0 := p.IV() + nx := x0.Len() if nx == 0 { return nil, errors.New("state vector length can not be equal to zero. Has ivp been set?") } - size := 8 * (nx + 1) * expectedLength - if size > maxAllocGB*1e9 { - return nil, fmt.Errorf("solution exceeds %dGB or not initialized (size is %dMB)", maxAllocGB, size/1e6) - } + // calculate expected size of results, these may differ + domainLength := tend - t0 + expectedLength := int(domainLength/stepsize) + 1 results = make([]result, 0, expectedLength) - for integrated < domainLength { - res := make([]float64, nx) + + // t stores current domain + t := t0 + for t < tend { + res := mat.NewVecDense(nx, nil) + + if t-tend > 1e-10 { + stepsize = math.Min(stepsize, (t-tend)*(1+1e-3)) + } stepsize, err = solver.Step(res, stepsize) if err != nil { return results, err @@ -89,8 +90,8 @@ func Solve(solver Integrator, stepsize, domainLength float64) (results []result, if stepsize <= 0 { return results, errors.New("got zero or negative step size from Integrator") } - integrated += stepsize - results = append(results, result{DomainStartOffset: integrated, X: res}) + t += stepsize + results = append(results, result{T: t, X: res}) } return results, nil } diff --git a/ode/ivp_test.go b/ode/ivp_test.go new file mode 100644 index 0000000..21d564f --- /dev/null +++ b/ode/ivp_test.go @@ -0,0 +1,172 @@ +// Copyright ©2021 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_test + +import ( + "fmt" + "log" + "math" + "testing" + + "gonum.org/v1/exp/ode" + "gonum.org/v1/gonum/mat" +) + +type TestModel struct { + *ode.Model + solution *ode.Model + err func(h, i float64) float64 +} + +func TestSolve(t *testing.T) { + // domain start + quad := quadTestModel(t) + solver := ode.NewDormandPrince5() + solver.Set(quad) + stepsize := 0.5 / 15. + end := 0.5 + results, err := ode.Solve(quad, solver, stepsize, end) + if err != nil { + t.Fatal(err) + } + // run the solver and corsscheck values + quadsol := quad.solution + _, x0 := quadsol.IV() + + sol := quadsol.Func + + solresults := mat.NewVecDense(x0.Len(), nil) + for i := range results { + sol(solresults, results[i].T, results[i].X) + for j := range solresults.RawVector().Data { + sol := solresults.AtVec(j) + res := results[i].X.AtVec(j) + got := math.Abs(solresults.AtVec(j) - results[i].X.AtVec(j)) + expect := quad.err(stepsize, float64(i)) + if got > expect { + t.Errorf("error %g is greater than permitted tolerance %g. solution:[%0.3g], result:[%0.3g]", got, expect, sol, res) + } + } + } + +} + +func Example_solve() { + const ( + g = -10. // gravity field [m.s^-2] + ) + // we declare our physical model. First argument is initial time, which is 0 seconds. + // Next is the initial state vector, which corresponds to 100 meters above the ground + // with 0 m/s velocity. + ballModel, err := ode.NewModel(0, mat.NewVecDense(2, []float64{100., 0.}), + func(yvec *mat.VecDense, _ float64, xvec mat.Vector) { + // this anonymous function defines the physics. + // The first variable xvec[0] corresponds to position + // second variable xvec[1] is velocity. + Dx := xvec.AtVec(1) + // yvec represents change in xvec, or derivative with respect to domain + // Change in position will be equal to velocity, which is the second variable: + // thus yvec[0] = xvec[1], which is the same as saying "change in xvec[0]" is equal to xvec[1] + yvec.SetVec(0, Dx) + // change in velocity is acceleration. We suppose our ball is on earth accelerating at `g` + yvec.SetVec(1, g) + }) + if err != nil { + log.Fatal(err) + } + // Here we choose our algorithm. Runge-Kutta 4th order is used + var solver ode.Integrator = ode.NewDormandPrince5() + // 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. + results, err := ode.Solve(ballModel, solver, 0.1, 10.) + fmt.Println(results) +} + +// Quadratic model may be used for future algorithms +func quadTestModel(t *testing.T) *TestModel { + t0 := 0.0 + Quadratic := new(TestModel) + quad, err := ode.NewModel(t0, mat.NewVecDense(2, []float64{0, 0}), + func(dst *mat.VecDense, t float64, x mat.Vector) { + dst.SetVec(0, x.AtVec(1)) + dst.SetVec(1, 1.) + }) + if err != nil { + t.Fatal(err) + } + Quadratic.Model = quad + quadsol, err := ode.NewModel(t0, mat.NewVecDense(2, []float64{0, 0}), + func(dst *mat.VecDense, t float64, x mat.Vector) { + dst.SetVec(0, t*t/2.) + dst.SetVec(1, t) + }) + if err != nil { + t.Fatal(err) + } + Quadratic.solution = quadsol + Quadratic.err = func(h, i float64) float64 { return math.Pow(h*i, 4) + 1e-10 } + return Quadratic +} + +// exponential unidimensional model may be used for future algorithms +// 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 + Quadratic := new(TestModel) + quad, err := ode.NewModel(t0, mat.NewVecDense(1, []float64{1.}), + func(dst *mat.VecDense, t float64, x mat.Vector) { + dst.SetVec(0, tau*x.AtVec(0)) + }) + if err != nil { + t.Fatal(err) + } + Quadratic.Model = quad + quadsol, err := ode.NewModel(t0, mat.NewVecDense(1, []float64{0}), + func(dst *mat.VecDense, t float64, x mat.Vector) { + dst.SetVec(0, math.Exp(tau*t)) + + }) + if err != nil { + t.Fatal(err) + } + Quadratic.solution = quadsol + Quadratic.err = func(h, i float64) float64 { return math.Pow(h*i, 4) + 1e-10 } + return Quadratic +} + +// func TestQuadratic(t *testing.T) { +// Quadratic := quadTestModel(t) +// solver := ivp.NewDormandPrince5() + +// solver.Set(Quadratic.Model) + +// _, x0 := Quadratic.Model.IV() +// steps := 10 +// dt := 0.1 + +// results := make([]float64, nx) + +// solmodel := Quadratic.solution +// soleq := solmodel.Equations() +// solDims, _ := solmodel.Dims() +// solution := make([]float64, solDims) +// for i := 1.; i < float64(steps+1); i++ { +// dom := i * dt +// solver.Step(results, dt) +// soleq(solution, dom, results) +// for j := range results { +// got := math.Abs(solution[j] - results[j]) +// expected := Quadratic.err(dt, i) +// if got > expected { +// t.Errorf("error %e greater than threshold %e", got, expected) +// } + +// } +// } +// } diff --git a/ode/model.go b/ode/model.go new file mode 100644 index 0000000..1208588 --- /dev/null +++ b/ode/model.go @@ -0,0 +1,41 @@ +// Copyright ©2021 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 ( + "errors" + "math" + + "gonum.org/v1/gonum/mat" +) + +// Model implements IVP interface. X vector and equations can not be nil or zero length. +type Model struct { + x0 mat.Vector + t0 float64 + xeq func(dst *mat.VecDense, dom float64, x mat.Vector) +} + +// IV returns initial values of the IVP. First returned parameter is the +// starting x vector and second parameter are inputs when solving non-autonomous +// ODEs. +func (m *Model) IV() (t0 float64, x0 mat.Vector) { return m.t0, m.x0 } + +// Equations returns differential equations relating to state vector x and input functions +// for non-autonomous ODEs. +// +// Input functions may be nil (ueq). +func (m *Model) Func(y *mat.VecDense, dom float64, x mat.Vector) { + m.xeq(y, dom, x) +} + +// NewModel returns a IVP given initial conditions (x0,u0), differential equations (xeq) and +// input functions for non-autonomous ODEs (ueq). +func NewModel(t0 float64, x0 mat.Vector, xeq func(y *mat.VecDense, dom float64, x mat.Vector)) (*Model, error) { + if x0 == nil || math.IsNaN(t0) || xeq == nil { + return nil, errors.New("bad model value") + } + return &Model{xeq: xeq, x0: x0, t0: t0}, nil +} From ef1f6e432158a64b22d35d48bedb11f81d0ccf18 Mon Sep 17 00:00:00 2001 From: soypat Date: Wed, 2 Jun 2021 23:15:12 -0300 Subject: [PATCH 08/10] api changes. kudos to @vladimir-ch --- ode/dopri5.go | 16 ++++++--- ode/ivp.go | 96 +++++++++++++------------------------------------ ode/ivp_test.go | 31 ++++++++-------- ode/model.go | 41 --------------------- ode/ode.go | 65 +++++++++++++++++++++++++++++++++ 5 files changed, 117 insertions(+), 132 deletions(-) delete mode 100644 ode/model.go create mode 100644 ode/ode.go diff --git a/ode/dopri5.go b/ode/dopri5.go index 16089cd..07724b2 100644 --- a/ode/dopri5.go +++ b/ode/dopri5.go @@ -11,7 +11,7 @@ type DoPri5 struct { adaptive bool x, aux *mat.VecDense // solutions - y5, err45 *mat.VecDense + y4, y5, err45 *mat.VecDense // integration coefficients (for Butcher Tableau) k1, k2, k3, k4, k5, k6, k7 *mat.VecDense dom float64 @@ -20,11 +20,11 @@ type DoPri5 struct { // Set implements the Integrator interface. Initializes a Dormand Prince method // for use as a solver -func (dp *DoPri5) Set(ivp IVP) { +func (dp *DoPri5) Init(ivp IVP) { dp.fx = ivp.Func - t0, x0 := ivp.IV() + t0, x0 := ivp.T0, mat.VecDenseCopyOf(ivp.Y0) dp.dom = t0 nx := x0.Len() dp.k1, dp.k2, dp.k3 = mat.NewVecDense(nx, nil), mat.NewVecDense(nx, nil), mat.NewVecDense(nx, nil) @@ -34,14 +34,20 @@ func (dp *DoPri5) Set(ivp IVP) { dp.x = mat.NewVecDense(nx, nil) dp.x.CloneFromVec(x0) + dp.y4 = mat.NewVecDense(nx, nil) dp.y5 = mat.NewVecDense(nx, nil) dp.err45 = mat.NewVecDense(nx, nil) } +func (dp *DoPri5) State(s *State) { + s.T = dp.dom + s.Y.CloneFromVec(dp.x) +} + // 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(y4 *mat.VecDense, h float64) (float64, error) { +func (dp *DoPri5) Step(h float64) (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. @@ -59,6 +65,7 @@ func (dp *DoPri5) Step(y4 *mat.VecDense, h float64) (float64, error) { aux := dp.aux F := dp.fx t := dp.dom + y4 := dp.y4 // if adaptive then this tag will be used until h==DoPri5.MinStep SOLVE: F(aux, t, x) @@ -94,6 +101,7 @@ SOLVE: F(aux, t+h*c60, dp.k6) dp.k6.ScaleVec(h, aux) + // TODO y4 could be dp.x, is there a disadvantage to this change? y4.AddScaledVec(x, b1, dp.k1) y4.AddScaledVec(y4, b3, dp.k3) y4.AddScaledVec(y4, b4, dp.k4) diff --git a/ode/ivp.go b/ode/ivp.go index 4ac174b..47689bc 100644 --- a/ode/ivp.go +++ b/ode/ivp.go @@ -11,87 +11,39 @@ import ( "gonum.org/v1/gonum/mat" ) -// IVP defines a multivariable, non-autonomous initial value problem. -// It is worth mentioning the system need not be non-autonomous. https://en.wikipedia.org/wiki/Autonomous_system_(mathematics) +// IVP defines a multivariable, initial value problem represented by a system of ordinary differential equations. // -// These problems have the form (capital letters are vectors in this example) -// X'(t) = F(t, X) -// X(0) = X_0 +// These problems have the form +// y'(t) = f(t, y(t)) +// y(0) = y_0 // // Where: // t is a scalar representing the integration domain, which is time for most physical problems. -// X' is the vector of first derivatives of the state vector X. -// F would be xequations as returned by Equations(). +// y is the state vector. +// y' is the vector of first derivatives of the state vector y. +// f are the differential equations represented by Func. // An initial value problem is characterized by the initial conditions imposed -// on the state vector X at the beginning of the integration domain. These +// on the state vector y at the beginning of the integration domain. These // initial conditions are returned by the IV() method for the state vector -// as x0 and for the input vector as u0. +// as y0. // // The term "state vector" and "state variables" are used interchangeably -// throughout the code and refer to X vector of independent variables. -type IVP interface { - // Initial values vector for state vector x. x0 defines - // the first values the state vector takes when integration begins. - IV() (t0 float64, x0 mat.Vector) - // Func is the coupled, non-linear algebraic differential - // equations for the state vector (x) - // Results are stored in dst which are the length of x. - // The scalar (float64) argument is the domain over which the - // problem is integrated, which is usually time for most physical problems. - Func(dst *mat.VecDense, t float64, x mat.Vector) +// throughout the code and refer to y vector of independent variables. +type IVP struct { + // Initial values for the state vector + Y0 mat.Vector + // Independent variable point at which Y0 is evaluated + T0 float64 + // Func are the differential equations f(t,y(t)). + // The result is y'(t) and is stored in dst. + Func func(dst *mat.VecDense, t float64, y mat.Vector) } -// Integrator abstracts algorithm specifics. For anyone looking to -// implement it, Set(ivp) should be called first to initialize the IVP with -// initial values. Step will calculate the next x values and store them in y -// u values should not be stored as they can easily be obtained if one has -// x values. Integrator should store 1 or more (depending on algorithm used) -// of previously calculated x values to be able to integrate. -type Integrator interface { - // Set initializes an initial value problem. First argument - // is the initial domain integration point, is usually zero. - Set(IVP) - // Step integrates IVP and stores result in dst. step is a suggested step - // for the algorithm to take. The algorithm may decide that it is not sufficiently - // small or big enough (these are adaptive algorithms) and take a different step. - // The resulting step is returned as the first parameter - Step(dst *mat.VecDense, step float64) (float64, error) -} - -type result = struct { - T float64 - X *mat.VecDense -} - -// Solve solves an already initialized Integrator returning state vector results. -func Solve(p IVP, solver Integrator, stepsize, tend float64) (results []result, err error) { - t0, x0 := p.IV() - nx := x0.Len() - if nx == 0 { - return nil, errors.New("state vector length can not be equal to zero. Has ivp been set?") - } - // calculate expected size of results, these may differ - domainLength := tend - t0 - expectedLength := int(domainLength/stepsize) + 1 - results = make([]result, 0, expectedLength) - - // t stores current domain - t := t0 - for t < tend { - res := mat.NewVecDense(nx, nil) - - if t-tend > 1e-10 { - stepsize = math.Min(stepsize, (t-tend)*(1+1e-3)) - } - stepsize, err = solver.Step(res, stepsize) - if err != nil { - return results, err - } - if stepsize <= 0 { - return results, errors.New("got zero or negative step size from Integrator") - } - t += stepsize - results = append(results, result{T: t, X: res}) +// 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 { + return IVP{}, errors.New("bad model value") } - return results, nil + return IVP{Func: f, Y0: y0, T0: t0}, nil } diff --git a/ode/ivp_test.go b/ode/ivp_test.go index 21d564f..15de75f 100644 --- a/ode/ivp_test.go +++ b/ode/ivp_test.go @@ -15,35 +15,36 @@ import ( ) type TestModel struct { - *ode.Model - solution *ode.Model + ode.IVP + solution ode.IVP err func(h, i float64) float64 } func TestSolve(t *testing.T) { // domain start quad := quadTestModel(t) + ivp := quad.IVP solver := ode.NewDormandPrince5() - solver.Set(quad) + solver.Init(ivp) stepsize := 0.5 / 15. end := 0.5 - results, err := ode.Solve(quad, solver, stepsize, end) + results, err := ode.Solve(ivp, solver, stepsize, end) if err != nil { t.Fatal(err) } // run the solver and corsscheck values quadsol := quad.solution - _, x0 := quadsol.IV() + x0 := quadsol.Y0 sol := quadsol.Func solresults := mat.NewVecDense(x0.Len(), nil) for i := range results { - sol(solresults, results[i].T, results[i].X) + sol(solresults, results[i].T, results[i].Y) for j := range solresults.RawVector().Data { sol := solresults.AtVec(j) - res := results[i].X.AtVec(j) - got := math.Abs(solresults.AtVec(j) - results[i].X.AtVec(j)) + res := results[i].Y.AtVec(j) + got := math.Abs(solresults.AtVec(j) - results[i].Y.AtVec(j)) expect := quad.err(stepsize, float64(i)) if got > expect { t.Errorf("error %g is greater than permitted tolerance %g. solution:[%0.3g], result:[%0.3g]", got, expect, sol, res) @@ -60,7 +61,7 @@ func Example_solve() { // we declare our physical model. First argument is initial time, which is 0 seconds. // Next is the initial state vector, which corresponds to 100 meters above the ground // with 0 m/s velocity. - ballModel, err := ode.NewModel(0, mat.NewVecDense(2, []float64{100., 0.}), + ballModel, err := ode.NewIVP(0, mat.NewVecDense(2, []float64{100., 0.}), func(yvec *mat.VecDense, _ float64, xvec mat.Vector) { // this anonymous function defines the physics. // The first variable xvec[0] corresponds to position @@ -89,7 +90,7 @@ func Example_solve() { func quadTestModel(t *testing.T) *TestModel { t0 := 0.0 Quadratic := new(TestModel) - quad, err := ode.NewModel(t0, mat.NewVecDense(2, []float64{0, 0}), + quad, err := ode.NewIVP(t0, mat.NewVecDense(2, []float64{0, 0}), func(dst *mat.VecDense, t float64, x mat.Vector) { dst.SetVec(0, x.AtVec(1)) dst.SetVec(1, 1.) @@ -97,8 +98,8 @@ func quadTestModel(t *testing.T) *TestModel { if err != nil { t.Fatal(err) } - Quadratic.Model = quad - quadsol, err := ode.NewModel(t0, mat.NewVecDense(2, []float64{0, 0}), + Quadratic.IVP = quad + quadsol, err := ode.NewIVP(t0, mat.NewVecDense(2, []float64{0, 0}), func(dst *mat.VecDense, t float64, x mat.Vector) { dst.SetVec(0, t*t/2.) dst.SetVec(1, t) @@ -119,15 +120,15 @@ func exp1DTestModel(t *testing.T) *TestModel { tau := -2. t0 := 0.0 Quadratic := new(TestModel) - quad, err := ode.NewModel(t0, mat.NewVecDense(1, []float64{1.}), + quad, err := ode.NewIVP(t0, mat.NewVecDense(1, []float64{1.}), func(dst *mat.VecDense, t float64, x mat.Vector) { dst.SetVec(0, tau*x.AtVec(0)) }) if err != nil { t.Fatal(err) } - Quadratic.Model = quad - quadsol, err := ode.NewModel(t0, mat.NewVecDense(1, []float64{0}), + Quadratic.IVP = quad + quadsol, err := ode.NewIVP(t0, mat.NewVecDense(1, []float64{0}), func(dst *mat.VecDense, t float64, x mat.Vector) { dst.SetVec(0, math.Exp(tau*t)) diff --git a/ode/model.go b/ode/model.go deleted file mode 100644 index 1208588..0000000 --- a/ode/model.go +++ /dev/null @@ -1,41 +0,0 @@ -// Copyright ©2021 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 ( - "errors" - "math" - - "gonum.org/v1/gonum/mat" -) - -// Model implements IVP interface. X vector and equations can not be nil or zero length. -type Model struct { - x0 mat.Vector - t0 float64 - xeq func(dst *mat.VecDense, dom float64, x mat.Vector) -} - -// IV returns initial values of the IVP. First returned parameter is the -// starting x vector and second parameter are inputs when solving non-autonomous -// ODEs. -func (m *Model) IV() (t0 float64, x0 mat.Vector) { return m.t0, m.x0 } - -// Equations returns differential equations relating to state vector x and input functions -// for non-autonomous ODEs. -// -// Input functions may be nil (ueq). -func (m *Model) Func(y *mat.VecDense, dom float64, x mat.Vector) { - m.xeq(y, dom, x) -} - -// NewModel returns a IVP given initial conditions (x0,u0), differential equations (xeq) and -// input functions for non-autonomous ODEs (ueq). -func NewModel(t0 float64, x0 mat.Vector, xeq func(y *mat.VecDense, dom float64, x mat.Vector)) (*Model, error) { - if x0 == nil || math.IsNaN(t0) || xeq == nil { - return nil, errors.New("bad model value") - } - return &Model{xeq: xeq, x0: x0, t0: t0}, nil -} diff --git a/ode/ode.go b/ode/ode.go new file mode 100644 index 0000000..7d6da5e --- /dev/null +++ b/ode/ode.go @@ -0,0 +1,65 @@ +// Copyright ©2021 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 ( + "errors" + "math" + + "gonum.org/v1/gonum/mat" +) + +// Integrator can integrate an initial-value problem (IVP) for a first-order +// system of ordinary differential equations (ODEs). +type Integrator interface { + // Init initializes the integrator and sets the initial condition. + 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) + + // State stores the current state of the integrator in-place in dst. + State(dst *State) +} + +type State struct { + T float64 + Y *mat.VecDense +} + +// Solve solves an already initialized Integrator returning state vector results. +func Solve(p IVP, solver Integrator, stepsize, tend float64) (results []State, err error) { + t0, x0 := p.T0, mat.VecDenseCopyOf(p.Y0) + nx := x0.Len() + if nx == 0 { + return nil, errors.New("state vector length can not be equal to zero. Has ivp been set?") + } + // calculate expected size of results, these may differ + domainLength := tend - t0 + expectedLength := int(domainLength/stepsize) + 1 + results = make([]State, 0, expectedLength) + + // t stores current domain + t := t0 + for t < tend { + res := State{Y: mat.NewVecDense(nx, nil)} + + if t-tend > 1e-10 { + stepsize = math.Min(stepsize, (t-tend)*(1+1e-3)) + } + stepsize, err = solver.Step(stepsize) + if err != nil { + return results, err + } + if stepsize <= 0 { + return results, errors.New("got zero or negative step size from Integrator") + } + solver.State(&res) + t += stepsize + results = append(results, res) + } + return results, nil +} From 24e818fbd986025563c7a6500bf8a07224f94c3a Mon Sep 17 00:00:00 2001 From: soypat Date: Thu, 3 Jun 2021 14:34:06 -0300 Subject: [PATCH 09/10] rename Solve->SolveIVP --- ode/ivp_test.go | 4 ++-- ode/ode.go | 4 ++-- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/ode/ivp_test.go b/ode/ivp_test.go index 15de75f..58d69b3 100644 --- a/ode/ivp_test.go +++ b/ode/ivp_test.go @@ -28,7 +28,7 @@ func TestSolve(t *testing.T) { solver.Init(ivp) stepsize := 0.5 / 15. end := 0.5 - results, err := ode.Solve(ivp, solver, stepsize, end) + results, err := ode.SolveIVP(ivp, solver, stepsize, end) if err != nil { t.Fatal(err) } @@ -82,7 +82,7 @@ func Example_solve() { // 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. - results, err := ode.Solve(ballModel, solver, 0.1, 10.) + results, err := ode.SolveIVP(ballModel, solver, 0.1, 10.) fmt.Println(results) } diff --git a/ode/ode.go b/ode/ode.go index 7d6da5e..953382f 100644 --- a/ode/ode.go +++ b/ode/ode.go @@ -30,8 +30,8 @@ type State struct { Y *mat.VecDense } -// Solve solves an already initialized Integrator returning state vector results. -func Solve(p IVP, solver Integrator, stepsize, tend float64) (results []State, err error) { +// SolveIVP solves an already initialized Integrator returning state vector results. +func SolveIVP(p IVP, solver Integrator, stepsize, tend float64) (results []State, err error) { t0, x0 := p.T0, mat.VecDenseCopyOf(p.Y0) nx := x0.Len() if nx == 0 { From 5a7413e9fd9ff849913db334b5c7e0e4baf3fab3 Mon Sep 17 00:00:00 2001 From: soypat Date: Fri, 18 Jun 2021 22:29:48 -0300 Subject: [PATCH 10/10] add first posit implementation --- ode/README.md | 3 - ode/config.go | 80 ------------------- ode/dopri5.go | 169 ----------------------------------------- ode/ivp.go | 49 ------------ ode/ivp_test.go | 173 ------------------------------------------ ode/ode.go | 65 ---------------- {ode => posit}/doc.go | 6 +- posit/posit.go | 143 ++++++++++++++++++++++++++++++++++ posit/posit_test.go | 121 +++++++++++++++++++++++++++++ 9 files changed, 267 insertions(+), 542 deletions(-) delete mode 100644 ode/README.md delete mode 100644 ode/config.go delete mode 100644 ode/dopri5.go delete mode 100644 ode/ivp.go delete mode 100644 ode/ivp_test.go delete mode 100644 ode/ode.go rename {ode => posit}/doc.go (52%) create mode 100644 posit/posit.go create mode 100644 posit/posit_test.go diff --git a/ode/README.md b/ode/README.md deleted file mode 100644 index 68459d5..0000000 --- a/ode/README.md +++ /dev/null @@ -1,3 +0,0 @@ -# Gonum ivp [![GoDoc](https://godoc.org/gonum.org/v1/gonum/ivp?status.svg)](https://godoc.org/gonum.org/v1/gonum/ivp) - -package ode provides numerical methods solving multivariable ODE initial value problems for the Go programming language. diff --git a/ode/config.go b/ode/config.go deleted file mode 100644 index 1975bd6..0000000 --- a/ode/config.go +++ /dev/null @@ -1,80 +0,0 @@ -package ode - -import "errors" - -// Configuration represents a change in default parameters in the solving -// of IVPs. Parameters are often needed to enable step size adaptivity and to -// overwrite solver default values. -type Configuration func(*parameters) error - -type parameters struct { - // step length - minStep, maxStep float64 - // Relative and absolute error tolerances - rtol, atol float64 - // parameters for step size selection. - // The new step size is subject to restriction minStepRatio <= hnew/hold <= maxStepRatio - minStepRatio, maxStepRatio float64 - // relaxation factor offers stability for implicit solvers - relax float64 -} - -// ConfigStepLimits sets minimum step length and max step length -// an integrator can take when using an adaptive method. -func ConfigStepLimits(minStep, maxStep float64) Configuration { - return func(p *parameters) error { - if minStep <= 0 || minStep > maxStep { - return errors.New("minimum step length too small or greater than max step length") - } - p.maxStep, p.minStep = maxStep, minStep - return nil - } -} - -// ConfigScalarTolerance Sets relative and absolute tolerances for controlling step error. -// If value passed is zero the corresponding tolerance is not modified. -func ConfigScalarTolerance(rtol, atol float64) Configuration { - return func(p *parameters) error { - if rtol < 0 || atol < 0 { - return errors.New("negative error tolerance") - } - if rtol > 0 { - p.rtol = rtol - } - if atol > 0 { - p.atol = atol - } - return nil - } -} - -// ConfigAdaptiveRatioLimits sets parameters for step size selection. -// The new step size is subject to restriction min <= hnew/hold <= max. -// If min or max are zero then value is not changed -func ConfigAdaptiveRatioLimits(min, max float64) Configuration { - return func(p *parameters) error { - if min < 0 || max < min { - return errors.New("bad step size selection ratio") - } - if min > 0 { - p.minStepRatio = min - } - if max > 0 { - p.maxStepRatio = max - } - return nil - } -} - -// ConfigRelaxation sets the relaxation factor for methods -// like Newton-Raphson which can be unstable in certain cases -// without it. -func ConfigRelaxation(alpha float64) Configuration { - return func(p *parameters) error { - if alpha < 0 || alpha >= 1. { - return errors.New("relaxation factor must be in [0,1)") - } - p.relax = alpha - return nil - } -} diff --git a/ode/dopri5.go b/ode/dopri5.go deleted file mode 100644 index 07724b2..0000000 --- a/ode/dopri5.go +++ /dev/null @@ -1,169 +0,0 @@ -package ode - -import ( - "math" - - "gonum.org/v1/gonum/mat" -) - -type DoPri5 struct { - param parameters - adaptive bool - x, aux *mat.VecDense - // solutions - y4, y5, err45 *mat.VecDense - // integration coefficients (for Butcher Tableau) - k1, k2, k3, k4, k5, k6, k7 *mat.VecDense - dom float64 - fx func(y *mat.VecDense, t float64, x mat.Vector) -} - -// Set implements the Integrator interface. Initializes a Dormand Prince method -// for use as a solver -func (dp *DoPri5) Init(ivp IVP) { - - dp.fx = ivp.Func - - t0, x0 := ivp.T0, mat.VecDenseCopyOf(ivp.Y0) - dp.dom = t0 - nx := x0.Len() - dp.k1, dp.k2, dp.k3 = mat.NewVecDense(nx, nil), mat.NewVecDense(nx, nil), mat.NewVecDense(nx, nil) - dp.k4, dp.k5, dp.k6 = mat.NewVecDense(nx, nil), mat.NewVecDense(nx, nil), mat.NewVecDense(nx, nil) - dp.k7 = mat.NewVecDense(nx, nil) - dp.aux = mat.NewVecDense(nx, nil) - dp.x = mat.NewVecDense(nx, nil) - dp.x.CloneFromVec(x0) - - dp.y4 = mat.NewVecDense(nx, nil) - dp.y5 = mat.NewVecDense(nx, nil) - dp.err45 = mat.NewVecDense(nx, nil) - -} - -func (dp *DoPri5) State(s *State) { - s.T = dp.dom - s.Y.CloneFromVec(dp.x) -} - -// 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) { - 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. - const c50, c51, c52, c53, c54 = 8. / 9., 19372. / 6561., -25360. / 2187., 64448. / 6561., -212. / 729. - const c60, c61, c62, c63, c64, c65 = 1., 9017. / 3168., -355. / 33., 46732. / 5247., 49. / 176., -5103. / 18656. - const c70, c71, c72, c73, c74, c75, c76 = 1., 35. / 384., 0., 500. / 1113., 125. / 192., -2187. / 6784., 11. / 84. - // Alternate solution for error calculation - const a1, a3, a4, a5, a6, a7 = 5179. / 57600., 7571. / 16695., 393. / 640., -92097. / 339200., 187. / 2100., 1. / 40. - // Fourth order solution is used to advance - const b1, b3, b4, b5, b6 = 35. / 384., 500. / 1113., 125. / 192., -2187. / 6784., 11. / 84. - - // prettier variable names, also someone once said pointers are equivalent to JMP - x := dp.x - // auxiliary vector for storing results of calling ivp.Func - aux := dp.aux - F := dp.fx - t := dp.dom - y4 := dp.y4 - // if adaptive then this tag will be used until h==DoPri5.MinStep -SOLVE: - F(aux, t, x) - dp.k1.ScaleVec(h, aux) - - dp.k2.AddScaledVec(x, c21, dp.k1) - F(aux, t+c20*h, dp.k2) - dp.k2.ScaleVec(h, aux) - - dp.k3.AddScaledVec(x, c31, dp.k1) - dp.k3.AddScaledVec(dp.k3, c32, dp.k2) - F(aux, t+h*c30, dp.k3) - dp.k3.ScaleVec(h, aux) - - dp.k4.AddScaledVec(x, c41, dp.k1) - dp.k4.AddScaledVec(dp.k4, c42, dp.k2) - dp.k4.AddScaledVec(dp.k4, c43, dp.k3) - F(aux, t+h*c40, dp.k4) - dp.k4.ScaleVec(h, aux) - - dp.k5.AddScaledVec(x, c51, dp.k1) - dp.k5.AddScaledVec(dp.k5, c52, dp.k2) - dp.k5.AddScaledVec(dp.k5, c53, dp.k3) - dp.k5.AddScaledVec(dp.k5, c54, dp.k4) - F(aux, t+h*c50, dp.k5) - dp.k5.ScaleVec(h, aux) - - dp.k6.AddScaledVec(x, c61, dp.k1) - dp.k6.AddScaledVec(dp.k6, c62, dp.k2) - dp.k6.AddScaledVec(dp.k6, c63, dp.k3) - dp.k6.AddScaledVec(dp.k6, c64, dp.k4) - dp.k6.AddScaledVec(dp.k6, c65, dp.k5) - F(aux, t+h*c60, dp.k6) - dp.k6.ScaleVec(h, aux) - - // TODO y4 could be dp.x, is there a disadvantage to this change? - y4.AddScaledVec(x, b1, dp.k1) - y4.AddScaledVec(y4, b3, dp.k3) - y4.AddScaledVec(y4, b4, dp.k4) - y4.AddScaledVec(y4, b5, dp.k5) - y4.AddScaledVec(y4, b6, dp.k6) - // fourth order approximation used to advance solution - - if dp.adaptive { - y5 := dp.y5 - dp.k7.AddScaledVec(x, c71, dp.k1) - dp.k7.AddScaledVec(dp.k7, c72, dp.k2) - dp.k7.AddScaledVec(dp.k7, c73, dp.k3) - dp.k7.AddScaledVec(dp.k7, c74, dp.k4) - dp.k7.AddScaledVec(dp.k7, c75, dp.k5) - dp.k7.AddScaledVec(dp.k7, c76, dp.k6) - // calculation of order 5 coefficients - - F(aux, t+h*c70, dp.k7) - dp.k7.ScaleVec(h, aux) - - y5.AddScaledVec(x, a1, dp.k1) - y5.AddScaledVec(y5, a3, dp.k3) - y5.AddScaledVec(y5, a4, dp.k4) - y5.AddScaledVec(y5, a5, dp.k5) - y5.AddScaledVec(y5, a6, dp.k6) - y5.AddScaledVec(y5, a7, dp.k7) - dp.err45.SubVec(y4, y5) - - // compute error between fifth order solution and fourth order solution - errRatio := dp.param.atol / mat.Norm(dp.err45, math.Inf(1)) - hnew := math.Min(math.Max(0.9*h*math.Pow(errRatio, .2), dp.param.minStep), dp.param.maxStep) - // given "bad" error ratio, algorithm will recompute with a smaller step, given it is not the minimum permissible step or less - if errRatio < 1 && h > dp.param.minStep { - h = hnew - goto SOLVE - } - } - // advance solution with fourth order solution - dp.x.CloneFromVec(y4) - dp.dom += h - return h, nil -} - -// 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)) -// -// If a invalid configuration is passed the function panics. -func NewDormandPrince5(configs ...Configuration) *DoPri5 { - dp := new(DoPri5) - for i := range configs { - err := configs[i](&dp.param) - if err != nil { - panic(err) - } - } - - if dp.param.atol > 0 && dp.param.minStep > 0 { - dp.adaptive = true - } - return dp -} diff --git a/ode/ivp.go b/ode/ivp.go deleted file mode 100644 index 47689bc..0000000 --- a/ode/ivp.go +++ /dev/null @@ -1,49 +0,0 @@ -// Copyright ©2021 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 ( - "errors" - "math" - - "gonum.org/v1/gonum/mat" -) - -// 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(0) = y_0 -// -// Where: -// t is a scalar representing the integration domain, which is time for most physical problems. -// y is the state vector. -// y' is the vector of first derivatives of the state vector y. -// f are the differential equations represented by Func. -// An initial value problem is characterized by the initial conditions imposed -// on the state vector y at the beginning of the integration domain. These -// initial conditions are returned by the IV() method for the state vector -// as y0. -// -// The term "state vector" and "state variables" are used interchangeably -// throughout the code and refer to y vector of independent variables. -type IVP struct { - // Initial values for the state vector - Y0 mat.Vector - // Independent variable point at which Y0 is evaluated - T0 float64 - // Func are the differential equations f(t,y(t)). - // The result is y'(t) and is stored in dst. - Func func(dst *mat.VecDense, t float64, y mat.Vector) -} - -// 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 { - return IVP{}, errors.New("bad model value") - } - return IVP{Func: f, Y0: y0, T0: t0}, nil -} diff --git a/ode/ivp_test.go b/ode/ivp_test.go deleted file mode 100644 index 58d69b3..0000000 --- a/ode/ivp_test.go +++ /dev/null @@ -1,173 +0,0 @@ -// Copyright ©2021 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_test - -import ( - "fmt" - "log" - "math" - "testing" - - "gonum.org/v1/exp/ode" - "gonum.org/v1/gonum/mat" -) - -type TestModel struct { - ode.IVP - solution ode.IVP - err func(h, i float64) float64 -} - -func TestSolve(t *testing.T) { - // domain start - quad := quadTestModel(t) - ivp := quad.IVP - solver := ode.NewDormandPrince5() - solver.Init(ivp) - stepsize := 0.5 / 15. - end := 0.5 - results, err := ode.SolveIVP(ivp, solver, stepsize, end) - if err != nil { - t.Fatal(err) - } - // run the solver and corsscheck values - quadsol := quad.solution - x0 := quadsol.Y0 - - sol := quadsol.Func - - solresults := mat.NewVecDense(x0.Len(), nil) - for i := range results { - sol(solresults, results[i].T, results[i].Y) - for j := range solresults.RawVector().Data { - sol := solresults.AtVec(j) - res := results[i].Y.AtVec(j) - got := math.Abs(solresults.AtVec(j) - results[i].Y.AtVec(j)) - expect := quad.err(stepsize, float64(i)) - if got > expect { - t.Errorf("error %g is greater than permitted tolerance %g. solution:[%0.3g], result:[%0.3g]", got, expect, sol, res) - } - } - } - -} - -func Example_solve() { - const ( - g = -10. // gravity field [m.s^-2] - ) - // we declare our physical model. First argument is initial time, which is 0 seconds. - // Next is the initial state vector, which corresponds to 100 meters above the ground - // with 0 m/s velocity. - ballModel, err := ode.NewIVP(0, mat.NewVecDense(2, []float64{100., 0.}), - func(yvec *mat.VecDense, _ float64, xvec mat.Vector) { - // this anonymous function defines the physics. - // The first variable xvec[0] corresponds to position - // second variable xvec[1] is velocity. - Dx := xvec.AtVec(1) - // yvec represents change in xvec, or derivative with respect to domain - // Change in position will be equal to velocity, which is the second variable: - // thus yvec[0] = xvec[1], which is the same as saying "change in xvec[0]" is equal to xvec[1] - yvec.SetVec(0, Dx) - // change in velocity is acceleration. We suppose our ball is on earth accelerating at `g` - yvec.SetVec(1, g) - }) - if err != nil { - log.Fatal(err) - } - // Here we choose our algorithm. Runge-Kutta 4th order is used - var solver ode.Integrator = ode.NewDormandPrince5() - // 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. - results, err := ode.SolveIVP(ballModel, solver, 0.1, 10.) - fmt.Println(results) -} - -// Quadratic model may be used for future algorithms -func quadTestModel(t *testing.T) *TestModel { - t0 := 0.0 - Quadratic := new(TestModel) - quad, err := ode.NewIVP(t0, mat.NewVecDense(2, []float64{0, 0}), - func(dst *mat.VecDense, t float64, x mat.Vector) { - dst.SetVec(0, x.AtVec(1)) - dst.SetVec(1, 1.) - }) - if err != nil { - t.Fatal(err) - } - Quadratic.IVP = quad - quadsol, err := ode.NewIVP(t0, mat.NewVecDense(2, []float64{0, 0}), - func(dst *mat.VecDense, t float64, x mat.Vector) { - dst.SetVec(0, t*t/2.) - dst.SetVec(1, t) - }) - if err != nil { - t.Fatal(err) - } - Quadratic.solution = quadsol - Quadratic.err = func(h, i float64) float64 { return math.Pow(h*i, 4) + 1e-10 } - return Quadratic -} - -// exponential unidimensional model may be used for future algorithms -// 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 - Quadratic := new(TestModel) - quad, err := ode.NewIVP(t0, mat.NewVecDense(1, []float64{1.}), - func(dst *mat.VecDense, t float64, x mat.Vector) { - dst.SetVec(0, tau*x.AtVec(0)) - }) - if err != nil { - t.Fatal(err) - } - Quadratic.IVP = quad - quadsol, err := ode.NewIVP(t0, mat.NewVecDense(1, []float64{0}), - func(dst *mat.VecDense, t float64, x mat.Vector) { - dst.SetVec(0, math.Exp(tau*t)) - - }) - if err != nil { - t.Fatal(err) - } - Quadratic.solution = quadsol - Quadratic.err = func(h, i float64) float64 { return math.Pow(h*i, 4) + 1e-10 } - return Quadratic -} - -// func TestQuadratic(t *testing.T) { -// Quadratic := quadTestModel(t) -// solver := ivp.NewDormandPrince5() - -// solver.Set(Quadratic.Model) - -// _, x0 := Quadratic.Model.IV() -// steps := 10 -// dt := 0.1 - -// results := make([]float64, nx) - -// solmodel := Quadratic.solution -// soleq := solmodel.Equations() -// solDims, _ := solmodel.Dims() -// solution := make([]float64, solDims) -// for i := 1.; i < float64(steps+1); i++ { -// dom := i * dt -// solver.Step(results, dt) -// soleq(solution, dom, results) -// for j := range results { -// got := math.Abs(solution[j] - results[j]) -// expected := Quadratic.err(dt, i) -// if got > expected { -// t.Errorf("error %e greater than threshold %e", got, expected) -// } - -// } -// } -// } diff --git a/ode/ode.go b/ode/ode.go deleted file mode 100644 index 953382f..0000000 --- a/ode/ode.go +++ /dev/null @@ -1,65 +0,0 @@ -// Copyright ©2021 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 ( - "errors" - "math" - - "gonum.org/v1/gonum/mat" -) - -// Integrator can integrate an initial-value problem (IVP) for a first-order -// system of ordinary differential equations (ODEs). -type Integrator interface { - // Init initializes the integrator and sets the initial condition. - 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) - - // State stores the current state of the integrator in-place in dst. - State(dst *State) -} - -type State struct { - T float64 - Y *mat.VecDense -} - -// SolveIVP solves an already initialized Integrator returning state vector results. -func SolveIVP(p IVP, solver Integrator, stepsize, tend float64) (results []State, err error) { - t0, x0 := p.T0, mat.VecDenseCopyOf(p.Y0) - nx := x0.Len() - if nx == 0 { - return nil, errors.New("state vector length can not be equal to zero. Has ivp been set?") - } - // calculate expected size of results, these may differ - domainLength := tend - t0 - expectedLength := int(domainLength/stepsize) + 1 - results = make([]State, 0, expectedLength) - - // t stores current domain - t := t0 - for t < tend { - res := State{Y: mat.NewVecDense(nx, nil)} - - if t-tend > 1e-10 { - stepsize = math.Min(stepsize, (t-tend)*(1+1e-3)) - } - stepsize, err = solver.Step(stepsize) - if err != nil { - return results, err - } - if stepsize <= 0 { - return results, errors.New("got zero or negative step size from Integrator") - } - solver.State(&res) - t += stepsize - results = append(results, res) - } - return results, nil -} diff --git a/ode/doc.go b/posit/doc.go similarity index 52% rename from ode/doc.go rename to posit/doc.go index b647ba7..6ebebfb 100644 --- a/ode/doc.go +++ b/posit/doc.go @@ -2,6 +2,6 @@ // Use of this source code is governed by a BSD-style // license that can be found in the LICENSE file. -// package ode provides interfaces and structures to solve -// multivariable ODE initial value problems. -package ode // import "gonum.org/v1/gonum/ode" +// Package posit implements a generalized unum type II implementation +// as a byte slice and provides primitives for working with posits. +package posit diff --git a/posit/posit.go b/posit/posit.go new file mode 100644 index 0000000..b263cfd --- /dev/null +++ b/posit/posit.go @@ -0,0 +1,143 @@ +// Copyright ©2021 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 posit + +import ( + "fmt" + "math" +) + +type Posit []byte + +func (p Posit) regime() (k int) { + // Regime lead bit. + const rbit = 0b0100_0000 + // Indicates if leading Regime bit is a one. + oneRun := rbit&p[0] != 0 +outer: + for i := range p { + for j := 0; j < 8; j++ { + // Ignore sign bit. + if i == 0 && j == 0 { + j++ + } + // If current bit is part of the run, add to r. + if (oneRun && p[i]&(1<<(7-j)) != 0) || (!oneRun && p[i]&(1<<(7-j)) == 0) { + k++ + } else { + break outer + } + } + } + if !oneRun { + k *= -1 + } else { + k-- + } + return k +} + +// sign returns true if posit negative +func (p Posit) sign() bool { return p[0]&(1<<7) != 0 } + +func (p Posit) bits() int { return 8 * len(p) } + +func (p Posit) String() string { return fmt.Sprintf("%08b", []byte(p)) } + +// es represents max amount of exponent bits that can be present in the posit. +// +// Values taken from http://www.johngustafson.net/pdfs/BeatingFloatingPoint.pdf +// (table 3) to match or exceed IEEE float dynamic range. +func (p Posit) es() (es int) { + switch p.bits() { + case 16: + es = 1 + case 32: + es = 3 + case 64: + es = 4 + case 128: + es = 7 + case 256: + es = 10 + default: + es = p.bits() / 16 // 8 bit posit has es == 0. + } + return es +} + +// exp returns the exponent part of the posit (2**e). +func (p Posit) exp() (exp int) { + // bits in front of exp. + flen := p.regimeLen() + 2 // sign and opposite bits included + es := p.es() + // Check if exp bits present for quick return. + if flen >= p.bits() || es == 0 { + return 0 + } + + expcount := 0 +outer: + for i := flen / 8; i < len(p); i++ { + for j := flen % 8; j < 8; j++ { + if expcount == es { + break outer + } + exp <<= 1 + exp |= (int(p[i]) & (1 << (7 - j))) >> (7 - j) + expcount++ + } + } + return exp +} + +// returns regime length for a given posit in number of bits. +func (p Posit) regimeLen() int { + r := p.regime() + if r < 0 { + return -r + } + return r + 1 +} + +// useed defines the midway point of accuracy from 1 to +inf and +// conversely from 1 to 0. It depends on es. +func (p Posit) useed() int { return 1 << (1 << (p.es())) } + +// fraction returns the numerator of the fraction part. +func (p Posit) fraction() (frac int) { + // bits in front of fraction. + flen := p.regimeLen() + p.es() + 2 // sign and opposite bits included + // Check if exp bits present for quick return. + if flen >= p.bits() { + return 0 + } + + for i := flen / 8; i < len(p); i++ { + for j := flen % 8; j < 8; j++ { + frac <<= 1 + frac |= (int(p[i]) & (1 << (7 - j))) >> (7 - j) + } + } + return frac +} + +func (p Posit) ToFloat64() float64 { + reg := float64(p.regime()) + useed := float64(p.useed()) + exp := 1 << p.exp() + return math.Pow(useed, reg) * float64(exp) * (1 + float64(p.fraction())/useed) +} + +// Format implements fmt.Formatter. +func (p Posit) Format(fs fmt.State, c rune) { + switch c { + case 'v', 'f': + fmt.Fprintf(fs, "%T{%f}", p, p.ToFloat64()) + default: + fmt.Fprintf(fs, "%%!%c(%T=%[2]v)", c, p) + return + } +} diff --git a/posit/posit_test.go b/posit/posit_test.go new file mode 100644 index 0000000..e8bb5da --- /dev/null +++ b/posit/posit_test.go @@ -0,0 +1,121 @@ +// Copyright ©2021 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 posit + +import ( + "math" + "testing" +) + +func TestRegime(t *testing.T) { + for _, test := range []struct { + p Posit + expect int + }{ + // posit8 + {Posit{0b_100_0000}, 0}, + {Posit{0b_111_0000}, 2}, + {Posit{0b_010_0000}, -1}, + {Posit{0b_111_1000}, 3}, + {Posit{0b_000_0100}, -4}, + + // posit16 + {Posit{0, 0b0010_0000}, -9}, + {Posit{0xff, 0b0010_0000}, 6}, + + // posit 32 + {Posit{0, 0, 0b1000_0001, 0}, -7 - 8}, + {Posit{0, 0, 0xff, 0}, -7 - 8}, + {Posit{0xff, 0xff, 0xff, 0}, 7 + 8 + 8 - 1}, + {Posit{0xff, 0, 0xff, 0}, 7 - 1}, + } { + got := test.p.regime() + if got != test.expect { + t.Errorf("posit %s expected regime %v, got %v", test.p, test.expect, got) + } + } +} + +func TestSign(t *testing.T) { + for _, test := range []struct { + p Posit + expect bool + }{ + {Posit{0b_1100_0000}, true}, + {Posit{0b_0111_0000}, false}, + {Posit{0b1010_0000}, true}, + {Posit{0b_0111_1000}, false}, + {Posit{0b_1000_0100}, true}, + {Posit{0, 0b0010_0000}, false}, + } { + got := test.p.sign() + if got != test.expect { + t.Errorf("posit %s expected sign %t, got %t", test.p, test.expect, got) + } + } +} + +func TestUseed(t *testing.T) { + for _, test := range []struct { + p Posit + }{ + {Posit{0b1100_0000}}, + + {Posit{0, 0b0010_0000}}, + {Posit{0, 0b0010_0000, 0}}, + } { + es := test.p.es() + expect := int(math.RoundToEven(math.Pow(2, math.Pow(2, float64(es))))) + got := test.p.useed() + if expect != got { + t.Errorf("posit %s expected useed %v, got %v", test.p, expect, got) + } + } +} + +func TestExp(t *testing.T) { + for _, test := range []struct { + p Posit + expect int + }{ + // posit16 has es == 1 + {Posit{0b_101_0000, 0}, 1}, + {Posit{0b_100_0000, 0}, 0}, + {Posit{0b_111_1001, 0}, 0}, + {Posit{0b_111_1011, 0}, 1}, + // posit 32 has es == 3 + {Posit{0, 0, 0b10111, 0}, 3}, + {Posit{0, 0, 0b11111, 0}, 7}, + } { + got := test.p.exp() + expect := test.expect + + if expect != got { + t.Errorf("posit %s expected exp %v, got %v", test.p, expect, got) + } + } +} +func TestFraction(t *testing.T) { + for _, test := range []struct { + p Posit + expect int + }{ + // posit16 has es == 1 + {Posit{0b_101_0000, 0}, 1}, + {Posit{0b_100_0000, 0}, 0}, + {Posit{0b_111_1001, 0}, 0}, + {Posit{0b_111_1011, 0}, 1}, + // posit 32 has es == 3 + {Posit{0, 0, 0b10111, 0}, 3}, + {Posit{0, 0, 0b11111, 0}, 7}, + } { + got := test.p.fraction() + expect := test.expect + + if expect != got { + t.Errorf("posit %s expected exp %v, got %v", test.p, expect, got) + } + } +}