Setting min value of variables within a function
1 view (last 30 days)
Show older comments
Kosgey Kip
on 15 Oct 2019
Commented: Walter Roberson
on 16 Oct 2019
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
0 Comments
Accepted Answer
Fabio Freschi
on 15 Oct 2019
Edited: Fabio Freschi
on 15 Oct 2019
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 Comments
More Answers (1)
Walter Roberson
on 15 Oct 2019
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 Comments
Walter Roberson
on 16 Oct 2019
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.
See Also
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!