matlab/octave random event ode45

Question

I have some troubles in understanding how to implement events in octave/matlab, in the resolution of differential equations.

Consider for example this simple code to solve the differential equation y' = -y:

function dy = odefun(t,y)
    dy = -y;
endfunction

options = odeset('RelTol',1e-4,'AbsTol',[1e-4]);
[T,Y] = ode45(@odefun,[0 12],[1],options);

Now I would like to introduce a random event. For example at some fixed time step I would like to randomly change the value of y and then continue the evolution according to the differential equation. How can I do this?

Answer 1

The way to do this is to integrate the system piecewise and append the resultant outputs from each run together:

% t = 0 to t = 6
tspan = [0 6];
y0 = 1;
options = odeset('RelTol',1e-4,'AbsTol',1e-4);
[T,Y] = ode45(@odefun,tspan,y0,options);

% t = 6 to t = 12
tspan = [6 12];
y0 = Y(end,:)+randn; % Pertubation
[t,y] = ode45(@odefun,tspan,y0,options);
T = [T;t(2:end)];  % Remove first value as it will be same as last of previous run
Y = [Y;y(2:end,:)];

The Matlab editor may complain about the array T and Y not being preallocated and/or growing, but it's fine in this case as they're growing in large chunks only a few times.

What you don't want to try to do (but many attempt) is add your perturbation inside the integration function. First, if you want it to occur exactly at a particular time or when particular conditions are met, this can't be done from within the ODE function. Second, inserting large discontinuities in your ODE can lead to less precision and longer computation times (especially with ode45 - ode15s might be a better option or at least make sure that your absolute and relative tolerances are suitable). You would have effectively produced a very stiff system.

If you do want your perturbation occur when particular conditions are met (rather than at a particular time) then you'll need to use an event function (see this question, for example). You'll need to use the events function to terminate integration when the conditions are met and then apply your perturbation as I've shown above.