How do I add an input in the middle of an ode45 integration?

I have a group of odes that I'm integrating with ode45. Let's say from time 0 to 100. I want to add an instantaneous (or nearly) impulse at time 50. How would I do that?

I tried adding an if loop like below in the derivative routine, but I got the following error.

if t >= 50
    original function + impulse
elseif t >= 50.1
    original function
else
    original function
end

Subscript indices must either be real positive integers or logicals.

[Walter Roberson]

Please show the actual code you attempting to add.

original function

would be the code for calling a routine named "original" with a single argument which was the string named 'function', the same thing as

original('function')

[Jordan Falgoust]

Here's the script.

    global Timp
    Timp = 50;
    Tend = 100;
    global mu
    mu = 3.986004418e14;
    global m
    m = 8211;
    r0 = [0,7e5,0];
    v0 = [1.1e4,0,0];
    global Fimp
    Fimp = [100,0,200];

    [tv,Yv] = ode45('state_int',[0 Tend],[r0,v0]);

Here's the function set I'm trying to integrate.

function Fv = state_int(t,Y)
%reshape inputs as column vectors
r0 = [Y(1);Y(2);Y(3)];
v0 = [Y(4);Y(5);Y(6)];

%calculate Force from position
global mu
global m
Fmag = mu*m/(norm(r0))^2;
Fx = -Y(1)/norm(r0)*Fmag;
Fy = -Y(2)/norm(r0)*Fmag;
Fz = -Y(3)/norm(r0)*Fmag;
F = [Fx; Fy; Fz];

r = v0;
v = 1/m*F;

%outputs
Fv([1,2,3],1) = r;
Fv([4,5,6],1) = v;

But I'd like to add an impulsive vector at Timp. So I replaced v with the following if statement.

global Timp
if t <= Timp
    v(t) = 1/m*F;
elseif t <= Timp+0.001
    v(t) = 1/m*(F+Fimp);
else
    v(t) = 1/m*F;
end

And I get the error mentioned in the original post.

Answer 1

[Joe]

You can use logical arithmetic to set conditions for your force. Taking the code snippet you posted above, you would assign "v" as follows:

v = 1/m*F + (t > Timp)*Fimp(:); % (t > Timp) evalutes to 0 if t is under Timp and 1 otherwise

Note that I made "Fimp" a column vector because you will get dimensional incompatibility otherwise.

[Jordan Falgoust]

If I just want Fimp to be applied at a single instant, how would I do that? If I used the following instead of what you suggested, it seems the integrator never recognizes the t value desired.

(t == Timp)

[Joe]

That is tougher, because you can't uniquely define "instant" with respect to ode45. It has an adaptive time-step routine that will most certainly never be exactly equal to "Timp." Worse, even if you specify an interval for the force to be active, i.e.

 interval = 1E-3;
 condition = (t > Timp) && (t < (Timp + interval));

If may still never trip, because the stepper might overstep that interval during the integration.

In this case, your best approach would be to look to Steven's recommendation and stop the integration just prior to your force application. Then, start a new run with the force applied over some interval. Make sure to conserve total impulse (Force*time) and that the interval is fast with respect to your system dynamics. When the time exceeds that interval, stop that run, and start up a third with no applied force.

The conditional approach which I initially suggested works fine when your force is applied over a long time interval, but for impulsive forces, your results can vary widely depending upon the timesteps which the integrator is taking, which you can't control. You can, however, guarantee a force application over a set period of time using the integration time command on ode45, so the separate integration runs is the way to go in that case.

Answer 2

[Steven Lord]

I recommend solving your system of ODEs from t = 0 to t = 50, apply the impulse to the final results from that initial call to the ODE solver to generate new initial conditions, then use those new initial conditions to solve the system of ODEs from t = 50 to t = 100.

The BALLODE example does a variant of this (using the speed of the ball when it hits the ground to set its speed upward after it bounces) but since it doesn't know when the ball is going to hit the ground it needs to use an events function. You DO know when your "event" will occur, so you can just set up the ODE solver call to stop its first run at that time for adjustment.