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

Hi,

Kindly, can someone assist to figure out what mistake i made on the attached code that is generating the error below?

'Warning: Rank deficient, rank = 2, tol =

7.141946e-13'

Thanks in advance

Answer 1

[Walter Roberson]

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

[KIPROTICH KOSGEY]

Thank you @Walter.

Basically what i am trying to do is to simulate the performance of a column biorector for a period of 535 days based on the input parameters provided. functions 1-6 represent the populations of different bacteria while functions 7-12 represent either substrate or products of metabolic activities. The equations representing these calculations are either ODEs or PDEs.

I have never coded for coupled ODE and PDE systems hence these difficulties and the obvious mistakes i am making.

1. if i can not use ODE sovers, what should i use?
2. i am still stack at n=94, please can you advice on what i should do?

Thanks in advace

Kosgey

[Walter Roberson]

I do not recommend coding your own pde solver.

In theory, you can use a pde solver with ode: they just do not happen to have coupling terms (or you could say that the coupling terms are 0)

I have not looked at the pde solvers much, but it would not surprise me if they required function continuity as well, so likely also requiring a split at each discontinuity... but I do not know for sure at the moment.

[KIPROTICH KOSGEY]

Ok thank you.

I have figured out how to code the reactor using a system of ODEs. However, i am still having a diffulty on how to include events in oder to take care of the discontinuities....

I attach the files.

Thanks again in advance.

[Walter Roberson]

Do not use events for this.

y0=[0.5*18.16;0.5*0.52;0.5*25.64;0.5*3.6;0.5*172.6;0.5*179.48;0.0005;2;0.0005;55;50;0.5];

tspans = [0, 46, 65, 67, 74, 80, 86, 87, 102, 108, 121, 134, 137, 151, 159, 176, 177, 180, 197, 208, 212, 214, 221, 225, 228, 237, 245, 250, 250.1, 254, 256, 264, 278, 288, 295, 302, 311, 313, 322, 327, 355, 425, 434, 475, 491, 535];
nspans = length(tspans) - 1;
tout = cell(nspans,1);
yout = cell(nspans,1);

hold off

for idx = 1 : nspans
    % Solve until the first terminal event.
    [t,y] = ode15s(@(t,y) MBBRFun5(t,y), linspace(tspans(idx),tspans(idx+1),25), y0);

    % Accumulate output.  This could be passed out as output arguments.
    
    if idx == 1
        tout{idx} = t;
        yout{idx} = y;
    else
        tout{idx} = t(2:end);
        yout{idx} = y(2:end,:);
    end
    for col = 1 : 12
        ax = subplot(3,4,col);
        plot(ax, t, y(:,col));
        title(ax, sprintf('y(:,%d)', col) );
        hold(ax, 'on');
    end
    drawnow
    
    %set new initial conditions
    y0 = y(end,:);
end

all_t = vertcat(tout{:});
all_y = vertcat(yout{:});

The slowest part is the plotting. If you were to postpone that to the end, then it would be faster.

[KIPROTICH KOSGEY]

Thank you soo much Walter for your efforts and time.

I realise that this codes are giving exactly the same results as i had before stopping at the discontinuities. However, my codes were difficult to execute because of high stiffness so i kept altering the tolerances. My previous script is attached.

As for the coupled PDEs and ODEs, i bet i will have to dig a lot deep in future.

I cannot thank you enough!