## C-integrators and specifying systems of eqs.

Help
Paw
2013-09-18
2013-11-04

• Paw
2013-09-18

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)
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.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()
```

• Paw
2013-09-18

Hmm, I did a poor job posting the code.
:)

• Paw
2013-09-18

I forgot to add that I'm running 32 bit linux and test  Dopri_backwards_test.py works fine.

• Rob Clewley
2013-09-18

Hi,

Correct me if I'm overlooking something, but I think the solution using the C integrators is relatively easy. You do not seem to need the meshgrid functionality in your final simulation. If that's the case, then it should be fairly simple to use some custom C code and a C math library that supports complex numbers and arithmetic. You can generate the C source code for your model, again inserting a call to your C vector field function (that calls your C vel function). Alternatively, you can make a dummy 4D RHS and edit it later (probably easiest unless you want to try to make it more automated).

In either case, just follow the example in tests/HH_model_Cintegrator.py. Use the nobuild=True option then edit the C source file to include your new functions and the edited vector field function (to redirect to yours). You just have to enter new values into the "return" array (it's call by reference so there's no explicit return). Then add include any header files for other libraries you used in the call to makeLib.

Let me know if you have more questions.
-Rob

• Paw
2013-09-19

Thanks for the quick answer. Will try it tomorrow.
In the mean time, I guess there is no way I can do something like

```    udot = '-(u-velx(x,y))'
vdot = '-(v-vely(x,y))'
xdot = 'c*u'
ydot = 'c*v'
aux = {
'velx': (['x','y'], 'real(1*(1-1**2/(x+1j*y)**2))' ),
'vely': (['x','y'], '-imag(1*(1-1**2/(x+1j*y)**2))' )
}
[\code]
?
And then a last question:
When I call dfdt(t,up,vp,xp,yp,C), can it be done so the call is dfdt(t,f,C) instead, where f = [up,vp,xp,yp]?
Thanks a lot
Paw
```

• Rob Clewley
2013-09-19

Yes, you can do that provided you have the functions real and imag defined in C. I don't think they are in the default C math library but I'm rusty on what's in there. If necessary, you can add a list of custom libraries to include when you call makeLib, and you could call that right after nobuild=True even if you don't go "outside" to edit anything by hand in between those two steps.

I'm not sure at which level your question about dfdt refers to. You can't specify a vector argument directly in the PyDSTool spec strings for functions or RHS's. You can make your custom python or C function accept any form of vector argument when you call it from the regular, original RHS vector field function.

HTH.

• Paw
2013-09-19

Got it to work with

```udot = '-(u-velx(x,y))'
vdot = '-(v-vely(x,y))'
xdot = 'c*u'
ydot = 'c*v'
aux = {
'velx': (['x','y'], 'creal(1*(1-1**2/( (x+1*I*y)*(x+1*I*y) )))' ),
'vely': (['x','y'], '-cimag(1*(1-1**2/( (x+1*I*y)*(x+1*I*y) )))' )
}
DSargs['nobuild']=True
DSargs.ignorespecial = ['cimag','creal','I']
...
ode.makeLib(include=['complex.h'])
```

Note that **2 cannot be used in vel. It compiles to pow() - as expected - and we need cpow().
It is lightening fast. Especially when you come from an matlab environment.

Thank you! This is nice :)

Btw.
Can you explain what you mean by a you can make a dummy 4D RHS and edit it later (probably easiest unless you want to try to make it more automated).. I do not quite get it.

Paw

• Rob Clewley
2013-09-22

The dummy RHS is for the C code situation, where you set each RHS = 0 or something like that. Then you edit the C code for the vector field and insert any function calls you want. I think you've already found the best solution to your problem here.