Saving time-dependent variable inside ODE

Specific Domains
,
1

Hi, I’d like to extract time-dependent variable from inside an ODE, but am having trouble doing so given that it is not in the state vector. I have my du function (called pendulum!) inside a system function. I use this setup so that I can pass parameters as arguments to the system function rather than hard coding them. See my code below:

using DifferentialEquations
using Plots

function system(l, m, g)
    function pendulum!(du, u, p, t)
        θ, ω = u
        l, m, g = p
        M = 0.1sin(ω*t)
        du[1] = ω
        du[2] = -3g/(2l)*sin(θ) + 3/(m*l^2)*M
    end

    θ_0 = 0.01
    ω_0 = 0.0
    u_0 = [θ_0, ω_0]
    tspan = (0.0, 10.0)
    params = [l, m, g]
    prob = ODEProblem(pendulum!, u_0, tspan, params)
    sol = solve(prob)
    return sol
end

mySol = system(1.0, 1.0, 9.81)
plot(mySol, linewidth = 2, xaxis = "t", label = ["θ" "ω"], layout = (2,1))

This all works great until I’d like to see how M = 0.1sin(ω*t) variable changes overtime (ideally plot it). I haven’t found a way to save that data outside of the pendulum function, so in a way it is trapped in there.

I have thought of appending the value of M to a local variable at each time step and returning it after the integration, however the local variable gets reset at each time step.

Additionally, I have considered making a global variable outside of the function but it cannot be recognized inside the pendulum function. Is there an elegant way to do this? Any ideas would be much appreciated. Thanks in advance!

2

For discrete variables, that are updated only in certain times, you can use callbacks and a vector inside the parameters p to store them. However, for continuous values as you want, I am not sure what is the best approach.

3

Hi, I added the following to your code. See if it helps and if this is what you want. mySol(t) can be treated as a function of time itself. I also removed your pendulum function from within the system function (perhaps it doesn’t matter?).

function M(t)
    ω = mySol(t)[1]
    M = 0.1*sin(ω*t)
    return M
end

plot(t->M(t), xlims=(0,10), label="M(t)")

JuliaDiscourseReply

Gabriel

1 Like
4

The easiest way is to use ModelingToolkit.

eqs = [
        M ~ 0.1sin(ω*t)
        D(θ) ~ ω
        D(ω) ~ -3g/(2l)*sin(θ) + 3/(m*l^2)*M
]

structural_simplify with that set of equations will eliminate the extra equation, but sol[M] will grab it from the observed, or plot(sol,vars=M) would plot it.

See the first MTK tutorial:

https://mtk.sciml.ai/dev/tutorials/ode_modeling/

The other way is a post analysis as @gnicolosi shows

5 Likes
6

Thanks for your reply! This works well in the example I’ve provided, but quickly becomes unideal when M is more than just one line (e.g. interpolating several variables, algebraically manipulating several variables in a long state vector, etc). I wanted to avoid having to re-write the definition of the desired variable M because theoretically it could take hundreds of lines to define.

7

Thank you! This is exactly what I was looking for

8

From my understanding ModelingToolkit cannot handle vector-valued variables. Is this correct? I would like a way to both grab sol[M] but also allow my variables to be vectors-- just wanted to see if this was possible in Julia

9

If you make M a vector variable that will work. It doesn’t internally use it, but it externally allows one to define equations with them.