aboutsummaryrefslogtreecommitdiff
path: root/examples
diff options
context:
space:
mode:
Diffstat (limited to 'examples')
-rw-r--r--examples/arenstorf.go74
-rw-r--r--examples/arenstorf.pngbin0 -> 62397 bytes
-rw-r--r--examples/arenstorf.txt1
-rw-r--r--examples/exp.go74
-rw-r--r--examples/exp.txt5
-rw-r--r--examples/hooke.go77
-rw-r--r--examples/hooke.pngbin0 -> 260239 bytes
-rw-r--r--examples/hooke_check.go91
-rw-r--r--examples/hooke_check.txt4
-rw-r--r--examples/parabolic.go110
-rw-r--r--examples/parabolic.pngbin0 -> 122494 bytes
-rw-r--r--examples/parabolic_check.go108
-rw-r--r--examples/parabolic_check.txt3
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
Binary files differ
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
Binary files differ
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
Binary files differ
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.
generated by cgit v1.2.3 (git 2.46.0) at 2025-10-26 09:58:39 +0000