# How do I fix negative values in a function to zero ?

8 views (last 30 days)
KIPROTICH KOSGEY on 9 Jul 2023
Commented: Steven Lord on 1 Sep 2023
I have the following script and function with a code to fix the negative values to zero but it is failing for some reason. What could be the issue?
%Script:
%initial conditions
y0=[18.16;0.52;25.64;172.6;179.48;0.0005;2;0.0005;55;50;0.5];
tSpan=[0 535];
M1=eye(11);
options = odeset('RelTol',1e-20,'AbsTol',1e-10,'Mass',M1,'MassSingular', 'maybe','MStateDependence','strong','InitialStep',1e-10, 'MStateDependence','strong');
[tSol, ySol]=ode15s(@(t,y) MBBR12(t,y), tSpan, y0, options);
%Funtion as follows:
function fval=MBBR12(t,y)
%Define the three variables
Xaob=y(1);
Xnb=y(2);
Xnsp=y(3);
Xamx=y(4);
Xhan=y(5);
Xs=y(6);
Sse=y(7);
Sno3=y(8);
Sno2=y(9);
Snh4=y(10);
So2=y(11);
T=303;
%aob
O1=68/(0.00831*293*(273+T)); u1=0.054*exp(O1*(T+293)); Ko2aob=0.6; Knh4aob=2.4; Yaob=0.15; Kno3aob=0.5;
%nb
O2=44/(0.00831*293*(273+T)); u2=0.047*exp(O2*(T+293)); Ko2nb=2.2; Kno2nb=0.5*1.375; Ynb=0.041;Kno3nb=0.5;
%nsp
O3=44/(0.00831*293*(273+T)); u3=0.047*exp(O3*(T+293)); Ko2nsp=2.2; Kno2nsp=5.5; Ynsp=0.041;Kno3nsp=0.5;
%amx
O5=71/(0.00831*293*(273+T)); u5=0.0022*exp(O5*(T+293)); Kno2amx=0.05; Knh4amx=0.07; Yamx=0.159; Ko2amx=0.01; Kno3amx=0.5;
%han
O6=0.069;u6=7.2*exp((T+293)*0.069); KSsehan=2; Kno2han=0.5;Kno3han=0.5; YNo2han=0.3; Ko2han=0.2;
YNo3han=0.43; Yhaer=0.54;
%other constants
baob=0.0054*exp(O6*(T+293)); bnb=0.0025*exp(O6*(T+293));bnsp=bnb; bamx=0.00013*exp(O6*(T+293)); bhan=0.008*exp(O6*(T+293));
INBM=0.07; fI=0.1; namx=0.5; nhan=0.6;naob=0.5; nnb=0.5; nnsp=0.5;
KX=1; KH=0.125*exp(0.069*(T+293)); INXI=0.02;
% input parameters
Sseo=300; Sno3o=0.2; Qo=75; Xso=250;
Snh4o=650; Sno2o=0.2; So2o=0;
%oxygen transfer parameters
KaL=80;
V=1400; SRT=40;
%fix output values to zero
Sno3(Sno3<0)=0;Xaob(Xaob<0)=0;Xnb(Xnb<0)=0;Xnsp(Xnsp<0)=0;Xamx(Xamx<0)=0;Xhan(Xhan<0)=0;Xs(Xs<0)=0;
Sse(Sse<0)=0;Sno2(Sno2<0)=0;Snh4(Snh4<0)=0;So2(So2<0)=0;
% main functions CSTR 1
fval=[((0-V*Xaob/SRT)/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)
((0-V*Xnb/SRT)/V + u2*(Sno2/(Kno2nb+Sno2))*(So2/(Ko2nb+So2))*Xnb - bnb*(So2/(Ko2nb+So2))*Xnb - bnb*nnb*((Ko2nb/(Ko2nb+So2))*(Sno2+Sno3)/(Kno3nb+Sno2+Sno3))*Xnb)
((0-V*Xnsp/SRT)/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)
((0-V*Xamx/SRT)/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)
((0-V*Xhan/SRT)/V + (u6*(Sse/(KSsehan+Sse))*(Sno2/(Kno2han+Sno2))*(Ko2han/(Ko2han+So2)) + u6*(Sse/(KSsehan+Sse))*(Sno3/(Kno3han+Sno3))*(Ko2han/(Ko2han+So2))+u6*((Sse/(KSsehan+Sse))*(So2/(Ko2han+So2)))*Xhan)- bhan*(So2/(Ko2han+So2))*Xhan - bhan*nhan*((Ko2han/(Ko2han+So2))*(Sno2+Sno3)/(Kno3han+Sno2+Sno3))*Xhan)
((Qo*Xso-Qo*Xs)/V - KH*((Xs/Xhan)/(KX+(Xs/Xhan)))*Xhan)
((Qo*Sseo-Qo*Sse)/V - (u6*nhan*(1/YNo2han)*(Sse/(KSsehan+Sse))*(Sno2/(Kno2han+Sno2))*(Ko2han/(Ko2han+So2)))*Xhan -(u6*nhan*(1/YNo3han)*(Sse/(KSsehan+Sse))*(Sno3/(Kno3han+Sno3))*(Ko2han/(Ko2han+So2)))*Xhan - ((u6/Yhaer)*(Sse/(KSsehan+Sse))*(So2/(Ko2han+So2))*Xhan) + KH*((Xs/Xhan)/(KX+(Xs/Xhan)))*Xhan)
((Qo*Sno3o-Qo*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) - ((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)*bamx*namx*(Ko2amx/(Ko2amx+So2))*(Sno2+Sno3)/(Kno3amx+Sno2+Sno3)*Xamx-((1-fI)/2.86)*bhan*nhan*(Ko2han/(Ko2han+So2))*(Sno2+Sno3)/(Kno3han+Sno2+Sno3)*Xhan)
((Qo*Sno2o-Qo*Sno2)/V - (((1-YNo2han)/(1.71*YNo2han))*u6*nhan*(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/(Kno2nsp+Sno2))*(So2/(Ko2nsp+So2)))*Xnsp)
((Qo*Snh4o-Qo*Snh4)/V - (u1*(INBM+1/Yaob)*(So2/(Ko2aob+So2))*(Snh4/(Knh4aob+Snh4)))*Xaob -(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*nhan*((Sse/(KSsehan+Sse))*(Sno2/(Kno2han+Sno2))*(Ko2han/(Ko2han+So2)))*Xhan - u6*INBM*nhan*((Sse/(KSsehan+Sse))*(Sno3/(Kno3han+Sno3))*(Ko2han/(Ko2han+So2)))*Xhan - u6*INBM*((Sse/(KSsehan+Sse))*(So2/(Ko2han+So2)))*Xhan + (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))*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))*baob*(So2/(Ko2aob+So2))*Xaob + (INBM-(fI*INXI))*bnb*(So2/(Ko2nb+So2))*Xnb +(INBM-(fI*INXI))*bnsp*(So2/(Ko2nsp+So2))*Xnsp + (INBM-(fI*INXI))*bamx*(So2/(Ko2amx+So2))*Xamx + (INBM-(fI*INXI))*bhan*(So2/(Ko2han+So2))*Xhan)
((Qo*So2o-Qo*So2)/V +KaL*(8-So2)-(((1-Yhaer)/Yhaer)*u6*(Sse/(KSsehan+Sse))*(So2/(Ko2han+So2)))*Xhan -(((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-fI)*baob*(So2/(Ko2aob+So2))*Xaob - (1-fI)*bnb*(So2/(Ko2nb+So2))*Xnb -(1-fI)*bnsp*(So2/(Ko2nsp+So2))*Xnsp - (1-fI)*bamx*(So2/(Ko2amx+So2))*Xamx - (1-fI)*bhan*(So2/(Ko2han+So2))*Xhan)];
end

John D'Errico on 9 Jul 2023
Have something like this be the last thing you do in the function.
x = max(x,0);
Or,
x(x<0) = 0;
##### 3 CommentsShow 1 older commentHide 1 older comment
Moved: Bruno Luong on 1 Sep 2023
Thank You.
Steven Lord on 1 Sep 2023
For this particular scenario, where you're trying to constrain the solution to a system of ODEs to be non-negative, rather than calling max at the end of the ODE function I'd use the NonNegative option for ODEs if it's supported for the ODE solver you're using. See this example in the documentation.

### Categories

Find more on Programming in Help Center and File Exchange

R2023a

### Community Treasure Hunt

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

Start Hunting!