diff options
Diffstat (limited to 'mstep')
| -rw-r--r-- | mstep/adams.go | 105 | ||||
| -rw-r--r-- | mstep/common.go | 53 | ||||
| -rw-r--r-- | mstep/implicit.go | 164 | ||||
| -rw-r--r-- | mstep/pc.go | 300 | ||||
| -rw-r--r-- | mstep/root.go | 71 |
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 + } +} |