编程时用到了ode中的events，用的人很少，贴出来跟大家讨论一下
既然是用ode.当然最先要贴出来的是我的微分方程组。这个应该没什么问题，大家都会，只不过我的自由度稍微多一点。function dy=odefun(t,y)
global w0 w v z f g p
dy=zeros(12,1);
dy(1)=y(2);
dy(2)=p(1)*sin(w0*t+z(1))-w(1)*w(1)*y(1)-2*f(1)*w(1)*y(2);
dy(3)=y(4);
dy(4)=-v(1)*v(1)*y(3)-2*g(1)*v(1)*y(4);
dy(5)=y(6);
dy(6)=p(2)*sin(w0*t+z(2))-w(2)*w(2)*y(5)-2*f(2)*w(2)*y(6);
dy(7)=y(8);
dy(8)=-v(2)*v(2)*y(7)-2*g(2)*v(2)*y(8);
dy(9)=y(10);
dy(10)=p(3)*sin(w0*t+z(3))-w(3)*w(3)*y(9)-2*f(3)*w(3)*y(10);
dy(11)=y(12);
dy(12)=-v(3)*v(3)*y(11)-2*g(3)*v(3)*y(12);
-------------------------------------------------------------------------------------------分割线
-------------------------------------------------------------------------------------------分割线
接下来就是我的主程序啦
initialization; %给上面odefun里面的矩阵赋值
global d1 d2
tspan=; %这些都是ode常用的代号，就是求解范围，初值
y0=zeros(12,1);
options=odeset('events',@events);
tout=tspan(1);
yout=y0';
teout=[];
yeout=[];
for i=1:10 %作十次运算，每次运算到满足events条件就停止
=ode45(@odefun,tspan,y0,options);
end
---------------------------------------------------------------------分割线
----------------------------------------------------------------------分割线
最后是events的m文件。这个就是出错的地方了
function =events(t,y) %当value等于1的时候，isterminal等于1，这时候ode停止计算。
global d1 d2 d k0 %下面两个d1,d2都是用来判断value什么时候等于1
d1=d-(k0*y(end,1)-y(end,3))/k0+(k0*y(end,5)-y(end,7))/k0; %direction为1的时候，表示value从负到正，否则value从正到负
d2=d-(k0*y(end,5)-y(end,7))/k0+(k0*y(end,9)-y(end,7))/k0;
if d1>0 & d2>0
value=1;
else
value=0;
end
isterminal=1;
direction=-1;
运行main的时候，报错为
Attempted to access y(12,3); index out of bounds because size(y)=.
这是为什么呢？ 已经解决了。。。
在events中不要直接用y(end,3)来运算。
首先，在odefun中，定义几个参数
function dy=odefun(t,y)
global w0 w v z f g p x1 x2 x3 y1 y2 y3 Vx1 Vx2 Vx3 Vy1 Vy2 Vy3
dy=zeros(12,1);
x1=y(1);
Vx1=y(2);
y1=y(3);
Vy1=y(4);
x2=y(5);
Vx2=y(6);
y2=y(7);
Vy2=y(8);
x3=y(9);
Vx3=y(10);
y3=y(11);
Vy3=y(12);dy(1)=y(2);
dy(2)=p(1)*sin(w0*t+z(1))-w(1)*w(1)*y(1)-2*f(1)*w(1)*y(2);
dy(3)=y(4);
dy(4)=-v(1)*v(1)*y(3)-2*g(1)*v(1)*y(4);
dy(5)=y(6);
dy(6)=p(2)*sin(w0*t+z(2))-w(2)*w(2)*y(5)-2*f(2)*w(2)*y(6);
dy(7)=y(8);
dy(8)=-v(2)*v(2)*y(7)-2*g(2)*v(2)*y(8);
dy(9)=y(10);
dy(10)=p(3)*sin(w0*t+z(3))-w(3)*w(3)*y(9)-2*f(3)*w(3)*y(10);
dy(11)=y(12);
dy(12)=-v(3)*v(3)*y(11)-2*g(3)*v(3)*y(12);
-----------------------------------------------------------------------------------------------分割线
-----------------------------------------------------------------------------------------------分割线
events.m也改一下
function =events(t,y)
global d1 d2 d k0 x1 x2 x3 y1 y2 y3
d1=d-(k0*x1-y1)/k0+(k0*x2-y2)/k0;
d2=d-(k0*x2-y2)/k0+(k0*x3-y3)/k0;
if d1>0 & d2>0
value=1;
else
value=0;
end
isterminal=1;
direction=-1;
虽然这样改了之后可以运行，但是不知道为什么要这样？
页:
[1]