From: Xiaobing Xu on 20 Jul 2010 02:27 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 20 Jul 2010 09:27 "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
|
Pages: 1 Prev: Help importing a Java class!!! Next: Invisible watermarking |