diff options
Diffstat (limited to 'method')
| -rw-r--r-- | method/euler.go | 79 | ||||
| -rw-r--r-- | method/fehlberg.go | 208 | ||||
| -rw-r--r-- | method/rk.go | 107 |
3 files changed, 394 insertions, 0 deletions
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 <https://www.gnu.org/licenses/>. + +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 <https://www.gnu.org/licenses/>. + +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 <https://www.gnu.org/licenses/>. + +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, + } +} |