Hi,

I'm having trouble specifying the systems of eqs. in a separate function and at the same time use Dopri/Radau integrators.

Using the python/VODE integrator works, but as I would like to get event detection, c-integrators seems to be the way to go.

I guess the problem at hand is that PyDSTool need the system of eqs. as strings for the generation of c-files.

But I really like the way of using a function for providing the eqs. Especially in my case, where I need to call a function when setting up the eqs.

Are there any way around this?

Regards, Paw

Ps. Thanks to amanderson for help with the 'custom code calls'

Sample code:

import numpy as np import matplotlib.pyplot as plt import time import PyDSTool as pdt domain = np.array([-10, 5, -3, 3])# % domain boundaries gspace = 20 * 0.05 # % grid spacing tspan = np.array([0, 5e2])# % time span vector for integration def vel(x,y): r = 1; U = 1; z = x+1j*y np.seterr(divide='ignore') vel = U*(1-r**2/z**2); u = vel.real v = -vel.imag return u,v ## Compute potential flow field x = np.arange(domain[0],domain[1]+gspace,gspace) y = np.arange(domain[2],domain[3]+gspace,gspace) [X,Y] = np.meshgrid(x,y) [u,v] = vel(X,Y) def dfdt(t,up,vp,xp,yp,C): # setup df = [dup dvp dxp dyp] df = np.zeros(4) # flow velocity [u,v] = vel(xp,yp) df[0] = -(up-u) df[1] = -(vp-v) df[2:4] = C*np.array([up,vp]) return df udot = 'df[0]';vdot = 'df[1]';xdot = 'df[2]';ydot = 'df[3]' ics = {'u':1,'v':0,'x':-32,'y':0.1} pars = {'c': 0.3} DSargs = pdt.args() DSargs.name = 'test' DSargs.ics = ics DSargs.pars = pars DSargs.tdata = tspan DSargs.varspecs = {'u':udot,'v':vdot,'x':xdot,'y':ydot} DSargs.vfcodeinsert_start = 'df = ds.df(t,u,v,x,y,c)' DSargs.ignorespecial = ['df'] DSargs.algparams = {'init_step': 0.1,'atol':1e-5, 'rtol':1e-5 } ode = pdt.Generator.Vode_ODEsystem(DSargs) # ode = pdt.Generator.Dopri_ODEsystem(DSargs) # attach function to ODE definition ode.df = dfdt t = time.time() result = ode.compute('test') print "time for integration: %f" %(time.time() - t) pts = result.sample() # PLOTS fig = plt.figure(1); plt.clf(); plt.hold(True) ax = fig.add_subplot(1, 1, 1) p = ax.plot(pts['x'],pts['y'],'g:',linewidth=2) ax.streamplot(X,Y,u,v, density=[1,0.4],color='cornflowerblue') # add a proxy artist, because streamplot doesn't support legend p2 = plt.Rectangle((0, 0), 1, 1, fc='cornflowerblue') circ = plt.Circle((0, 0), radius=1, color='r') ax.add_patch(circ) ax.legend([p,p2],["particle path, c=%0.1f"%0.3,"Potential flow field"]) plt.legend(); plt.hold(False); plt.grid(True) plt.xlabel('x'); plt.ylabel('y') plt.axis('equal'); plt.axis((domain)) # run matplotlib in interactive mode plt.ion(); plt.show()