Basic problem creating autonomous external input

[AM]

Hi, I am a newcomer to Python and PyDSTool having previously used matlab, xpp etc.

I am working in the Ipython notebook using the latest version of PyDSTool from git.

So far it seems like a really great toolbox :) Unfortunately, I am having difficulty getting the autonomous external inputs working correctly. I have managed to run and modify both test_hybrid_extinputs.py and test_hybrid_extinputs_Cintegrator.py from the included examples.

The system I am working on does not really require any sort of hybrid features at the moment, I am really just interested in having an external input varying as a function of time in a small system of nonlinear ODEs. I tried to create a minimal example of a one-dimensional system with a single external input (below), as a proof of concept. The problem is that I can observe the external input when it is recorded as an auxiliary variable, but it does not seem to be recognised as part of the RHS. The example I have implemented looks like this:

:::python
import numpy as np
import matplotlib.pyplot as plt
import PyDSTool as dst

DSargs = dst.args()
DSargs.name = 'control_minimal'
DSargs.ics = {'x': 0}
DSargs.pars = {'k': 10, 'alpha': 2} 

# Case #1 - this does not seem to work:
DSargs.varspecs = {'x': 'k*(in - x) - alpha*x', 'inval': 'in'}

# Case #2 - an example of what I would like to happen:
# DSargs.varspecs = {'x': 'k*(if(t > 10, 1.0, 0.0) - x) - alpha*x', 'inval': 'in'}

DSargs.vars = ['x']

tdomain = [0, 20]
timeData = np.linspace(tdomain[0], tdomain[1], 100)
inputData = 0.0 + 1.0*(timeData > 10)
xData = {'in': inputData}

# After trying to generate 'Variable'/'Trajectory' etc as recommended
# in test_hybrid_extinputs.py, I have resorted to directly copying
# the method used in that file
my_input = dst.InterpolateTable({'tdata': timeData,
                              'ics': xData,
                              'name': 'interp1d',
                              'method': 'linear',
                              'checklevel': 1,
                              'abseps': 1e-6
                              }).compute('interp')

DSargs.tdata = tdomain
DSargs.inputs = {'in': my_input.variables['in']}

DS = dst.Generator.Vode_ODEsystem(DSargs)

traj = DS.compute('test_traj')
pts = traj.sample()

# This doesn't work in case #1 above
plt.figure()
plt.plot(pts['t'], pts['x'], label = 'x')

# However this seems to indicate the correct input is being applied
plt.figure()
plt.plot(pts['t'], pts['inval'])

I can't see where I am going wrong so would appreciate any help or tips with implementing this. Thanks in advance!

[Rob Clewley]

Hi, thanks for posting. Very interesting situation. For reasons that are still unclear to me, the imported Vode integrator will not behave properly in the situation of an instantaneous signal change when not given an explicit max_step. As soon as I provide a max_step to alg_params, your code works. It also works without an explicit max_step if the input signal varies smoothly. Now, I'm not super surprised that there is a problem given how ill-behaved things can get with discrete changes, but I don't specifically see where the code messes up. In any case, it's good to have a record of this problem and I'll try to figure it out. In the meantime, add an explicit max time step and/or provide an explicit event (non-terminal) to the system that corresponds to any known sudden changes in the input signal. Another alternative is to slightly smooth out your input signal (e.g. using a very tight tanh()) so that there are more "continuous" changes at the scale of the initial time step.

[AM]

Hi, thanks for the rapid response!

Yes, I guess that it is unsurprising given that the solver (possibly) has no idea there is a discontinuity in this case.

Setting the max_step parameter works fine with both the Vode and Radau integrators for a discontinuous input - thank you.

provide an explicit event (non-terminal) to the system that corresponds to any known sudden changes in the input signal

So, am I correct in thinking you can use the events mechanism to 'hint' where some discontinuity, or other 'meaningful' point will occur?