matlab微分方程ODE求解器的事件(Event)属性
在特定的微分方程求解过程中,比如碰撞、车辆刹车,这种特殊运动时间简单的时序求解不够完善,故需要用到一个ode求解器的事件(Event)属性
首先假定一个微分方程
dy1=y2
dy2=y1+1
其中y1不能超过4
求解改微分方程
event时间定义:
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);
figure
plot(t,y)
结果为:
events函数解析:
%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中的最后一列就是我想要的值;
X(end)
ans =