Solvned project from 2020

Note that Go didn't have generics back then.

13 files changed, 547 insertions, 0 deletions

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 <https://www.gnu.org/licenses/>.
+
+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
--- /dev/null
+++ b/examples/arenstorf.png
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 <https://www.gnu.org/licenses/>.
+
+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 <https://www.gnu.org/licenses/>.
+
+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
--- /dev/null
+++ b/examples/hooke.png
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 <https://www.gnu.org/licenses/>.
+
+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 <https://www.gnu.org/licenses/>.
+
+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
--- /dev/null
+++ b/examples/parabolic.png
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 <https://www.gnu.org/licenses/>.
+
+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.