Solving and ODE system starting from different initial points


I have a simple system of which i would like to build the phase portrait.
I thought of creating a mesh of interesting points and solve the problem for each of them.
This is what i naively tried until now but it does not work.

xs = -0.5:0.1:1.5
ys = -0.5:0.1:3
u0s = [[x,y] for x in xs, y in ys]
tspan = (0.0, 200.0)
λ = 0
prob = ODEProblem(f, u0s, tspan, λ)
sol = solve(prob, progress=true)

I could not find anything specific in the documentation on how to approach this problem, i know about remake but i am not sure if it is useful or a good idea to use it with so many new initialisations.

Could a for loop work where every time i store the solution for every initial point?
something on the likes of

for u0 in u0s
    sol = solve(remake(prob; u0=u0))

and then store sol somewhere for later access and/or plotting.
(in this case, how do i store the solutions correctly?)

Any smarter ideas, this feels really dumb.


EDIT: added the fix in the second answer

Two possibilities:

  • Shouldn’t for u0 in eachindex(u0s) instead be for u0 in u0s?
  • The remake documentation suggests you should use keyword assignment; this may be unnecessary, though?? Thus sol = solve(remake(prob; u0 = u0))

You are right about both points, i fixed the code and now at least runs.
But now if i try to do


with the new computed solution it does not get updated.

Anyhow, i was mainly wondering if the ODEProblem interface accepted some kind of “multiple initial point” kind of thing, the current approach feels rather hacky.

The following works – I’ve used a boring pendulum model, since you don’t give your model:

function pendulum(du,u,p,t)
    g = 9.81
    L = 2
    alpha = u[1]
    omega = u[2]
    du[1] = omega
    du[2] = -g/L*sin(alpha)
as = -2pi:pi/6:2pi
ws = -0.5:0.1:0.5
us = [[a,w] for a in as, w in ws]
tspan = (0,5.)
prob = ODEProblem(pendulum,us[1],tspan)
sol = solve(prob)
for u in us
    sol = solve(remake(prob,u0 = u))

The result is as follows (not very interesting…):

Thanks a bunch!
I had tried before without the final plot!() and it did not work since i forgot the display inside the for loop! is for that.

1 Like

I had trouble finding the right terminology in the documentations, thanks for pointing me to it!