CODEWITHSUNDEEP
pythonscipy
SourBitter
SourBitter

Reputation: 157

Scipy ode solver

Using scipy version 0.18.0 and really scratching my head at the moment. I reduced the code to the following minimal example:

try to solve the easiest differential equation possible

def phase(t, y):
    c1 = y
    dydt = - c1 
    return dydt

c1 = 1.0
y0 = c1
t = np.linspace(0,1,100)
ode_obj = sp.integrate.ode(phase)
ode_obj.set_initial_value(y0)
sol = ode_obj.integrate(t)

The object tells me it was successful, but the return value is a numpy array of length 1.

The documentation is extremely sparse and I'm not sure if I'm using it incorrectly.

Thanks for the help guys.

Greetings.

Upvotes: 1

Views: 823

Answers (1)

SourBitter
SourBitter

Reputation: 157

Found this issue on github https://github.com/scipy/scipy/issues/1976 This basically explains how it works and when you think how these initial value solvers work ( stepwise into direction of the final point ) it gets clearer why this is done this way. My code from above would look like this:

import scipy as sp
import pylab as pb

def phase(t, y):
    c1 = y
    dydt = - c1 
    return dydt

c1 = 1.0
t0 = 0.0 # not necessary here
y0 = c1
t1 = 5.0

t = []
sol = []
ode_obj = sp.integrate.ode(phase)
ode_obj.set_initial_value(y0, t=t0)
ode_obj.set_integrator('vode',method='bdf',rtol=1.e-12)
sol = ode_obj.integrate(t)
while ode_obj.successful() and ode_obj.t < t1:
    ode_obj.integrate(t1,step=100)
    t.append(ode_obj.t)
    sol.append(ode_obj.y)

pb.plot(t,sol)
pb.show()

This will produce the following output: enter image description here

Upvotes: 1

Related Questions