Using desolve_odeint with external stimulus

I'm trying to numerically solve a two-variable system of differential equations with a periodic external driver (square wave). Right now, the code I have is:

from scipy.signal import square
def sqstim(t, amp=1):
return (amp/2)*float(square(t))+amp/2

var("V,w")
trange=srange(0,100,0.1)

def eqns(t):
return [1/100*(-w+(V^3-V)+sqstim(t)), V-0.2*w]

sol=desolve_odeint(eqns, times=trange, ics=[0.1,0.1], dvars=[V,w])

This gives "TypeError: 'function' object is not iterable". What should I do instead?

edit retag close merge delete

Sort by » oldest newest most voted

I've had difficulty getting Sage to handle Python functions, like your sqstim, in differential equations. Sage wants to execute sqstim(t) immediately when it is called in the desolve_odeint command.

Here is an option using numpy and scipy.

import numpy as np
from scipy.integrate import odeint

def eqns(y,t):
V=y
w=y
dydt = [1/100.*(-w+(V^3-V)+sqstim(t)), V-0.2*w]
return(dydt)

tspan=np.linspace(0,100,1000)

sol = odeint(eqns, [0.1,0.1], tspan)

tV = zip(tspan,[s for s in sol])
tw = zip(tspan,[s for s in sol])

Then, you can plot easily using Sage commands.

line(tV,color='red')

and

line(tw)
more

Thanks. Is there a way to incorporate something like a square wave into an ODE in Sage without using Python functions?

This is an alternative approach, using GSL numerical library for ODE's which can be used in Sage via ode_solver().

from scipy.signal import square

def sqstim(t, amp=1):
return (amp/2)*float(square(t))+amp/2

# note that we pass the vector y as well
def eqns(t, y):
V, w = y, y
return [1/100*(-w+(V^3-V)+sqstim(t)), V-0.2*w]

# instantiate a solver object
T = ode_solver()

# set the system of equations and time interval
T.function = eqns
trange=srange(0,100,0.1)

# solve for given initial condition. the default algorithm is "runga-kutta-felhberg", see the help for more options
T.ode_solve(y_0=[0.1,0.1], t_span=[trange, trange[-1]], num_points=len(trange))

We can plot the result with these commands:

# plot solutions
V_t = T.interpolate_solution(0)
w_t = T.interpolate_solution(1)
Fig = Graphics()
Fig += plot(V_t, trange, trange[-1], color='red', legend_label='V')
Fig += plot(w_t, trange, trange[-1], legend_label='w')
show(Fig)
more

As a side-note, the approach by @calc314 is better in computational time since it took me around 250ms, against 800ms using ode_solver(). Moreover, I thought the time would improve using fast_float for eqns and sqstim, but it didn't.