Revision history [back]

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

%%time
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)

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

%%time
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)