Solvned project from 2020

Note that Go didn't have generics back then.

5 files changed, 693 insertions, 0 deletions

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