Bug in Vode integrator with external data?

[Paw]

Dear all,
I want to do an integration with external data, much like here
So for x_ddot + d*x_dot + k*x = a*sin(omega*t), instead of the sine excitation, I want to supply my own array. I have a time array and corresponding forces.

From examples/interp_vode_test.py I made the posted example. According to Rob's reply max_step is set.

The interpolated force is returning sane values, but when used in the ode, there seems to be a constant phase lage for the interpolated force. This is only the case for the vode integrator. For the Dopri integrator there is no phase lag. See the code where the abs. error is plotted.

But even with the Dopri integrator, the two solutions differs. No matter if I decrease max_step and increase max_pts, the error between the two solutions seems constant. The interpolated force and sine returns the same value, ie. both RHS' seems to be equal.
I would very much like to two solutions to be identical. What can I do to ensure that? Is there some settings in the docs I should use?

Another thing: How do I:
Better now is to create a Pointset of data then use pointsettotraj, and pass the variable from that in to the creation of the Generator instead of using the InterpolateTable?

Best regards,
Paw

import numpy as np
import matplotlib.pyplot as plt
import PyDSTool as dst
from copy import copy

fs = 10
timeData = np.arange(501)/fs
tdomain = [timeData[0], timeData[-1]]
inputData = np.sin(timeData)
xData = {'force': inputData}

my_input = dst.InterpolateTable({'tdata': timeData,
                                 'ics': xData,
                                 'name': 'interp1d',
                                 'method': 'linear', # next 3 not necessary
                                 'checklevel': 1,
                                 'abseps': 1e-6,
                              }).compute('interp1d')
print(my_input(5.9, ['force']), np.sin(5.9))
DSargs = dst.args()
DSargs.tdata = tdomain
DSargs.pars = {'m':0.1, 'c':1 ,'k':0.5}
DSargs.inputs = my_input.variables['force']
DSargs.varspecs = {'y': 'v', 'v': '-c*v/m - k*y/m + force/m ', 'inval': 'force'}
DSargs.vars = ['y', 'v']
DSargs.algparams = {'init_step' :0.01, 'max_step': 0.5}
DSargs.checklevel = 2
DSargs.ics = {'y': 0, 'v': 0}
DSargs.name = 'ODEtest'

python = True
if python:
    DS = dst.Generator.Vode_ODEsystem(DSargs)
else:
    DS = dst.Generator.Dopri_ODEsystem(DSargs)
traj = DS.compute('in-table')
pts = traj.sample(dt=1/fs, precise=True)

DSargs2=copy(DSargs)
DSargs2.name = 'sin'
DSargs2.varspecs = {'y': 'v', 'v': '-c*v/m - k*y/m + sin(t)/m ' }
DSargs.pop('inputs', None)
DS2 = dst.Generator.Dopri_ODEsystem(DSargs2)
DS2.makeLib()
traj2 = DS2.compute('sin')
pts2 = traj2.sample(dt=1/fs, precise=True)

# Both RHS seems to be the same at t:5.7
print("DS: y_ddot, y_dot:", DS.Rhs(5.7, {'y':1., 'v':0.5}, DSargs.pars))
print("DS2: y_ddot, y_dot:", DS2.Rhs(5.7, {'y':1., 'v':0.5}, DSargs.pars))

error = (pts['inval']-np.sin(pts['t']))
error2 = (pts['y']-pts2['y'])
plt.figure()
plt.plot(pts['t'], pts['y'], '-ok',label = 'in-table')
plt.plot(pts2['t'], pts2['y'],'--xr', label = 'sin')
plt.legend()

plt.figure()
plt.plot(pts['t'], error, label='error: in-table/sin')
plt.plot(pts['t'], error2,'--k', label='error: y-vals')
plt.legend()
plt.show()

[Rob Clewley]

This is definitely a problem, thanks for posting. I've had some success with addressing it here (https://github.com/robclewley/pydstool/pull/128) but it's not complete and I don't know how to proceed right now. Dopri is definitely the superior way in all cases, though ...