Использование функции события в Matlab ode45 для многомерного вектора состояния

Во-первых, я полагаю, вы знаете, что можете использовать экспоненциальную матрицу (expm в Matlab), чтобы напрямую решить вашу систему линейных дифференциальных уравнений.

Есть много способов добиться того, что вы пытаетесь сделать. Все они немного зависят от вашей системы, ее поведения и конкретного события, которое вы хотите зафиксировать. Вот небольшой пример для системы линейных дифференциальных уравнений 2 на 2:

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;

Интеграция останавливается, когда оба состояния находятся в пределах tol из thresh. Это достигается путем настройки isterminal вывода функции событий. Обратите внимание, что отдельные переменные допуска и порога на самом деле не нужны — вам просто нужно определить значение пересечения.

Если ваша система колеблется по мере приближения к устойчивому состоянию (если A имеет комплексные собственные значения), вам нужно будет проделать дополнительную работу. Но ваши вопросы не указывают на это. И опять же, численное интегрирование может быть не самым простым/лучшим способом решения вашей проблемы, связанной с такой системой. Вот как вы могли бы использовать expm в сочетании с небольшим количеством символической математики:

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));

Однако это будет работать только для довольно маленьких A из-за символической математики. Численные подходы с использованием expm также возможны и, вероятно, более масштабируемы.