diff options
Diffstat (limited to 'ivp')
| -rw-r--r-- | ivp/arenstorf.go | 68 | ||||
| -rw-r--r-- | ivp/spring.go | 79 | ||||
| -rw-r--r-- | ivp/throw.go | 68 |
3 files changed, 215 insertions, 0 deletions
diff --git a/ivp/arenstorf.go b/ivp/arenstorf.go new file mode 100644 index 0000000..05b1527 --- /dev/null +++ b/ivp/arenstorf.go @@ -0,0 +1,68 @@ +// Arenstorf orbit problem. +// +// 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 ivp + +import ( + "math" + "github.com/JwanMan/mned" +) + +var arenstorfOrbit = mned.IVP{ + Derivative: func(p mned.Point) ([]float64, bool) { + const μ = 0.012277471 + const μ2 = 1 - μ + xμ := p.Value[0] + μ + xμ2 := p.Value[0] - μ2 + D1 := xμ*xμ + p.Value[1]*p.Value[1] + D1 *= math.Sqrt(D1) + D2 := xμ2*xμ2 + p.Value[1]*p.Value[1] + D2 *= math.Sqrt(D2) + return []float64{ + p.Value[2], + p.Value[3], + p.Value[0] + 2*p.Value[3] - μ2*xμ/D1 - μ*xμ2/D2, + p.Value[1]*(1-μ2/D1-μ/D2) - 2*p.Value[2], + }, true + }, + Start: mned.Point{ + Time: 0, + Value: []float64{0.994, 0, 0, -2.001585106}, + }, +} + +// The Arenstorf orbit is a stable orbit of a satellite around the Earth taking +// into account the gravity of the Earth and the Moon and considering the +// satellite to have a negligible mass compared to both. The orbit was +// discovered by Richard Arenstorf and is given by the following formulae: +// +// ``` +// μ = 0.012277471, μ' = 1 - μ; +// D1 = ((x+μ)^2 + y^2)^(3/2), D2 = ((x-μ')^2 + y^2)^(3/2); +// x" = x + 2y' - μ'(x+μ)/D1 - μ(x-μ')/D2; +// y" = y - 2x' - μ'y/D1 - μy/D2; +// x(0) = 0.994, y(0) = 0, x'(0) = 0, y'(0) = -2.001585106. +// ``` +// +// The orbit has three big loops at one side and a small one at the other, and +// varying μ changes the number of loops. Because of the unstability, this +// problem is often used to test different ODE solving methods. +func ArenstorfOrbit() mned.IVP { + return arenstorfOrbit.Clone() +} diff --git a/ivp/spring.go b/ivp/spring.go new file mode 100644 index 0000000..4c8798d --- /dev/null +++ b/ivp/spring.go @@ -0,0 +1,79 @@ +// Spring movement problem based on 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 ivp + +import ( + "math" + "github.com/JwanMan/mned" +) + +// A HookeSpring is a spring moving following Hooke's law. +// +// The spring is laid out horizontally over a table and has one end attached to +// a body which receives a certain friction from the table and such that the +// mass of the spring is negligible compared to the mass of the body. The other +// end of the spring is attached to a harmonic oscillator. The position of the +// body is measured as the (signed) distance to the rest position of the +// oscillator, and the velocity is measured in terms of such position. +// +// Only Mass and Spring are required, although the result of only specifying +// that would be really boring. Note also that, since we're following Hooke's +// law, the length of the spring can get negative. +type HookeSpring struct { + Mass float64 // Mass of the body (kg). + Spring float64 // Hooke's constant of the spring (kg/s^2). + X0 float64 // Initial position of the body (m). + V0 float64 // Initial velocity of the body (m/s). + Length float64 // Length of the spring at rest (m). + Friction float64 // Friction between the body and the table (kg/s). + Amplitude float64 // Amplitude of the oscillator's force (N). + AngleSpeed float64 // Angle speed of the oscillator (rad/s). +} + +// Get the IVP associated to the problem. The solution values have the form +// (length, speed), time starts at 0, and the harmonic oscillator starts at +// its natural position. +// +// The ODE is `m * x''(t) = A*sin(wt) - k*(x(t)-l) - rx'(t)`, where m is the +// mass of the body, x is its position, A is the amplitude of the oscillator, +// w is its angular velocity, k is the spring's constant, l is its length at +// rest, and r is the friction coefficient. +func (h *HookeSpring) ToIVP() mned.IVP { + adjSpring := h.Spring / h.Mass + adjFriction := h.Friction / h.Mass + adjAmplitude := h.Amplitude / h.Mass + w := h.AngleSpeed + length := h.Length + + return mned.IVP{ + Derivative: func(p mned.Point) ([]float64, bool) { + return []float64{ + p.Value[1], + adjAmplitude*math.Sin(w*p.Time) - + adjSpring*(p.Value[0]-length) - + adjFriction*p.Value[1], + }, true + }, + Start: mned.Point{ + Time: 0, + Value: []float64{h.X0, h.V0}, + }, + } +} diff --git a/ivp/throw.go b/ivp/throw.go new file mode 100644 index 0000000..e91d989 --- /dev/null +++ b/ivp/throw.go @@ -0,0 +1,68 @@ +// Parabolic throw problem 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 ivp + +import ( + "math" + "github.com/JwanMan/mned" +) + +// Gravity acceleration in the Earth surface +const EARTH_GRAVITY float64 = 9.806 + +// A ParabolicThrow is the result of throwing an object, which we consider to +// be punctual, from a given height and with some initial velocity. We assume +// that Newton mechanics apply and that air resistance has a force proportional +// to the square of the velocity, and opposite to it. +// +// The resulting ODE is `x'' = -gj - (k/m)*x'*|x'|`, where `x` is the position, +// `g` is the gravity acceleration, `j` is the vertical unit vector (upwards), +// `k` is the air resistance, and `m` is the mass of the object. +// +// The comments on fields assume standard SI units, but other units may be used. +type ParabolicThrow struct { + Height float64 // Initial height over the floor (m). + Mass float64 // Mass of the object (kg). *Must* be positive. + V0 [2]float64 // Initial velocity (x,y) (m/s). + Resistance float64 // Air resistance (kg/m). + Gravity float64 // Gravity acceleration (m/s^2). +} + +// Create an IVP from the data on this problem. +// +// Result values have the form (x,y,vx,vy), where (x,y) is the position and +// (vx,vy) is the velocity. The initial x and time are 0. +func (p *ParabolicThrow) ToIVP() mned.IVP { + ratio := p.Resistance / p.Mass + gravity := p.Gravity + return mned.IVP{ + Derivative: func(p mned.Point) ([]float64, bool) { + x := p.Value + air := -ratio * math.Sqrt(x[2]*x[2]+x[3]*x[3]) + return []float64{ + x[2], x[3], air * x[2], air*x[3] - gravity, + }, true + }, + Start: mned.Point{ + Time: 0, + Value: []float64{0, p.Height, p.V0[0], p.V0[1]}, + }, + } +} |