在matlab ode45中使用事件函数进行多维状态向量 [英] Using event function in matlab ode45 for multi-dimensional state vector

问题描述

我有一套以矩阵形式写成的odes $ X'= AX $;我还有一个所需的状态$ X_des $的值。 $ X $是一个五维向量。所有国家达到预期的价值后，我想停止整合（或至少接近他们$ 1e {-3} $）。如何在matlab中使用事件函数呢？ （我所看到的所有帮助都是关于1维状态）

PS：我知道所有的状态都能长时间接近他们想要的值。我只想在所需值中为$ 1e {-3} $时停止整合。

解决方案

首先，我假设您知道您可以使用矩阵指数（ expm 在Matlab中）解决您的系统线性微分方程直接。

有很多方法来完成你想要做的事情。它们都取决于您的系统，它的行为方式，以及您要捕获的特定事件。这是一个2-by-2系统的线性微分方程的一个小例子：

  function multipleeventsdemo 
 A = [ -1 1; 1 -2]; ％示例A矩阵
 tspan = [0 50]; ％初始和最终时间
 x0 = [1; 1]; ％初始条件
 f = @（t，y）A * y; ％ODE函数
 thresh = 0; ％阈值
 tol = 1e-3; ％Tolerance on threshold 
 
 opts = odeset（'Events'，@（t，y）events（t，y，thresh，tol））; ％创建事件函数
 [t，y] = ode45（f，tspan，x0，opts）; ％集成选项
数字; 
 plot（t，y）; 
 
函数[value，isterminal，direction] = events（t，y，thresh，tol）
 value = y-thresh-tol; 
 isterminal = all（y-thresh-tol <= 0）+零（大小（y））; ％变更终止条件
 direction = -1; 
  

当两个状态都在 tol thresh 。这是通过调整事件函数的 isterminal 输出来实现的。请注意，单独的公差和阈值变量并不是必需的 - 您只需要定义交叉值。

如果您的系统在接近稳定状态时发生振荡（如果 A 具有复杂的特征值），那么你需要做更多的工作。但你的问题并不表示这一点。而且，数值整合可能不是解决你这个系统的问题的最简单/最好的方法。以下是您如何使用 expm 与一点符号数学：

  A = [-1 1; 1 -2]; 
 x0 = [1; 1]; 
 tol = 1e-3; 
 syms t_sym 
 y = simplified（expm（A * t_sym）* x0）％Y作为t 
 t0 = NaN（1，length（x0））的函数; 
 i = 1：length（x0）
 sol = double（solve（y（i）== tol，t_sym））％当y（i）等于tol 
时，求解t如果〜isempty（sol）％可以不是解，那么NaN 
 t0（i）= max（sol）; ％或多个解决方案，最大
 end 
 end 
 f = matlabFunction（y）; ％创建向量化函数t 
 t_vec = linspace（0，max（t0），1e2）; ％时间向量
数字; 
 plot（t_vec，f（t_vec））; 
  

这只适用于相当小的 A 然而，因为象征性的数学。使用 expm 的数值方法也是可能的，可能更具可扩展性。

I have a set of odes written in matrix form as $X' = AX$; I also have a desired value of the states $X_des$. $X$ is a five dimensional vector. I want to stop the integration after all the states reach their desired values (or atleast close to them by $1e{-3}$). How do I use event function in matlab to do this? (All the help I have seen are about 1 dimensional states)

PS: I know for sure that all the states approach their desired values after long time. I just want to stop the integration once they are $1e{-3}$ within the desired values.

解决方案

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.

这篇关于在matlab ode45中使用事件函数进行多维状态向量的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！