1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
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))
}
}
|
