ODE - It's a question of options

Recently, I have modeled - very summarily - a dampened spring driven by a given signal on one end and free on the other with Octave. The first results baffled me: when the signal was not starting at the beginning of my time range and when it returned to its original value under a certain time, the output would be completely flat, as shown in the following figure.

This is clearly wrong: there is no way that kind of signal won't result in a residual oscillations that will dampen over time. Interestingly enough, when the driving signal repeats, the output is different. But do you notice something weird?

Yup! The oscillations start only at the third descending edge: there is no way this can be the correct solution.

So where did it go wrong?

Let's change the start time for the driving oscillation.

When the first falling edge is at t=1, t=3 or t=5, the simulation shows no oscillation. When the first falling edge occurs at t=2 or t=4, the model seems to behave correctly, and when the first falling edge is at t=6, the first two falling edges are missed. This looks like the time step is variable, and if it falls on an "unchanging" output, it doesn't bother computing the intermediate values.

Octave's lsode_options() allows to get or set the options used by lsode(). These options include the tolerances, the integration method used, but also the minimum and maximum step size. By default the minimum is set to 0 (no minimum) and the maximum is set to -1 (no maximum): the step size can take any value and skip has many points as needed.

Here is the output from lsode_options() right after Octave was started.

Options for LSODE include:

  keyword                                             value
  -------                                             -----
  absolute tolerance                                  1.4901e-08
  relative tolerance                                  1.4901e-08
  integration method                                  stiff
  initial step size                                   -1
  maximum order                                       -1
  maximum step size                                   -1
  minimum step size                                   0
  step limit                                          100000

What if we change the maximum step size to be half of the pulse duration? Let's try that with the command lsode_options("maximum step size",0.5). The options are now

Options for LSODE include:

  keyword                                             value
  -------                                             -----
  absolute tolerance                                  1.4901e-08
  relative tolerance                                  1.4901e-08
  integration method                                  stiff
  initial step size                                   -1
  maximum order                                       -1
  maximum step size                                   0.5
  minimum step size                                   0

  step limit                                          100000

And my test set looks like this.

Which is way better: this is consistent with the spring oscillations starting at the first falling edge. Redoing the first example with the new options, I have:

And that  is consistent with the real life experience.