Solvned project from 2020

Note that Go didn't have generics back then.

1 files changed, 324 insertions, 0 deletions

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 <https://www.gnu.org/licenses/>.
+
+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)
+}