MATLAB: Setting min value of variables within a function

min valuesvariables

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

Best Answer

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.

Related Solutions

MATLAB: Warning: Rank deficient, rank = 2, tol = 7.141946e-13

Your function needed a lot of cleanup, and I had to guess about what was desired so you should review it.

With the current code, when you run, at time 0 you will hit a breakpoint in calculating dSno3dt at n = 94, at the new temp2 partial calculation (the lines were too long to debug... copying and pasting is messy in the editor when it exceeds your screen width.)

The reason for the breakpoint is that the calculation will overflow to -infinity . That -infinity will then mix with other calculations leading to an attempt to calculate infinity - infinity which gives nan. Therefore the later fval all return as NaN, leading to an integration failure at time 0.

Please have a careful look at the changes I made. I had to index a number of arrays, and you need to review whether I indexed properly for your calculation purposes. You should also double-check that I got the sign right on adding / subtracting the sub-expressions.

When you work on the code, I recommend that you break the remaining code even further into sub-expressions, to make it easier to debug. Chasing down which operator is causing a problem in a long line of operations is a pain.

By the way:

    %% conditions
    if le (t,46)

Your functions are discontinuous in time. That is not permitted for the ode* solvers: the mathematics of the Runge-Kutta solvers is only valid for continuous functions. You must terminate the ode15s call at each one of those time boundaries and re-start integration.

function fval=UASBFun6PDE(t,y)
    
    dbstop if naninf
    
    %% Define the variables
    Xaob=y(1);
    Xnb=y(2);
    Xnsp=y(3);
    Xcmx=y(4);
    Xamx=y(5);
    Xhan=y(6);
    Xs=y(7);
    Sse=y(8);
    Sno3=y(9);
    Sno2=y(10);
    Snh4=y(11);
    So2=y(12);
    
    %% parameters
    Xo=0; Xe=Xo; Xso=0;
    
    %aob
    O1=68/(0.00831*293*303); u1=1.443*exp(O1*15); Ko2aob=0.85*0.594; Knh4aob=0.55*2.4; Yaob=3*0.15; Kno3aob=0.5;
    %nb
    O2=44/(0.00831*293*303); u2=1.8*0.48*exp(O2*15); Ko2nb=0.98*0.435; Kno2nb=0.6*1.375; Ynb=1.1*0.13;Kno3nb=0.5;
    %nsp
    O3=44/(0.00831*293*303); u3=1*0.69*exp(O3*15); Ko2nsp=0.95*0.435; Kno2nsp=1.0*0.76; Ynsp=1.0*0.14;Kno3nsp=0.5;
    %cmx
    O4=44/(0.00831*293*303); u4=0.5*0.72*exp(O4*15); Ko2cmx=0.98*0.43; Kno2cmx=0.15*6.29;Knh4cmx=70*0.011; Ycmx=0.1*0.47;YcmxNo2=2*Ycmx;Kno3cmx=0.5;
    %amx
    O5=71/(0.00831*293*303); u5=1.2*0.0664*exp(O5*15); Kno2amx=1.2*0.175; Kno3amx=0.5; Knh4amx=2*0.185; Yamx=0.95*0.159; Ko2amx=0.01; Kno3amx=0.5;
    %han
    u6=7.2*exp(15*0.069); KSsehan=2; Kno2han=0.5; YNo2han=1.8*0.49; Ko2han=0.2;
    YNo3han=0.79; Kno3han=0.32;Yhaer=1.8*0.37;
    
    %other constants
    baob=0.1296*exp(0.11*15); bnb=0.069*exp(0.11*15);bnsp=0.06*exp(0.11*15);
    bcmx=0.06*exp(0.11*15); bamx=0.00312*exp(0.11*15); bhaer=0.192*exp(15*0.11);
    bhan=0.192*exp(0.11*15);
    INBM=0.0583; fI=0.08; namx=0.5; nhan=0.6;naob=0.5; nnb=0.5;nhaer=0.5;nnsp=0.5;ncmx=0.5;
    KX=0.03; KH=3*exp(0.069*15); INXI=0.02;
    
    % input parameters
    Sseo=2; Sno3o=0;
    
    %oxygen transfer
    KaL=222*(1.02^(-5));
    
    %% conditions
    if le (t,46)
        Snh4o=50; Sno2o=55; So2G=1;So2o=2; V=0.0038/5; Qo=0.000864; Qe=Qo;
    elseif and (gt (t,46),le (t,65))
        Snh4o=60.4; Sno2o=25.3; So2G=1;So2o=8.5; V=0.0038/5; Qo=0.000864;  Qe=Qo;
    elseif and (gt (t,65) , le (t,67))
        Snh4o=60.4; Sno2o=25.3; So2G=1; So2o=8.5; V=0.0038/5; Qo=0.000864;  Qe=Qo;
    elseif and (gt (t,67) , le (t,74))
        Snh4o=76.72; Sno2o=45.25; So2G=1;So2o=8.5; V=0.0038/5; Qo=0.002016; Qe=Qo;
    elseif and (gt (t,74) , le (t,79))
        Snh4o=76.72; Sno2o=45.25; So2G=1; So2o=8.5; V=0.0038/5; Qo=0.002592; Qe=Qo;
    elseif and (gt (t,80) , le (t,86))
        Snh4o=76.72; Sno2o=45.25; So2o=8.5; V=0.005/5; Qo=0.002592; Qe=Qo; So2G=6.7521;  %So2G=14.65-0.41*35+7.99*(10^-3)*(35^2)-7.78*(10^-5)*35^3;
    elseif and (gt (t,86) , le (t,87))
        Snh4o=76.72; Sno2o=45.25; So2G=0.1;So2o=0.2; V=0.005/5; Qo= 0.002736; Qe=Qo;
    elseif and (gt (t,87) , le (t,102))
        Snh4o=84.9; Sno2o=50.9; So2G=0.1;So2o=0.2; V=0.005/5; Qo=0.002736; Qe=Qo;
    elseif and (gt (t,102) , le (t,108))
        Snh4o=76.72; Sno2o=45.25; So2G=0.1;So2o=0.2; V=0.005/5; Qo=0.003168; Qe=Qo;
    elseif and (gt (t,108) , le (t,120))
        Snh4o=76.72; Sno2o=45.25; So2G=0.1;So2o=0.2; V=0.005/5; Qo=0.004176; Qe=Qo;
    elseif and (gt (t,121) , le (t,134))
        Snh4o=75.1; Sno2o=62.8; So2G=0.1;So2o=0.2; V=0.005/5; Qo=0.0054; Qe=Qo;
    elseif and (gt (t,134) , le (t,137))
        Snh4o=61.2; Sno2o=45.7; So2G=0.1;So2o=0.2; V=0.005/5; Qo=0.0054; Qe=Qo;
    elseif and (gt (t,137) , le (t,151))
        Snh4o=61.2; Sno2o=45.7; So2G=0.1;So2o=0.2; V=0.005/5; Qo=0.006120; Qe=Qo;
    elseif and (gt (t,151) , le (t,159))
        Snh4o=61.2; Sno2o=45.7; So2G=0.1;So2o=0.2; V=0.005/5; Qo=0.007056; Qe=Qo;
    elseif and (gt (t,159) , le (t,176))
        Snh4o=68.7; Sno2o=49.7; So2G=0.1;So2o=0.2; V=0.005/5; Qo=0.007920; Qe=Qo;
    elseif and (gt (t,176) , le (t,177))
        Snh4o=68.7; Sno2o=49.7; So2G=0.1;So2o=0.2; V=0.005/5; Qo=0.006120; Qe=Qo;
    elseif and (gt (t,177) , le (t,180))
        Snh4o=68.7; Sno2o=49.7; So2G=0.1;So2o=0.2; V=0.005/5; Qo=0.006120; Qe=Qo;
    elseif and (gt (t,180) , le (t,197))
        Snh4o=50.4; Sno2o=42; So2G=0.1;So2o=0.2; V=0.0038/5; Qo=0.006120; Qe=Qo;
    elseif and (gt (t,197) , le (t,208))
        Snh4o=59.3; Sno2o=58.2; So2G=0.1;So2o=0.2; V=0.0038/5; Qo=0.006120; Qe=Qo;
    elseif and (gt (t,208) , le (t,212))
        Snh4o=59.3; Sno2o=58.2; So2G=0.1;So2o=0.2; V=0.0038/5; Qo=0.007920; Qe=Qo;
    elseif and (gt (t,212) , le (t,214))
        Snh4o=59.3; Sno2o=58.2; So2G=0.1;So2o=0.2; V=0.0038/5; Qo=0.007920; Qe=Qo;
    elseif and (gt (t,214) , le (t,221))
        Snh4o=65; Sno2o=75; So2G=0.1;So2o=0.2; V=0.0038/5; Qo=0.007920; Qe=Qo;
    elseif and (gt (t,221) , le (t,225))
        Snh4o=65; Sno2o=75; So2G=0.1;So2o=0.2; V=0.0038/5; Qo=0.007056; Qe=Qo;
    elseif and (gt (t,225) , le (t,228))
        Snh4o=65; Sno2o=75; So2G=0.1;So2o=0.2; V=0.0038/5; Qo=0.007920; Qe=Qo;
    elseif and (gt (t,228) , le (t,237))
        Snh4o=72.8; Sno2o=89.6; So2G=0.1;So2o=0.2; V=0.0038/5; Qo=0.007920; Qe=Qo;
    elseif and (gt (t,237) , le (t,245))
        Snh4o=99.7; Sno2o=116.3; So2G=0.1;So2o=0.2; V=0.0038/5; Qo=0.007920; Qe=Qo;
    elseif and (gt (t,245) , lt (t,250))
        Snh4o=99.7; Sno2o=116.3; So2G=0.1;So2o=0.2; V=0.0038/5; Qo=0.007920; Qe=Qo;
    elseif eq(t,250)
        Snh4o=99.7; Sno2o=116.3; So2G=0.1;So2o=0.2; V=0.0038/5; Qo=0.007920; Qe=Qo;
    elseif and (gt (t,250) , le (t,254))
        Snh4o=112.19; Sno2o=145.38; So2G=0.1;So2o=0.2; V=0.0038/5; Qo=0.01872; Qe=Qo;
    elseif and (gt (t,254) , le (t,256))
        Snh4o=112.19; Sno2o=145.38; So2G=0.1;So2o=0.2; V=0.0038/5; Qo=0.007920; Qe=Qo;
    elseif and (gt (t,256) , le (t,264))
        Snh4o=112.19; Sno2o=145.38; So2G=0.1;So2o=0.2; V=0.0038/5; Qo=0.006588; Qe=Qo;
    elseif and (gt (t,264) , le (t,278))
        Snh4o=112.19; Sno2o=145.38; So2G=0.1;So2o=0.2; V=0.0038/5; Qo=0.007056; Qe=Qo;
    elseif and (gt (t,278) , le (t,288))
        Snh4o=108.43; Sno2o=140.786; So2G=0.1;So2o=0.2; V=0.0038/5; Qo=0.007056; Qe=Qo;
    elseif and (gt (t,288) , le (t,295))
        Snh4o=108.43; Sno2o=140.786; So2G=0.1;So2o=0.2; V=0.005/5; Qo=0.007920; Qe=Qo;
    elseif and (gt (t,295) , le (t,302))
        Snh4o=140.43; Sno2o=196.4; So2G=0.1;So2o=0.2; V=0.005/5; Qo=0.007920; Qe=Qo;
    elseif and (gt (t,302) , le (t,311))
        Snh4o=160; Sno2o=223.85; So2G=0.1;So2o=0.2; V=0.005/5; Qo=0.007920; Qe=Qo;
    elseif and (gt (t,311) , le (t,313))
        Snh4o=180.64; Sno2o=252.08; So2G=0.1;So2o=0.2; V=0.005/5; Qo=0.007920; Qe=Qo;
    elseif and (gt (t,313) , le (t,322))
        Snh4o=180.64; Sno2o=252.08; So2G=0.1;So2o=0.2; V=0.005/5; Qo=0.01; Qe=Qo;
    elseif and (gt (t,322) , le (t,327))
        Snh4o=200; Sno2o=280; So2G=0.1;So2o=0.2; V=0.005/5; Qo=0.01; Qe=Qo;
    elseif and (gt (t,327) , le (t,355))
        Snh4o=220; Sno2o=308; So2G=0.1;So2o=0.2; V=0.005/5; Qo=0.01; Qe=Qo;
    elseif and (gt (t,355) , le (t,425))
        Snh4o=200; Sno2o=264; So2G=0.1;So2o=0.2; V=0.005/5; Qo=0.01; Qe=Qo;
    elseif and (gt (t,425) , le (t,434))
        Snh4o=110; Sno2o=145; So2G=0.1;So2o=0.2; V=0.005/5; Qo=0.0166667; Qe=Qo;
    elseif and (gt (t,434) , le (t,475))
        Snh4o=200; Sno2o=264; So2G=0.1;So2o=0.2; V=0.005/5; Qo=0.005; Qe=Qo;
    elseif and (gt (t,475) , le (t,491))
        Snh4o=200; Sno2o=264; So2G=0.1;So2o=0.2; V=0.005/5; Qo=0.00667; Qe=Qo;
    else
        Snh4o=200; Sno2o=264; So2G=0.1;So2o=0.2; V=0.005/5; Qo=0.01; Qe=Qo;
    end
    
    %% ODE
    fval(1)=((0-V*Xaob/250)/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)=((0-V*Xnb/250)/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);
    
    fval(3)=((0-V*Xnsp/250)/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)=((0-V*Xcmx/250)/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)=((0-V*Xamx/250)/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)=((0-V*Xhan/250)/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);
    
    fval(7)=((0-V*Xs/250)/V - KH*((Xs/Xhan)/(KX+(Xs/Xhan)))*Xhan + (1-fI)*(baob*Xaob+bnb*Xnb+bnsp*Xnsp+bcmx*Xcmx+bamx*Xamx+bhan*Xhan));
    
    %% PDE
    L=0.31;
    N=535/1;
    J=5;
    hz=L/J;
    for j=1:J+1; Sse(1,j)=0.0005; end
    for n=1:N+1; Sse(n,1)=Sseo; end
    
    for j=1:J+1; Sno3(1,j)=0.0005; end
    for n=1:N+1; Sno3(n,1)=Sno3o; end
    
    for j=1:J+1; Sno2(1,j)=55; end
    for n=1:N+1; Sno2(n,1)=Sno2o; end
    
    for j=1:J+1; Snh4(1,j)=50; end
    for n=1:N+1; Snh4(n,1)=Snh4o; end
    
    for j=1:J+1; So2(1,j)=0.5; end
    for n=1:N+1; So2(n,1)=So2o; end
    
    dSsedt=zeros(N+1,J+1);
    dSno3dt=zeros(N+1,J+1);
    dSno2dt=zeros(N+1,J+1);
    dSnh4dt=zeros(N+1,J+1);
    dSo2dt=zeros(N+1,J+1);
    Dma=3;
    U=1;
    
    for n=1:N
        for j=2:J
            temp1 = -U*(dSsedt(n+1,j)-dSsedt(n,j));
            temp2 = Dma*(dSsedt(n,j+1)-2*dSsedt(n,j)+(dSsedt(n,j-1)))/hz^2;
            temp3 = (u6*nhan*(1/YNo2han)*(Sse(n,J)./(KSsehan+Sse(n,j)))*(Sno2(n,j)./(Kno2han+Sno2(n,j)))*(Ko2han/(Ko2han+So2(n,j))))*Xhan;
            temp4 = (u6*nhan*(1/YNo3han)*(Sse(n,j)/(KSsehan+Sse(n,j)))*(Sno3(n,j)/(Kno3han+Sno3(n,j)))*(Ko2han/(Ko2han+So2(n,j))))*Xhan;
            dSsedt(n+1,j)= temp1 + temp2 - temp3 - temp4 - ((u6/Yhaer)*(Sse(n,j)/(KSsehan+Sse(n,j)))*(So2(n,j)/(Ko2han+So2(n,j)))*Xhan) + KH*(Xs/(Xhan))/(KX+(Xs/(Xhan)))*(Xhan)+0.0*(Xaob+Xnb+Xnsp+Xcmx+Xamx+Xhan);    
            
            temp1 = -U*(dSno3dt(n+1,j)-dSno3dt(n,j));
            temp2 = Dma*(dSno3dt(n,j+1)-2*dSno3dt(n,j)+dSno3dt(n,j-1))/hz^2;
            temp3 = (u6*nhan*((1-YNo3han)/(2.86*YNo3han))*(Sse(n,j)/(KSsehan+Sse(n,j)))*(Sno3(n,j)/(Kno3han+Sno3(n,j)))*(Ko2han/(Ko2han+So2(n,j))))*Xhan;
            temp4 = (u5/1.14)*(Snh4(n,j)/(Knh4amx+Snh4(n,j)))*(Sno2(n,j)/(Kno2amx+Sno2(n,j)))*(Ko2amx/(Ko2amx+So2(n,j)))*Xamx;
            dSno3dt(n+1,j) = temp1 + temp2 - temp3 + temp4 + (u2/Ynb*(Sno2(n,j)/(Kno2nb+Sno2(n,j)))*(So2(n,j)./(Ko2nb+So2(n,j)))*Xnb) + ((u3/Ynsp)*(Sno2(n,j)/(Kno2nsp+Sno2(n,j)))*(So2(n,j)./(Ko2nsp+So2(n,j)))*Xnsp) + ((u4/YcmxNo2)*(Sno2(n,j)/(Kno2cmx+Sno2(n,j)))*(So2(n,j)/(Ko2cmx+So2(n,j)))*Xcmx)+((u4/Ycmx)*(Snh4(n,j)/(Knh4cmx+Snh4(n,j)))*(So2(n,j)./(Ko2cmx+So2(n,j)))*Xcmx) - ((1-fI)/2.86)*baob*naob*(Ko2aob/(Ko2aob+So2(n,j)))*(Sno2(n,j)+Sno3(n,j))/(Kno3aob+Sno2(n,j)+Sno3(n,j))*Xaob - ((1-fI)/2.86)*bnb*nnb*(Ko2nb/(Ko2nb+So2(n,j)))*(Sno2(n,j)+Sno3(n,j))/(Kno3nb+Sno2(n,j)+Sno3(n,j))*Xnb - ((1-fI)/2.86)*bnsp*nnsp*(Ko2nsp/(Ko2nsp+So2(n,j)))*(Sno2(n,j)+Sno3(n,j))/(Kno3nsp+Sno2(n,j)+Sno3(n,j))*Xnsp - ((1-fI)/2.86)*bcmx*ncmx*(Ko2cmx./(Ko2cmx+So2(n,j)))*(Sno2(n,j)+Sno3(n,j))/(Kno3cmx+Sno2(n,j)+Sno3(n,j))*Xcmx -  ((1-fI)/2.86)*bamx*namx*(Ko2amx/(Ko2amx+So2(n,j)))*(Sno2(n,j)+Sno3(n,j))/(Kno3amx+Sno2(n,j)+Sno3(n,j))*Xamx-((1-fI)/2.86)*bhan*nhan*(Ko2han/(Ko2han+So2(n,j)))*(Sno2(n,j)+Sno3(n,j))/(Kno3han+Sno2(n,j)+Sno3(n,j))*Xhan;
            
            temp1 = -U*(dSno2dt(n+1,j)-dSno2dt(n,j));
            temp2 = Dma*(dSno2dt(n,j+1)-2*dSno2dt(n,j)+dSno2dt(n,j-1))/hz^2;
            temp3 = (((1-YNo2han)/(1.71*YNo2han))*u6*nhan*(Sse(n,j)/(KSsehan+Sse(n,j)))*(Sno2(n,j)/(Kno2han+Sno2(n,j)))*(Ko2han/(Ko2han+So2(n,j)))*Xhan);
            temp4 = (((u5/Yamx)+u5/1.14)*(Snh4(n,j)/(Knh4amx+Snh4(n,j)))*(Sno2(n,j)/(Kno2amx+Sno2(n,j)))*(Ko2amx/(Ko2amx+So2(n,j))))*Xamx;
            dSno2dt(n+1,j) = temp1 + temp2 - temp3 - temp4 + ((u1/Yaob)*(So2(n,j)/(Ko2aob+So2(n,j)))*(Snh4(n,j)/(Knh4aob+Snh4(n,j))))*Xaob - ((u2/Ynb)*(Sno2(n,j)/(Kno2nb+Sno2(n,j)))*(So2(n,j)/(Ko2nb+So2(n,j)))*Xnb) - ((u3/Ynsp)*(Sno2(n,j)/(Kno2nb+Sno2(n,j)))*(So2(n,j)/(Ko2nb+So2(n,j))))*Xnsp - ((u4/YcmxNo2)*(Sno2(n,j)/(Kno2nb+Sno2(n,j)))*(So2(n,j)/(Ko2nb+So2(n,j))))*Xcmx;
            temp1 = -U*(dSnh4dt(n+1,j)-dSnh4dt(n,j));
            temp2 = Dma*(dSnh4dt(n,j+1)-2*dSnh4dt(n,j)+dSnh4dt(n,j-1))/hz^2;
            temp3 = (u1*(INBM+1/Yaob)*(So2(n,j)/(Ko2aob+So2(n,j)))*(Snh4(n,j)/(Knh4aob+Snh4(n,j))))*Xaob;
            temp4 = (u4*(INBM+1/Ycmx)*(So2(n,j)/(Ko2cmx+So2(n,j)))*(Snh4(n,j)/(Knh4cmx+Snh4(n,j))))*Xcmx;
            dSnh4dt(n+1,j) = temp1 + temp2 - temp3 - temp4 -(u5*(INBM+1/Yamx)*(Snh4(n,j)/(Knh4amx+Snh4(n,j)))*(Sno2(n,j)/(Kno2amx+Sno2(n,j)))*(Ko2amx/(Ko2amx+So2(n,j))))*Xamx - (u2*INBM*(Sno2(n,j)/(Kno2nb+Sno2(n,j)))*(So2(n,j)/(Ko2nb+So2(n,j))))*Xnb -(u3*INBM*(Sno2(n,j)/(Kno2nsp+Sno2(n,j)))*(So2(n,j)/(Ko2nsp+So2(n,j))))*Xnsp - INBM*u6*nhan*((Sse(n,j)/(KSsehan+Sse(n,j)))*(Sno2(n,j)/(Kno2han+Sno2(n,j)))*(Ko2han/(Ko2han+So2(n,j))))*Xhan - u6*INBM*nhan*((Sse(n,j)/(KSsehan+Sse(n,j)))*(Sno3(n,j)/(Kno3han+Sno3(n,j)))*(Ko2han/(Ko2han+So2(n,j))))*Xhan - u6*INBM*((Sse(n,j)/(KSsehan+Sse(n,j)))*(So2(n,j)/(Ko2han+So2(n,j))))*Xhan + (INBM-(fI*INXI))*baob*naob*((Ko2aob/(Ko2aob+So2(n,j)))*(Sno2(n,j)+Sno3(n,j))/(Kno3aob+Sno2(n,j)+Sno3(n,j)))*Xaob + (INBM-(fI*INXI))*bnb*nnb*((Ko2nb/(Ko2nb+So2(n,j)))*(Sno2(n,j)+Sno3(n,j))/(Kno3nsp+Sno2(n,j)+Sno3(n,j)))*Xnb +(INBM-(fI*INXI))*bnsp*nnsp*((Ko2nsp/(Ko2nsp+So2(n,j)))*(Sno2(n,j)+Sno3(n,j))/(Kno3nsp+Sno2(n,j)+Sno3(n,j)))*Xnsp + (INBM-(fI*INXI))*(bcmx*ncmx*(Ko2cmx/(Ko2cmx+So2(n,j)))*(Sno2(n,j)+Sno3(n,j))/(Kno3cmx+Sno2(n,j)+Sno3(n,j)))*Xcmx +(INBM-(fI*INXI))*bamx*namx*((Ko2amx/(Ko2amx+So2(n,j)))*(Sno2(n,j)+Sno3(n,j))/(Kno3amx+Sno2(n,j)+Sno3(n,j)))*Xamx + (INBM-(fI*INXI))*bhan*nhan*((Ko2han/(Ko2han+So2(n,j)))*(Sno2(n,j)+Sno3(n,j))/(Kno3han+Sno2(n,j)+Sno3(n,j)))*Xhan +(INBM-(fI*INXI))*baob*(So2(n,j)/(Ko2aob+So2(n,j)))*Xaob + (INBM-(fI*INXI))*bnb*(So2(n,j)/(Ko2nb+So2(n,j)))*Xnb +(INBM-(fI*INXI))*bnsp*(So2(n,j)/(Ko2nsp+So2(n,j)))*Xnsp + (INBM-(fI*INXI))*bcmx*(So2(n,j)/(Ko2cmx+So2(n,j)))*Xcmx + (INBM-(fI*INXI))*bamx*(So2(n,j)/(Ko2amx+So2(n,j)))*Xamx + (INBM-(fI*INXI))*bhan*(So2(n,j)/(Ko2han+So2(n,j)))*Xhan;
            
            temp1 = -U*(dSo2dt(n+1,j)-dSo2dt(n,j));
            temp2 = Dma*(dSo2dt(n,j+1)-2*dSo2dt(n,j)+dSo2dt(n,j-1))/hz^2;
            temp3 = (KaL*(So2G-So2(n,j)))-(((1-Yhaer)/Yhaer)*u6*(Sse(n,j)/(KSsehan+Sse(n,j)))*(So2(n,j)/(Ko2han+So2(n,j))))*Xhan;
            temp4 = (((3.43-Yaob)/Yaob)*u1*(So2(n,j)/(Ko2aob+So2(n,j)))*(Snh4(n,j)/(Knh4aob+Snh4(n,j))))*Xaob;
            dSo2dt(n+1,j) = temp1 + temp2 + temp3 - temp4 - (((1.14-Ynb)/Ynb)*u2*(Sno2(n,j)/(Kno2nb+Sno2(n,j)))*(So2(n,j)/(Ko2nb+So2(n,j))))*Xnb - (((1.14-Ynsp)/Ynsp)*u3*(Sno2(n,j)/(Kno2nsp+Sno2(n,j)))*(So2(n,j)/(Ko2nsp+So2(n,j))))*Xnsp - (((4.57-Ycmx)/Ycmx)*u4*(So2(n,j)/(Ko2cmx+So2(n,j)))*(Snh4(n,j)/(Knh4cmx+Snh4(n,j))))*Xcmx -(((1.14-YcmxNo2)/YcmxNo2)*u4*(Sno2(n,j)/(Kno2cmx+Sno2(n,j)))*(So2(n,j)/(Ko2cmx+So2(n,j))))*Xcmx- (1-fI)*baob*(So2(n,j)/(Ko2aob+So2(n,j)))*Xaob - (1-fI)*bnb*(So2(n,j)/(Ko2nb+So2(n,j)))*Xnb -(1-fI)*bnsp*(So2(n,j)/(Ko2nsp+So2(n,j)))*Xnsp - (1-fI)*bcmx*(So2(n,j)/(Ko2cmx+So2(n,j)))*Xcmx - (1-fI)*bamx*(So2(n,j)/(Ko2amx+So2(n,j)))*Xamx - (1-fI)*bhan*(So2(n,j)/(Ko2han+So2(n,j)))*Xhan;
        end
    end
    fval(8)=dSsedt(n+1,j);
    fval(9)=dSno3dt(n+1,j);
    fval(10)=dSno2dt(n+1,j);
    fval(11)=dSnh4dt(n+1,j);
    fval(12)=dSo2dt(n+1,j);
    
    fval = fval(:);
    
    dbclear if naninf
end

MATLAB: How to get the missing peak points

I don't see anything missing: I see it plotting all of the points that your program defines for it to plot.

If you see something as being missing, you are going to have to describe for us what you want to be plotted and what it produces instead. And keep in mind that we don't have your input file, and you have (ahem) omitted posting copies of the output, so we cannot see what you are seeing.

Best Answer

Related Solutions

MATLAB: Warning: Rank deficient, rank = 2, tol = 7.141946e-13

MATLAB: How to get the missing peak points

Related Question