Ask Your Question
1

Using desolve_odeint with external stimulus

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

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?

edit retag flag offensive close merge delete

2 Answers

Sort by ยป oldest newest most voted
2

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

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)
edit flag offensive delete link more

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 ( 2017-01-01 22:32:36 +0200 )edit
1

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

mforets gravatar image

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

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)
edit flag offensive delete link more

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 ( 2017-01-01 14:53:46 +0200 )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

Stats

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

Seen: 522 times

Last updated: Jan 01 '17