CODEWITHSUNDEEP
pythonpython-3.xscipyode
Ezbob
Ezbob

Reputation: 449

Using adaptive time step for scipy.integrate.ode when solving ODE systems

I have to just read Using adaptive step sizes with scipy.integrate.ode and the accepted solution to that problem, and have even reproduced the results by copy-and-paste in my Python interpreter.

My problem is that when I try and adapt the solution code to my own code I only get flat lines.

My code is as follows:

from scipy.integrate import ode
from matplotlib.pyplot import plot, show

initials = [1,1,1,1,1]
integration_range = (0, 100)

f = lambda t,y: [1.0*y[0]*y[1], -1.0*y[0]*y[1], 1.0*y[2]*y[3] - 1.0*y[2], -1.0*y[2]*y[3], 1.0*y[2], ]

y_solutions = []
t_solutions = []
def solution_getter(t,y): 
   t_solutions.append(t)
   y_solutions.append(y) 


backend = "dopri5"
ode_solver = ode(f).set_integrator(backend)
ode_solver.set_solout(solution_getter)
ode_solver.set_initial_value(y=initials, t=0)

ode_solver.integrate(integration_range[1])

plot(t_solutions,y_solutions)
show()

And the plot it yields: enter image description here

Upvotes: 2

Views: 1474

Answers (1)

Lutz Lehmann
Lutz Lehmann

Reputation: 25992

In the line 

   y_solutions.append(y)

you think that you are appending the current vector. What actally happens is that you are appending the object reference to y. Since apparently the integrator reuses the vector y during the integration loop, you are always appending the same object reference. Thus at the end, each position of the list is filled by the same reference pointing to the vector of the last state of y.

Long story short: replace with

    y_solutions.append(y.copy())

and everything is fine.

Upvotes: 4

Related Questions