From: Xiaobing Xu on
I have a set of ODEs as following:

function xb=generation(t,x)
global vl mds1 mds2 mds3 c1 c2 c3 fh nh cm fm km nm
vl=0.333;
mds1=667*0.59*0.64;
mds2=667*0.10*0.40;
mds3=667*0.10*0.10;
c1=0.00023;
c2=0.00013;
c3=0.00003;
fh=1.0;
kvfa=0.2;
nh=1.0;
cm=0.5;
fm=1.0;
km=0.023;
nm=1.0;
kdm=0.02;
xb=[nh*(mds1*c1*fh*exp(-kvfa*(x(1)-x(2)))*exp(-c1*fh*exp(-kvfa*(x(1)-x(2)))*t)+...
mds2*c2*fh*exp(-kvfa*(x(1)-x(2)))*exp(-c2*fh*exp(-kvfa*(x(1)-x(2)))*t)+...
mds3*c3*fh*exp(-kvfa*(x(1)-x(2)))*exp(-c3*fh*exp(-kvfa*(x(1)-x(2)))*t))/vl;...
cm*fm*(x(1)-x(2))*x(3)/(km+x(1)-x(2));...
(cm*(x(1)-x(2))/(km+x(1)-x(2))-kdm)*x(3)];
end

ODE45 was used to solve above questions:

[t,x]=ode45('generation',[0 36500],[0 0 0.3]);
plot(t,x(:,1),'-',t,x(:,2),'--',t,x(:,3),'-.');
xlabel('Time/day');
ylabel('Value');
legend('Gh','Gm','Gb');

However, the results showed that both the values of x(2) (or Gm) and x(3) (or Gb) become zero after several timesteps. Actually, there exists a solution which has positive values for x(2) and x(3). How can I get positive results? In matlab, 'odeset' only has 'nonnegative' option.

Sincere thanks for your help.

Xiaobing
From: Steven_Lord on


"Xiaobing Xu" <xxb_2002(a)163.com> wrote in message
news:i23fkb$mmj$1(a)fred.mathworks.com...
> I have a set of ODEs as following:
>
> function xb=generation(t,x)
> global vl mds1 mds2 mds3 c1 c2 c3 fh nh cm fm km nm
> vl=0.333; mds1=667*0.59*0.64;
> mds2=667*0.10*0.40;
> mds3=667*0.10*0.10; c1=0.00023; c2=0.00013; c3=0.00003;
> fh=1.0; kvfa=0.2; nh=1.0; cm=0.5;
> fm=1.0; km=0.023; nm=1.0;
> kdm=0.02;
> xb=[nh*(mds1*c1*fh*exp(-kvfa*(x(1)-x(2)))*exp(-c1*fh*exp(-kvfa*(x(1)-x(2)))*t)+...
>
> mds2*c2*fh*exp(-kvfa*(x(1)-x(2)))*exp(-c2*fh*exp(-kvfa*(x(1)-x(2)))*t)+...
>
> mds3*c3*fh*exp(-kvfa*(x(1)-x(2)))*exp(-c3*fh*exp(-kvfa*(x(1)-x(2)))*t))/vl;...
> cm*fm*(x(1)-x(2))*x(3)/(km+x(1)-x(2));...
> (cm*(x(1)-x(2))/(km+x(1)-x(2))-kdm)*x(3)];
> end
>
> ODE45 was used to solve above questions:
>
> [t,x]=ode45('generation',[0 36500],[0 0 0.3]);
> plot(t,x(:,1),'-',t,x(:,2),'--',t,x(:,3),'-.');
> xlabel('Time/day');
> ylabel('Value');
> legend('Gh','Gm','Gb');
>
> However, the results showed that both the values of x(2) (or Gm) and x(3)
> (or Gb) become zero after several timesteps. Actually, there exists a
> solution which has positive values for x(2) and x(3). How can I get
> positive results? In matlab, 'odeset' only has 'nonnegative' option.

Look at the BALLODE demo. It contains an example that uses an events
function to detect when one of the variables becomes 0, and handles that
appropriately by stopping the solver, modifying the initial conditions, and
restarting the solver from where it left off. You could probably do the
same.

--
Steve Lord
slord(a)mathworks.com
comp.soft-sys.matlab (CSSM) FAQ: http://matlabwiki.mathworks.com/MATLAB_FAQ
To contact Technical Support use the Contact Us link on
http://www.mathworks.com