This post does not present a specific question or problem. I have working code—a tiny simulation of a dynamic system involving some springs in 2d space—and I am hoping for general feedback regarding code idiomaticity, performance, and correctness of the physical model.
My main goals at present are to learn the DifferentialEquations.jl ecosystem and to dust off my knowledge of high school physics.
The system I’m simulating is as follows:
- There are two springs, each with one end anchored at a known point (given by
as), and with known spring constants (ks). The spring’s obey Hooke’s law, and their resting lengths (Ls) are known. - The free ends of the spring are attached to a single point mass (
m). - There is gravity (
g). - We know the initial position of the point mass (
u0).
Here is a diagram (graphic design is my passion, etc.):
Here is my code:
using DifferentialEquations
using Plots
# Picture: We have a point mass affixed to a number of springs
# The springs themselves are anchored to various points in space
# We want to visualize how the point mass moves around in time
mutable struct Params
as::Vector{Vector{Float64}} # Points to which springs are anchored
ks::Vector{Float64} # Spring coefficients
Ls::Vector{Float64} # Spring length at rest
m::Float64 # Mass of the point mass
g::Float64 # Gravity in positive y direction
end
# Euclidean distance between two points
function dist(x1::T, y1::T, x2::U, y2::U) where {T<:Number, U<:Number}
return sqrt((x1 - x2)^2 + (y1 - y2)^2)
end
function f!(du, u, p, t)
x, dx, y, dy = u
du[1] = dx
du[2] = 0
du[3] = dy
du[4] = p.g # Gravity
for (a, k, L) in zip(p.as, p.ks, p.Ls)
L_t = dist(x, y, a[1], a[2]) # Current length of the spring
F = k * (L_t - L) # Hooke's law
# F is the magnitude of the force the spring exerts on the mass
# It is directed from the mass toward the anchor, so:
# Its x component is F cos θ = - F * (x - a[1]) / L_t
du[2] -= F * (x - a[1]) / (L_t * p.m) # = d^2 x
# Its y component is F sin θ = - F * (y - a[2]) / L_t
du[4] -= F * (y - a[2]) / (L_t * p.m) # = d^2 y
end
end
function solveexampleproblem()
as = [[0.0, 0.0], [1.0, 1.0]]
ks = [100.0, 50.0]
Ls = [0.5, 0.5]
m = 10.0
g = -9.81
p = Params(as, ks, Ls, m, g)
u0 = [0.3, 0.0, 0.4, 0.0]
tspan = (0.0, 50.0)
prob = ODEProblem(f!, u0, tspan, p)
@info "Solving example problem ..."
solve(prob, AutoTsit5(Rosenbrock23()))
end
function makeanimation(
sol::OrdinaryDiffEq.ODECompositeSolution,
filename::String
)
@info "Making animation ..."
animate(
sol,
filename,
idxs=(1, 3),
c=:crimson,
label=nothing,
figsize=(500, 800)
)
end
let sol = solveexampleproblem()
makeanimation(sol, "2dHookeAnim.gif")
@info "Done"
end
And here is the gif it outputs:
A general code review, as well as suggestions for follow-up projects, would be more than welcome. Thank you in advance for any feedback.