First time here? Check out the FAQ!

Ask Your Question
1

Using desolve_odeint with external stimulus

asked 8 years ago

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

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?

Preview: (hide)

2 Answers

Sort by » oldest newest most voted
2

answered 8 years ago

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):
    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)
Preview: (hide)
link

Comments

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

jaia gravatar imagejaia ( 8 years ago )
1

answered 8 years ago

mforets gravatar image

updated 8 years ago

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)
Preview: (hide)
link

Comments

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 ( 8 years ago )

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

Stats

Asked: 8 years ago

Seen: 582 times

Last updated: Jan 01 '17