Setting min value of variables within a function

1 次查看（过去 30 天）
Kosgey Kip
Kosgey Kip 2019-10-15
评论： Walter Roberson 2019-10-16
Hi,
Please, how can i set a minimum within a function to zero?
I realise matlab is calculation with negative values and in the process returning some negative figures, so i would like to fix the min values of the variables to zero. this is a function of many variables that i am trying to solve with ode45.
function fval=MBBRFun4(t,y)
%Define the three variables
Xaob=y(1);
Xnb=y(2);
Xnsp=y(3);
Xcmx=y(4);
Xamx=y(5);
Xhan=y(6);
Xhaer=y(7);
Xs=y(8);
Sse=y(9);
Sno3=y(10);
Sno2=y(11);
Snh4=y(12);
So2=y(13);
% main functions
fval(1,1)=(0-V*Xaob/100)/V + u1*(So2/(Ko2aob+So2))*(Snh4/(Knh4aob+Snh4))*Xaob - baob*(So2/(Ko2aob+So2))*Xaob - baob*naob*(Ko2aob/(Ko2aob+So2))*(Sno2+Sno3)/(Kno3aob+Sno2+Sno3)*Xaob;
fval(2,1)=(0-V*Xnb/100)/V + u2*(Sno2/(Kno2nb+Sno2))*(So2/(Ko2nb+So2))*Xnb - bnb*(So2/(Ko2nb+So2))*Xnb - bnb*nnb*(Ko2nb/(Ko2nb+So2))*(Sno2+Sno3)/(Kno3nsp+Sno2+Sno3)*Xnb;
fval(3,1)=(0-V*Xnsp/100)/V + u3*(Sno2/(Kno2nsp+Sno2))*(So2/(Ko2nsp+So2))*Xnsp - bnsp*(So2/(Ko2nsp+So2))*Xnsp - bnsp*nnsp*(Ko2nsp/(Ko2nsp+So2))*(Sno2+Sno3)/(Kno3nsp+Sno2+Sno3)*Xnsp;
fval(4,1)=(0-V*Xcmx/100)/V + u4*((Sno2/(Kno2cmx+Sno2))*(So2/(Ko2cmx+So2)) + u4*(So2/(Ko2cmx+So2))*(Snh4/(Knh4cmx+Snh4)))*Xcmx - bcmx*(So2/(Ko2cmx+So2))*Xcmx - bcmx*ncmx*(Ko2cmx/(Ko2cmx+So2))*(Sno2+Sno3)/(Kno3cmx+Sno2+Sno3)*Xcmx;
fval(5,1)=(0-V*Xamx/100)/V + u5*(Snh4/(Knh4amx+Snh4))*(Sno2/(Kno2amx+Sno2))*(Ko2amx/(Ko2amx+So2))*Xamx-bamx*(So2/(Ko2amx+So2))*Xamx - bamx*namx*(Ko2amx/(Ko2amx+So2))*(Sno2+Sno3)/(Kno3amx+Sno2+Sno3)*Xamx;
fval(6,1)=(0-V*Xhan/100)/V + (u6*nhan*(Sse/(KSsehan+Sse))*(Sno2/(Kno2han+Sno2))*(Ko2han/(Ko2han+So2)) + u7*nhan*(Sse/(KSsehan+Sse))*(Sno3/(Kno3han+Sno3))*(Ko2han/(Ko2han+So2)))*Xhan - bhan*(So2/(Ko2han+So2))*Xhan - bhan*nhan*(Ko2han/(Ko2han+So2))*(Sno2+Sno3)/(Kno3han+Sno2+Sno3)*Xhan;
fval(7,1)=(0-V*Xhaer/100)/V + u8*(Sse/(KSsehaer+Sse))*(So2/(Ko2haer+So2))*Xhaer - bhaer*(So2/(Ko2haer+So2))*Xhaer - bhaer*nhaer*(Ko2haer/(Ko2haer+So2))*(Sno2+Sno3)/(Kno3haer+Sno2+Sno3)*Xhaer;
fval(8,1)=(0-V*Xs*100)/V - KH*(Xs/(Xhan+Xhaer))/(KX+(Xs/(Xhan+Xhaer)))*(Xhan+Xhaer) + (1-fI)*(baob*Xaob+bnb*Xnb+bnsp*Xnsp+bcmx*Xcmx+bamx*Xamx+bhan*Xhan+bhaer*Xhaer);
fval(9,1)=(Qo*Sseo-Qe*Sse)/V - (u6*nhan*(1/YNo2han)*(Sse/(KSsehan+Sse))*(Sno2/(Kno2han+Sno2))*(Ko2han/(Ko2han+So2)))*Xhan - (u7*nhan*1/YNo3han*(Sse/(KSsehan+Sse))*(Sno3/(Kno3han+Sno3))*(Ko2han/(Ko2han+So2)))*Xhan - (u8/Yhaer*(Sse/(KSsehaer+Sse))*(So2/(Ko2haer+So2))*Xhaer) + KH*(Xs/(Xhan+Xhaer))/(KX+(Xs/(Xhan+Xhaer)))*(Xhan+Xhaer);
fval(10,1)=(Qo*Sno3o-Qe*Sno3)/V - (u6*nhan*(1-YNo3han)/(2.86*YNo3han)*(Sse/(KSsehan+Sse))*(Sno3/(Kno3han+Sno3))*(Ko2han/(Ko2han+So2)))*Xhan + u5/1.14*(Snh4/(Knh4amx+Snh4))*(Sno2/(Kno2amx+Sno2))*(Ko2amx/(Ko2amx+So2))*Xamx + (u2/Ynb*(Sno2/(Kno2nb+Sno2))*(So2/(Ko2nb+So2))*Xnb) + (u3/Ynsp*(Sno2/(Kno2nsp+Sno2))*(So2/(Ko2nsp+So2))*Xnsp) + (u4/Ycmx*(Sno2/(Kno2cmx+Sno2))*(So2/(Ko2cmx+So2))*Xcmx) - ((1-fI)/2.86)*baob*naob*(Ko2aob/(Ko2aob+So2))*(Sno2+Sno3)/(Kno3aob+Sno2+Sno3)*Xaob - ((1-fI)/2.86)*bnb*nnb*(Ko2nb/(Ko2nb+So2))*(Sno2+Sno3)/(Kno3nb+Sno2+Sno3)*Xnb - ((1-fI)/2.86)*bnsp*nnsp*(Ko2nsp/(Ko2nsp+So2))*(Sno2+Sno3)/(Kno3nsp+Sno2+Sno3)*Xnsp - ((1-fI)/2.86)*bcmx*ncmx*(Ko2cmx/(Ko2cmx+So2))*(Sno2+Sno3)/(Kno3cmx+Sno2+Sno3)*Xcmx - ((1-fI)/2.86)*bamx*namx*(Ko2amx/(Ko2amx+So2))*(Sno2+Sno3)/(Kno3amx+Sno2+Sno3)*Xamx - ((1-fI)/2.86)*bhaer*nhaer*(Ko2haer/(Ko2haer+So2))*(Sno2+Sno3)/(Kno3han+Sno2+Sno3)*Xhaer - ((1-fI)/2.86)*bhan*nhan*(Ko2han/(Ko2han+So2))*(Sno2+Sno3)/(Kno3han+Sno2+Sno3)*Xhan;
fval(11,1)=(Qo*Sno2o-Qe*Sno2)/V - (((1-YNo2han)/(1.71*YNo2han))*u6*(Sse/(KSsehan+Sse))*(Sno2/(Kno2han+Sno2))*(Ko2han/(Ko2han+So2))*Xhan) - (((u5/Yamx)+u5/1.14)*(Snh4/(Knh4amx+Snh4))*(Sno2/(Kno2amx+Sno2))*(Ko2amx/(Ko2amx+So2))*Xamx) + (u1/Yaob*(So2/(Ko2aob+So2))*(Snh4/(Knh4aob+Snh4))*Xaob) - (u2/Ynb*(Sno2/(Kno2nb+Sno2))*(So2/(Ko2nb+So2))*Xnb) - (u3/Ynsp*(Sno2/(Kno2nb+Sno2))*(So2/(Ko2nb+So2))*Xnsp) - (u4/Ycmx*(Sno2/(Kno2nb+Sno2))*(So2/(Ko2nb+So2))*Xcmx);
fval(12,1)=(Qo*Snh4o-Qe*Snh4)/V - (u1*(INBM+1/Yaob)*(So2/(Ko2aob+So2))*(Snh4/(Knh4aob+Snh4))*Xaob) - (u4*(INBM+1/Ycmx)*(So2/(Ko2cmx+So2))*(Snh4/(Knh4cmx+Snh4))*Xcmx) - (u5*(INBM+1/Yamx)*(Snh4/(Knh4amx+Snh4))*(Sno2/(Kno2amx+Sno2))*(Ko2amx/(Ko2amx+So2))*Xamx) - (u2*INBM*(Sno2/(Kno2nb+Sno2))*(So2/(Ko2nb+So2))*Xnb) -(u3*INBM*(Sno2/(Kno2nsp+Sno2))*(So2/(Ko2nsp+So2))*Xnsp) - INBM*u6*(Sse/(KSsehan+Sse))*(Sno2/(Kno2han+Sno2))*(Ko2han/(Ko2han+So2))*Xhan - u7*INBM*(Sse/(KSsehan+Sse))*(Sno3/(Kno3han+Sno3))*(Ko2han/(Ko2han+So2))*Xhan - u8*INBM*(Sse/(KSsehaer+Sse))*(So2/(Ko2haer+So2))*Xhaer + (INBM-(fI*INXI))*baob*naob*(Ko2aob/(Ko2aob+So2))*(Sno2+Sno3)/(Kno3aob+Sno2+Sno3)*Xaob + (INBM-(fI*INXI))*bnb*nnb*(Ko2nb/(Ko2nb+So2))*(Sno2+Sno3)/(Kno3nsp+Sno2+Sno3)*Xnb + (INBM-(fI*INXI))*bnsp*nnsp*(Ko2nsp/(Ko2nsp+So2))*(Sno2+Sno3)/(Kno3nsp+Sno2+Sno3)*Xnsp + (INBM-(fI*INXI))*bcmx*ncmx*(Ko2cmx/(Ko2cmx+So2))*(Sno2+Sno3)/(Kno3cmx+Sno2+Sno3)*Xcmx + (INBM-(fI*INXI))*bamx*namx*(Ko2amx/(Ko2amx+So2))*(Sno2+Sno3)/(Kno3amx+Sno2+Sno3)*Xamx + (INBM-(fI*INXI))*bhan*nhan*(Ko2han/(Ko2han+So2))*(Sno2+Sno3)/(Kno3han+Sno2+Sno3)*Xhan + (INBM-(fI*INXI))*bhaer*nhaer*(Ko2haer/(Ko2haer+So2))*(Sno2+Sno3)/(Kno3haer+Sno2+Sno3)*Xhaer+(INBM-(fI*INXI))*baob*(So2/(Ko2aob+So2))*Xaob + (INBM-(fI*INXI))*bnb*(So2/(Ko2nb+So2))*Xnb + (INBM-(fI*INXI))*bnsp*(So2/(Ko2nsp+So2))*Xnsp + (INBM-(fI*INXI))*bcmx*(So2/(Ko2cmx+So2))*Xcmx + (INBM-(fI*INXI))*bamx*(So2/(Ko2amx+So2))*Xamx + (INBM-(fI*INXI))*bhan*(So2/(Ko2han+So2))*Xhan + (INBM-(fI*INXI))*bhaer*(So2/(Ko2haer+So2))*Xhaer;
fval(13,1)=(Qo*So2o-Qe*So2)/V + (KaL*(So2G-So2))-((1-Yhaer)/Yhaer)*u8*(Sse/(KSsehaer+Sse))*(So2/(Ko2haer+So2))*Xhaer-((3.43-Yaob)/Yaob)*u1*(So2/(Ko2aob+So2))*(Snh4/(Knh4aob+Snh4))*Xaob - ((1.14-Ynb)/Ynb)*u2*(Sno2/(Kno2nb+Sno2))*(So2/(Ko2nb+So2))*Xnb - ((1.14-Ynsp)/Ynsp)*u3*(Sno2/(Kno2nsp+Sno2))*(So2/(Ko2nsp+So2))*Xnsp - ((1.14-Ycmx)/Ycmx)*u4*((Sno2/(Kno2cmx+Sno2))*(So2/(Ko2cmx+So2)) - ((4.57-Ycmx)/Ycmx)*u4*(So2/(Ko2cmx+So2))*(Snh4/(Knh4cmx+Snh4)))*Xcmx - (1-fI)*baob*(So2/(Ko2aob+So2))*Xaob - (1-fI)*bnb*(So2/(Ko2nb+So2))*Xnb - (1-fI)*bnsp*(So2/(Ko2nsp+So2))*Xnsp - (1-fI)*bcmx*(So2/(Ko2cmx+So2))*Xcmx - (1-fI)*bamx*(So2/(Ko2amx+So2))*Xamx - (1-fI)*bhan*(So2/(Ko2han+So2))*Xhan - (1-fI)*bhaer*(So2/(Ko2haer+So2))*Xhaer;
script
%initial conditions
y0=[0.02352;0.00048;0.0288;0.00432;0.216;0.1104;0.1104;2;2;2;55;50;0.5];
h=0.0006944444;
tSpan=[0 535];
[tSol, ySol]=ode15s(@(t, y) MBBRFun4(t,y), tSpan, y0);
Please people assist
Thanks

请先登录，再回答此问题。

采纳的回答

Fabio Freschi
Fabio Freschi 2019-10-15
编辑：Fabio Freschi 2019-10-15
Maybe I don't understand correctly: if you want to set all negative entries of fval to 0 you can add
fval(fval < 0) = 0;
at the end of your function. Otherwise follow Walter answer.
  4 个评论
显示 2更早的评论
Kosgey Kip
Kosgey Kip 2019-10-16
Thanks but this is also giving some ambiguously high values

请先登录，再进行评论。

更多回答（1 个）

Walter Roberson
Walter Roberson 2019-10-15
https://www.mathworks.com/help/matlab/ref/odeset.html#d117e993600
You can set the NonNegative option to the list of variable indices that must not be negative.
However, the intention for this is strictly for variables that might become negative due to arithmetic round-off, and which would not give negative in infinite precision calculations. If that is not your situation, then instead of using NonNegative you should be using event functions to detect the component going negative and restarting the calculation as if you had "bounced" off of the negative section; see the ballode example.
  3 个评论
显示 1更早的评论
Kosgey Kip
Kosgey Kip 2019-10-16
@Walter, please advice where i have gone wrong here, i have created an event_function besides the main function, and my script now looks like this:
y0=[0.02352;0.00048;0.0288;0.00432;0.216;0.1104;0.1104;2;2;2;55;50;0.5];
tSpan=[0 535];
options = odeset('Events',@eventFun);
[tSol, ySol]=ode15s(@(t, y) MBBRFun4(t,y), tSpan, y0, options);
events_function
function [value,isterminal,direction] = eventFun(t,y)
%Define the three variables
Xaob=y(1);
Xnb=y(2);
Xnsp=y(3);
Xcmx=y(4);
Xamx=y(5);
Xhan=y(6);
Xhaer=y(7);
Xs=y(8);
Sse=y(9);
Sno3=y(10);
Sno2=y(11);
Snh4=y(12);
So2=y(13);
fval(1,1)=(0-V*Xaob/100)/V + u1*(So2/(Ko2aob+So2))*(Snh4/(Knh4aob+Snh4))*Xaob - baob*(So2/(Ko2aob+So2))*Xaob - baob*naob*(Ko2aob/(Ko2aob+So2))*(Sno2+Sno3)/(Kno3aob+Sno2+Sno3)*Xaob;
fval(2,1)=(0-V*Xnb/100)/V + u2*(Sno2/(Kno2nb+Sno2))*(So2/(Ko2nb+So2))*Xnb - bnb*(So2/(Ko2nb+So2))*Xnb - bnb*nnb*(Ko2nb/(Ko2nb+So2))*(Sno2+Sno3)/(Kno3nsp+Sno2+Sno3)*Xnb;
fval(3,1)=(0-V*Xnsp/100)/V + u3*(Sno2/(Kno2nsp+Sno2))*(So2/(Ko2nsp+So2))*Xnsp - bnsp*(So2/(Ko2nsp+So2))*Xnsp - bnsp*nnsp*(Ko2nsp/(Ko2nsp+So2))*(Sno2+Sno3)/(Kno3nsp+Sno2+Sno3)*Xnsp;
fval(4,1)=(0-V*Xnsp/100)/V+ u4*((Sno2/(Kno2cmx+Sno2))*(So2/(Ko2cmx+So2)) + u4*(So2/(Ko2cmx+So2))*(Snh4/(Knh4cmx+Snh4)))*Xcmx - bcmx*(So2/(Ko2cmx+So2))*Xcmx - bcmx*ncmx*(Ko2cmx/(Ko2cmx+So2))*(Sno2+Sno3)/(Kno3cmx+Sno2+Sno3)*Xcmx;
fval(5,1)= (0-V*Xcmx/100)/V + u5*(Snh4/(Knh4amx+Snh4))*(Sno2/(Kno2amx+Sno2))*(Ko2amx/(Ko2amx+So2))*Xamx-bamx*(So2/(Ko2amx+So2))*Xamx - bamx*namx*(Ko2amx/(Ko2amx+So2))*(Sno2+Sno3)/(Kno3amx+Sno2+Sno3)*Xamx;
fval(6,1)=(0-V*Xhan/100)/V + (u6*nhan*(Sse/(KSsehan+Sse))*(Sno2/(Kno2han+Sno2))*(Ko2han/(Ko2han+So2)) + u7*nhan*(Sse/(KSsehan+Sse))*(Sno3/(Kno3han+Sno3))*(Ko2han/(Ko2han+So2)))*Xhan - bhan*(So2/(Ko2han+So2))*Xhan - bhan*nhan*(Ko2han/(Ko2han+So2))*(Sno2+Sno3)/(Kno3han+Sno2+Sno3)*Xhan;
fval(7,1)=(0-V*Xhaer/100)/V + u8*(Sse/(KSsehaer+Sse))*(So2/(Ko2haer+So2))*Xhaer - bhaer*(So2/(Ko2haer+So2))*Xhaer - bhaer*nhaer*(Ko2haer/(Ko2haer+So2))*(Sno2+Sno3)/(Kno3haer+Sno2+Sno3)*Xhaer;
fval(8,1)=(0-V*Xs)/V - KH*(Xs/(Xhan+Xhaer))/(KX+(Xs/(Xhan+Xhaer)))*(Xhan+Xhaer) + (1-fI)*(baob*Xaob+bnb*Xnb+bnsp*Xnsp+bcmx*Xcmx+bamx*Xamx+bhan*Xhan+bhaer*Xhaer);
fval(9,1)=(Qo*Sseo-Qe*Sse)/V - (u6*nhan*(1/YNo2han)*(Sse/(KSsehan+Sse))*(Sno2/(Kno2han+Sno2))*(Ko2han/(Ko2han+So2)))*Xhan - (u7*nhan*1/YNo3han*(Sse/(KSsehan+Sse))*(Sno3/(Kno3han+Sno3))*(Ko2han/(Ko2han+So2)))*Xhan - (u8/Yhaer*(Sse/(KSsehaer+Sse))*(So2/(Ko2haer+So2))*Xhaer) + KH*(Xs/(Xhan+Xhaer))/(KX+(Xs/(Xhan+Xhaer)))*(Xhan+Xhaer);
fval(10,1)=(Qo*Sno3o-Qe*Sno3)/V - (u6*nhan*(1-YNo3han)/(2.86*YNo3han)*(Sse/(KSsehan+Sse))*(Sno3/(Kno3han+Sno3))*(Ko2han/(Ko2han+So2)))*Xhan + u5/1.14*(Snh4/(Knh4amx+Snh4))*(Sno2/(Kno2amx+Sno2))*(Ko2amx/(Ko2amx+So2))*Xamx + (u2/Ynb*(Sno2/(Kno2nb+Sno2))*(So2/(Ko2nb+So2))*Xnb) + (u3/Ynsp*(Sno2/(Kno2nsp+Sno2))*(So2/(Ko2nsp+So2))*Xnsp) + (u4/Ycmx*(Sno2/(Kno2cmx+Sno2))*(So2/(Ko2cmx+So2))*Xcmx) - ((1-fI)/2.86)*baob*naob*(Ko2aob/(Ko2aob+So2))*(Sno2+Sno3)/(Kno3aob+Sno2+Sno3)*Xaob - ((1-fI)/2.86)*bnb*nnb*(Ko2nb/(Ko2nb+So2))*(Sno2+Sno3)/(Kno3nb+Sno2+Sno3)*Xnb - ((1-fI)/2.86)*bnsp*nnsp*(Ko2nsp/(Ko2nsp+So2))*(Sno2+Sno3)/(Kno3nsp+Sno2+Sno3)*Xnsp - ((1-fI)/2.86)*bcmx*ncmx*(Ko2cmx/(Ko2cmx+So2))*(Sno2+Sno3)/(Kno3cmx+Sno2+Sno3)*Xcmx - ((1-fI)/2.86)*bamx*namx*(Ko2amx/(Ko2amx+So2))*(Sno2+Sno3)/(Kno3amx+Sno2+Sno3)*Xamx - ((1-fI)/2.86)*bhaer*nhaer*(Ko2haer/(Ko2haer+So2))*(Sno2+Sno3)/(Kno3han+Sno2+Sno3)*Xhaer - ((1-fI)/2.86)*bhan*nhan*(Ko2han/(Ko2han+So2))*(Sno2+Sno3)/(Kno3han+Sno2+Sno3)*Xhan;
fval(11,1)=(Qo*Sno2o-Qe*Sno2)/V - (((1-YNo2han)/(1.71*YNo2han))*u6*(Sse/(KSsehan+Sse))*(Sno2/(Kno2han+Sno2))*(Ko2han/(Ko2han+So2))*Xhan) - (((u5/Yamx)+u5/1.14)*(Snh4/(Knh4amx+Snh4))*(Sno2/(Kno2amx+Sno2))*(Ko2amx/(Ko2amx+So2))*Xamx) + (u1/Yaob*(So2/(Ko2aob+So2))*(Snh4/(Knh4aob+Snh4))*Xaob) - (u2/Ynb*(Sno2/(Kno2nb+Sno2))*(So2/(Ko2nb+So2))*Xnb) - (u3/Ynsp*(Sno2/(Kno2nb+Sno2))*(So2/(Ko2nb+So2))*Xnsp) - (u4/Ycmx*(Sno2/(Kno2nb+Sno2))*(So2/(Ko2nb+So2))*Xcmx);
fval(12,1)=(Qo*Snh4o-Qe*Snh4)/V - (u1*(INBM+1/Yaob)*(So2/(Ko2aob+So2))*(Snh4/(Knh4aob+Snh4))*Xaob) - (u4*(INBM+1/Ycmx)*(So2/(Ko2cmx+So2))*(Snh4/(Knh4cmx+Snh4))*Xcmx) - (u5*(INBM+1/Yamx)*(Snh4/(Knh4amx+Snh4))*(Sno2/(Kno2amx+Sno2))*(Ko2amx/(Ko2amx+So2))*Xamx) - (u2*INBM*(Sno2/(Kno2nb+Sno2))*(So2/(Ko2nb+So2))*Xnb) -(u3*INBM*(Sno2/(Kno2nsp+Sno2))*(So2/(Ko2nsp+So2))*Xnsp) - INBM*u6*(Sse/(KSsehan+Sse))*(Sno2/(Kno2han+Sno2))*(Ko2han/(Ko2han+So2))*Xhan - u7*INBM*(Sse/(KSsehan+Sse))*(Sno3/(Kno3han+Sno3))*(Ko2han/(Ko2han+So2))*Xhan - u8*INBM*(Sse/(KSsehaer+Sse))*(So2/(Ko2haer+So2))*Xhaer + (INBM-(fI*INXI))*baob*naob*(Ko2aob/(Ko2aob+So2))*(Sno2+Sno3)/(Kno3aob+Sno2+Sno3)*Xaob + (INBM-(fI*INXI))*bnb*nnb*(Ko2nb/(Ko2nb+So2))*(Sno2+Sno3)/(Kno3nsp+Sno2+Sno3)*Xnb + (INBM-(fI*INXI))*bnsp*nnsp*(Ko2nsp/(Ko2nsp+So2))*(Sno2+Sno3)/(Kno3nsp+Sno2+Sno3)*Xnsp + (INBM-(fI*INXI))*bcmx*ncmx*(Ko2cmx/(Ko2cmx+So2))*(Sno2+Sno3)/(Kno3cmx+Sno2+Sno3)*Xcmx + (INBM-(fI*INXI))*bamx*namx*(Ko2amx/(Ko2amx+So2))*(Sno2+Sno3)/(Kno3amx+Sno2+Sno3)*Xamx + (INBM-(fI*INXI))*bhan*nhan*(Ko2han/(Ko2han+So2))*(Sno2+Sno3)/(Kno3han+Sno2+Sno3)*Xhan + (INBM-(fI*INXI))*bhaer*nhaer*(Ko2haer/(Ko2haer+So2))*(Sno2+Sno3)/(Kno3haer+Sno2+Sno3)*Xhaer+(INBM-(fI*INXI))*baob*(So2/(Ko2aob+So2))*Xaob + (INBM-(fI*INXI))*bnb*(So2/(Ko2nb+So2))*Xnb + (INBM-(fI*INXI))*bnsp*(So2/(Ko2nsp+So2))*Xnsp + (INBM-(fI*INXI))*bcmx*(So2/(Ko2cmx+So2))*Xcmx + (INBM-(fI*INXI))*bamx*(So2/(Ko2amx+So2))*Xamx + (INBM-(fI*INXI))*bhan*(So2/(Ko2han+So2))*Xhan + (INBM-(fI*INXI))*bhaer*(So2/(Ko2haer+So2))*Xhaer;
fval(13,1)=(Qo*So2o-Qe*So2)/V + (KaL*(So2G-So2))-((1-Yhaer)/Yhaer)*u8*(Sse/(KSsehaer+Sse))*(So2/(Ko2haer+So2))*Xhaer-((3.43-Yaob)/Yaob)*u1*(So2/(Ko2aob+So2))*(Snh4/(Knh4aob+Snh4))*Xaob - ((1.14-Ynb)/Ynb)*u2*(Sno2/(Kno2nb+Sno2))*(So2/(Ko2nb+So2))*Xnb - ((1.14-Ynsp)/Ynsp)*u3*(Sno2/(Kno2nsp+Sno2))*(So2/(Ko2nsp+So2))*Xnsp - ((1.14-Ycmx)/Ycmx)*u4*((Sno2/(Kno2cmx+Sno2))*(So2/(Ko2cmx+So2)) - ((4.57-Ycmx)/Ycmx)*u4*(So2/(Ko2cmx+So2))*(Snh4/(Knh4cmx+Snh4)))*Xcmx - (1-fI)*baob*(So2/(Ko2aob+So2))*Xaob - (1-fI)*bnb*(So2/(Ko2nb+So2))*Xnb - (1-fI)*bnsp*(So2/(Ko2nsp+So2))*Xnsp - (1-fI)*bcmx*(So2/(Ko2cmx+So2))*Xcmx - (1-fI)*bamx*(So2/(Ko2amx+So2))*Xamx - (1-fI)*bhan*(So2/(Ko2han+So2))*Xhan - (1-fI)*bhaer*(So2/(Ko2haer+So2))*Xhaer;
value = y(:);
isterminal = 13;
direction = 0;
end
However, this throws out the following error:
Attempted to access isterminal(5); index out of bounds because numel(isterminal)=1.
Error in odezero (line 142) if any(isterminal(indzc))
Error in ode15s (line 815)
[te,ye,ie,valt,stop] =
odezero(@ntrp15s,eventFcn,eventArgs,valt,...
Error in MBBR2 (line 7)
[tSol, ySol]=ode15s(@(t, y) MBBRFun4(t,y),
tSpan, y0, options);
Please advce
Walter Roberson
Walter Roberson 2019-10-16
Do not set isterminal to the index of an event that is terminal. Instead set it to a vector that indicates whether each value is terminal.
Note that for your purposes, you probably only care about terminal events, so you probably do not even need to compute the other ones.

请先登录，再进行评论。

类别

在 Help Center 和 File Exchange 中查找有关 Solver Outputs and Iterative Display 的更多信息

标签

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by