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

1 次查看（过去 30 天）
KIPROTICH KOSGEY
KIPROTICH KOSGEY 2023-7-9
评论： Steven Lord 2023-9-1
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
John D'Errico 2023-7-9
Have something like this be the last thing you do in the function.
x = max(x,0);
Or,
x(x<0) = 0;
  3 个评论
显示 1更早的评论
Steven Lord
Steven Lord 2023-9-1
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.

请先登录，再进行评论。

更多回答（0 个）

类别

在 Help Center 和 File Exchange 中查找有关 Ordinary Differential Equations 的更多信息

标签

产品


版本

R2023a

Community Treasure Hunt

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

Start Hunting!

Translated by