How can I simplify the creation of a system of equations in a loop:
using ModelingToolkit, OrdinaryDiffEq, LinearAlgebra
G_EARTH = Float64[0.0, 0.0, -9.81] # gravitational acceleration [m/s²]
L0::Float64 = 5.0 # initial segment length [m]
V0::Float64 = 2 # initial velocity of lowest mass [m/s]
M0::Float64 = 0.5 # mass per particle [kg]
C_SPRING::Float64 = 50 # spring constant
segments::Int64 = 5 # number of tether segments [-]
α0 = π/10 # initial tether angle [rad]
duration = 30.0 # duration of the simulation [s]
POS0 = zeros(3, segments+1)
VEL0 = zeros(3, segments+1)
ACC0 = zeros(3, segments+1)
SEGMENTS0 = zeros(3, segments)
UNIT_VECTORS0 = zeros(3, segments)
for i in 1:segments+1
local l0
l0 = -(i-1)*L0
v0 = (i-1)*V0/segments
POS0[:, i] .= [sin(α0) * l0, 0, cos(α0) * l0]
VEL0[:, i] .= [sin(α0) * v0, 0, cos(α0) * v0]
end
for i in 2:segments+1
ACC0[:, i] .= G_EARTH
end
for i in 1:segments
UNIT_VECTORS0[:, i] .= [0, 0, 1.0]
SEGMENTS0[:, i] .= POS0[:, i+1] - POS0[:, i]
end
# defining the model, Z component upwards
@parameters mass=M0 c_spring0=C_SPRING damping=0.5 l_seg=L0
@variables t
@variables pos(t)[1:3, 1:segments+1] = POS0
@variables vel(t)[1:3, 1:segments+1] = VEL0
@variables acc(t)[1:3, 1:segments+1] = ACC0
@variables segment(t)[1:3, 1:segments] = SEGMENTS0
@variables unit_vector(t)[1:3, 1:segments] = UNIT_VECTORS0
@variables norm1(t)[1:segments] = l_seg * ones(segments)
@variables rel_vel(t)[1:3, 1:segments] = zeros(3, segments)
@variables spring_vel(t)[1:segments] = zeros(segments)
@variables c_spring(t)[1:segments] = c_spring0 * ones(segments)
@variables spring_force(t)[1:3, 1:segments] = zeros(3, segments)
@variables total_force(t)[1:3, 1:segments] = zeros(3, segments)
D = Differential(t)
eqs1 = vcat(D.(pos) ~ vel,
D.(vel) ~ acc)
eqs2 = []
for i in segments:-1:1
global eqs2
eqs2 = vcat(eqs2, segment[:, i] ~ pos[:, i+1] - pos[:, i])
eqs2 = vcat(eqs2, norm1[i] ~ norm(segment[:, i]))
eqs2 = vcat(eqs2, unit_vector[:, i] ~ -segment[:, i]/norm1[i])
eqs2 = vcat(eqs2, rel_vel[:, i] ~ vel[:, i+1] - vel[:, i])
eqs2 = vcat(eqs2, spring_vel[i] ~ -unit_vector[:, i] ⋅ rel_vel[:, i])
eqs2 = vcat(eqs2, c_spring[i] ~ c_spring0 * (norm1[i] > l_seg))
eqs2 = vcat(eqs2, spring_force[:, i] ~ (c_spring[i] * (norm1[i] - l_seg) + damping * spring_vel[i]) * unit_vector[:, i])
if i == segments
eqs2 = vcat(eqs2, total_force[:, i] ~ spring_force[:, i])
else
eqs2 = vcat(eqs2, total_force[:, i] ~ spring_force[:, i]- spring_force[:, i+1])
end
eqs2 = vcat(eqs2, acc[:, i+1] .~ G_EARTH + total_force[:, i] / mass)
end
eqs2 = vcat(eqs2, acc[:, 1] .~ zeros(3))
eqs = vcat(eqs1..., eqs2)
@named sys = ODESystem(eqs, t)
simple_sys = structural_simplify(sys)
For example:
Can I use only one vcat call for these two lines?
Or can I use some other function than vcat?
eqs2 = vcat(eqs2, unit_vector[:, i] ~ -segment[:, i]/norm1[i])
eqs2 = vcat(eqs2, rel_vel[:, i] ~ vel[:, i+1] - vel[:, i])
I tried push! and append! with no success…
I found out, some of the equations in the array eqs2 have the type Equation and some the type BasicSymbol. Perhaps that causes the issues?
Why is this not of type Equation:
julia> eqs2[1]
(broadcast(~, (segment(t))[Colon(), 5], broadcast(-, (pos(t))[Colon(), 6], (pos(t))[Colon(), 5])))[1]
julia> typeof(eqs2[1])
SymbolicUtils.BasicSymbolic{Any}
But this is:
julia> eqs2[end]
damping(t) ~ 473.0 / ((1//5)*length(t))
julia> typeof(eqs2[end])
Equation