Runge-Kutta with impulse

I'm trying to study a simple predator-prey type ODE system of two variables, but I'd like to analyze the impulse response. Really, I have two versions of the same dynamical system, so two copies of each variable, where all the variables interact according to the same predator-prey model (so the herbivores in both systems can be eaten by carnivores in both systems), but where one of the systems is always on and where the other is off until some time when I subject it to an impulse. I'm using Runge-Kutta for integration, and I'd like to study the impulse response of that parallel system.

Note that I want the total sum of the copies of the two variables to be constant (when I impulse the one system I 'negative' impulse the other, to keep this consistent). Really I'm trying to study a single system, but separate out the contribution of the driving signal in a single time step so I can track how long it contributes to system dynamics.

I noticed that if I start the one system in steady state, so nothing is changing, and then impulse the other, the first system is kicked out of steady state and starts oscillating (while the other decays to 0, as expected). This oscillation suggested something was wrong, since it should return to steady state when the other system dies out. I've since learned that Runge-Kutta was not meant to deal with impulse.

I've tried two things: the first is to only have the impulse in the first of the four elements of the runge-kutta scheme at impulse time. This seems to work, and gives me back my steady state after the other system decays, but there is some very weird transience in between (when the sum of variables from the two systems is no longer constant). I have also tried explicitly using only Euler's method at the impulse time, and this gives a similar result.

Is there something else I can try or a better solution to this problem?

This is actually a fairly complicated system. Let me see if I am interpreting what you have properly.

You begin with some system with initial conditions, and a forcing function, but there is no "input" term, i.e. $$y' = f(t,y).$$

At some point $t = T$, you want to add an impulse of some magnitude. So what you really have now is

$$y' = f(t,y)+u(t)$$

where, with $K$ being the amplitude of your impulse,

$$u(t) = \left\{ \begin{array}{cc} 0, & t \neq T, \\ K, & t = T, \end{array}\right.$$

Presumably, you can compute with some accuracy a solution of your system prior to $t=T$.

I would then examine solving the finite impulse response of the modified system

$$y'(t-T) = f(t-T,y)+K\delta(0)$$

where $\delta$ is the impulse "function". You can then use your solution at time $t=T$ as the initial conditions for this system. There might be some better ways to handle the impulse response at non-zero $t$, but effectively this is what it's doing.

You can analyze this system with a Laplace tranform, but what you're really looking for is the finite impulse response of your ODE system. I might also recommend looking at SignalProcessing.SE for some advice.

In general, however, numerical schemes do very poorly with discontinuous input terms. The problem is that most numerical methods "look ahead" either to predict where the system will be (in a PECE technique) and then adjust, or use some moving estimate of the prior slope of the solution to project the solution forward. In either case, the approach has a really hard time with the discontinuity.

A PECE method might give you better results, but it really all depends on many factors.

• Just to be slightly more explicit about what Arkamis is saying (if I understand correctly): use Runge-Kutta or whatever to solve your system up to slightly before the impulse, giving you $y(T-)$ where $T$ is the time of the impulse and $T-$ is "just before"; then add in the impulse, giving you $y(T+) = y(T-) + K$; and then use Runge-Kutta again with that as initial condition. Mar 8 '13 at 9:07