From 46afb6bf501f4001c4bdb7bd2f2db7c466f95554 Mon Sep 17 00:00:00 2001 From: Juan Marin Noguera Date: Tue, 25 Jul 2023 18:22:04 +0200 Subject: Solvned project from 2020 Note that Go didn't have generics back then. --- COPYING.LESSER | 166 ++++++++++++++++++++++ README.md | 10 ++ caching.go | 301 ++++++++++++++++++++++++++++++++++++++++ common.go | 75 ++++++++++ dense.go | 174 +++++++++++++++++++++++ event.go | 95 +++++++++++++ examples/arenstorf.go | 74 ++++++++++ examples/arenstorf.png | Bin 0 -> 62397 bytes examples/arenstorf.txt | 1 + examples/exp.go | 74 ++++++++++ examples/exp.txt | 5 + examples/hooke.go | 77 ++++++++++ examples/hooke.png | Bin 0 -> 260239 bytes examples/hooke_check.go | 91 ++++++++++++ examples/hooke_check.txt | 4 + examples/parabolic.go | 110 +++++++++++++++ examples/parabolic.png | Bin 0 -> 122494 bytes examples/parabolic_check.go | 108 +++++++++++++++ examples/parabolic_check.txt | 3 + interpolate.go | 97 +++++++++++++ ivp/arenstorf.go | 68 +++++++++ ivp/spring.go | 79 +++++++++++ ivp/throw.go | 68 +++++++++ method/euler.go | 79 +++++++++++ method/fehlberg.go | 208 +++++++++++++++++++++++++++ method/rk.go | 107 ++++++++++++++ mstep/adams.go | 105 ++++++++++++++ mstep/common.go | 53 +++++++ mstep/implicit.go | 164 ++++++++++++++++++++++ mstep/pc.go | 300 +++++++++++++++++++++++++++++++++++++++ mstep/root.go | 71 ++++++++++ step.go | 324 +++++++++++++++++++++++++++++++++++++++++++ 32 files changed, 3091 insertions(+) create mode 100644 COPYING.LESSER create mode 100644 README.md create mode 100644 caching.go create mode 100644 common.go create mode 100644 dense.go create mode 100644 event.go create mode 100644 examples/arenstorf.go create mode 100644 examples/arenstorf.png create mode 100644 examples/arenstorf.txt create mode 100644 examples/exp.go create mode 100644 examples/exp.txt create mode 100644 examples/hooke.go create mode 100644 examples/hooke.png create mode 100644 examples/hooke_check.go create mode 100644 examples/hooke_check.txt create mode 100644 examples/parabolic.go create mode 100644 examples/parabolic.png create mode 100644 examples/parabolic_check.go create mode 100644 examples/parabolic_check.txt create mode 100644 interpolate.go create mode 100644 ivp/arenstorf.go create mode 100644 ivp/spring.go create mode 100644 ivp/throw.go create mode 100644 method/euler.go create mode 100644 method/fehlberg.go create mode 100644 method/rk.go create mode 100644 mstep/adams.go create mode 100644 mstep/common.go create mode 100644 mstep/implicit.go create mode 100644 mstep/pc.go create mode 100644 mstep/root.go create mode 100644 step.go diff --git a/COPYING.LESSER b/COPYING.LESSER new file mode 100644 index 0000000..5357f69 --- /dev/null +++ b/COPYING.LESSER @@ -0,0 +1,166 @@ + GNU LESSER GENERAL PUBLIC LICENSE + Version 3, 29 June 2007 + + Copyright (C) 2007 Free Software Foundation, Inc. + Everyone is permitted to copy and distribute verbatim copies + of this license document, but changing it is not allowed. + + + This version of the GNU Lesser General Public License incorporates +the terms and conditions of version 3 of the GNU General Public +License, supplemented by the additional permissions listed below. + + 0. Additional Definitions. + + As used herein, "this License" refers to version 3 of the GNU Lesser +General Public License, and the "GNU GPL" refers to version 3 of the GNU +General Public License. + + "The Library" refers to a covered work governed by this License, +other than an Application or a Combined Work as defined below. + + An "Application" is any work that makes use of an interface provided +by the Library, but which is not otherwise based on the Library. +Defining a subclass of a class defined by the Library is deemed a mode +of using an interface provided by the Library. + + A "Combined Work" is a work produced by combining or linking an +Application with the Library. The particular version of the Library +with which the Combined Work was made is also called the "Linked +Version". + + The "Minimal Corresponding Source" for a Combined Work means the +Corresponding Source for the Combined Work, excluding any source code +for portions of the Combined Work that, considered in isolation, are +based on the Application, and not on the Linked Version. + + The "Corresponding Application Code" for a Combined Work means the +object code and/or source code for the Application, including any data +and utility programs needed for reproducing the Combined Work from the +Application, but excluding the System Libraries of the Combined Work. + + 1. Exception to Section 3 of the GNU GPL. + + You may convey a covered work under sections 3 and 4 of this License +without being bound by section 3 of the GNU GPL. + + 2. Conveying Modified Versions. + + If you modify a copy of the Library, and, in your modifications, a +facility refers to a function or data to be supplied by an Application +that uses the facility (other than as an argument passed when the +facility is invoked), then you may convey a copy of the modified +version: + + a) under this License, provided that you make a good faith effort to + ensure that, in the event an Application does not supply the + function or data, the facility still operates, and performs + whatever part of its purpose remains meaningful, or + + b) under the GNU GPL, with none of the additional permissions of + this License applicable to that copy. + + 3. Object Code Incorporating Material from Library Header Files. + + The object code form of an Application may incorporate material from +a header file that is part of the Library. You may convey such object +code under terms of your choice, provided that, if the incorporated +material is not limited to numerical parameters, data structure +layouts and accessors, or small macros, inline functions and templates +(ten or fewer lines in length), you do both of the following: + + a) Give prominent notice with each copy of the object code that the + Library is used in it and that the Library and its use are + covered by this License. + + b) Accompany the object code with a copy of the GNU GPL and this license + document. + + 4. Combined Works. + + You may convey a Combined Work under terms of your choice that, +taken together, effectively do not restrict modification of the +portions of the Library contained in the Combined Work and reverse +engineering for debugging such modifications, if you also do each of +the following: + + a) Give prominent notice with each copy of the Combined Work that + the Library is used in it and that the Library and its use are + covered by this License. + + b) Accompany the Combined Work with a copy of the GNU GPL and this license + document. + + c) For a Combined Work that displays copyright notices during + execution, include the copyright notice for the Library among + these notices, as well as a reference directing the user to the + copies of the GNU GPL and this license document. + + d) Do one of the following: + + 0) Convey the Minimal Corresponding Source under the terms of this + License, and the Corresponding Application Code in a form + suitable for, and under terms that permit, the user to + recombine or relink the Application with a modified version of + the Linked Version to produce a modified Combined Work, in the + manner specified by section 6 of the GNU GPL for conveying + Corresponding Source. + + 1) Use a suitable shared library mechanism for linking with the + Library. A suitable mechanism is one that (a) uses at run time + a copy of the Library already present on the user's computer + system, and (b) will operate properly with a modified version + of the Library that is interface-compatible with the Linked + Version. + + e) Provide Installation Information, but only if you would otherwise + be required to provide such information under section 6 of the + GNU GPL, and only to the extent that such information is + necessary to install and execute a modified version of the + Combined Work produced by recombining or relinking the + Application with a modified version of the Linked Version. (If + you use option 4d0, the Installation Information must accompany + the Minimal Corresponding Source and Corresponding Application + Code. If you use option 4d1, you must provide the Installation + Information in the manner specified by section 6 of the GNU GPL + for conveying Corresponding Source.) + + 5. Combined Libraries. + + You may place library facilities that are a work based on the +Library side by side in a single library together with other library +facilities that are not Applications and are not covered by this +License, and convey such a combined library under terms of your +choice, if you do both of the following: + + a) Accompany the combined library with a copy of the same work based + on the Library, uncombined with any other library facilities, + conveyed under the terms of this License. + + b) Give prominent notice with the combined library that part of it + is a work based on the Library, and explaining where to find the + accompanying uncombined form of the same work. + + 6. Revised Versions of the GNU Lesser General Public License. + + The Free Software Foundation may publish revised and/or new versions +of the GNU Lesser General Public License from time to time. Such new +versions will be similar in spirit to the present version, but may +differ in detail to address new problems or concerns. + + Each version is given a distinguishing version number. If the +Library as you received it specifies that a certain numbered version +of the GNU Lesser General Public License "or any later version" +applies to it, you have the option of following the terms and +conditions either of that published version or of any later version +published by the Free Software Foundation. If the Library as you +received it does not specify a version number of the GNU Lesser +General Public License, you may choose any version of the GNU Lesser +General Public License ever published by the Free Software Foundation. + + If the Library as you received it specifies that a proxy can decide +whether future versions of the GNU Lesser General Public License shall +apply, that proxy's public statement of acceptance of any version is +permanent authorization for you to choose that version for the +Library. + diff --git a/README.md b/README.md new file mode 100644 index 0000000..3e32a67 --- /dev/null +++ b/README.md @@ -0,0 +1,10 @@ +# Solvned: A numerical ODE framework for Go + +Solvned is a framework for solving ordinary differential equations in Go. It +provides a variety of methods to solve these equations, from the simple Euler +method to the sophisticated Runge-Kutta-Fehlberg method, as well as the means +to define new methods and operate with them in a general way, with +interpolation and an event mechanism for iteration until a given point. + +Check out the visualization examples at the `examples` directory and the +[pkg.go.dev documentation](https://pkg.go.dev/github.com/JwanMan/mned). diff --git a/caching.go b/caching.go new file mode 100644 index 0000000..1d31b94 --- /dev/null +++ b/caching.go @@ -0,0 +1,301 @@ +// Caching of solution points. +// +// Copyright (C) 2020 Juan Marín Noguera +// +// This file is part of Solvned. +// +// Solvned is free software: you can redistribute it and/or modify it under the +// terms of the GNU Lesser General Public License as published by the Free +// Software Foundation, either version 3 of the License, or (at your option) any +// later version. +// +// Solvned is distributed in the hope that it will be useful, but WITHOUT ANY +// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR +// A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more +// details. +// +// You should have received a copy of the GNU Lesser General Public License +// along with Solvned. If not, see . + +package mned + +import "sort" + +// A CacheSolver is like a dynamic version of a DenseSolution. It stores +// points of the solution of a problem in a conceptually maximal interval that +// is evaluated lazily; that is, it starts with the solutions in a range that +// only includes the initial point and, when a point is asked for a Time value +// outside the range, the range is expanded to cover that point. +type CacheSolver struct { + backward []Point // Invariant: Decreasing times, times before forward's. + forward []Point // Invariant: Non-empty, increasing times. + backStep Stepper // Invariant: Decreasing times from end of backward. + forStep Stepper // Invariant: Increasing times from end of forward. + evs []Event + interp Interpolator +} + +// Create a CacheSolver to solve the given IVP with the given solving Method, +// with the given interpolator to calculate points between steps and taking into +// account the events provided. +func CacheSolve( + m Method, ivp *IVP, interp Interpolator, evs ...Event, +) CacheSolver { + return CacheSolver{ + forward: []Point{ivp.Start.Clone()}, + backward: []Point{}, + forStep: m.Forward(ivp), + backStep: m.Backward(ivp), + evs: evs, + interp: interp, + } +} + +// Get the time of the stored solution point with the lowest time. +func (c *CacheSolver) Start() float64 { + if len(c.backward) > 0 { + return c.backward[len(c.backward)-1].Time + } else { + return c.forward[0].Time + } +} + +// Get the time of the stored solution point with the greatest time. +func (c *CacheSolver) End() float64 { + return c.forward[len(c.forward)-1].Time +} + +func (c *CacheSolver) getActions( + current *Point, next *Point, vals []float64, +) requiredActions { + var actions requiredActions = make([]requiredAction, 0) + for i := 0; i < len(vals); i++ { + newVal := c.evs[i].Cross(next) + if differentSign(vals[i], newVal) { + point := c.evs[i].FindPoint(c.interp, current, next) + actions = append( + actions, + requiredAction{index: i, point: point}, + ) + } + vals[i] = newVal + } + sort.Sort(actions) + return actions +} + +func (c *CacheSolver) expandBackward(t float64) bool { + if c.backStep == nil { + return false + } + + var current Point + if len(c.backward) == 0 { + current = c.forward[0] + } else { + current = c.backward[len(c.backward)-1] + } + + vals := make([]float64, len(c.evs)) + for i := 0; i < len(vals); i++ { + vals[i] = c.evs[i].Cross(¤t) + } + + for t < current.Time { + next, ok := c.backStep.Next() + if !ok { + // TODO Consider bisection for reasonable last point + c.backStep = nil + return false + } + actions := c.getActions(¤t, next, vals) + for i := len(actions) - 1; i >= 0; i-- { + if !c.evs[actions[i].index].Action(&actions[i].point) { + return false + } + } + current = next.Clone() + c.backward = append(c.backward, current) + } + return true +} + +func (c *CacheSolver) expandForward(t float64) bool { + if c.forStep == nil { + return false + } + + current := c.forward[len(c.forward)-1] + vals := make([]float64, len(c.evs)) + for i := 0; i < len(vals); i++ { + vals[i] = c.evs[i].Cross(¤t) + } + + for t < current.Time { + next, ok := c.forStep.Next() + if !ok { + // TODO Consider bisection for reasonable last point + c.forStep = nil + return false + } + actions := c.getActions(¤t, next, vals) + for i := 0; i < len(actions); i++ { + if !c.evs[actions[i].index].Action(&actions[i].point) { + return false + } + } + current = next.Clone() + c.forward = append(c.forward, current) + } + return true +} + +// Get the value of the solution for a given time. +// +// If the value is outside of the bounds of what's stored, values are added +// until reaching the given time. If ok is false, the value is out of bounds for +// the problem or some event and points have only been added as far as it was +// possible. +func (c *CacheSolver) Get(t float64) (x []float64, ok bool) { + if t < c.Start() { + if ok := c.expandBackward(t); !ok { + return nil, false + } + } + if t > c.End() { + if ok := c.expandForward(t); !ok { + return nil, false + } + } + if t >= c.forward[0].Time { + i := 1 + for t > c.forward[i].Time { + i++ + } + return c.interp.FindValue( + &c.forward[i-1], &c.forward[i], t, + ), true + } + + i := 0 + for t > c.backward[i].Time { + i++ + } + if i == 0 { + return c.interp.FindValue( + &c.backward[0], &c.forward[0], t, + ), true + } + return c.interp.FindValue(&c.backward[i], &c.backward[i-1], t), true +} + +// Get the value at `c.Start() - h` as in `c.Get` except that, if the underlying +// backward Stepper is a ConfigurableStepper, its method `NextStep(-h)` is used. +// The result value should *NOT* be mutated. +// +// This is most useful with fixed step methods to specify the step. +func (c *CacheSolver) StepBackward(h float64) ([]float64, bool) { + if cs, ok := c.backStep.(ConfigurableStepper); ok { + if point, ok := cs.NextStep(-h); ok { + if h > 0 { + c.backward = append(c.backward, point.Clone()) + } + return point.Value, true + } else { + return nil, false + } + } + return c.Get(c.Start() - h) +} + +// Get the value at `c.End() + h` as in `c.Get` except that, if the underlying +// forward Stepper is a ConfigurableStepper, its method `NextStep(h)` is used. +// The result value should *NOT* be mutated. +// +// This is most useful with fixed step methods to specify the step. +func (c *CacheSolver) StepForward(h float64) ([]float64, bool) { + if cs, ok := c.forStep.(ConfigurableStepper); ok { + if point, ok := cs.NextStep(h); ok { + if h > 0 { + c.forward = append(c.forward, point.Clone()) + } + return point.Value, true + } else { + return nil, false + } + } + return c.Get(c.End() + h) +} + +// Get the points calculated from the initial values to earlier times. The +// points **MUST NOT** be modified. +func (c *CacheSolver) BackwardPoints() []Point { + return c.backward +} + +// Get the points calculated from the initial values to later times. The points +// **MUST NOT** be modified. +func (c *CacheSolver) ForwardPoints() []Point { + return c.forward +} + +// Rearrange the stored solution points in a result such that the i-th point has +// time `result[0][i]` and value `(result[1][i],...,result[size][i])`, where +// `size` is the number of dimensions of a value. This is useful for plotting. +func (c *CacheSolver) PointCoords() [][]float64 { + size := len(c.forward[0].Value) + result := make([][]float64, size+1) + backs := len(c.backward) + elems := backs + len(c.forward) + + for i := 0; i <= size; i++ { + result[i] = make([]float64, elems) + } + for i := 0; i < backs; i++ { + result[0][i] = c.backward[backs-i-1].Time + for j := 0; j < size; j++ { + result[j+1][i] = c.backward[backs-i-1].Value[j] + } + } + for i := backs; i < elems; i++ { + result[0][i] = c.forward[i-backs].Time + for j := 0; j < size; j++ { + result[j+1][i] = c.forward[i-backs].Value[j] + } + } + return result +} + +// Force the generation of values to the left until the domain ends or an event +// tells the stepping to stop. +func (c *CacheSolver) StepToBeginning() { + var init Point + if len(c.backward) > 0 { + init = c.backward[len(c.backward)-1] + } else { + init = c.forward[0] + } + StepUntil( + c.backStep, + init, + c.interp, + func(point *Point) { + c.backward = append(c.backward, point.Clone()) + }, + c.evs..., + ) +} + +// Force the generation of values to the right until the domain ends or an event +// tells the stepping to stop. +func (c *CacheSolver) StepToEnd() { + StepUntil( + c.forStep, + c.forward[len(c.forward)-1], + c.interp, + func(point *Point) { + c.forward = append(c.forward, point.Clone()) + }, + c.evs..., + ) +} diff --git a/common.go b/common.go new file mode 100644 index 0000000..9b41888 --- /dev/null +++ b/common.go @@ -0,0 +1,75 @@ +// Core data type definitions. +// +// Copyright (C) 2020 Juan Marín Noguera +// +// This file is part of Solvned. +// +// Solvned is free software: you can redistribute it and/or modify it under the +// terms of the GNU Lesser General Public License as published by the Free +// Software Foundation, either version 3 of the License, or (at your option) any +// later version. +// +// Solvned is distributed in the hope that it will be useful, but WITHOUT ANY +// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR +// A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more +// details. +// +// You should have received a copy of the GNU Lesser General Public License +// along with Solvned. If not, see . + +package mned + +/* +This package implements several methods to compute approximations of +solutions of initial value problems (IVPs) with first-order ODEs. These problems +have the form +``` +x'(t) = f(t, x(t)), +x(t0) = x0, +``` +where `t0 : R` and `x0 : R^n` are the **initial values**, +`f : Ω ⊆ R x R^n -> R^n` is the **derivative function**, and `x : I ⊆ R -> R^n` +is the unknown of the equation. Here `I` is an open interval of `R` containing +`t0` and `Ω` is an open subset of `R x R^n` containing `(t0, x0)`. + +We refer to `t` as the **independent variable** and to `x(t)` as the `dependent +variable`. We know that `x` is uniquely defined in a neighbourhood of `t0`. +Given a point `t` in such a neighbourhood, we say that `(t,x(t))` is a point of +the solution of the initial value problem. + +The results calculated will be approximations, and judgement is required to get +approximations suitable for a particular problem; nevertheless, for the purpose +of explanation, we will use the same terminology for the approximations than for +the actual functions. +*/ + +// A Point is an element of the solution, given by the values of the independent +// and dependent variables. +type Point struct { + Time float64 // The independent variable. + Value []float64 // The dependent variable. +} + +// Deep-copy a point. +func (p *Point) Clone() Point { + value := make([]float64, len(p.Value)) + copy(value, p.Value) + return Point{Time: p.Time, Value: value} +} + +// An IVP is an initial value problem given by a first-order ODE. It's given by +// the Derivative function, which is a pure function, and the initial values. +// +// The second return value of Derivative indicates if the Point was in the +// domain of the function. If it wasn't the first return value should not be +// used. Ideally, the domain should be restricted to the connected component of +// the initial value, as otherwise a solving method could "jump" to another +// component and the results from there on would be invalid. +type IVP struct { + Derivative func(Point) ([]float64, bool) + Start Point +} + +func (ivp *IVP) Clone() IVP { + return IVP{Derivative: ivp.Derivative, Start: ivp.Start.Clone()} +} diff --git a/dense.go b/dense.go new file mode 100644 index 0000000..8bfe9c5 --- /dev/null +++ b/dense.go @@ -0,0 +1,174 @@ +// Dense solutions. +// +// Copyright (C) 2020 Juan Marín Noguera +// +// This file is part of Solvned. +// +// Solvned is free software: you can redistribute it and/or modify it under the +// terms of the GNU Lesser General Public License as published by the Free +// Software Foundation, either version 3 of the License, or (at your option) any +// later version. +// +// Solvned is distributed in the hope that it will be useful, but WITHOUT ANY +// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR +// A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more +// details. +// +// You should have received a copy of the GNU Lesser General Public License +// along with Solvned. If not, see . + +package mned + +// A DenseSolution is a structure with precalculated points of a solution within +// a given interval. It allows getting the value of any point within that +// interval via interpolation. The zero value is meaningless. +type DenseSolution struct { + points []Point // This can't be empty. Points are ordered increasingly. + interp Interpolator +} + +// Get the earliest time of the points in the solution. +func (d *DenseSolution) Start() float64 { + return d.points[0].Time +} + +// Get the latest time of the points in the solution. +func (d *DenseSolution) End() float64 { + return d.points[len(d.points)-1].Time +} + +func (d *DenseSolution) search(time float64, min int, max int) []float64 { + if min == max { + return d.points[min].Value + } + if max-min == 1 { + switch time { + case d.points[min].Time: + return d.points[min].Clone().Value + case d.points[max].Time: + return d.points[max].Clone().Value + } + return d.interp.FindValue(&d.points[min], &d.points[max], time) + } + mid := (max - min) / 2 + midtime := d.points[mid].Time + switch { + case time == midtime: + return d.points[mid].Value + case time > midtime: + return d.search(time, mid, max) + default: + return d.search(time, min, mid) + } +} + +// Get the value for a point in the solution with a given time between +// `d.Start()` and `d.End()`. If the second return value is false, the point +// was not in the interval and the first return value is meaningless. +func (d *DenseSolution) Get(time float64) ([]float64, bool) { + if time < d.Start() { + return d.points[0].Value, false + } + if time > d.End() { + return d.points[len(d.points)-1].Value, false + } + return d.search(time, 0, len(d.points)-1), true +} + +// Return the list of all explicitly-calculated points. This slice *MUST NOT* +// be modified, nor any of the points it contains. +func (d *DenseSolution) InnerPoints() []Point { + return d.points +} + +// Like CacheSolver.PointCoords but for a DenseSolution. +func (d *DenseSolution) PointCoords() [][]float64 { + size := len(d.points[0].Value) + result := make([][]float64, size+1) + for i := range result { + result[i] = make([]float64, len(d.points)) + } + for i, p := range d.points { + result[0][i] = p.Time + for j, v := range p.Value { + result[j+1][i] = v + } + } + return result +} + +// Solve an initial value problem in an interval `[start, end]` with the given +// method. The points are calculated eagerly and, when a point is requested not +// given by the method, interpolation is used with the given interpolator. +// +// If the second return value is false, the method didn't get to generate any +// point and the first return value is meaningless. If it's true, the first +// return value contains a solution for an interval that contains at most two +// explicitly-stored points outside the interval requested (the two extremes) +// but which might not contain the whole interval requested if the method +// stopped. +func DenseSolve( + method Method, + ivp *IVP, + start float64, + end float64, + interp Interpolator, +) (DenseSolution, bool) { + forward := method.Forward(ivp) + backward := method.Backward(ivp) + points := make([]Point, 0) + current := &ivp.Start + pendingAddInitial := true + var prev Point + var ok bool + + if end < current.Time { // Special case: end is on the left + for end < current.Time { + prev = current.Clone() + if current, ok = backward.Next(); !ok { + return DenseSolution{}, false + } + } + if current.Time < end { + points = append(points, prev) + } + points = append(points, current.Clone()) + pendingAddInitial = false + } + for start < current.Time { // Move backwards until the start + if current, ok = backward.Next(); !ok { + break + } + points = append(points, current.Clone()) + } + size := len(points) // Reverse backwards points + for i := 0; i < size/2; i++ { + points[i], points[size-i] = points[size-i], points[i] + } + + current = &ivp.Start + if start > current.Time { // Special case: start is on the right + for start > current.Time { + prev = current.Clone() + if current, ok = forward.Next(); !ok { + return DenseSolution{}, false + } + } + if current.Time > start { + points = append(points, prev) + } + points = append(points, current.Clone()) + pendingAddInitial = false + } + if pendingAddInitial { + points = append(points, ivp.Start) + } + for end > current.Time { // Move forwards until the end + if current, ok = forward.Next(); !ok { + break + } + points = append(points, current.Clone()) + } + + return DenseSolution{points: points, interp: interp}, true +} diff --git a/event.go b/event.go new file mode 100644 index 0000000..c4a687b --- /dev/null +++ b/event.go @@ -0,0 +1,95 @@ +// Events mechanism for early iteration ending. +// +// Copyright (C) 2020 Juan Marín Noguera +// +// This file is part of Solvned. +// +// Solvned is free software: you can redistribute it and/or modify it under the +// terms of the GNU Lesser General Public License as published by the Free +// Software Foundation, either version 3 of the License, or (at your option) any +// later version. +// +// Solvned is distributed in the hope that it will be useful, but WITHOUT ANY +// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR +// A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more +// details. +// +// You should have received a copy of the GNU Lesser General Public License +// along with Solvned. If not, see . + +package mned + +import "math" + +// An Event is a condition that can happen in a point in the solution of an +// initial value problem and an associated action to take. The event happens +// at the point of the solution where a given function changes its sign. +// +// The Cross function is a continuous function whose zeroes are the points +// where the event happens. The Tolerance is the margin of error allowed; the +// maximum absolute value of `Cross(p)` such that `p` is considered to be close +// enough to an event to be passed to Action. The Action is to be called when +// an occurrence of the event is found; it can use the point in some way and it +// returns a boolean indicating whether the calculation should continue. +type Event struct { + Cross func(*Point) float64 + Tolerance float64 + Action func(*Point) bool +} + +// If e.Cross has a zero between p1 and p2, find such a zero or a point `p` +// close enough to the zero that `|e.Cross(p)| < e.Tolerance`. To get the +// intermediante points, the given interpolator is used. +func (e *Event) FindPoint(i Interpolator, p1 *Point, p2 *Point) Point { + var min, max, mid Point + if p1.Time < p2.Time { + min, max = *p1, *p2 + } else { + max, min = *p1, *p2 + } + downwards := e.Cross(p1) > 0 + for { + mid.Time = (min.Time + max.Time) / 2 + mid.Value = i.FindValue(p1, p2, mid.Time) + value := e.Cross(&mid) + if math.Abs(value) < e.Tolerance { + return mid + } + if downwards { + if value > 0 { + min = mid + } else { + max = mid + } + } else { + if value > 0 { + max = mid + } else { + min = mid + } + } + } +} + +type requiredAction struct { + index int + point Point +} + +type requiredActions []requiredAction + +func (r requiredActions) Len() int { + return len(r) +} + +func (r requiredActions) Less(i, j int) bool { + return r[i].point.Time < r[j].point.Time +} + +func (r requiredActions) Swap(i, j int) { + r[i], r[j] = r[j], r[i] +} + +func differentSign(f1 float64, f2 float64) bool { + return (f1 <= 0 && f2 >= 0) || (f1 >= 0 && f2 <= 0) +} diff --git a/examples/arenstorf.go b/examples/arenstorf.go new file mode 100644 index 0000000..6506a6a --- /dev/null +++ b/examples/arenstorf.go @@ -0,0 +1,74 @@ +// Example visualization of the Arenstorf orbit. +// +// Copyright (C) 2020 Juan Marín Noguera +// +// This file is part of Solvned. +// +// Solvned is free software: you can redistribute it and/or modify it under the +// terms of the GNU Lesser General Public License as published by the Free +// Software Foundation, either version 3 of the License, or (at your option) any +// later version. +// +// Solvned is distributed in the hope that it will be useful, but WITHOUT ANY +// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR +// A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more +// details. +// +// You should have received a copy of the GNU Lesser General Public License +// along with Solvned. If not, see . + +package main + +import ( + "fmt" + "github.com/Arafatk/glot" + "github.com/JwanMan/mned" + "github.com/JwanMan/mned/ivp" + "github.com/JwanMan/mned/method" + "os" +) + +func main() { + var backed bool = false + problem := ivp.ArenstorfOrbit() + solution := mned.CacheSolve( + method.RKFehlberg(0.01, 0.01, 1e-7, 0.1), + &problem, + mned.HermiteForIVP(&problem), + mned.Event{ + Cross: func(p *mned.Point) float64 { + return p.Value[1] + }, + Tolerance: 0.01, + Action: func(p *mned.Point) bool { + if p.Value[0] < 0 { + backed = true + } + return true + }, + }, + mned.Event{ + Cross: func(p *mned.Point) float64 { + return p.Value[1] + }, + Tolerance: 0.01, + Action: func(p *mned.Point) bool { + return p.Value[0] <= 0.9 || !backed + }, + }, + ) + solution.StepToEnd() + coords := solution.PointCoords() + fmt.Printf("Took %v steps.\n", len(coords[0])) + plot, err := glot.NewPlot(2, true, false) + if err != nil { + fmt.Printf("Couldn't create plot: %v\n", err) + os.Exit(1) + } + if err := plot.AddPointGroup("Orbit", "lines", [][]float64{ + coords[1], coords[2], + }); err != nil { + fmt.Printf("Couldn't draw the points: %v\n", err) + os.Exit(1) + } +} diff --git a/examples/arenstorf.png b/examples/arenstorf.png new file mode 100644 index 0000000..22cf99b Binary files /dev/null and b/examples/arenstorf.png differ diff --git a/examples/arenstorf.txt b/examples/arenstorf.txt new file mode 100644 index 0000000..b9d3337 --- /dev/null +++ b/examples/arenstorf.txt @@ -0,0 +1 @@ +Took 357 steps. diff --git a/examples/exp.go b/examples/exp.go new file mode 100644 index 0000000..1021be6 --- /dev/null +++ b/examples/exp.go @@ -0,0 +1,74 @@ +// Example of trapezium method for exponential solutions. +// +// Copyright (C) 2020 Juan Marín Noguera +// +// This file is part of Solvned. +// +// Solvned is free software: you can redistribute it and/or modify it under the +// terms of the GNU Lesser General Public License as published by the Free +// Software Foundation, either version 3 of the License, or (at your option) any +// later version. +// +// Solvned is distributed in the hope that it will be useful, but WITHOUT ANY +// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR +// A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more +// details. +// +// You should have received a copy of the GNU Lesser General Public License +// along with Solvned. If not, see . + +package main + +import ( + "fmt" + "math" + "github.com/JwanMan/mned" + "github.com/JwanMan/mned/method" + "github.com/JwanMan/mned/mstep" +) + +type methodDescriptor struct { + name string + m mned.Method +} + +const step = 0.01 +const tol = 0.01 + +var df mstep.YDerivative = func(t float64, v float64) float64 { return -1 } +var methods = [5]methodDescriptor{ + {name: "Euler", m: method.Euler(0.01)}, + {name: "RK4", m: method.RK4(0.1)}, + {name: "Implicit Euler", m: mstep.ImplicitEuler(0.01, 0.01, df)}, + {name: "Trapezium", m: mstep.Trapezium(0.01, 0.01, df)}, + {name: "RKF", m: method.RKFehlberg(0.01, 0.01, 1e-7, 0.1)}, +} + +func main() { + problem := mned.IVP{ + Derivative: func(p mned.Point) ([]float64, bool) { + return []float64{-p.Value[0]}, true + }, + Start: mned.Point{Time: 0, Value: []float64{1}}, + } + for _, m := range methods { + solution, _ := mned.DenseSolve( + m.m, &problem, 0, 10, mned.HermiteForIVP(&problem), + ) + points := solution.InnerPoints() + maxerr := float64(0) + for _, p := range points { + realVal := math.Exp(-p.Time) + err := math.Abs(p.Value[0] / realVal) + if err < 1 { + err = 1 / err + } + err -= 1 + if err > maxerr { + maxerr = err + } + } + fmt.Printf("%v method got %v steps, max rel err = %v.\n", + m.name, len(points), maxerr) + } +} diff --git a/examples/exp.txt b/examples/exp.txt new file mode 100644 index 0000000..3c3eddf --- /dev/null +++ b/examples/exp.txt @@ -0,0 +1,5 @@ +Euler method got 1002 steps, max rel err = 0.05167716448732684. +RK4 method got 102 steps, max rel err = 9.149053072698976e-06. +Implicit Euler method got 1002 steps, max rel err = 0.05097553729673554. +Trapezium method got 1002 steps, max rel err = 8.342139748207522e-05. +RKF method got 103 steps, max rel err = 1.5683599816629368e-06. diff --git a/examples/hooke.go b/examples/hooke.go new file mode 100644 index 0000000..36229c1 --- /dev/null +++ b/examples/hooke.go @@ -0,0 +1,77 @@ +// Example visualization of Hooke's law. +// +// Copyright (C) 2020 Juan Marín Noguera +// +// This file is part of Solvned. +// +// Solvned is free software: you can redistribute it and/or modify it under the +// terms of the GNU Lesser General Public License as published by the Free +// Software Foundation, either version 3 of the License, or (at your option) any +// later version. +// +// Solvned is distributed in the hope that it will be useful, but WITHOUT ANY +// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR +// A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more +// details. +// +// You should have received a copy of the GNU Lesser General Public License +// along with Solvned. If not, see . + +package main + +import ( + "fmt" + "github.com/Arafatk/glot" + "mned" + "mned/ivp" + "mned/method" +) + +func spring1(speed float64) { + spring := ivp.HookeSpring{ + Mass: 1, + Spring: 1.5, + Length: 0.7, + Friction: 0.3, + Amplitude: 0.4, + AngleSpeed: speed, + } + problem := spring.ToIVP() + plot, err := glot.NewPlot(2, true, false) + if err != nil { + fmt.Printf("Creating the plot: %v\n", err) + return + } + if err := plot.SetXLabel( + fmt.Sprintf("Time (s) -- w = %v", speed), + ); err != nil { + fmt.Printf("Drawing the X label: %v\n", err) + return + } + + solution, _ := mned.DenseSolve( + method.RKFehlberg(0.0001, 0.01, 1e-7, 0.5), + &problem, + 0, + 40, + mned.HermiteForIVP(&problem), + ) + coord := solution.PointCoords() + if err := plot.AddPointGroup("Position (m)", "lines", [][]float64{ + coord[0], coord[1], + }); err != nil { + fmt.Printf("Drawing the position: %v\n", err) + } + if err := plot.AddPointGroup("Velocity (m/s)", "lines", [][]float64{ + coord[0], coord[2], + }); err != nil { + fmt.Printf("Drawing the velocity: %v\n", err) + } +} + +func main() { + spring1(0) + spring1(1.3) + spring1(2.4) + spring1(3.5) +} diff --git a/examples/hooke.png b/examples/hooke.png new file mode 100644 index 0000000..f8bbd32 Binary files /dev/null and b/examples/hooke.png differ diff --git a/examples/hooke_check.go b/examples/hooke_check.go new file mode 100644 index 0000000..c51339c --- /dev/null +++ b/examples/hooke_check.go @@ -0,0 +1,91 @@ +// Example comparing the formula and numeric approximation for Hooke's law +// +// Copyright (C) 2020 Juan Marín Noguera +// +// This file is part of Solvned. +// +// Solvned is free software: you can redistribute it and/or modify it under the +// terms of the GNU Lesser General Public License as published by the Free +// Software Foundation, either version 3 of the License, or (at your option) any +// later version. +// +// Solvned is distributed in the hope that it will be useful, but WITHOUT ANY +// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR +// A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more +// details. +// +// You should have received a copy of the GNU Lesser General Public License +// along with Solvned. If not, see . + +package main + +import ( + "fmt" + "math" + "github.com/JwanMan/mned" + "github.com/JwanMan/mned/ivp" + "github.com/JwanMan/mned/method" +) + +type methodDescriptor struct { + name string + m mned.Method +} + +const step = 0.01 +const tol = 0.01 + +var methods = [4]methodDescriptor{ + {name: "Euler", m: method.Euler(step)}, + {name: "RK4", m: method.RK4(step)}, + {name: "Adaptive RK4", m: method.AdaptiveRK4(tol, step, 1e-6, 1)}, + {name: "RK-Fehlberg", m: method.RKFehlberg(tol, step, 1e-6, 1)}, +} + +func normdiff(p1 []float64, p2 []float64) float64 { + var norm float64 = 0 + for i, v := range p1 { + norm += (v - p2[i]) * (v - p2[i]) + } + return norm +} + +func maxerr(points []mned.Point, sol func(float64) []float64) float64 { + var max float64 = 0 + for _, p := range points { + actual := sol(p.Time) + norm := normdiff(actual, p.Value) + if norm > max { + max = norm + } + } + return max +} + +func main() { + spring := ivp.HookeSpring{ + Mass: 1, + Spring: 0.7, + X0: 1, + Length: 0.7, + } + problem := spring.ToIVP() + amplitude := spring.X0 - spring.Length + sqratio := math.Sqrt(spring.Spring / spring.Mass) + realSol := func(t float64) []float64 { + return []float64{ + amplitude*math.Cos(sqratio*t) + spring.Length, + -amplitude * sqratio * math.Sin(sqratio*t), + } + } + interp := mned.HermiteForIVP(&problem) + + for _, m := range methods { + solution, _ := mned.DenseSolve( + m.m, &problem, 0, 40, interp, + ) + points := solution.InnerPoints() + fmt.Printf("%v: Error of %v in %v steps.\n", + m.name, maxerr(points, realSol), len(points)) + } +} diff --git a/examples/hooke_check.txt b/examples/hooke_check.txt new file mode 100644 index 0000000..8081664 --- /dev/null +++ b/examples/hooke_check.txt @@ -0,0 +1,4 @@ +Euler: Error of 0.0017828128093108954 in 4001 steps. +RK4: Error of 1.6399864728279775e-19 in 4001 steps. +Adaptive RK4: Error of 7.496803458130451e-07 in 45 steps. +RK-Fehlberg: Error of 7.835831771324868e-05 in 45 steps. diff --git a/examples/parabolic.go b/examples/parabolic.go new file mode 100644 index 0000000..f45872b --- /dev/null +++ b/examples/parabolic.go @@ -0,0 +1,110 @@ +// Example visualization of parabolic throw with air friction. +// +// Copyright (C) 2020 Juan Marín Noguera +// +// This file is part of Solvned. +// +// Solvned is free software: you can redistribute it and/or modify it under the +// terms of the GNU Lesser General Public License as published by the Free +// Software Foundation, either version 3 of the License, or (at your option) any +// later version. +// +// Solvned is distributed in the hope that it will be useful, but WITHOUT ANY +// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR +// A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more +// details. +// +// You should have received a copy of the GNU Lesser General Public License +// along with Solvned. If not, see . + +package main + +import ( + "fmt" + "github.com/Arafatk/glot" + "github.com/JwanMan/mned" + "github.com/JwanMan/mned/ivp" + "github.com/JwanMan/mned/method" + "os" +) + +func addPlot( + plot *glot.Plot, name string, arrays [][]float64, idx int, idy int, +) { + if err := plot.AddPointGroup( + name, + "lines", + [][]float64{arrays[idx], arrays[idy]}); err != nil { + fmt.Printf("Error adding %v: %w\n", name, err) + os.Exit(1) + } +} + +func main() { + var end mned.Point + var err error + + throw := ivp.ParabolicThrow{ + Height: 300, + V0: [2]float64{100, 0}, + Mass: 1, + Resistance: 0.01, + Gravity: ivp.EARTH_GRAVITY, + } + problem := throw.ToIVP() + solution := mned.CacheSolve( + method.Euler(0.01), + &problem, + mned.LinearInterpolator{}, + mned.Event{ + Cross: func(p *mned.Point) float64 { + return p.Value[1] + }, + Tolerance: 0.01, + Action: func(p *mned.Point) bool { + end = *p + return false + }, + }) + solution.StepToEnd() + fmt.Printf( + "Reached ground after %v s at %v m, with (%v, %v) m/s.\n", + end.Time, + end.Value[0], + end.Value[2], + end.Value[3], + ) + + arrays := solution.PointCoords() + byTime, err := glot.NewPlot(2, true, false) + if err != nil { + fmt.Printf("Creating plot window: %w\n") + os.Exit(1) + } + defer byTime.Close() + if err := byTime.SetXLabel("Time (s)"); err != nil { + fmt.Printf("Adding X label: %w\n") + os.Exit(1) + } + addPlot(byTime, "X position (m)", arrays, 0, 1) + addPlot(byTime, "Y position (m)", arrays, 0, 2) + addPlot(byTime, "X speed (m/s)", arrays, 0, 3) + addPlot(byTime, "Y speed (m/s)", arrays, 0, 4) + + byTraj, err := glot.NewPlot(2, true, false) + if err != nil { + fmt.Printf("Creating plot window: %w\n") + os.Exit(1) + } + defer byTraj.Close() + if err := byTraj.SetXLabel("X coordinate"); err != nil { + fmt.Printf("Adding X label: %w\n") + os.Exit(1) + } + if err := byTraj.SetYLabel("Y coordinate"); err != nil { + fmt.Printf("Adding Y label: %w\n") + os.Exit(1) + } + addPlot(byTraj, "Trajectory (m)", arrays, 1, 2) + addPlot(byTraj, "Speed trajectory (m/s)", arrays, 3, 4) +} diff --git a/examples/parabolic.png b/examples/parabolic.png new file mode 100644 index 0000000..a8dffb6 Binary files /dev/null and b/examples/parabolic.png differ diff --git a/examples/parabolic_check.go b/examples/parabolic_check.go new file mode 100644 index 0000000..4189fe9 --- /dev/null +++ b/examples/parabolic_check.go @@ -0,0 +1,108 @@ +// Example comparing the formula and numeric approximation for parabolic throw. +// +// Copyright (C) 2020 Juan Marín Noguera +// +// This file is part of Solvned. +// +// Solvned is free software: you can redistribute it and/or modify it under the +// terms of the GNU Lesser General Public License as published by the Free +// Software Foundation, either version 3 of the License, or (at your option) any +// later version. +// +// Solvned is distributed in the hope that it will be useful, but WITHOUT ANY +// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR +// A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more +// details. +// +// You should have received a copy of the GNU Lesser General Public License +// along with Solvned. If not, see . + +package main + +import ( + "fmt" + "math" + "github.com/JwanMan/mned" + "github.com/JwanMan/mned/ivp" + "github.com/JwanMan/mned/method" +) + +type methodDescriptor struct { + name string + m mned.Method +} + +const step = 0.01 + +var methods = [3]methodDescriptor{ + {name: "Euler", m: method.Euler(step)}, + {name: "Modified Euler", m: method.ModifiedEuler(step)}, + {name: "RK4", m: method.RK4(step)}, +} + +func normdiff(p1 []float64, p2 []float64) float64 { + var norm float64 = 0 + for i, v := range p1 { + norm += (v - p2[i]) * (v - p2[i]) + } + return norm +} + +func maxerr(points []mned.Point, sol func(float64) []float64) float64 { + var max float64 = 0 + for _, p := range points { + actual := sol(p.Time) + norm := normdiff(actual, p.Value) + if norm > max { + max = norm + } + } + return max +} + +func main() { + // Here we test the precision of different methods against the known + // problem of parabolic throw without air resistance. The problem is + // x''(t) = -gj, so x(t) = -g(t^2/2)j + v0*t + hj, thus + // x(t)[1] = 0 when t = (v0[1] + sqrt(v0[1]^2 + 2hg)) / g + var fallen float64 + throw := ivp.ParabolicThrow{ + Height: 300, + V0: [2]float64{100, 0}, + Mass: 1, + Gravity: ivp.EARTH_GRAVITY, + } + problem := throw.ToIVP() + realSolution := func(t float64) []float64 { + return []float64{ + throw.V0[0] * t, + throw.Height + throw.V0[1]*t - throw.Gravity*t*t/2, + throw.V0[0], + throw.V0[1] - throw.Gravity*t, + } + } + realEnd := (throw.V0[1] + math.Sqrt( + throw.V0[1]*throw.V0[1]+2*throw.Height*throw.Gravity, + )) / throw.Gravity * throw.V0[0] + + for _, method := range methods { + solution := mned.CacheSolve( + method.m, &problem, mned.LinearInterpolator{}, + mned.Event{ + Cross: func(p *mned.Point) float64 { + return p.Value[1] + }, + Tolerance: 0.01, + Action: func(p *mned.Point) bool { + fallen = p.Value[0] + return false + }, + }) + + solution.StepToEnd() + soldiff := maxerr(solution.ForwardPoints(), realSolution) + enddiff := math.Abs(fallen - realEnd) + fmt.Printf("%v: Point diff @ %v, end diff @ %v.\n", + method.name, soldiff, enddiff) + } +} diff --git a/examples/parabolic_check.txt b/examples/parabolic_check.txt new file mode 100644 index 0000000..373b723 --- /dev/null +++ b/examples/parabolic_check.txt @@ -0,0 +1,3 @@ +Euler: Point diff @ 0.1470067554878007, end diff @ 0.49737244764885347. +Modified Euler: Point diff @ 1.9920977079185658e-22, end diff @ 0.0026275523512602206. +RK4: Point diff @ 1.9926361270567946e-22, end diff @ 0.0026275523512602206. diff --git a/interpolate.go b/interpolate.go new file mode 100644 index 0000000..075a534 --- /dev/null +++ b/interpolate.go @@ -0,0 +1,97 @@ +// Interpolation routines. +// +// Copyright (C) 2020 Juan Marín Noguera +// +// This file is part of Solvned. +// +// Solvned is free software: you can redistribute it and/or modify it under the +// terms of the GNU Lesser General Public License as published by the Free +// Software Foundation, either version 3 of the License, or (at your option) any +// later version. +// +// Solvned is distributed in the hope that it will be useful, but WITHOUT ANY +// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR +// A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more +// details. +// +// You should have received a copy of the GNU Lesser General Public License +// along with Solvned. If not, see . + +package mned + +// An Interpolator tries to approximate a function based on a finite set of +// known points in it. We'll focus just on interpolation with two points since +// this is precise enough for our purposes, although the interpolator might be +// constructed with some knowledge about the underlying function. +type Interpolator interface { + // Approximate the value of a function in a point `t` given two + // **different** points `p1` and `p2` of the function such that `t` + // is between `p1.Time` and `p2.Time`. + FindValue(p1 *Point, p2 *Point, t float64) []float64 +} + +// A LinearInterpolator is an interpolator that assumes that there's a straight +// line between the two points given. It can be constructed with `new`. +type LinearInterpolator struct{} + +func (li LinearInterpolator) FindValue( + p1 *Point, p2 *Point, t float64, +) []float64 { + // The line between (t1,x1) and (t2,x2) is x1 + (x2-x1)*(t-t1)/(t2-t1). + // If ratio:=(t-t1)/(t2-t1), this is x2*ratio + x1*(1-ratio) + size := len(p1.Value) + result := make([]float64, size) + ratio := (t - p1.Time) / (p2.Time - p1.Time) + for i := 0; i < size; i++ { + result[i] = p1.Value[i]*(1-ratio) + p2.Value[i]*ratio + } + return result +} + +// A HermiteInterpolator is an interpolator that takes into account the +// derivative of the solution function in the end points. +type HermiteInterpolator struct { + F func(Point) ([]float64, bool) // The derivative function. +} + +func (h HermiteInterpolator) FindValue( + p1 *Point, p2 *Point, t float64, +) []float64 { + diff := p2.Time - p1.Time + off := t - p1.Time + d01, ok := h.F(*p1) + if !ok { // Can't derive: fallback + return LinearInterpolator{}.FindValue(p1, p2, t) + } + d23, ok := h.F(*p2) + if !ok { + return LinearInterpolator{}.FindValue(p1, p2, t) + } + d12 := make([]float64, len(d01)) + for i := range d12 { + d12[i] = (p2.Value[i] - p1.Value[i]) / diff + } + d02 := make([]float64, len(d01)) + for i := range d02 { + d02[i] = (d12[i] - d01[i]) / diff + } + d13 := make([]float64, len(d01)) + for i := range d13 { + d13[i] = (d23[i] - d12[i]) / diff + } + d03 := make([]float64, len(d01)) + for i := range d03 { + d03[i] = (d13[i] - d02[i]) / diff + } + result := make([]float64, len(d01)) + for i := range result { + result[i] = p1.Value[i] + off* + (d01[i]+off*(d02[i]+(off-diff)*d03[i])) + } + return result +} + +// Create a Hermite interpolator suitable for the given IVP. +func HermiteForIVP(ivp *IVP) HermiteInterpolator { + return HermiteInterpolator{F: ivp.Derivative} +} diff --git a/ivp/arenstorf.go b/ivp/arenstorf.go new file mode 100644 index 0000000..05b1527 --- /dev/null +++ b/ivp/arenstorf.go @@ -0,0 +1,68 @@ +// Arenstorf orbit problem. +// +// Copyright (C) 2020 Juan Marín Noguera +// +// This file is part of Solvned. +// +// Solvned is free software: you can redistribute it and/or modify it under the +// terms of the GNU Lesser General Public License as published by the Free +// Software Foundation, either version 3 of the License, or (at your option) any +// later version. +// +// Solvned is distributed in the hope that it will be useful, but WITHOUT ANY +// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR +// A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more +// details. +// +// You should have received a copy of the GNU Lesser General Public License +// along with Solvned. If not, see . + +package ivp + +import ( + "math" + "github.com/JwanMan/mned" +) + +var arenstorfOrbit = mned.IVP{ + Derivative: func(p mned.Point) ([]float64, bool) { + const μ = 0.012277471 + const μ2 = 1 - μ + xμ := p.Value[0] + μ + xμ2 := p.Value[0] - μ2 + D1 := xμ*xμ + p.Value[1]*p.Value[1] + D1 *= math.Sqrt(D1) + D2 := xμ2*xμ2 + p.Value[1]*p.Value[1] + D2 *= math.Sqrt(D2) + return []float64{ + p.Value[2], + p.Value[3], + p.Value[0] + 2*p.Value[3] - μ2*xμ/D1 - μ*xμ2/D2, + p.Value[1]*(1-μ2/D1-μ/D2) - 2*p.Value[2], + }, true + }, + Start: mned.Point{ + Time: 0, + Value: []float64{0.994, 0, 0, -2.001585106}, + }, +} + +// The Arenstorf orbit is a stable orbit of a satellite around the Earth taking +// into account the gravity of the Earth and the Moon and considering the +// satellite to have a negligible mass compared to both. The orbit was +// discovered by Richard Arenstorf and is given by the following formulae: +// +// ``` +// μ = 0.012277471, μ' = 1 - μ; +// D1 = ((x+μ)^2 + y^2)^(3/2), D2 = ((x-μ')^2 + y^2)^(3/2); +// x" = x + 2y' - μ'(x+μ)/D1 - μ(x-μ')/D2; +// y" = y - 2x' - μ'y/D1 - μy/D2; +// x(0) = 0.994, y(0) = 0, x'(0) = 0, y'(0) = -2.001585106. +// ``` +// +// The orbit has three big loops at one side and a small one at the other, and +// varying μ changes the number of loops. Because of the unstability, this +// problem is often used to test different ODE solving methods. +func ArenstorfOrbit() mned.IVP { + return arenstorfOrbit.Clone() +} diff --git a/ivp/spring.go b/ivp/spring.go new file mode 100644 index 0000000..4c8798d --- /dev/null +++ b/ivp/spring.go @@ -0,0 +1,79 @@ +// Spring movement problem based on Hooke's law. +// +// Copyright (C) 2020 Juan Marín Noguera +// +// This file is part of Solvned. +// +// Solvned is free software: you can redistribute it and/or modify it under the +// terms of the GNU Lesser General Public License as published by the Free +// Software Foundation, either version 3 of the License, or (at your option) any +// later version. +// +// Solvned is distributed in the hope that it will be useful, but WITHOUT ANY +// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR +// A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more +// details. +// +// You should have received a copy of the GNU Lesser General Public License +// along with Solvned. If not, see . + +package ivp + +import ( + "math" + "github.com/JwanMan/mned" +) + +// A HookeSpring is a spring moving following Hooke's law. +// +// The spring is laid out horizontally over a table and has one end attached to +// a body which receives a certain friction from the table and such that the +// mass of the spring is negligible compared to the mass of the body. The other +// end of the spring is attached to a harmonic oscillator. The position of the +// body is measured as the (signed) distance to the rest position of the +// oscillator, and the velocity is measured in terms of such position. +// +// Only Mass and Spring are required, although the result of only specifying +// that would be really boring. Note also that, since we're following Hooke's +// law, the length of the spring can get negative. +type HookeSpring struct { + Mass float64 // Mass of the body (kg). + Spring float64 // Hooke's constant of the spring (kg/s^2). + X0 float64 // Initial position of the body (m). + V0 float64 // Initial velocity of the body (m/s). + Length float64 // Length of the spring at rest (m). + Friction float64 // Friction between the body and the table (kg/s). + Amplitude float64 // Amplitude of the oscillator's force (N). + AngleSpeed float64 // Angle speed of the oscillator (rad/s). +} + +// Get the IVP associated to the problem. The solution values have the form +// (length, speed), time starts at 0, and the harmonic oscillator starts at +// its natural position. +// +// The ODE is `m * x''(t) = A*sin(wt) - k*(x(t)-l) - rx'(t)`, where m is the +// mass of the body, x is its position, A is the amplitude of the oscillator, +// w is its angular velocity, k is the spring's constant, l is its length at +// rest, and r is the friction coefficient. +func (h *HookeSpring) ToIVP() mned.IVP { + adjSpring := h.Spring / h.Mass + adjFriction := h.Friction / h.Mass + adjAmplitude := h.Amplitude / h.Mass + w := h.AngleSpeed + length := h.Length + + return mned.IVP{ + Derivative: func(p mned.Point) ([]float64, bool) { + return []float64{ + p.Value[1], + adjAmplitude*math.Sin(w*p.Time) - + adjSpring*(p.Value[0]-length) - + adjFriction*p.Value[1], + }, true + }, + Start: mned.Point{ + Time: 0, + Value: []float64{h.X0, h.V0}, + }, + } +} diff --git a/ivp/throw.go b/ivp/throw.go new file mode 100644 index 0000000..e91d989 --- /dev/null +++ b/ivp/throw.go @@ -0,0 +1,68 @@ +// Parabolic throw problem with air friction. +// +// Copyright (C) 2020 Juan Marín Noguera +// +// This file is part of Solvned. +// +// Solvned is free software: you can redistribute it and/or modify it under the +// terms of the GNU Lesser General Public License as published by the Free +// Software Foundation, either version 3 of the License, or (at your option) any +// later version. +// +// Solvned is distributed in the hope that it will be useful, but WITHOUT ANY +// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR +// A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more +// details. +// +// You should have received a copy of the GNU Lesser General Public License +// along with Solvned. If not, see . + +package ivp + +import ( + "math" + "github.com/JwanMan/mned" +) + +// Gravity acceleration in the Earth surface +const EARTH_GRAVITY float64 = 9.806 + +// A ParabolicThrow is the result of throwing an object, which we consider to +// be punctual, from a given height and with some initial velocity. We assume +// that Newton mechanics apply and that air resistance has a force proportional +// to the square of the velocity, and opposite to it. +// +// The resulting ODE is `x'' = -gj - (k/m)*x'*|x'|`, where `x` is the position, +// `g` is the gravity acceleration, `j` is the vertical unit vector (upwards), +// `k` is the air resistance, and `m` is the mass of the object. +// +// The comments on fields assume standard SI units, but other units may be used. +type ParabolicThrow struct { + Height float64 // Initial height over the floor (m). + Mass float64 // Mass of the object (kg). *Must* be positive. + V0 [2]float64 // Initial velocity (x,y) (m/s). + Resistance float64 // Air resistance (kg/m). + Gravity float64 // Gravity acceleration (m/s^2). +} + +// Create an IVP from the data on this problem. +// +// Result values have the form (x,y,vx,vy), where (x,y) is the position and +// (vx,vy) is the velocity. The initial x and time are 0. +func (p *ParabolicThrow) ToIVP() mned.IVP { + ratio := p.Resistance / p.Mass + gravity := p.Gravity + return mned.IVP{ + Derivative: func(p mned.Point) ([]float64, bool) { + x := p.Value + air := -ratio * math.Sqrt(x[2]*x[2]+x[3]*x[3]) + return []float64{ + x[2], x[3], air * x[2], air*x[3] - gravity, + }, true + }, + Start: mned.Point{ + Time: 0, + Value: []float64{0, p.Height, p.V0[0], p.V0[1]}, + }, + } +} diff --git a/method/euler.go b/method/euler.go new file mode 100644 index 0000000..8370ae9 --- /dev/null +++ b/method/euler.go @@ -0,0 +1,79 @@ +// Euler method and variants. +// +// Copyright (C) 2020 Juan Marín Noguera +// +// This file is part of Solvned. +// +// Solvned is free software: you can redistribute it and/or modify it under the +// terms of the GNU Lesser General Public License as published by the Free +// Software Foundation, either version 3 of the License, or (at your option) any +// later version. +// +// Solvned is distributed in the hope that it will be useful, but WITHOUT ANY +// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR +// A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more +// details. +// +// You should have received a copy of the GNU Lesser General Public License +// along with Solvned. If not, see . + +package method + +import "github.com/JwanMan/mned" + +func eulerStep(ivp *mned.IVP, h float64) bool { + deriv, ok := ivp.Derivative(ivp.Start) + if !ok { + return false + } + for i := 0; i < len(deriv); i++ { + ivp.Start.Value[i] += h * deriv[i] + } + return true +} + +// Make an instance of the Euler method with the given step, which must be +// positive. For a problem `x'(t)=f(t,x(t)), x(t0)=x0`, a step of the Euler +// method is given by `x(t0+step) = x(t0) + step*f(t0,x(t0))`. +func Euler(step float64) mned.Method { + return mned.FixedStepMethod{Step: step, Method: eulerStep} +} + +// Make an instance of the adaptive version of the Euler method given the +// tolerance, the initial step, and the minimum and maximum steps. +func AdaptiveEuler( + tolerance float64, step float64, hmin float64, hmax float64, +) mned.Method { + return &mned.AdaptiveStepMethod{ + Step: step, + Hmin: hmin, + Hmax: hmax, + Tolerance: tolerance, + Method: eulerStep, + Order: 1, + } +} + +func modEulerStep(ivp *mned.IVP, h float64) bool { + deriv, ok := ivp.Derivative(ivp.Start) + if !ok { + return false + } + point := ivp.Start.Clone() + point.Time += h + for i, v := range deriv { + point.Value[i] += h * v + } + deriv2, ok := ivp.Derivative(point) + for i, _ := range deriv { + ivp.Start.Value[i] += h * (deriv[i] + deriv2[i]) / 2 + } + return true +} + +// Make an instance of the modified Euler method with the given step, which must +// be positive. This is a FixedStepMethod given by +// `x(t+h)=x(t) + h/2 * (f(t,x(t)) + f(t+h,x(t)+h*f(t,x(t))))`. +func ModifiedEuler(step float64) mned.Method { + return mned.FixedStepMethod{Step: step, Method: modEulerStep} +} diff --git a/method/fehlberg.go b/method/fehlberg.go new file mode 100644 index 0000000..3ea14d5 --- /dev/null +++ b/method/fehlberg.go @@ -0,0 +1,208 @@ +// Runge-Kutta-Fehlberg method. +// +// Copyright (C) 2020 Juan Marín Noguera +// +// This file is part of Solvned. +// +// Solvned is free software: you can redistribute it and/or modify it under the +// terms of the GNU Lesser General Public License as published by the Free +// Software Foundation, either version 3 of the License, or (at your option) any +// later version. +// +// Solvned is distributed in the hope that it will be useful, but WITHOUT ANY +// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR +// A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more +// details. +// +// You should have received a copy of the GNU Lesser General Public License +// along with Solvned. If not, see . + +package method + +import ( + "math" + "github.com/JwanMan/mned" +) + +type rkFehlbergStepper struct { + current mned.Point + deriv func(mned.Point) ([]float64, bool) + step float64 + tol float64 + hmin float64 + hmax float64 +} + +func (s *rkFehlbergStepper) Next() (*mned.Point, bool) { + for s.step >= s.hmin { + k1, ok := s.deriv(s.current) + if !ok { + s.step *= 0.5 + continue + } + for i := range k1 { + k1[i] *= s.step + } + + p2 := s.current.Clone() + p2.Time += s.step / 4 + for i := range p2.Value { + p2.Value[i] += k1[i] / 4 + } + k2, ok := s.deriv(p2) + if !ok { + s.step *= 0.5 + continue + } + for i := range k2 { + k2[i] *= s.step + } + + p3 := s.current.Clone() + p3.Time += s.step * 3 / 8 + for i := range p3.Value { + p3.Value[i] += (3*k1[i] + 9*k2[i]) / 32 + } + k3, ok := s.deriv(p3) + if !ok { + s.step *= 0.5 + continue + } + for i := range k3 { + k3[i] *= s.step + } + + p4 := s.current.Clone() + p4.Time += s.step * 12 / 13 + for i := range p4.Value { + p4.Value[i] += + (1932*k1[i] - 7200*k2[i] + 7296*k3[i]) / 2197 + } + k4, ok := s.deriv(p4) + if !ok { + s.step *= 0.5 + continue + } + for i := range k4 { + k4[i] *= s.step + } + + p5 := s.current.Clone() + p5.Time += s.step + for i := range p5.Value { + p5.Value[i] += 439*k1[i]/216 - 8*k2[i] + + 3680*k3[i]/513 - 845*k4[i]/4104 + } + k5, ok := s.deriv(p5) + if !ok { + s.step *= 0.5 + continue + } + for i := range k5 { + k5[i] *= s.step + } + + p6 := s.current.Clone() + p6.Time += s.step / 2 + for i := range p6.Value { + p6.Value[i] += + -8*k1[i]/27 + 2*k2[i] - 3544*k3[i]/2656 + + 1859*k4[i]/4104 - 11*k5[i]/40 + } + k6, ok := s.deriv(p6) + if !ok { + s.step *= 0.5 + continue + } + for i := range k6 { + k6[i] *= s.step + } + + sqerr := float64(0) + // Get \tilde{w}_{i+1}-w_{i+1} directly: + // * (- 16/135 25/216) + // 1/360 + // * (- 6656/12825 1408/2565) + // -128/4275 + // * (- 28561/56430 2197/4104) + // -2197/75240 + // * (- -9/50 -1/5) + // 1/50 + for i := range s.current.Value { + v := k1[i]/360 - 128*k3[i]/4275 - + 2197*k4[i]/75240 + k5[i]/50 + 2*k6[i]/55 + sqerr += v * v + } + err := math.Sqrt(sqerr) + + if err > s.tol { + s.step = mned.AdaptiveUpdateStep( + s.step, s.tol, err, s.hmax, 4, + ) + continue + } + + for i := range s.current.Value { + s.current.Value[i] += 25*k1[i]/216 + 1408*k3[i]/2565 + + 2197*k4[i]/4104 - k5[i]/5 + } + s.current.Time += s.step + s.step = mned.AdaptiveUpdateStep(s.step, s.tol, err, s.hmax, 4) + return &s.current, true + } + return nil, false +} + +type rkFehlbergMethod struct { + step float64 + tol float64 + hmin float64 + hmax float64 +} + +func (m *rkFehlbergMethod) withStep(step float64, ivp *mned.IVP) mned.Stepper { + return &rkFehlbergStepper{ + step: step, + tol: m.tol, + hmin: m.hmin, + hmax: m.hmax, + current: ivp.Start.Clone(), + deriv: ivp.Derivative, + } +} + +func (m *rkFehlbergMethod) Forward(ivp *mned.IVP) mned.Stepper { + return m.withStep(m.step, ivp) +} + +func (m *rkFehlbergMethod) Backward(ivp *mned.IVP) mned.Stepper { + return m.withStep(-m.step, ivp) +} + +// Create an instance of the Runge-Kutta-Fehlberg method, an adaptive, +// fourth-order method, with the given parameter. +// +// This method is given by the following parameters: +// ``` +// k1 = hf(t, x(t)) +// k2 = hf(t+h/4, x(t)+k1/4) +// k3 = hf(t+3h/8, x(t)+3k1/32+9k2/32) +// k4 = hf(t+12h/13, x(t)+1932k1/2197-7200k2/2197+7296k3/2197) +// k5 = hf(t+h, x(t)+439k1/216-8k2+3680k3/513-845k4/4104) +// k6 = hf(t+h/2, x(t)-8k1/27+2k2-3544k3/2656+1859k4/4104-11k5/40) +// e = k1/360-128k3/4275-2197k4/75240+k5/50+2k6/55 +// x(t+h) ~ x(t)+25k1/216+1408k3/2565+2197k4/4104-k5/5 +// ``` +// If `|e|` is greater than `tol`, we let `q:=(h*tol/(2*|e|))^(1/4)`, saturated +// into `[0.1,4]`, and redo the operations with `h*q`, saturating over `hmax` +// and stopping below `hmin`. +func RKFehlberg( + tolerance float64, step float64, hmin float64, hmax float64, +) mned.Method { + return &rkFehlbergMethod{ + step: step, + tol: tolerance, + hmin: hmin, + hmax: hmax, + } +} diff --git a/method/rk.go b/method/rk.go new file mode 100644 index 0000000..391fcb8 --- /dev/null +++ b/method/rk.go @@ -0,0 +1,107 @@ +// Runge-Kutta method and adaptive variant. +// +// Copyright (C) 2020 Juan Marín Noguera +// +// This file is part of Solvned. +// +// Solvned is free software: you can redistribute it and/or modify it under the +// terms of the GNU Lesser General Public License as published by the Free +// Software Foundation, either version 3 of the License, or (at your option) any +// later version. +// +// Solvned is distributed in the hope that it will be useful, but WITHOUT ANY +// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR +// A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more +// details. +// +// You should have received a copy of the GNU Lesser General Public License +// along with Solvned. If not, see . + +package method + +import "github.com/JwanMan/mned" + +func RK4Step(ivp *mned.IVP, h float64) ([]float64, bool) { + k1, ok := ivp.Derivative(ivp.Start) + if !ok { + return nil, false + } + for i, _ := range k1 { + k1[i] *= h + } + + p1 := ivp.Start.Clone() + p1.Time += h / 2 + for i, v := range k1 { + p1.Value[i] += v / 2 + } + k2, ok := ivp.Derivative(p1) + if !ok { + return nil, false + } + for i, _ := range k2 { + k2[i] *= h + } + + p2 := ivp.Start.Clone() + p2.Time += h / 2 + for i, v := range k2 { + p2.Value[i] += v / 2 + } + k3, ok := ivp.Derivative(p2) + if !ok { + return nil, false + } + for i, _ := range k3 { + k3[i] *= h + } + + p3 := ivp.Start.Clone() + p3.Time += h + for i, v := range k3 { + p3.Value[i] += v + } + k4, ok := ivp.Derivative(p3) + if !ok { + return nil, false + } + for i, _ := range k4 { + k4[i] *= h + } + + for i, _ := range ivp.Start.Value { + ivp.Start.Value[i] += (k1[i] + 2*k2[i] + 2*k3[i] + k4[i]) / 6 + } + return k1, true +} + +func rk4Step(ivp *mned.IVP, h float64) bool { + _, ok := RK4Step(ivp, h) + return ok +} + +// Initialize the 4th-order Runge-Kutta method, a.k.a. RK4. This is a fixed step +// method given by the following formulae: +// ``` +// k1 = h*f(t, x(t)) +// k2 = h*f(t+h/2, x(t)+k1/2) +// k3 = h*f(t+h/2, x(t)+k2/2) +// k4 = h*f(t+h, x(t)+k3) +// x(t+h) = x(t) + (k1 + 2*k2 + 2*k3 + k4) / 6 +// ``` +func RK4(step float64) mned.Method { + return mned.FixedStepMethod{Step: step, Method: rk4Step} +} + +func AdaptiveRK4( + tolerance float64, step float64, hmin float64, hmax float64, +) mned.Method { + return &mned.AdaptiveStepMethod{ + Step: step, + Hmin: hmin, + Hmax: hmax, + Tolerance: tolerance, + Method: rk4Step, + Order: 4, + } +} diff --git a/mstep/adams.go b/mstep/adams.go new file mode 100644 index 0000000..e74269c --- /dev/null +++ b/mstep/adams.go @@ -0,0 +1,105 @@ +// Adams stepping methods. +// +// Copyright (C) 2020 Juan Marín Noguera +// +// This file is part of Solvned. +// +// Solvned is free software: you can redistribute it and/or modify it under the +// terms of the GNU Lesser General Public License as published by the Free +// Software Foundation, either version 3 of the License, or (at your option) any +// later version. +// +// Solvned is distributed in the hope that it will be useful, but WITHOUT ANY +// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR +// A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more +// details. +// +// You should have received a copy of the GNU Lesser General Public License +// along with Solvned. If not, see . + +package mstep + +import ( + "github.com/JwanMan/mned" + "github.com/JwanMan/mned/method" +) + +func adamsBashfordStep(steps []Step, h float64) []float64 { + return Implicit(steps, []float64{1}, []float64{ + float64(55) / 24, float64(-59) / 24, + float64(37) / 24, float64(-9) / 24, + }, h) +} + +func adamsMoultonStep(steps []Step, h float64) []float64 { + return Implicit(steps, []float64{0, 1}, []float64{ + float64(9) / 24, float64(19) / 24, + float64(-5) / 24, float64(1) / 24, + }, h) +} + +type adamsBashfordStepper struct { + steps [4]Step + h float64 + problem mned.IVP + init uint8 +} + +func (s *adamsBashfordStepper) NextStep(h float64) (*mned.Point, bool) { + var ok bool + if s.h == 0 { // Error indicator + return nil, false + } + if s.init != 0 { + deriv, ok := method.RK4Step(&s.problem, h) + if !ok { + return nil, false + } + s.problem.Start.Time += h + s.steps[s.init].Deriv = deriv + s.init-- + s.steps[s.init].Value = make([]float64, len(deriv)) + copy(s.steps[s.init].Value, s.problem.Start.Value) + return &s.problem.Start, true + } + + s.steps[0].Deriv, ok = s.problem.Derivative(s.problem.Start) + if !ok { + s.h = 0 + return nil, false + } + s.problem.Start.Value = adamsBashfordStep(s.steps[:], h) + s.problem.Start.Time += h + Shift(s.steps[:]) + s.steps[0].Value = s.problem.Start.Value + return &s.problem.Start, true +} + +func (s *adamsBashfordStepper) Next() (*mned.Point, bool) { + return s.NextStep(s.h) +} + +type adamsBashfordMethod float64 + +func (m adamsBashfordMethod) withStep(h float64, ivp *mned.IVP) mned.Stepper { + result := adamsBashfordStepper{ + h: h, + problem: ivp.Clone(), + init: 3, + } + result.steps[3].Value = make([]float64, len(ivp.Start.Value)) + copy(result.steps[3].Value, ivp.Start.Value) + return &result +} + +func (m adamsBashfordMethod) Forward(ivp *mned.IVP) mned.Stepper { + return m.withStep(float64(m), ivp) +} + +func (m adamsBashfordMethod) Backward(ivp *mned.IVP) mned.Stepper { + return m.withStep(-float64(m), ivp) +} + +func AdamsBashford(step float64) mned.Method { + return adamsBashfordMethod(step) +} diff --git a/mstep/common.go b/mstep/common.go new file mode 100644 index 0000000..2ecff09 --- /dev/null +++ b/mstep/common.go @@ -0,0 +1,53 @@ +// Common code to define stepping methods. +// +// Copyright (C) 2020 Juan Marín Noguera +// +// This file is part of Solvned. +// +// Solvned is free software: you can redistribute it and/or modify it under the +// terms of the GNU Lesser General Public License as published by the Free +// Software Foundation, either version 3 of the License, or (at your option) any +// later version. +// +// Solvned is distributed in the hope that it will be useful, but WITHOUT ANY +// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR +// A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more +// details. +// +// You should have received a copy of the GNU Lesser General Public License +// along with Solvned. If not, see . + +package mstep + +type Step struct { + Value []float64 + Deriv []float64 +} + +func Implicit(steps []Step, a []float64, b []float64, step float64) []float64 { + if len(steps) == 0 { + panic("There must be at least one step!\n") + } + if len(steps) < len(a) || len(steps) < len(b) { + panic("There must be at least as many steps as coefficients.\n") + } + + result := make([]float64, len(steps[0].Value)) + for i, v := range a { + for j := range result { + result[j] += v * steps[i].Value[j] + } + } + for i, v := range b { + for j := range result { + result[j] += step * v * steps[i].Deriv[j] + } + } + return result +} + +func Shift(steps []Step) { + for i := len(steps) - 1; i > 0; i-- { + steps[i] = steps[i-1] + } +} diff --git a/mstep/implicit.go b/mstep/implicit.go new file mode 100644 index 0000000..228a3b8 --- /dev/null +++ b/mstep/implicit.go @@ -0,0 +1,164 @@ +// Implicit and trapezium stepping methods. +// +// Copyright (C) 2020 Juan Marín Noguera +// +// This file is part of Solvned. +// +// Solvned is free software: you can redistribute it and/or modify it under the +// terms of the GNU Lesser General Public License as published by the Free +// Software Foundation, either version 3 of the License, or (at your option) any +// later version. +// +// Solvned is distributed in the hope that it will be useful, but WITHOUT ANY +// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR +// A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more +// details. +// +// You should have received a copy of the GNU Lesser General Public License +// along with Solvned. If not, see . + +package mstep + +import "github.com/JwanMan/mned" + +type YDerivative func(float64, float64) float64 + +type implicitEulerStepper struct { + f func(mned.Point) ([]float64, bool) + dfy YDerivative + step float64 + tolerance float64 + current mned.Point +} + +func (s *implicitEulerStepper) NextStep(h float64) (*mned.Point, bool) { + inner := s.f + time := s.current.Time + s.step + value := s.current.Value[0] + deriv := s.dfy + fn := func(w float64) (float64, bool) { + v, oki := inner(mned.Point{Time: time, Value: []float64{w}}) + if oki { + return w - value - h*v[0], true + } else { + return 0, false + } + } + result, ok := NewtonRoot1D( + s.current.Value[0], s.tolerance, fn, + func(w float64) (float64, bool) { + return 1 - h*deriv(time, w), true + }) + if !ok { + return nil, false + } + s.current.Value[0] = result + s.current.Time += s.step + return &s.current, true +} + +func (s *implicitEulerStepper) Next() (*mned.Point, bool) { + return s.NextStep(s.step) +} + +type implicitEulerMethod struct { + dfy YDerivative + h float64 + tolerance float64 +} + +func (m *implicitEulerMethod) withStep(h float64, ivp *mned.IVP) mned.Stepper { + return &implicitEulerStepper{ + f: ivp.Derivative, + current: ivp.Start.Clone(), + step: h, + dfy: m.dfy, + tolerance: m.tolerance, + } +} + +func (m *implicitEulerMethod) Forward(ivp *mned.IVP) mned.Stepper { + return m.withStep(m.h, ivp) +} + +func (m *implicitEulerMethod) Backward(ivp *mned.IVP) mned.Stepper { + return m.withStep(-m.h, ivp) +} + +func ImplicitEuler(step, tolerance float64, df YDerivative) mned.Method { + return &implicitEulerMethod{h: step, tolerance: tolerance, dfy: df} +} + +type trapeziumStepper struct { + f func(mned.Point) ([]float64, bool) + dfy YDerivative + step float64 + tolerance float64 + current mned.Point + prevf float64 +} + +func (s *trapeziumStepper) NextStep(h float64) (*mned.Point, bool) { + inner := s.f + time := s.current.Time + s.step + pv, ok := s.f(mned.Point{ + Time: s.current.Time, + Value: []float64{s.current.Value[0]}, + }) + if !ok { + return nil, false + } + value := s.current.Value[0] + h/2*pv[0] + fn := func(w float64) (float64, bool) { + v, oki := inner(mned.Point{Time: time, Value: []float64{w}}) + if oki { + return w - value - h/2*v[0], true + } else { + return 0, false + } + } + deriv := s.dfy + result, ok := NewtonRoot1D( + s.current.Value[0], s.tolerance, fn, + func(w float64) (float64, bool) { + return 1 - h/2*deriv(time, w), true + }) + if !ok { + return nil, false + } + s.current.Value[0] = result + s.current.Time += s.step + return &s.current, true +} + +func (s *trapeziumStepper) Next() (*mned.Point, bool) { + return s.NextStep(s.step) +} + +type trapeziumMethod struct { + dfy YDerivative + h float64 + tolerance float64 +} + +func (m *trapeziumMethod) withStep(h float64, ivp *mned.IVP) mned.Stepper { + return &trapeziumStepper{ + f: ivp.Derivative, + current: ivp.Start.Clone(), + step: h, + dfy: m.dfy, + tolerance: m.tolerance, + } +} + +func (m *trapeziumMethod) Forward(ivp *mned.IVP) mned.Stepper { + return m.withStep(m.h, ivp) +} + +func (m *trapeziumMethod) Backward(ivp *mned.IVP) mned.Stepper { + return m.withStep(-m.h, ivp) +} + +func Trapezium(step, tolerance float64, df YDerivative) mned.Method { + return &trapeziumMethod{h: step, tolerance: tolerance, dfy: df} +} diff --git a/mstep/pc.go b/mstep/pc.go new file mode 100644 index 0000000..b6147bb --- /dev/null +++ b/mstep/pc.go @@ -0,0 +1,300 @@ +// Correction routines for optimal step estimation. +// +// Copyright (C) 2020 Juan Marín Noguera +// +// This file is part of Solvned. +// +// Solvned is free software: you can redistribute it and/or modify it under the +// terms of the GNU Lesser General Public License as published by the Free +// Software Foundation, either version 3 of the License, or (at your option) any +// later version. +// +// Solvned is distributed in the hope that it will be useful, but WITHOUT ANY +// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR +// A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more +// details. +// +// You should have received a copy of the GNU Lesser General Public License +// along with Solvned. If not, see . + +package mstep + +import ( + "math" + "github.com/JwanMan/mned" + "github.com/JwanMan/mned/method" +) + +type pcStepper struct { + steps [5]Step + problem mned.IVP + h float64 + init uint8 +} + +func (s *pcStepper) NextStep(h float64) (*mned.Point, bool) { + var ok bool + if s.h == 0 { + return nil, false + } + if s.init != 0 { + deriv, ok := method.RK4Step(&s.problem, h) + if !ok { + return nil, false + } + s.problem.Start.Time += h + s.steps[s.init+1].Deriv = deriv + s.steps[s.init].Value = make([]float64, len(deriv)) + copy(s.steps[s.init].Value, s.problem.Start.Value) + return &s.problem.Start, true + } + + s.steps[1].Deriv, ok = s.problem.Derivative(s.problem.Start) + if !ok { + s.h = 0 + return nil, false + } + s.steps[0].Value = adamsBashfordStep(s.steps[1:], h) + s.problem.Start.Time += h + s.problem.Start.Value = s.steps[0].Value + s.steps[0].Deriv, ok = s.problem.Derivative(s.problem.Start) + if !ok { + s.problem.Start.Time -= h + s.problem.Start.Value = s.steps[1].Value + return nil, false + } + s.steps[0].Value = adamsMoultonStep(s.steps[:], h) + s.problem.Start.Value = s.steps[0].Value + Shift(s.steps[:]) + s.steps[0].Value = nil + s.steps[0].Deriv = nil + s.steps[1].Deriv = nil + return &s.problem.Start, true +} + +func (s *pcStepper) Next() (*mned.Point, bool) { + return s.NextStep(s.h) +} + +type pcMethod float64 + +func (m pcMethod) withStep(h float64, ivp *mned.IVP) mned.Stepper { + result := pcStepper{ + problem: ivp.Clone(), + h: h, + init: 3, + } + result.steps[4].Value = make([]float64, len(ivp.Start.Value)) + copy(result.steps[4].Value, ivp.Start.Value) + return &result +} + +func (m pcMethod) Forward(ivp *mned.IVP) mned.Stepper { + return m.withStep(float64(m), ivp) +} + +func (m pcMethod) Backward(ivp *mned.IVP) mned.Stepper { + return m.withStep(-float64(m), ivp) +} + +func PredictorCorrector(step float64) mned.Method { + return pcMethod(step) +} + +type adaptivePCStepper struct { + steps [5]Step + lasttime float64 + deriv func(mned.Point) ([]float64, bool) + h float64 + hmin float64 + hmax float64 + tolerance float64 + pos uint8 +} + +func (s *adaptivePCStepper) tryStep() (errNorm float64, ok bool) { + s.steps[1].Deriv, ok = s.deriv(mned.Point{ + Time: s.lasttime, + Value: s.steps[1].Value, + }) + if !ok { + return 0, false + } + s.steps[0].Value = adamsBashfordStep(s.steps[1:], s.h) + + s.steps[0].Deriv, ok = s.deriv(mned.Point{ + Time: s.lasttime + s.h, + Value: s.steps[0].Value, + }) + if !ok { + return 0, false + } + bash := s.steps[0].Value + moulton := adamsMoultonStep(s.steps[:], s.h) + s.steps[0].Value = moulton + s.steps[0].Deriv = nil + + sqerr := float64(0) + for i, v := range s.steps[0].Value { + diff := v - bash[i] + sqerr += diff * diff + } + err := (math.Sqrt(sqerr) * 19) / 270 + return err, true +} + +func (s *adaptivePCStepper) adjustStep(err float64) { + q := math.Pow((s.tolerance*math.Abs(s.h))/(2*err), 0.25) + if q < 0.1 { + q = 0.1 + } + if q > 4 { + q = 4 + } + s.h *= q + if math.Abs(s.h) > s.hmax { + if s.h > 0 { + s.h = s.hmax + } else { + s.h = -s.hmax + } + } +} + +func (s *adaptivePCStepper) init() bool { + var ok bool + size := len(s.steps[4].Value) + ivp := mned.IVP{Derivative: s.deriv} + for math.Abs(s.h) >= s.hmin { + ivp.Start.Time = s.lasttime + ivp.Start.Value = make([]float64, size) + copy(ivp.Start.Value, s.steps[4].Value) + + for i := 4; i > 1; i-- { + s.steps[i].Deriv, ok = method.RK4Step(&ivp, s.h) + if !ok { + return false + } + ivp.Start.Time += s.h + s.steps[i-1].Value = make([]float64, size) + copy(s.steps[i-1].Value, ivp.Start.Value) + } + s.lasttime = ivp.Start.Time + errn, ok := s.tryStep() + s.lasttime += s.h + if !ok { + s.lasttime -= 4 * s.h + return false + } + tol := s.tolerance * s.h + if errn > tol { + s.lasttime -= 4 * s.h + s.adjustStep(errn) + continue + } + s.pos = 4 + return true + } + return false +} + +func (s *adaptivePCStepper) tryInit() bool { + for math.Abs(s.h) >= s.hmin { + if s.init() { + Shift(s.steps[:]) + return true + } + for i := 0; i < 4; i++ { + s.steps[i].Value = nil + s.steps[i].Deriv = nil + } + s.steps[4].Deriv = nil + s.h /= 2 + } + return false +} + +func (s *adaptivePCStepper) Next() (*mned.Point, bool) { + var result mned.Point + + if s.pos == 5 { + if !s.tryInit() { + return nil, false + } + } + for s.h >= s.hmin { + if s.pos > 0 { + result.Time = s.lasttime - s.h*float64(s.pos-1) + result.Value = s.steps[s.pos].Value + s.pos-- + return &result, true + } + + errn, ok := s.tryStep() + if !ok { + s.h /= 2 + continue + } + tol := s.tolerance * s.h + if errn > tol { + s.steps[4] = s.steps[1] + s.adjustStep(errn) + if !s.tryInit() { + return nil, false + } + continue + } + s.lasttime += s.h + Shift(s.steps[:]) + result.Time = s.lasttime + result.Value = s.steps[1].Value + if 10*errn < tol { + s.steps[4] = s.steps[1] + s.adjustStep(errn) + s.pos = 5 + } + return &result, true + } + return nil, false +} + +type adaptivePCMethod struct { + h float64 + hmax float64 + hmin float64 + tolerance float64 +} + +func (m *adaptivePCMethod) withStep(h float64, ivp *mned.IVP) mned.Stepper { + result := adaptivePCStepper{ + h: h, + hmax: m.hmax, + hmin: m.hmin, + tolerance: m.tolerance, + lasttime: ivp.Start.Time, + deriv: ivp.Derivative, + pos: 5, + } + result.steps[4].Value = make([]float64, len(ivp.Start.Value)) + copy(result.steps[4].Value, ivp.Start.Value) + return &result +} + +func (m *adaptivePCMethod) Forward(ivp *mned.IVP) mned.Stepper { + return m.withStep(m.h, ivp) +} + +func (m *adaptivePCMethod) Backward(ivp *mned.IVP) mned.Stepper { + return m.withStep(-m.h, ivp) +} + +func AdaptivePredictorCorrector( + tolerance, step, hmin, hmax float64, +) mned.Method { + return &adaptivePCMethod{ + h: step, + hmax: hmax, + hmin: hmin, + tolerance: tolerance, + } +} diff --git a/mstep/root.go b/mstep/root.go new file mode 100644 index 0000000..20758f1 --- /dev/null +++ b/mstep/root.go @@ -0,0 +1,71 @@ +// Root-finding algorithms. +// +// Copyright (C) 2020 Juan Marín Noguera +// +// This file is part of Solvned. +// +// Solvned is free software: you can redistribute it and/or modify it under the +// terms of the GNU Lesser General Public License as published by the Free +// Software Foundation, either version 3 of the License, or (at your option) any +// later version. +// +// Solvned is distributed in the hope that it will be useful, but WITHOUT ANY +// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR +// A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more +// details. +// +// You should have received a copy of the GNU Lesser General Public License +// along with Solvned. If not, see . + +package mstep + +import "math" + +func NewtonRoot1D( + x0, tolerance float64, f, df func(float64) (float64, bool), +) (float64, bool) { + prev := x0 + for { + fprev, ok := f(prev) + if !ok { + return 0, false + } + dfprev, ok := df(prev) + if !ok { + return 0, false + } + off := fprev / dfprev + next := prev - off + if math.Abs(off) <= tolerance { + return next, true + } + prev = next + } +} + +func SecantRoot1D( + x0, x1, tolerance float64, f func(float64) (float64, bool), +) (float64, bool) { + prev := x0 + fprev, ok := f(prev) + if !ok { + return 0, false + } + cur := x1 + fcur, ok := f(cur) + if !ok { + return 0, false + } + for { + off := fcur * (prev - cur) / (fprev - fcur) + next := cur - off + if math.Abs(off) <= tolerance { + return next, true + } + fnext, ok := f(next) + if !ok { + return 0, false + } + prev, fprev, cur, fcur = cur, fcur, next, fnext + } +} diff --git a/step.go b/step.go new file mode 100644 index 0000000..b4a3657 --- /dev/null +++ b/step.go @@ -0,0 +1,324 @@ +// Stepping routines and method interface. +// +// Copyright (C) 2020 Juan Marín Noguera +// +// This file is part of Solvned. +// +// Solvned is free software: you can redistribute it and/or modify it under the +// terms of the GNU Lesser General Public License as published by the Free +// Software Foundation, either version 3 of the License, or (at your option) any +// later version. +// +// Solvned is distributed in the hope that it will be useful, but WITHOUT ANY +// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR +// A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more +// details. +// +// You should have received a copy of the GNU Lesser General Public License +// along with Solvned. If not, see . + +package mned + +import ( + "math" + "sort" +) + +// A Stepper is an iterator over points of the solution of an IVP. ODE solver +// methods generally produce steppers that go from the initial values of the +// problems to lower or greater values. +type Stepper interface { + // Go to the next point and return a pointer to it. If the second value + // is false, there are no more points in the direction the Stepper is + // following and the value of Point should not be trusted. + // + // The point returned might change after the next call to next() on the + // same Stepper. + Next() (*Point, bool) +} + +// Starting in the point init, step until the domain of the function ends or +// one of the events says that the stepping should stop, calculating the +// intermediante points with the given interpolator and, if `callback` is not +// nil, calling the given callback with each new point. If ok is false, the +// iteration ended because the stepper returned false; otherwise index contains +// the index of the event that made the iteration stop. +func StepUntil( + s Stepper, + init Point, + interp Interpolator, + callback func(*Point), + evs ...Event, +) (index int, ok bool) { + var prev Point = init + vs := make([]float64, len(evs)) + for i := 0; i < len(evs); i++ { + vs[i] = evs[i].Cross(&prev) + } + for { + point, ok := s.Next() + if !ok { + return 0, false + } + actions := make(requiredActions, 0) + for i := 0; i < len(evs); i++ { + newVal := evs[i].Cross(point) + if differentSign(vs[i], newVal) { + actions = append(actions, requiredAction{ + index: i, + point: evs[i].FindPoint( + interp, &prev, point, + )}) + } + vs[i] = newVal + } + + sort.Sort(actions) + size := len(actions) + if point.Time < prev.Time { + for i := 0; i < size/2; i++ { + actions[i], actions[size-i] = + actions[size-i], actions[i] + } + } + for _, a := range actions { + if !evs[a.index].Action(&a.point) { + return a.index, true + } + } + + callback(point) + prev = point.Clone() + } +} + +// A ConfigurableStepper is a Stepper which allows specifying the delta of the +// next step; that is, the difference between the time of a point and that of +// the previous point. +type ConfigurableStepper interface { + Stepper + // Like Stepper.Next but the time of the next point is specified to be + // the time of the latest point plus the value given. + NextStep(float64) (*Point, bool) +} + +// A Method is a way of approximating points of the solution to an IVP. It works +// as a factory of Steppers that step forward and backward from the initial +// values of the problem. +type Method interface { + // Make a stepper over the solution of the given IVP that starts at the + // initial values of the problem and goes to ever increasing values for + // the independent variable. + Forward(*IVP) Stepper + // Make a stepper over the solution of the given IVP that starts at the + // initial values of the problem and goes to ever decreasing values for + // the independent variable. + Backward(*IVP) Stepper +} + +type fixedStepper struct { + state IVP + step float64 + method func(*IVP, float64) bool +} + +func (s *fixedStepper) NextStep(h float64) (*Point, bool) { + ok := s.method(&s.state, h) + if !ok { + return nil, false + } + s.state.Start.Time += h + return &s.state.Start, true +} + +func (s *fixedStepper) Next() (*Point, bool) { + return s.NextStep(s.step) +} + +// Create a fixed step method, one where all steps increment or decrement the +// time by the same amount, with the added condition that it shouldn't depend +// on extra state. The Step is the fixed increment/decrement of time, which +// must be positive. The steppers returned by this method are +// `ConfigurableStepper`s. +// +// The Method is a function that takes an IVP and a step and stores in +// IVP.Start.Value the value for IVP.Start.Time+step, while leaving +// IVP.Start.Time untouched. The return value of `method` is true if the +// operation completed successfully or false if there was a problem, such as +// going out of the domain of the IVP.Derivative, in which case the resulting +// state of the IVP is unspecified. +type FixedStepMethod struct { + Step float64 + Method func(*IVP, float64) bool +} + +func (m FixedStepMethod) Forward(ivp *IVP) Stepper { + return &fixedStepper{ + state: ivp.Clone(), + step: m.Step, + method: m.Method} +} + +func (m FixedStepMethod) Backward(ivp *IVP) Stepper { + return &fixedStepper{ + state: ivp.Clone(), + step: -m.Step, + method: m.Method} +} + +type adaptiveStepper struct { + state IVP + step float64 + hmin float64 + hmax float64 + tolerance float64 + method func(*IVP, float64) bool + order uint8 +} + +// Return the new step that should be taken by an AdaptiveStepMethod of the +// given order and with a given tolerance took a step of a given size yielding +// the given amount of error, with max as the maximum size of a step. +func AdaptiveUpdateStep( + step float64, tol float64, err float64, max float64, order uint8, +) float64 { + var q float64 + adjust := tol * math.Abs(step) / (2 * err) + if order == 1 { // Fast path + q = adjust + } else { + q = math.Pow(adjust, 1/float64(order)) + } + switch { + case q < 0.1: + q = 0.1 + case q > 4: + q = 4 + } + newStep := q * step + if math.Abs(newStep) > max { + if newStep > 0 { + newStep = max + } else { + newStep = -max + } + } + return newStep +} + +func (s *adaptiveStepper) adaptStep(err float64) { + s.step = AdaptiveUpdateStep(s.step, s.tolerance, err, s.hmax, s.order) +} + +func (s *adaptiveStepper) getSteps(full []float64, half []float64) bool { + orig := s.state.Start.Value + copy(full, orig) + s.state.Start.Value = full + if !s.method(&s.state, s.step) { + return false + } + copy(half, orig) + s.state.Start.Value = half + if !s.method(&s.state, s.step/2) { + return false + } + if !s.method(&s.state, s.step/2) { + return false + } + s.state.Start.Value = orig + return true +} + +func (s *adaptiveStepper) Next() (*Point, bool) { + size := len(s.state.Start.Value) + full := make([]float64, size) + half := make([]float64, size) + + for s.step >= s.hmin { + for !s.getSteps(full, half) { + if s.step *= 0.5; s.step < s.hmin { + return nil, false + } + } + + orderm1 := 1 / (math.Ldexp(1, int(s.order)) - 1) + multiplier := 1 + orderm1 + sqerr := float64(0) + for i, v := range full { + diff := v - half[i] + sqerr += diff * diff + } + err := math.Sqrt(sqerr) * multiplier + tol := s.tolerance * math.Abs(s.step) + if err > tol { + s.adaptStep(err) + continue + } + + for i, v := range full { + s.state.Start.Value[i] = + multiplier*half[i] - orderm1*v + } + + s.state.Start.Time += s.step + s.adaptStep(err) + return &s.state.Start, true + } + + return nil, false // s.step < s.hmin +} + +// An AdaptiveStepMethod is a kind of adaptive method (that is, one that changes +// the step over time for greater precision and performance) that can be adapted +// from a fixed-step method. +// +// Step, Hmin, and Hmax are the initial, minimum, and maximum step. If the +// step must go lower than Hmin to satisfy the tolerance, the stepper ends; +// if it could get greater than Hmax, it saturates to Hmax. The Tolerance is +// such that the Euclidean norm of the estimated error of any step of size `h` +// must be lower than `h*Tolerance`. +// +// Method works the same way as in `FixedStepMethod` except that its error for a +// given interval of time must be on the order of `O(h^Order)`, with the +// variable `h` being the step. Order must be positive. +// +// The algorithm used is the following, where `(t,x(t))` are the initial values, +// `y(h,t,x(t))` is the result of applying the inner method with step `h` from +// `(t,x(t))`, `k` is the order of the method and `h` is the step. +// 1. First, `F:=y(h,t,x(t))` and `H:=y(h/2,t+h/2,y(h/2,t,x(t)))` are obtained, +// and the error is estimated as `e:=|(1+1/(2^k-1))*(H-F)|`. +// 2. If that error is lower than `|h|*Tolerance`, the step is accepted. +// Otherwise we set `q:=(|h|*Tolerance)/(2*e)`, saturate so that `q` is in +// `[0.1, 4]`, set `h:=q*h`, check bounds on `h` and return to the start. +// 3. The next value is taken to be `(t+h, (2^k*H - H)/(2^k - 1))`, and the +// next step size is calculated by changing `h` as in step 2. +// +// The steppers are *not* `ConfigurableStepper`s. +type AdaptiveStepMethod struct { + Step float64 + Hmin float64 + Hmax float64 + Tolerance float64 + Method func(*IVP, float64) bool + Order uint8 +} + +func (m *AdaptiveStepMethod) toStepper(ivp *IVP, step float64) Stepper { + return &adaptiveStepper{ + state: ivp.Clone(), + step: step, + hmin: m.Hmin, + hmax: m.Hmax, + tolerance: m.Tolerance, + method: m.Method, + order: m.Order, + } +} + +func (m *AdaptiveStepMethod) Forward(ivp *IVP) Stepper { + return m.toStepper(ivp, m.Step) +} + +func (m *AdaptiveStepMethod) Backward(ivp *IVP) Stepper { + return m.toStepper(ivp, -m.Step) +} -- cgit v1.2.3