# Feedback on a 2d spring simulation

**URL:** <https://discourse.julialang.org/t/feedback-on-a-2d-spring-simulation/88350>\
**Category:** Numerics\
**Created:** [October 6, 2022, 2:34pm UTC](https://discourse.julialang.org/t/feedback-on-a-2d-spring-simulation/88350 "2022-10-06T14:34:41Z")\
**Posts on this page:** 11\
**Page:** 1

<div class="post-metadata">

**Author:** ![maxkapur](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/maxkapur/32/21208_2.png) [@maxkapur](https://discourse.julialang.org/u/maxkapur)\
**Post date:** [October 6, 2022, 2:34pm UTC](https://discourse.julialang.org/t/feedback-on-a-2d-spring-simulation/88350/1 "2022-10-06T14:34:41Z")

</div>

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.):

 ![Schematic](https://global.discourse-cdn.com/julialang/original/3X/9/2/92c45dad71cad5f9339e9163785f71a7f299d89c.png)

Here is my code:

```julia
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:

![2dHookeAnim](https://global.discourse-cdn.com/julialang/original/3X/4/8/48d208defbd04b1055e34c8c27e7860c33baec87.gif)

A general code review, as well as suggestions for follow-up projects, would be more than welcome. Thank you in advance for any feedback.

---

<div class="post-metadata">

**Author:** ![lmiq](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/lmiq/32/18314_2.png) [@lmiq](https://discourse.julialang.org/u/lmiq)\
**Post date:** [October 6, 2022, 7:18pm UTC](https://discourse.julialang.org/t/feedback-on-a-2d-spring-simulation/88350/2 "2022-10-06T19:18:39Z")

</div>

At first sight the code is fine IMO. Some type-annotations seem odd `... where {T<:Number, U<:Number}`.

What can make it perhaps faster and also more idiomatic in some lines is to use `StaticArrays` to represent all types of coordinates. For instance, that would allow you to write:

```julia
using LinearAlgebra: norm
dist(x,y) = norm(x-y)

#and
du -= F * (x - a) / L_t * p.m

```

where `x`, `a` would be static arrays with both coordinates, and the result is also an array. That change would require some reworking of the code, but the equation would be closer to what they look on paper.

Also, you could simplify a bit the code using things like:

```julia
du .= (dx, 0, dy, p.g)

```

instead of assigning each element in a different line, index by index.

---

<div class="post-metadata">

**Author:** ![maxkapur](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/maxkapur/32/21208_2.png) [@maxkapur](https://discourse.julialang.org/u/maxkapur)\
**Post date:** [October 6, 2022, 9:02pm UTC](https://discourse.julialang.org/t/feedback-on-a-2d-spring-simulation/88350/3 "2022-10-06T21:02:11Z")

</div>

Thank you so much! All these tips are helpful.

> [@lmiq](#):
>
> Some type-annotations seem odd `... where {T<:Number, U<:Number}`.

Originally, I had `dist(x1::T, y1::T, x2::T, y2::T) where {T<:Number}` but this caused problems (for some solvers) in autodiff where one of the points was in `Float64`s and the other was in `Dual{Float64}`s. But yeah, a smarter point type from StaticArrays would obviate the need for this workaround.

---

<div class="post-metadata">

**Author:** ![lmiq](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/lmiq/32/18314_2.png) [@lmiq](https://discourse.julialang.org/u/lmiq)\
**Post date:** [October 6, 2022, 9:15pm UTC](https://discourse.julialang.org/t/feedback-on-a-2d-spring-simulation/88350/4 "2022-10-06T21:15:42Z")

</div>

It would be more idiomatic to remove the annotation completely, I think.

---

<div class="post-metadata">

**Author:** ![maxkapur](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/maxkapur/32/21208_2.png) [@maxkapur](https://discourse.julialang.org/u/maxkapur)\
**Post date:** [October 7, 2022, 4:48pm UTC](https://discourse.julialang.org/t/feedback-on-a-2d-spring-simulation/88350/5 "2022-10-07T16:48:34Z")

</div>

OK, I am working on redoing this with a more legible `Point` type, but I am running into missing method for `oneunit(::Type{Any})` error deep inside the solver algorithm.

Here is a minimal example of the problem I am running into:

```julia
using DifferentialEquations

struct Point{T<:Number}
    x::T
    y::T
end

function f!(du::Point, u::Point, _, _)
    du.x = -u.y
    du.y = u.x
end

let
    u0 = Point(1.0, 1.0)
    prob = ODEProblem(f!, u0, (0.0, 1.0))
    sol = solve(prob, Tsit5())
end

```

Error:

```julia
ERROR: MethodError: no method matching oneunit(::Type{Any})
Closest candidates are:
  oneunit(::Type{Union{Missing, T}}) where T at missing.jl:105
  oneunit(::Type{T}) where T at number.jl:358
  oneunit(::T) where T at number.jl:357
  ...
Stacktrace:
  [1] oneunit(#unused#::Type{Any})
    @ Base .\missing.jl:106
  ...

```

(Here I have used a `Point` type with x and y fields for simplicity, but in my actual implementation I am using StaticArrays as suggested and having the same issue.)

Am I using the wrong approach by trying to use a struct as my state variable instead of a vector?

It seems like being able to use a custom struct for `u` and `du` would enable a very expressive coding style, because your state variable could have legible fields corresponding to physical “objects” in the scene, etc.

---

<div class="post-metadata">

**Author:** ![lmiq](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/lmiq/32/18314_2.png) [@lmiq](https://discourse.julialang.org/u/lmiq)\
**Post date:** [October 7, 2022, 6:04pm UTC](https://discourse.julialang.org/t/feedback-on-a-2d-spring-simulation/88350/6 "2022-10-07T18:04:57Z")

</div>

> [@maxkapur](#):
>
> ```julia
> function f!(du::Point, u::Point, _, _)
> du.x = -y
> du.y = x
> end
> 
> ```

`x` and `y` are not defined here, maybe you meant `u.x` and `u.y`?

Apart from that it feels strange to me that `f!` has four parameters while your initial point has two. But I don’t know the interface of `ODEProblem` from the top of my head, so this comment may be irrelevant.

---

<div class="post-metadata">

**Author:** ![maxkapur](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/maxkapur/32/21208_2.png) [@maxkapur](https://discourse.julialang.org/u/maxkapur)\
**Post date:** [October 7, 2022, 6:07pm UTC](https://discourse.julialang.org/t/feedback-on-a-2d-spring-simulation/88350/7 "2022-10-07T18:07:12Z")

</div>

Yes, will edit.

---

<div class="post-metadata">

**Author:** ![maxkapur](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/maxkapur/32/21208_2.png) [@maxkapur](https://discourse.julialang.org/u/maxkapur)\
**Post date:** [October 7, 2022, 6:13pm UTC](https://discourse.julialang.org/t/feedback-on-a-2d-spring-simulation/88350/8 "2022-10-07T18:13:52Z")

</div>

This is required (I think?) by the interface of `ODEProblem`. The remaining two arguments of `f!` are for parameters (e.g. the spring coefficients in my OP) and time (in case there is a forcing function). Neither of these are necessary in my minimal example.

---

<div class="post-metadata">

**Author:** ![lmiq](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/lmiq/32/18314_2.png) [@lmiq](https://discourse.julialang.org/u/lmiq)\
**Post date:** [October 7, 2022, 6:13pm UTC](https://discourse.julialang.org/t/feedback-on-a-2d-spring-simulation/88350/9 "2022-10-07T18:13:56Z")

</div>

Note that `Point` is immutable, so your `f!` will not work as you expect. To use it like that you need to define `mutable struct`. (but that is not related to that specific issue).

Yes, seems that the interface of ODEProblem is fine as it is.

---

<div class="post-metadata">

**Author:** ![lmiq](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/lmiq/32/18314_2.png) [@lmiq](https://discourse.julialang.org/u/lmiq)\
**Post date:** [October 7, 2022, 6:17pm UTC](https://discourse.julialang.org/t/feedback-on-a-2d-spring-simulation/88350/10 "2022-10-07T18:17:01Z")

</div>

The example runs if you define:

```julia
julia> using StaticArrays

julia> mutable struct Point{T} <: FieldVector{2,T}
           x::T
           y::T
       end

```

Which will give you all the algebra of arrays to your point. The error of missing `oneunit(::Any)` is pretty bad, though.

---

<div class="post-metadata">

**Author:** ![maxkapur](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/maxkapur/32/21208_2.png) [@maxkapur](https://discourse.julialang.org/u/maxkapur)\
**Post date:** [October 7, 2022, 6:26pm UTC](https://discourse.julialang.org/t/feedback-on-a-2d-spring-simulation/88350/11 "2022-10-07T18:26:05Z")

</div>

Yikes, so how do I initialize a `FieldVector`? Not in the [docstring](https://juliaarrays.github.io/StaticArrays.jl/stable/pages/api/#StaticArraysCore.FieldVector) AFAICT.

```julia
julia> a = FieldVector{2}(1.0, 1.0)
ERROR: The constructor for FieldVector{2, Float64}(::Float64, ::Float64) is missing!
...

julia> FieldVector{2, Float64}([2.0, 4.0])
ERROR: The constructor for FieldVector{2, Float64}(::Float64, ::Float64) is missing!
...

julia> FieldVector{2, Float64}(2.0, 4.0)
ERROR: The constructor for FieldVector{2, Float64}(::Float64, ::Float64) is missing!
...

julia> FieldVector{2, Float64}((2.0, 4.0))
ERROR: The constructor for FieldVector{2, Float64}(::Float64, ::Float64) is missing!
...

julia> FieldVector(2.0, 4.0)
ERROR: FieldVector has no static size!

```

Edit: Nevermind, FieldVector is an abstract type, so you have to define your own concrete type like `Point`.