Ask Your Question

Using desolve_odeint with external stimulus

asked 2016-12-31 22:54:09 +0100

jaia gravatar image

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


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 flag offensive close merge delete

2 Answers

Sort by ยป oldest newest most voted

answered 2017-01-01 05:02:00 +0100

calc314 gravatar image

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):
    dydt = [1/100.*(-w+(V^3-V)+sqstim(t)), V-0.2*w]


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.



edit flag offensive delete link more


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

jaia gravatar imagejaia ( 2017-01-01 22:32:36 +0100 )edit

answered 2017-01-01 14:52:40 +0100

mforets gravatar image

updated 2017-01-01 15:10:18 +0100

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    

# 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')
edit flag offensive delete link 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.

mforets gravatar imagemforets ( 2017-01-01 14:53:46 +0100 )edit

Your Answer

Please start posting anonymously - your entry will be published after you log in or create a new account.

Add Answer

Question Tools

1 follower


Asked: 2016-12-31 22:54:09 +0100

Seen: 433 times

Last updated: Jan 01 '17