Clear Filters
Clear Filters

Index exceeds the number of array elements. Index must not exceed 1.

4 views (last 30 days)
I keep getting the above error 'Index exceeds the number of array elements. Index must not exceed 1.' when I run the below program. What could be the problem?
Thanks in advance...
tstart = 0;
tend = 180*24;
nt = 14400000;
xstart = 0;
xend = 0.012;
nx = 2000;
t = linspace(tstart,tend,nt);
x = linspace(xstart,xend,nx);
% Method-of-lines solution
x = x.';
% Enlarge x and xhalf by a ghostcell
x = [x;x(end)+(x(end)-x(end-1))];
a=arrayfun(@(i)oscic(x(i)),1:nx,'UniformOutput',false);
b=cell2mat(a);
y01 = 0.5; y02 = 0.5; y03 = 0.5;y04 = 0.5; y05 = 0.5; y06 = 0.5;
y07 = b(7,:).'; y08 = b(8,:).'; y09 = b(9,:).'; y10 = b(10,:).'; y11 = b(11,:).'; y12 = b(12,:).';
y13 = b(13,:).'; y14 = b(14,:).'; y15 = b(15,:).'; y16 = b(16,:).'; y17 = 1;
y1=[y01; y02;y03;y04; y05;y06;y07; y08;y09;y10; y11;y12;y13; y14;y15;y16; y17];
%% constants
T=301;
%aob
O1=68/(0.00831*293*(273+T)); u1=1.443*exp(O1*8); Ko2aob=0.594; Knh4aob=1.5*2.4; Yaob=0.15; Kno3aob=0.5;
%nb
O2=44/(0.00831*293*(273+T)); u2=0.01*0.8*0.48*exp(O2*8); Ko2nb=0.98*0.435; Kno2nb=0.5*1.375; Ynb=1.0*0.13;Kno3nb=0.5;
%nsp
O3=44/(0.00831*293*(273+T)); u3=0.01*0.69*exp(O3*8); Ko2nsp=0.435; Kno2nsp=1.5*0.76; Ynsp=1.0*0.14;Kno3nsp=0.5;
%amx
O5=71/(0.00831*293*(273+T)); u5=5*0.0664*exp(O5*8); Kno2amx=0.2*0.175; Kno3amx=0.5; Knh4amx=0.2*0.185; Yamx=0.159; Ko2amx=0.1;
%han
O6=0.069;u6=7.2*exp(8*0.069); KSsehan=0.5*2; Kno2han=0*0.5;Kno3han=0.5*0.32; YNo2han=2*0.49; Ko2han=0.2;
YNo3han=0.5*0.79; Yhaer=0.37;
%other constants
baob=0.0054*exp(O6*(T-293)); bnb=3*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; Xso=355; Sno2o=0.2; So2o=0;
VR=350;
sigma=0; L=0.0001;
DLSse=4.8*10^-5; DLSno3=1.8*10^-4; DLSno2=1.8*10^-4; DLSnh4=1.8*10^-4; DLSo2=1.8*10^-4;
AL=120400; LL=2.0*10^-5;
DSse=5.8*10^-5; DSno3=1.4*10^-5; DSno2=1.4*10^-4; DSnh4=1.5*10^-4; DSo2=2.2*10^-4;
Qo=10; Snh4o=5; KaL=80;
%% execution
options = odeset('RelTol',1e-20,'AbsTol',1e-10,'InitialStep',1e-10, 'NonNegative',1);
[tSol, ySol]=ode15s(@(t,y) fun(t,y,Qo,nx,u1,Ko2aob,Knh4aob,Yaob,Kno3aob,u2,Ko2nb,Kno2nb,Ynb,Kno3nb,u3,Ko2nsp,Kno2nsp,Ynsp,Kno3nsp,u5,Ko2amx,Knh4amx,Yamx,Kno3amx,Kno2amx,u6,KSsehan,Kno2han,Kno3han,YNo2han,Ko2han,YNo3han,Yhaer,Sseo,Sno3o,Xso,Sno2o,So2o,Snh4o,VR,sigma,DLSse,DLSno3,DLSno2,DLSnh4,DLSo2,AL,LL,DSse,DSno3,DSno2,DSnh4,DSo2,baob,bnb,bnsp,bamx,bhan,INBM,fI,namx,nhan,naob,nnb,nnsp,KX,KH,INXI,KaL),t,y1, options);
%% feed parameters
function dydt=fun(t,y,Qo,nx,u1,Ko2aob,Knh4aob,Yaob,Kno3aob,u2,Ko2nb,Kno2nb,Ynb,Kno3nb,u3,Ko2nsp,Kno2nsp,Ynsp,Kno3nsp,u5,Ko2amx,Knh4amx,Yamx,Kno3amx,Kno2amx,u6,KSsehan,Kno2han,Kno3han,YNo2han,Ko2han,YNo3han,Yhaer,Sseo,Sno3o,Xso,Sno2o,So2o,Snh4o,VR,sigma,DLSse,DLSno3,DLSno2,DLSnh4,DLSo2,AL,LL,DSse,DSno3,DSno2,DSnh4,DSo2,baob,bnb,bnsp,bamx,bhan,INBM,fI,namx,nhan,naob,nnb,nnsp,KX,KH,INXI,KaL)
%Define the variables
dydt = zeros(length(y),1);
SseBL=y(1,:);
Sno3BL=y(2,:);
Sno2BL=y(3,:);
Snh4BL=y(4,:);
So2BL=y(5,:);
Xs=y(6,:);
SseBF=y(7:nx+6,:);
Sno3BF=y(nx+7:2*nx+6,:);
Sno2BF=y(2*nx+7:3*nx+6,:);
Snh4BF=y(3*nx+7:4*nx+6,:);
So2BF=y(4*nx+7:5*nx+6,:);
faob=y(5*nx+7:6*nx+6,:);
fnb=y(6*nx+7:7*nx+6,:);
fnsp=y(7*nx+7:8*nx+6,:);
famx=y(8*nx+7:9*nx+6,:);
fhan=y(9*nx+7:10*nx+6,:);
L=y(10*nx+7,:);
muaob=u1.*(So2BF./(Ko2aob+So2BF)).*(Snh4BF./(Knh4aob+Snh4BF)).*faob - baob.*(So2BF./(Ko2aob+So2BF)).*faob - baob*naob.*((Ko2aob./(Ko2aob+So2BF)).*(Sno2BF+Sno3BF)./(Kno3aob+Sno2BF+Sno3BF)).*faob;
munb=u2.*(Sno2BF./(Kno2nb+Sno2BF)).*(So2BF/(Ko2nb+So2BF)).*fnb - bnb.*(So2BF./(Ko2nb+So2BF)).*fnb - bnb*nnb.*((Ko2nb./(Ko2nb+So2BF)).*(Sno2BF+Sno3BF)./(Kno3nb+Sno2BF+Sno3BF)).*fnb;
munsp=u3.*(Sno2BF./(Kno2nsp+Sno2BF)).*(So2BF/(Ko2nsp+So2BF)).*fnsp - bnsp.*(So2BF./(Ko2nsp+So2BF)).*fnsp - bnsp*nnsp.*((Ko2nsp./(Ko2nsp+So2BF)).*(Sno2BF+Sno3BF)./(Kno3nsp+Sno2BF+Sno3BF)).*fnsp;
muamx=u5.*((Snh4BF./(Knh4amx+Snh4BF)).*(Sno2BF./(Kno2amx+Sno2BF)).*(Ko2amx./(Ko2amx+So2BF))).*famx - bamx.*(So2BF./(Ko2amx+So2BF)).*famx - bamx*namx.*((Ko2amx/(Ko2amx+So2BF)).*(Sno2BF+Sno3BF)./(Kno3amx+Sno2BF+Sno3BF)).*famx;
muhan=(u6.*(SseBF./(KSsehan+SseBF)).*(Sno2BF./(Kno2han+Sno2BF)).*(Ko2han./(Ko2han+So2BF)) + u6.*(SseBF/(KSsehan+SseBF)).*(Sno3BF/(Kno3han+Sno3BF)).*(Ko2han/(Ko2han+So2BF))+u6.*((SseBF/(KSsehan+SseBF)).*(So2BF/(Ko2han+So2BF))).*fhan)- bhan.*(So2BF/(Ko2han+So2BF)).*fhan - bhan*nhan.*((Ko2han./(Ko2han+So2BF)).*(Sno2BF+Sno3BF)./(Kno3han+Sno2BF+Sno3BF)).*fhan;
muaver=faob.*muaob + fnb.*munb + fnsp.*munsp + famx.*muamx + fhan.*muhan;
UL=L.*muaver+sigma;
VL=VR-AL.*(L + LL);
%BC-1
SseBFmin=SseBF(1); Sno3BFmin=Sno3BF(1); Sno2BFmin=Sno2BF(1); Snh4BFmin=Snh4BF(1); So2BFmin=So2BF(1);
faobmin=exp((muaob(1)-muaver(1))*0.0003);fnbmin=exp((munb(1)-muaver(1))*0.0003);fnspmin=exp((munsp(1)-muaver(1))*0.0003);famxmin=exp((muamx(1)-muaver(1))*0.0003);fhanmin=exp((muhan(1)-muaver(1))*0.0003);
%BC-2
SseBFmax=SseBL; Sno3BFmax=Sno3BL; Sno2BFmax=Sno2BL; Snh4BFmax=Snh4BL;
So2BFmax=So2BL;
SseBF=[SseBFmin;SseBFmax];Sno3BF=[Sno3BFmin;Sno3BFmax];Sno2BF=[Sno2BFmin;Sno2BFmax];Snh4BF=[Snh4BFmin;Snh4BFmax];So2BF=[So2BFmin;So2BFmax];
faob=[faobmin;faob];fnb=[fnbmin;fnb];fnsp=[fnspmin;fnsp];famx=[famxmin;famx];fhan=[fhanmin;fhan];
%dSseBLdt=zeros(nx,1);dSno3BLdt=zeros(nx,1);dSno2BLdt=zeros(nx,1);dSnh4BLdt=zeros(nx,1);dSo2BLdt=zeros(nx,1);dXsdt=zeros(nx,1);
dSseBFdt=zeros(nx,1);dSno3BFdt=zeros(nx,1);dSno2BFdt=zeros(nx,1);
dSnh4BFdt=zeros(nx,1);dSo2BFdt=zeros(nx,1);dfaobdt=zeros(nx,1);dfnbdt=zeros(nx,1);dfnspdt=zeros(nx,1);dfamxdt=zeros(nx,1);dfhandt=zeros(nx,1);
%dLdt=zeros(nx,1);
%% PDEs (7-16) and ODEs (1-6&17)
for ix=2:nx
hz=nx-(nx-1);
z=0.012; zeta=z./L;
U=L.*muaver*zeta;
dSseBLdt=Qo*(Sseo-SseBL)/VL(nx)-AL*((DSse/LL)-UL(nx))*(SseBL-SseBF(nx));
dSno3BLdt=Qo.*(Sno3o-Sno3BL)./VL(nx)-AL.*((DSno3/LL)-UL(nx)).*(Sno3BL-Sno3BF(nx));
dSno2BLdt=Qo.*(Sno2o-Sno2BL)./VL(nx)-AL.*((DSno2/LL)-UL(nx)).*(Sno2BL-Sno2BF(nx));
dSnh4BLdt=Qo.*(Snh4o-Snh4BL)./VL(nx)-AL.*((DSnh4/LL)-UL(nx)).*(Snh4BL-Snh4BF(nx));
dSo2BLdt=Qo.*(So2o-So2BL)./VL(nx) + KaL.*(8-So2BL) - AL.*((DSo2)-UL(nx)).*(So2BL-So2BF(nx));
dXsdt=Qo.*(Xso-Xs)./VL(nx) - KH.*((Xs./fhan(nx))./(KX+(Xs/fhan(nx)))).*fhan(nx);
dSseBFdt(ix)=DSse.*(SseBF(ix+1)-2.*SseBF(ix)+(SseBF(ix-1)))./(hz^2.*L(ix)^2) +zeta.*UL.*(SseBF(ix+1)-SseBF(ix))./L - (u6*nhan.*(1/YNo2han).*(SseBF(ix)./(KSsehan+SseBF(ix))).*(Sno2BF(ix)./(Kno2han+Sno2BF(ix))).*(Ko2han./(Ko2han+So2BF(ix)))).*fhan(ix) -(u6*nhan.*(1/YNo3han).*(SseBF(ix)./(KSsehan+SseBF(ix))).*(Sno3BF(ix)./(Kno3han+Sno3BF(ix))).*(Ko2han./(Ko2han+So2BF(ix)))).*fhan(ix) - (u6/Yhaer).*(SseBF(ix)./(KSsehan+SseBF(ix))).*(So2BF(ix)./(Ko2han+So2BF(ix))).*fhan(ix);
dSno3BFdt(ix)=DSno3.*(Sno3BF(ix+1)-2.*Sno3BF(ix)+(Sno3BF(ix-1)))./(hz.*L)^2 +zeta.*UL.*(Sno3BF(ix+1)-Sno3BF(ix))./L - (u6*nhan*((1-YNo3han)/(2.86*YNo3han)).*(SseBF(ix)./(KSsehan+SseBF(ix))).*(Sno3BF(ix)./(Kno3han+Sno3BF(ix))).*(Ko2han./(Ko2han+So2BF(ix)))).*fhan(ix) + (u5/1.14).*(Snh4BF(ix)./(Knh4amx+Snh4BF(ix))).*(Sno2BF(ix)./(Kno2amx+Sno2BF(ix))).*(Ko2amx./(Ko2amx+So2BF(ix))).*famx(ix) + (u2/Ynb*(Sno2BF(ix)./(Kno2nb+Sno2BF(ix))).*(So2BF(ix)./(Ko2nb+So2BF(ix))).*fnb(ix)) + ((u3/Ynsp).*(Sno2BF(ix)./(Kno2nsp+Sno2BF(ix))).*(So2BF(ix)./(Ko2nsp+So2BF(ix))).*fnsp(ix)) - ((1-fI)/2.86)*baob*naob.*(Ko2aob./(Ko2aob+So2BF(ix))).*(Sno2BF(ix)+Sno3BF(ix))./(Kno3aob+Sno2BF(ix)+Sno3BF(ix)).*faob(ix) - ((1-fI)/2.86)*bnb*nnb.*(Ko2nb/(Ko2nb+So2BF(ix))).*(Sno2BF(ix)+Sno3BF(ix))./(Kno3nb+Sno2BF(ix)+Sno3BF(ix)).*fnb(ix) - ((1-fI)/2.86)*bnsp*nnsp.*(Ko2nsp/(Ko2nsp+So2BF(ix))).*(Sno2BF(ix)+Sno3BF(ix))./(Kno3nsp+Sno2BF(ix)+Sno3BF(ix)).*fnsp(ix) - ((1-fI)/2.86)*bamx*namx.*(Ko2amx/(Ko2amx+So2BF(ix))).*(Sno2BF(ix)+Sno3(ix))./(Kno3amx+Sno2BF(ix)+Sno3BF(ix)).*famx(ix)-((1-fI)/2.86)*bhan*nhan.*(Ko2han/(Ko2han+So2BF(ix))).*(Sno2BF(ix)+Sno3BF(ix))./(Kno3han+Sno2BF(ix)+Sno3BF(ix)).*fhan(ix);
dSno2BFdt(ix)=DSno2.*(Sno2BF(ix+1)-2.*Sno2BF(ix)+(Sno2BF(ix-1)))./(hz.*L)^2 +zeta.*UL.*(Sno2BF(ix+1)-Sno2BF(ix))./L - (((1-YNo2han)/(1.71*YNo2han))*u6*nhan.*(SseBF(ix)./(KSsehan+SseBF(ix))).*(Sno2BF(ix)./(Kno2han+Sno2BF(ix))).*(Ko2han/(Ko2han+So2BF(ix))).*fhan(ix)) - (((u5/Yamx)+(u5/1.14)).*(Snh4BF(ix)./(Knh4amx+Snh4BF(ix))).*(Sno2BF(ix)/(Kno2amx+Sno2BF(ix))).*(Ko2amx/(Ko2amx+So2BF(ix)))).*famx(ix) + ((u1/Yaob).*(So2BF(ix)/(Ko2aob+So2(ix))).*(Snh4BF(ix)./(Knh4aob+Snh4BF(ix)))).*faob(ix) - ((u2/Ynb).*(Sno2BF(ix)./(Kno2nb+Sno2BF(ix))).*(So2BF(ix)./(Ko2nb+So2BF(ix))).*fnb(ix)) - ((u3/Ynsp).*(Sno2BF(ix)./(Kno2nsp+Sno2BF(ix))).*(So2BF(ix)./(Ko2nsp+So2BF(ix)))).*fnsp(ix);
dSnh4BFdt(ix)=DSnh4.*(Snh4BF(ix+1)-2.*Snh4BF(ix)+(Snh4BF(ix-1)))./(hz.*L)^2 +zeta.*UL.*(Snh4BF(ix+1)-Snh4BF(ix))./L - (u1*(INBM+1/Yaob).*(So2BF(ix)./(Ko2aob+So2BF(ix))).*(Snh4BF(ix)./(Knh4aob+Snh4BF(ix)))).*faob(ix) -(u5*(INBM+1/Yamx).*(Snh4BF(ix)./(Knh4amx+Snh4BF(ix))).*(Sno2BF(ix)./(Kno2amx+Sno2BF(ix))).*(Ko2amx./(Ko2amx+So2BF(ix)))).*famx(ix) - (u2*INBM.*(Sno2BF(ix)./(Kno2nb+Sno2BF(ix))).*(So2BF(ix)./(Ko2nb+So2BF(ix)))).*(ix) -(u3*INBM.*(Sno2BF(ix)./(Kno2nsp+Sno2BF(ix))).*(So2BF(ix)./(Ko2nsp+So2BF(ix)))).*fnsp(ix) - INBM*u6*nhan.*((SseBF(ix)./(KSsehan+SseBF(ix))).*(Sno2BF(ix)./(Kno2han+Sno2BF(ix))).*(Ko2han/(Ko2han+So2BF(ix)))).*fhan(ix) - u6*INBM*nhan.*((SseBF(ix)./(KSsehan+SseBF(ix))).*(Sno3BF(ix)./(Kno3han+Sno3BF(ix))).*(Ko2han./(Ko2han+So2BF(ix)))).*fhan(ix) - u6*INBM.*((SseBF(ix)./(KSsehan+SseBF(ix))).*(So2BF(ix)./(Ko2han+So2BF(ix)))).*fhan(ix) + (INBM-(fI*INXI))*baob*naob.*((Ko2aob/(Ko2aob+So2BF(ix))).*(Sno2BF(ix)+Sno3BF(ix))./(Kno3aob+Sno2BF(ix)+Sno3BF(ix))).*faob(ix) + (INBM-(fI*INXI))*bnb*nnb.*((Ko2nb/(Ko2nb+So2BF(ix))).*(Sno2BF(ix)+Sno3BF(ix))./(Kno3nsp+Sno2BF(ix)+Sno3BF(ix))).*fnb(ix) +(INBM-(fI*INXI))*bnsp*nnsp.*((Ko2nsp./(Ko2nsp+So2BF(ix))).*(Sno2BF(ix)+Sno3BF(ix))./(Kno3nsp+Sno2BF(ix)+Sno3BF(ix))).*fnsp(ix) + (INBM-(fI*INXI))*bamx*namx.*((Ko2amx./(Ko2amx+So2BF(ix))).*(Sno2BF(ix)+Sno3BF(ix))./(Kno3amx+Sno2BF(ix)+Sno3BF(ix))).*famx(ix) + (INBM-(fI*INXI))*bhan*nhan.*((Ko2han/(Ko2han+So2BF(ix))).*(Sno2BF(ix)+Sno3BF(ix))./(Kno3han+Sno2BF(ix)+Sno3BF(ix))).*fhan(ix) +(INBM-(fI*INXI))*baob.*(So2BF(ix)./(Ko2aob+So2BF(ix))).*faob(ix) + (INBM-(fI*INXI))*bnb.*(So2BF(ix)./(Ko2nb+So2BF(ix))).*fnb(ix) +(INBM-(fI*INXI))*bnsp.*(So2BF(ix)./(Ko2nsp+So2BF(ix))).*fnsp(ix) + (INBM-(fI*INXI))*bamx.*(So2BF(ix)./(Ko2amx+So2BF(ix))).*famx(ix) + (INBM-(fI*INXI))*bhan.*(So2BF(ix)./(Ko2han+So2BF(ix))).*fhan(ix);
dSo2BFdt(ix)=DSo2.*(So2BF(ix+1)-2.*So2BF(ix)+(So2BF(ix-1)))./(hz.*L)^2 +zeta.*UL.*(So2BF(ix+1)-So2BF(ix))./L - -(((1-Yhaer)/Yhaer)*u6.*(SseBF(ix)./(KSsehan+SseBF(ix))).*(So2BF(ix)./(Ko2han+So2BF(ix)))).*fhan(ix) -(((3.43-Yaob)/Yaob)*u1.*(So2BF(ix)./(Ko2aob+So2BF(ix))).*(Snh4BF(ix)./(Knh4aob+Snh4BF(ix)))).*faob(ix) - (((1.14-Ynb)/Ynb)*u2.*(Sno2BF(ix)/(Kno2nb+Sno2BF(ix))).*(So2BF(ix)./(Ko2nb+So2BF(ix)))).*fnb(ix) - (((1.14-Ynsp)/Ynsp)*u3.*(Sno2BF(ix)/(Kno2nsp+Sno2BF(ix))).*(So2BF(ix)./(Ko2nsp+So2BF(ix)))).*fnsp(ix) - (1-fI)*baob.*(So2BF(ix)./(Ko2aob+So2BF(ix))).*faob(ix) - (1-fI)*bnb.*(So2BF(ix)./(Ko2nb+So2BF(ix))).*fnb -(1-fI)*bnsp.*(So2BF(ix)./(Ko2nsp+So2BF(ix))).*fnsp(ix) - (1-fI)*bamx.*(So2BF(ix)./(Ko2amx+So2BF(ix))).*famx(ix) - (1-fI)*bhan.*(So2BF(ix)./(Ko2han+So2BF(ix))).*fhan(ix);
dfaobdt(ix)=(muaob-muaver).*faob(ix) - (U-zeta.*UL).*(faob(ix+1)-faob(ix))./hz;
dfnbdt(ix)=(munb-muaver).*fnb(ix) - (U-zeta.*UL).*(fnb(ix+1)-fnb(ix))./hz;
dfnspdt(ix)=(munsp-muaver).*fnsp(ix) - (U-zeta.*UL).*(fnsp(ix+1)-fnsp(ix))./hz;
dfamxdt(ix)=(muamx-muaver).*famx(ix) - (U-zeta.*UL).*(famx(ix+1)-famx(ix))./hz;
dfhandt(ix)=(muhan-muaver).*fhan(ix) - (U-zeta.*UL).*(fhan(ix+1)-fahan(ix))./hz;
dLdt(ix)=L(ix).*muaver.*zeta+sigma;
end
dydt=[dSseBLdt;dSno3BLdt;dSno2BLdt;dSnh4BLdt;dSo2BLdt;dXsdt;dSseBFdt;dSno3BFdt;dSno2BFdt;dSnh4BFdt;dSo2BFdt;dfaobdt;dfnbdt;dfnspdt;dfamxdt;dfhandt;dLdt];
end
%% initial conditions
function u0 = oscic(x)
u0 = [1;1;1;1;1;1;1;1;1;1;1;0.4;0.02;0.02;0.3;0.26;0.001];
end
  52 Comments
Walter Roberson
Walter Roberson on 4 Sep 2023
Index in position 1 exceeds array bounds. Index must not exceed 106.
Error in MBBR_aerobic (line 49)
Sno3BF=y(nx+6:2*nx+5,:);
y is length 106, nx is 100 so nx+6:2*nx+5 would be 106:205 which is mostly past the end of y
KIPROTICH KOSGEY
KIPROTICH KOSGEY on 5 Sep 2023
Edited: Walter Roberson on 5 Sep 2023
Dear Walter and Torsten,
Thank you so much for the useful comments.
Including odeset('NonNegative',1:6+10*nx) has prevented the solutions from going to negative but the computation keeps failing. Removing this conditions is not helping either as the computation still fail before it reaches tend.
Either way I am not winning.
This is the script and the main files are atatched:
close all; clear; clc;
dbstop if error
t1=[0 180*24];
xstart = 0;
xend = 0.012;
nx=10;
x = linspace(xstart,xend,nx);
% Method-of-lines solution
x = x.';
% Enlarge x and xhalf by a ghostcell
x = [x;x(end)+(x(end)-x(end-1))];
a=arrayfun(@(i)oscic(x(i)),1:nx,'UniformOutput',false);
b=cell2mat(a);
y0 = 0.5*ones(1,5); y16 = 1*10^-6;
y06 = b(6,:).'; y07 = b(7,:).'; y08 = b(8,:).'; y09 = b(9,:).'; y10 = b(10,:).'; y11 = b(11,:).';
y12 = b(12,:).'; y13 = b(13,:).'; y14 = b(14,:).'; y15 = b(15,:).';
y1=[y0(1); y0(2);y0(3);y0(4); y0(5);y06;y07; y08;y09;y10; y11;y12;y13; y14;y15;y16];
%% execution
options1 = odeset('RelTol',1e-14,'AbsTol',1e-16,'events',@eventsfcn_aerobic,'NonNegative',1:6+10*nx);
options2 = odeset('RelTol',1e-14,'AbsTol',1e-16,'events',@eventsfcn_anoxic,'NonNegative',1:6+10*nx);
%main
[tSol, ySol]=ode15s(@(t,y) MBBR_aerobic(t,y),t1,y1,options1);
tplot=tSol;
yplot=ySol;
for i=1:10
y2=[yplot(end,1);yplot(end,2);yplot(end,3);yplot(end,4);yplot(end,5);yplot(end,6:15)';yplot(end,16:25)';yplot(end,26:35)';yplot(end,36:45)';...
yplot(end,46:55)';yplot(end,56:65)';yplot(end,66:75)';yplot(end,76:85)';yplot(end,86:95)';yplot(end,96:105)';yplot(end,106)];
[tSol,ySol]=ode15s(@MBBR_anoxic,[tplot(end) 24*180],y2,options2);
tplot=[tplot;tSol];
yplot=[yplot;ySol];
y3=[yplot(end,1);yplot(end,2);yplot(end,3);yplot(end,4);yplot(end,5);yplot(end,6:15)';yplot(end,16:25)';yplot(end,26:35)';yplot(end,36:45)';...
yplot(end,46:55)';yplot(end,56:65)';yplot(end,66:75)';yplot(end,76:85)';yplot(end,86:95)';yplot(end,96:105)';yplot(end,106)];
[tSol, ySol]=ode15s(@(t,y) MBBR_aerobic(t,y),[tplot(end) 24*180],y3,options1);
tplot=[tplot;tSol];
yplot=[yplot;ySol];
end
%% initial conditions
function u0 = oscic(x)
u0 = [10;10;10;10;0.5;10;10;10;10;10;0.5;0.4;0.05;0.05;0.3;0.2;0.001];
end
%% events aerobic phase
function [value,isterminal,direction]=eventsfcn_aerobic(t,y)
z=y(5);
value=z - 1.5;
isterminal =1; % stop at 0
direction = 1; % increasing
end
function [value,isterminal,direction]=eventsfcn_anoxic(t,y)
z=y(5);
value=0.5 - z;
isterminal = 1 ; % stop at 0
direction = 1; % increasing
end

Sign in to comment.

Answers (0)

Tags

Products


Release

R2023a

Community Treasure Hunt

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

Start Hunting!