Using desolve_odeint with external stimulus

savecancel

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[0]
    w=y[1]
    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[0] for s in sol])
tw = zip(tspan,[s[1] for s in sol])

Then, you can plot easily using Sage commands.

line(tV,color='red')

and

line(tw)

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[0], y[1]
    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[0], 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[0], trange[-1], color='red', legend_label='V')
Fig += plot(w_t, trange[0], trange[-1], legend_label='w')
show(Fig)