请高手帮忙看看我的程序错在什么地方？帮忙修改，谢谢了。
本帖最后由 xueyongzhou 于 2010-8-19 23:18 编辑请高手帮忙看看我的程序错在什么地方？帮忙修改，谢谢了。
events.m
function =events(t,y)
global a1 b1 d1 a2 b2 d2
value=(a2*y(2)*y(3)/(1.0+b2*y(2)) - d2*y(3));
isterminal=0;
direction=0;
process.m
clear all;
global a1 b1 d1 a2 b2 d2
a1=5.; a2=0.1; b1=3.; b2=2.; d1=0.25;
% d2=0.0125;%bif par
t0=0.;
tf=20000.;
%y0=
y0=
datout=;
trans=800.;
%for d2=0.01434:-0.00002:0.0129
for d2=0.01:-0.00002:0.008
nsol=size(datout);
y0=datout(nsol(1),2:4);
options = odeset('Events',@events,'RelTol',1e-7,'AbsTol',);
=ode45(@rm,,y0,options);
npar=size(TE);
par=d2*ones(npar(1),1);
format long
funcout=;
datoutj=;
datouti=datoutj(trans:npar(1),:);
datout=cat(1,datout,datouti);
end
%save 'funcout.dat' funcout -ASCII
save 'datout.dat' datout -ASCII
disp('ready')
rm.m
function f = rm(t,y)
global a1 b1 d1 a2 b2 d2
f = zeros(3,1); % a column vector
f(1)=(y(1)*(1.0-y(1)) - a1*y(1)*y(2)/(1.0+b1*y(1)));
f(2)=(a1*y(1)*y(2)/(1.0+b1*y(1)) - d1*y(2) - a2*y(2)*y(3)/(1.0+b2*y(2)));
f(3)=(a2*y(2)*y(3)/(1.0+b2*y(2)) - d2*y(3));
页:
[1]