How to apply a symbolic equation on a vector?

Specific Domains Modelling & Simulations
, ,
1

I have the following code:

using ModelingToolkit, OrdinaryDiffEq, DataInterpolations, Plots

Cp=[0.024384377390186646,0.03811829845979822,0.05731175742000688,0.0803120051369546,0.10526323380995135,0.1312349486430945,
    0.1588590278697544,0.1882995422205645,0.21937359147848456,0.25269290703920616,0.2882162631147021,0.3240355349501659,
    0.3607290350346774,0.39295261085207084,0.4141821251393675,0.4274331262027119,0.43636422894857274,0.4414767530303278,
    0.443412806515125,0.44365385445693295,0.44317448744425714,0.44218641629727234,0.4407405317795181,0.43888039054963845,
    0.4362875461540461,0.4325702155989231,0.4278606704093836,0.4224263790651815,0.4164272616252002,0.40987027258773945,
    0.4027832291503407,0.39518785578723936,0.38719847687832065]
TSR = 2.0:0.25:10.0

cp = CubicSpline(Cp, TSR)

@variables t ω(t) Pr(t) λ(t)
@parameters J R ρ A U Pg
D = Differential(t)

eqs = [J * D(ω) * ω ~ Pr - Pg,
       λ            ~ ω * R/U,
       Pr           ~ 0.5 * ρ * A * U^3 * cp(λ)
      ]

@named sys = ODESystem(eqs, t)
sys = structural_simplify(sys)

u0 = [ω    => 0.9, # this is a fixed initial value
      D(ω) => 0.0] # the value for the derivative is an inital estimate

p = [ J => 4.0470e+07,
      R => 63,
      ρ => 1.225,
      A => 1.2469e+04,
      U => 8.1,
      Pg => 1771e3]

tspan = (0.0, 1000.0)
prob = ODEProblem(sys, u0, tspan, p)
sol = solve(prob, Rosenbrock23())

#### plot the result ####
X = sol.t
Y1 = [vec[1] for vec in sol.u]
Y2 = [vec[2] for vec in sol.u]
println("ω at 1000s: $(Y1[end])")
display(plot(X, Y1, label="ω"))

### How to plot λ ? I have already the equation:
eqs[2]

Because I simplify the system of equations the variable \lambda does not appear in the solution. But I want to plot \lambda(t).

How can I do this without writing the second equation again?

2

plot(sol, idxs=λ)

You can also index the solution with a symbolic index, I think the syntax is something like this

sol(t, idxs=λ)
2 Likes
3

This works:

But it has not the desired effect because I cannot use the parameter reuse=false to get two diagrams on the screen. So I need to extract the vector of values of \lambda(t) from the solution.
This doesn’t work:

sol(t, idxs=λ)

UPDATE:
This works:

TSR = sol(sol.t, idxs=λ)

Thank you!

1 Like
4

@baggepinnen Do you know when this was added, and what solution types it works with? I have an outstanding PR implementing this for ODESolution