Solvned project from 2020

Note that Go didn't have generics back then.

1 files changed, 108 insertions, 0 deletions

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)
+	}
+}