ode求解器的事件(Event)属性

检测事件

matlab微分方程如何设定变量的范围 
如
dy1=y2
dy2=y1+1
其中y1的范围为0

【解】

M文件：
function [value,isterminal,direction] = events1(t,y)
value = y(1)-4;
isterminal= 1;
direction = 0;

命令窗口：
dy = @(t,y) [y(2);y(1) + 1];
options=odeset('events',@events1);
[t,y] = ode45(dy,[0 12],[0 1],options);
plot(t,y)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%function [value,isterminal,direction]=events(t,x)

% 事件检查函数，此时需要做的是过零点检测
% ode45函数自动检查当value=0是否成立
% 如果我们要求检测Y=0的点，设置value=Y
% 这里我们要检测Y=4，那么就设置value=Y-4

% isterminal检测到指定条件时，是否终止ode45函数的运行
% 1表示终止，0表示继续
% 在我们这个问题上，我们只要检测到零点时就停止程序

% direction：value过零点检测的方向
% -1表示由正到负，+1表示由负到正
% 对于我们这个问题，当然是由正到负

%.........................

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

一个例子和说明：

这是一个用到简单微分方程的物理情景
      一个质量m=100kg的物体从高处竖直落下，加速度会受到空气阻力的影响，这里简单的认为重力加速度g=9.8不变，空气阻力f=k*v^2  ，简单起见，k=1。只考虑竖直方向速度v，且速度，加速度，竖直位移都以向下为正。初速度，初位移都为0；那么有以下微分方程：
dy/dt=v                   
dv/dt=9.8-1*v^2/m
m=100,v0=y0=0

然后我用MATLAB的ode45函数求这个微分方程的数值解
先编写函数

function dx=fun(t,x)
% x(1)表示下落的距离y(向下为正)，x(2)表示下落速度v(向下为正)
k=1;     % k=1为表示空气阻力的一个常量，这里简化空气阻力f=k*v^2
m=100;% m为质量=100kg
dx=zeros(2,1);
dx(1)=x(2);             % 下落距离对时间的导数=速度
a=9.8-(k*x(2)^2)/m;% a加速度(向下为正)=重力加速度 - 空气阻力产生的加速度
dx(2)=a;                 % 速度对时间的导数=加速度
end

现在我想要得到t=15s时的位移和速度
那么输入
[T,X]=ode45('fun',[0,15],[0 0]); 

返回的X中的最后一列就是我想要的值；

但假如我想知道当竖直向下的位移刚好=100米时的时间和速度，那该怎么办？现在我的做法是先将解一个充分大的时间，然后在里面找位移在100两侧的时间和速度，再通过插值得到位移刚好=100时的时间和速度。但这样很麻烦，也不见得准确，MATLAB有什么自带的语句能实现这个功能吗？或是有什么更好的方法？

在不知道结果时间的时候是需要先设定一个比较大的时间范围计算的
但是并不需要将整个范围的结果都算出来再插值
这个时候可以设定触发事件函数在一定条件下停止计算
用odeset可以为ode45求解器设定触发事件的函数
详细的用法要仔细查看matlab的帮助文件，这里我以你的题目，举个例子

微分方程还是用你的函数fun
在用ode45解方程之前，再写一个函数：事件触发函数eventfun,
它的格式固定要返回三个值，这三个值的意思是
当第一个值vaule的值到达0时，时间会触发
而根据第二个值isterminal设置，触发时间会否终止求解
第三个值是设置触发过0的方向

function [value,isterminal,direction] = eventfun(t,x)
value=x(1)-100; %触发时间，当其值为0的时候，时间会触发
isterminal=1; %设为1时会，触发时间会停止求解器，设0时触发不影响工作
direction=1; %触发方向设1时是上升触发，设-1是下降触发，设0是双向触发
end

写好fun和eventfun之后，你就可以调用ode45解方程了
但用ode45之前记得先用odeset，将触发函数加入哦

op=odeset(‘Events’,@eventfun);
[T,X,Tend,Xend,evennum]=ode45(@fun,[0,15],[0 0],op);

这样你看看，到达100米时，求解器就停住了
细心的你注意，ode45多返回了Tend,Xend,evennum三个参数
第一个Tend是触发事件发生的时间
第二个Xend是触发时间发生时刻的X
第三个evenum是标识触发事件的编号
由于这里只设置了一个触发事件，所以编号肯定是1

其实你也可以将该时间isterminal设为0
那么求解器不会因为触发事件而停止15秒内的解都会算出
你仍旧可以根据Tend 和 Xend得到到达100米时的时间和状态

额外提示，有时候你不知道到底取多大的时间范围才能够等到你要的触发时间
那么你可以用循环判断的方法，先设一个时间范围，然后求解
到最后都没有等到触发事件（Tend和Xend都是空矩阵）
那么适当延长求解时间区间，将上次的最后时刻和状态作为初始状态
再一次求解时间范围更大的解，如此直到找到触发事件才停止