First, I presume that you're aware that you can use the matrix exponential (expm
in Matlab) to solve your system of linear differential equations directly.
There are many ways to accomplish what you're trying to do. They all depend a bit on your system, how it behaves, and the particular event you want to capture. Here's a small example for a 2-by-2 system of linear differential equations:
function multipleeventsdemo
A = [-1 1;1 -2]; % Example A matrix
tspan = [0 50]; % Initial and final time
x0 = [1;1]; % Initial conditions
f = @(t,y)A*y; % ODE function
thresh = 0; % Threshold value
tol = 1e-3; % Tolerance on threshold
opts = odeset('Events',@(t,y)events(t,y,thresh,tol)); % Create events function
[t,y] = ode45(f,tspan,x0,opts); % Integrate with options
figure;
plot(t,y);
function [value,isterminal,direction] = events(t,y,thresh,tol)
value = y-thresh-tol;
isterminal = all(y-thresh-tol<=0)+zeros(size(y)); % Change termination condition
direction = -1;
Integration is stopped when both states are within tol
of thresh
. This is accomplished by adjusting the isterminal
output of the events function. Note that separate tolerance and threshold variables isn't really necessary – you simply need to define the crossing value.
If your system oscillates as it approaches it's steady state (if A
has complex eigenvalues), then you'll need to do more work. But you questions doesn't indicate this. And again, numerical integration may not be the easiest/best way to solve your problem which such a system. Here is how you could use expm
in conjunction with a bit of symbolic math:
A = [-1 1;1 -2];
x0 = [1;1];
tol = 1e-3;
syms t_sym
y = simplify(expm(A*t_sym)*x0) % Y as a function of t
t0 = NaN(1,length(x0));
for i = 1:length(x0)
sol = double(solve(y(i)==tol,t_sym)) % Solve for t when y(i) equal to tol
if ~isempty(sol) % Could be no solution, then NaN
t0(i) = max(sol); % Or more than one solution, take largest
end
end
f = matlabFunction(y); % Create vectorized function of t
t_vec = linspace(0,max(t0),1e2); % Time vector
figure;
plot(t_vec,f(t_vec));
This will only work for fairly small A
, however, because of the symbolic math. Numerical approaches using expm
are also possible and likely more scalable.