I am using DifferentialEquations.jl to solve an ODE. After solving it, I plot the different variables.
The problem is:
I want to add a parameter (p) to the same plot. Is it possible? I only see different ways of plotting the variables, but no way to add parameters to the plot.
The reason I want to see the parameter, is because it is changed a few times using iterative callbacks (DiffEqCallbacks).
Could you explain what you mean by “see the parameter”? How would it appear on the plot - is it a point on the plot, a line, some text somewhere?
Do you want to display the parameter as a text ( annotation) on your plot? Are you plotting a solution coordinates, t\to x_i(t), for each i, or the trajectory t\to (x_1(t),\ldots ,x_n(t)), n=2 or 3?
A minimal example and a more precise formulation of your intention to display the parameter for some interval of time, t, would be useful.
My plot is two variable plotted over time (which is the solution of the differential equation). While the parameter is a constant value given as an input to the differential equation.
The parameter would be a line on the plot. It would have been a straight horizontal line, but since the callback function changes the value at specified times, it should look like steps when plotted over time.
Here is my code, if it might help:
using Pkg, ModelingToolkit, BlockSystems, OrdinaryDiffEq, Plots, Interpolations, Symbolics function interp(t) time = [0.05, 1.0, 2.0, 3.0, 6.0, 7.0, 8.0] P_elec = [1.0, 1.0, 0.6, 1.0, 1.0, 1.6, 1.6] interpoleer = linear_interpolation(time, P_elec, extrapolation_bc=Line()) return interpoleer(t) end @parameters t M D P_m(t) P_e(t) @variables ω(t) dt = Differential(t) swing = IOBlock([dt(ω) ~ 1/M * (P_m - D*ω - P_e)], # equations [P_m], # inputs [ω]) # outputs swing = set_p(swing, :D=>1, :M=>1) @variables P_m(t) pfix = IOBlock([P_m ~ 1], , [P_m]) swing_with_P_m_P_e = IOSystem([pfix.P_m => swing.P_m], [pfix, swing]) |> connect_system gen = generate_io_function(swing_with_P_m_P_e) u0 = zeros(length(gen.states)) odef = ODEFunction((du, u, p, t) -> gen.f_ip(du, u, nothing, p, t); mass_matrix=gen.massm, syms=Symbol.(gen.states)) tspan = (0.0, 10.0) p = [0.99] #P_e as parameter prob = ODEProblem(odef, u0, tspan, p) time_choice(integrator) = 0.05 + integrator.t affect!(integrator) = integrator.p .= [interp(integrator.t)] cb = IterativeCallback(time_choice, affect!) sol = solve(prob, Rodas4(autodiff=false), callback = cb) size = (500, 300) linewidth=2 palette=:tab10 grid=false xlabel = "Time" plot(sol; linewidth, legend=:right, palette, grid, xlabel, ylabel="Frequency / Power", size, title="fixed P_m")
A horizontal line is:
hline!([parameter_value], label = "parameter")
but if you want
P_elec at different points in time it’s:
plot!(time, P_elec, label = "parameter")
In general the differential equations ecosystem uses Plots.jl recipes so all the plots once created can just be used like any other Plots.jl plot, i.e. they can be modified by the relevant commands with a bang
Though it didn’t really solve the problem. I am trying to plot “p” over “t”, not “P_elec” over “time”.
affect!(integrator) = integrator.p .= [interp(integrator.t)], calls values of P_elec at different times, but plotting P_elec over Time will not give me the same as “p” over “t”.
You’d have to save it yourself.
p isn’t logged.
I’ll go and try that. Do you have any references/links where I can learn more about logging parameters?
Just add a push in the callback
Thanks. It works!
For others that might want to see the push!():
function simulate(sys) gen = generate_io_function(sys) u0 = zeros(length(gen.states)) odef = ODEFunction((du, u, p, t) -> gen.f_ip(du, u, nothing, p, t); mass_matrix=gen.massm, syms=Symbol.(gen.states)) #@info "Gen state" gen.states tspan = (0, 5) p = [0.8] #P_e as parameter >>>Is this only P_e, or is this a combination of the D, M, and P_e parameters? #@info p prob = ODEProblem(odef, u0, tspan, p) # Logging callback logged_p =  logged_t =  time_choice(integrator) = 0.001 + integrator.t affect!(integrator) = begin integrator.p .= [interp(integrator.t)] push!(logged_p, integrator.p) push!(logged_t, integrator.t) end cb = IterativeCallback(time_choice, affect!) sol = solve(prob, Rodas4(autodiff=false), callback=cb) return sol, logged_p, logged_t end