2

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

Community
  • 1
  • 1
Ezbob
  • 329
  • 1
  • 9

1 Answers1

4

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.

Lutz Lehmann
  • 25,219
  • 2
  • 22
  • 51
  • Thank you very much! References still manages to sometimes confuse me... Btw, how did you come to this solution? – Ezbob Oct 03 '16 at 09:53
  • 2
    By debugging. The q'n'd solution is to add print statements inside the solution_getter and before the plot to see what the actual data is. – Lutz Lehmann Oct 03 '16 at 11:08