You are now following this question
- You will see updates in your followed content feed.
- You may receive emails, depending on your communication preferences.
Index exceeds the number of array elements. Index must not exceed 1.
11 views (last 30 days)
Show older comments
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
Dyuman Joshi
on 10 Aug 2023
1 - Please do not use dynamically named variables, it is not a good programming practice.
y01 = 0.5; y02 = 0.5; y03 = 0.5;y04 = 0.5; y05 = 0.5; y06 = 0.5;
As Steve Vai says (or rather plays), For the love of god, do not do this, simply define an array
y0 = 0.5*ones(1,6)
and use indexing.
2 - VL is a scalar i.e. it has only elements, and you are trying to access its nxth element, which is not possible. That's why you get the error.
% vvvvvv
dSseBLdt=Qo*(Sseo-SseBL)/VL(nx)-AL*((DSse/LL)-UL(nx))*(SseBL-SseBF(nx));
KIPROTICH KOSGEY
on 11 Aug 2023
I have edited but still the error persists.
New codes are:
y0 = 0.5*ones(1,6); y017=10^-6*ones(1,1);
y01=y0(1);y02=y0(2);y03=y0(3);y04=y0(4);y05=y0(5);y06=y0(6);y17=y017(1);
dSseBLdt=Qo*(Sseo-SseBL)/VL-AL*((DSse/LL)-UL)*(SseBL-SseBF(nx));
dSno3BLdt=Qo.*(Sno3o-Sno3BL)./VL-AL.*((DSno3/LL)-UL).*(Sno3BL-Sno3BF(nx)); dSno2BLdt=Qo.*(Sno2o-Sno2BL)./VL-AL.*((DSno2/LL)-UL).*(Sno2BL-Sno2BF(nx));
dSnh4BLdt=Qo.*(Snh4o-Snh4BL)./VL-AL.*((DSnh4/LL)-UL).*(Snh4BL-Snh4BF(nx));
dSo2BLdt=Qo.*(So2o-So2BL)./VL + KaL.*(8-So2BL) - AL.*((DSo2)-UL).*(So2BL-So2BF(nx));
dXsdt=Qo.*(Xso-Xs)./VL - KH.*((Xs./fhan(nx))./(KX+(Xs/fhan(nx)))).*fhan(nx);
Dyuman Joshi
on 11 Aug 2023
It's the same thing as above, Ssebf has only 2 elements, and you are trying to access the nxth element, which is not possible.
You should check the size of all the variables before running your code.
And, this is still not good -
y01=y0(1);y02=y0(2);y03=y0(3);y04=y0(4);y05=y0(5);y06=y0(6);y17=y017(1);
Do this -
%Define y0 as a column vector
y0 = 0.5*ones(6,1);
y17 = 10^-6*ones(1,1);
%USE INDEXING!!!
y1=[y0; b(7:16,:).'; y17];
KIPROTICH KOSGEY
on 11 Aug 2023
I am struggling to code the BC for pdes.
%BC-1
SseBF(1)=SseBF(2); Sno3BF(1)=Sno3BF(2); Sno2BF(1)=Sno2BF(2); Snh4BF(1)=Snh4BF(2); So2BF(1)=So2BF(2);
faob(1)=exp((muaob(1)-muaver(1))*0.0003);fnb(1)=exp((munb(1)-muaver(1))*0.0003);fnsp(1)=exp((munsp(1)-muaver(1))*0.0003);famx(1)=exp((muamx(1)-muaver(1))*0.0003);fhan(1)=exp((muhan(1)-muaver(1))*0.0003);
%BC-2
SseBF(nx)=SseBL; Sno3BF(nx)=Sno3BL; Sno2BF(nx)=Sno2BL; Snh4BF(nx)=Snh4BL;
So2BF(nx)=So2BL;
KIPROTICH KOSGEY
on 11 Aug 2023
Thank you Torsten, you have been so helpful in the course of my project. I didn't quite understand what you mean by 'set the boundary conditions separately.'
I have tried introducing a for loop for BCs but it aint helping:
for ix=1:nx
%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;
end
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];
Torsten
on 11 Aug 2023
The output of your function "fun" are time derivatives.
Thus if you want to define boundary conditions and if the function values in the boundary points are solution variables, you have to transfer their time derivatives to "ode15s".
If you want to prescribe a value for a variable and the value remains the same over time, you simply have to set d(variable)/dt = 0 for ix = 1 or ix = nx.
If you want to prescribe a value for a variable and the value changes over time (say like a function boundary(t)), you have to set d(variable)/dt = d(boundary)/dt for ix = 1 or ix = nx.
If you want to set d(variable)/dx = something as boundary condition, look at the test code we developed here:
KIPROTICH KOSGEY
on 12 Aug 2023
Edited: Walter Roberson
on 16 Aug 2023
Thank you so much.
I have nested another function to solve for the ODE for the left boundary conditions. However, I am getting an error regarding the indexing within the BC code:
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,:);
dydt=@fval;
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(nx-(nx-1)); Sno3BFmin=Sno3BF(nx-(nx-1)); Sno2BFmin=Sno2BF(nx-(nx-1)); Snh4BFmin=Snh4BF(nx-(nx-1));
So2BFmin=So2BF(nx-(nx-1));
%BC-2
SseBFmax=SseBF(nx-1)+L.*DLSse.*(SseBL(nx)-SseBF(nx))./LL; Sno3BFmax=Sno3BF(nx-1)+L.*DLSno3.*(Sno3BL(nx)-Sno3BF(nx))./LL;
Sno2BFmax=Sno2BF(nx-1)+L.*DLSno2.*(Sno2BL(nx)-Sno2BF(nx))/LL;
Snh4BFmax=Snh4BF(nx-1)+L.*DLSnh4.*(Snh4BL(nx)-Snh4BF(nx))/LL;
So2BFmax=So2BF(nx-1)+L.*DLSo2.*(So2BL(nx)-So2BF(nx))./LL;
faobBC(1)=0.5; fnbBC(1)=0.5;fnspBC(1)=0.5;famxBC(1)=0.5; fhanBC(1)=0.5;
SseBF=[SseBFmin;SseBFmax];Sno3BF=[Sno3BFmin;Sno3BFmax];Sno2BF=[Sno2BFmin;Sno2BFmax];Snh4BF=[Snh4BFmin;Snh4BFmax];So2BF=[So2BFmin;So2BFmax];
faob=[faobBC;faob];fnb=[fnbBC;fnb];fnsp=[fnspBC;fnsp];famx=[famxBC;famx];fhan=[fhanBC;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-1
hz=nx-(nx-1);
z=0.012; zeta=z./L;
U=L.*muaver*zeta;
dSseBLdt=Qo.*(Sseo-SseBL)/VL-AL*((DSse/LL)-UL)*(SseBL-SseBF(nx));
dSno3BLdt=Qo.*(Sno3o-Sno3BL)./VL-AL.*((DSno3/LL)-UL).*(Sno3BL-Sno3BF(nx));
dSno2BLdt=Qo.*(Sno2o-Sno2BL)./VL-AL.*((DSno2/LL)-UL).*(Sno2BL-Sno2BF(nx));
dSnh4BLdt=Qo.*(Snh4o-Snh4BL)./VL-AL.*((DSnh4/LL)-UL).*(Snh4BL-Snh4BF(nx));
dSo2BLdt=Qo.*(So2o-So2BL)./VL + KaL.*(8-So2BL) - AL.*((DSo2)-UL).*(So2BL-So2BF(nx));
dXsdt=Qo.*(Xso-Xs)./VL - 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(ix) - (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(ix)-zeta.*UL).*(faob(ix+1)-faob(ix))./hz;
dfnbdt(ix)=(munb-muaver).*fnb(ix) - (U(ix)-zeta.*UL).*(fnb(ix+1)-fnb(ix))./hz;
dfnspdt(ix)=(munsp-muaver).*fnsp(ix) - (U(ix)-zeta.*UL).*(fnsp(ix+1)-fnsp(ix))./hz;
dfamxdt(ix)=(muamx-muaver).*famx(ix) - (U(ix)-zeta.*UL).*(famx(ix+1)-famx(ix))./hz;
dfhandt(ix)=(muhan-muaver).*fhan(ix) - (U(ix)-zeta.*UL).*(fhan(ix+1)-fahan(ix))./hz;
dLdt=L(ix).*muaver.*zeta+sigma;
end
function fval=BC(t,yBC,muaob,munb,munsp,muamx,muhan,muaver)
fval=zeros(length(yBC),1);
faobBC=yBC(1,:);
fnbBC=yBC(2,:);
fnspBC=yBC(3,:);
famxBC=yBC(4,:);
fhanBC=yBC(5,:);
fval=[(muaob-muaver)*faobBC
(munb-muaver)*fnbBC
(munsp-muaver)*fnspBC
(muamx-muaver)*famxBC
(muhan-muaver)*fhanBC];
fval=[faobBC;fnbBC;fnspBC;famxBC;fhanBC];
end
dydt=[dSseBLdt;dSno3BLdt;dSno2BLdt;dSnh4BLdt;dSo2BLdt;dXsdt;dSseBFdt;dSno3BFdt;dSno2BFdt;dSnh4BFdt;dSo2BFdt;dfaobdt;dfnbdt;dfnspdt;dfamxdt;dfhandt;dLdt];
end
Torsten
on 12 Aug 2023
I cannot run your code since all constants and the call to ode15s is missing.
KIPROTICH KOSGEY
on 15 Aug 2023
Edited: Walter Roberson
on 16 Aug 2023
sorry for late response but my computer had crashed. Below is the code:
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 = b(1,:).'; y02=b(2,:).'; y03=b(3,:).';y04=b(4,:).'; y05 = b(5,:).'; 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,:).'; y16 = b(16,:).'; y17 = b(17,:).';
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);
%main
[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);
Index exceeds the number of array elements. Index must not exceed 2.
Error in solution>fun (line 100)
dSseBLdt=Qo.*(Sseo-SseBL)/VL-AL*((DSse/LL)-UL)*(SseBL-SseBF(nx));
Error in solution>@(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) (line 48)
[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);
Error in odearguments (line 92)
f0 = ode(t0,y0,args{:}); % ODE15I sets args{1} to yp0.
Error in ode15s (line 153)
odearguments(odeIsFuncHandle, odeTreatAsMFile, solver_name, ode, tspan, y0, options, varargin);
%% 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,:);
dydt=@fval;
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);
%BC-2
SseBFmax=SseBF(nx)+2*(nx-(nx-1))*L.*DLSse.*(SseBL-SseBF(nx))./LL;
Sno3BFmax=Sno3BF(nx)+2*(nx-(nx-1))*L.*DLSno3.*(Sno3BL-Sno3BF(nx))./LL;
Sno2BFmax=Sno2BF(nx)+2*(nx-(nx-1))*L.*DLSno2.*(Sno2BL-Sno2BF(nx))/LL;
Snh4BFmax=Snh4BF(nx)+2*(nx-(nx-1))*L.*DLSnh4.*(Snh4BL-Snh4BF(nx))/LL;
So2BFmax=So2BF(nx)+2*(nx-(nx-1))*L.*DLSo2.*(So2BL-So2BF(nx))./LL;
faobBC(1)=0.5; fnbBC(1)=0.5;fnspBC(1)=0.5;famxBC(1)=0.5; fhanBC(1)=0.5;
SseBF=[SseBFmin;SseBFmax];Sno3BF=[Sno3BFmin;Sno3BFmax];Sno2BF=[Sno2BFmin;Sno2BFmax];Snh4BF=[Snh4BFmin;Snh4BFmax];So2BF=[So2BFmin;So2BFmax];
faob=[faobBC;faob];fnb=[fnbBC;fnb];fnsp=[fnspBC;fnsp];famx=[famxBC;famx];fhan=[fhanBC;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-AL*((DSse/LL)-UL)*(SseBL-SseBF(nx));
dSno3BLdt=Qo.*(Sno3o-Sno3BL)./VL-AL.*((DSno3/LL)-UL).*(Sno3BL-Sno3BF(nx));
dSno2BLdt=Qo.*(Sno2o-Sno2BL)./VL-AL.*((DSno2/LL)-UL).*(Sno2BL-Sno2BF(nx));
dSnh4BLdt=Qo.*(Snh4o-Snh4BL)./VL-AL.*((DSnh4/LL)-UL).*(Snh4BL-Snh4BF(nx));
dSo2BLdt=Qo.*(So2o-So2BL)./VL + KaL.*(8-So2BL) - AL.*((DSo2)-UL).*(So2BL-So2BF(nx));
dXsdt=Qo.*(Xso-Xs)./VL - 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(ix).*UL(ix).*(SseBF(ix+1)-SseBF(ix))./L(ix) - (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^2.*L(ix)^2) +zeta(ix).*UL(ix).*(Sno3BF(ix+1)-Sno3BF(ix))./L(ix) - (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)+Sno3BF(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^2.*L(ix)^2) +zeta(ix).*UL(ix).*(Sno2BF(ix+1)-Sno2BF(ix))./L(ix) - (((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+So2BF(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^2.*L(ix)^2) +zeta(ix).*UL(ix).*(Snh4BF(ix+1)-Snh4BF(ix))./L(ix) - (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^2.*L(ix)^2) +zeta(ix).*UL(ix).*(So2BF(ix+1)-So2BF(ix))./L(ix) -(((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(ix) -(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(ix)-muaver(ix)).*faob(ix) - (U(ix)-zeta(ix).*UL(ix)).*(faob(ix+1)-faob(ix))./hz;
dfnbdt(ix)=(munb(ix)-muaver(ix)).*fnb(ix) - (U(ix)-zeta(ix).*UL(ix)).*(fnb(ix+1)-fnb(ix))./hz;
dfnspdt(ix)=(munsp(ix)-muaver(ix)).*fnsp(ix) - (U(ix)-zeta(ix).*UL(ix)).*(fnsp(ix+1)-fnsp(ix))./hz;
dfamxdt(ix)=(muamx(ix)-muaver(ix)).*famx(ix) - (U(ix)-zeta(ix).*UL(ix)).*(famx(ix+1)-famx(ix))./hz;
dfhandt(ix)=(muhan(ix)-muaver(ix)).*fhan(ix) - (U(ix)-zeta(ix).*UL(ix)).*(fhan(ix+1)-fhan(ix))./hz;
dLdt=L(ix).*muaver(ix).*zeta(ix)+sigma;
end
function fval=BC(t,yBC,nx,muaob,munb,munsp,muamx,muhan,muaver)
fval=zeros(length(yBC),1);
faobBC=yBC(1:nx,:);
fnbBC=yBC(nx+1:2*nx,:);
fnspBC=yBC(2*nx+1:3*nx,:);
famxBC=yBC(3*nx+1:4*nx,:);
fhanBC=yBC(4*nx+1:5*nx,:);
for i=1:nx
fval=[(muaob-muaver)*faobBC
(munb-muaver)*fnbBC
(munsp-muaver)*fnspBC
(muamx-muaver)*famxBC
(muhan-muaver)*fhanBC];
end
fval=[faobBC;fnbBC;fnspBC;famxBC;fhanBC];
end
dydt=[dSseBLdt;dSno3BLdt;dSno2BLdt;dSnh4BLdt;dSo2BLdt;dXsdt;dSseBFdt;dSno3BFdt;dSno2BFdt;dSnh4BFdt;dSo2BFdt;dfaobdt;dfnbdt;dfnspdt;dfamxdt;dfhandt;dLdt;faobBC;fnbBC;fnspBC;famxBC;fhanBC];
end
%% initial conditions
function u0 = oscic(x)
u0 = [0.5;0.5;0.5;0.5;0.5;0.5;1;1;1;1;1;0.4;0.02;0.02;0.3;0.26;0.001];
end
KIPROTICH KOSGEY
on 16 Aug 2023
Dear Torsten,
could you please advise what the problem is in the above code, it is still throwing the index error.
Walter Roberson
on 16 Aug 2023
SseBF=y(7:nx+6,:);
You create SSeBF as a vector (y only has one column)
SseBF=[SseBFmin;SseBFmax];
and there you overwrite the column vector with a different column vector of length 2.
KIPROTICH KOSGEY
on 16 Aug 2023
Thank you Walter!
I have corrected this as :
SseBF=[SseBF;SseBFmax];Sno3BF=[Sno3BF;Sno3BFmax];Sno2BF=[Sno2BF;Sno2BFmax];Snh4BF=[Snh4BF;Snh4BFmax];So2BF=[So2BF;So2BFmax];
However, further down I am getting the error:
Error using vertcat
Dimensions of arrays being concatenated are not consistent.
Error in MBBR>fun (line 155)
Please advise.
KIPROTICH KOSGEY
on 16 Aug 2023
Edited: Walter Roberson
on 16 Aug 2023
Dear Walter and Torsten,
May I thank you for pointing out errors in my below set of codes. Please I have one last error that I can't figure ot why its coming up yet everything seems fine in my coding: Error using vertcat Dimensions of arrays being concatenated are not consistent. Error in MBBR>fun (line 154)
dydt=[dSseBLdt;dSno3BLdt;......
Edited program:
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);
y0 = 0.5*ones(1,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,:).'; y16 = b(16,:).'; y17 = 1*10^-6;
y1=[y0(1); y0(2);y0(3);y0(4); y0(5);y0(6);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);
%main
[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);
Error using vertcat
Dimensions of arrays being concatenated are not consistent.
Dimensions of arrays being concatenated are not consistent.
Error in solution>fun (line 134)
dydt=[dSseBLdt;dSno3BLdt;dSno2BLdt;dSnh4BLdt;dSo2BLdt;dXsdt;dSseBFdt;dSno3BFdt;dSno2BFdt;dSnh4BFdt;dSo2BFdt;dfaobdt;dfnbdt;dfnspdt;dfamxdt;dfhandt;dLdt];
Error in solution>@(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) (line 48)
[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);
Error in odearguments (line 92)
f0 = ode(t0,y0,args{:}); % ODE15I sets args{1} to yp0.
Error in ode15s (line 153)
odearguments(odeIsFuncHandle, odeTreatAsMFile, solver_name, ode, tspan, y0, options, varargin);
%% 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,:);
dydt=@fval;
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);
%BC-2
SseBFmax=SseBF(nx)+2*(nx-(nx-1))*L.*DLSse.*(SseBL-SseBF(nx))./LL;
Sno3BFmax=Sno3BF(nx)+2*(nx-(nx-1))*L.*DLSno3.*(Sno3BL-Sno3BF(nx))./LL;
Sno2BFmax=Sno2BF(nx)+2*(nx-(nx-1))*L.*DLSno2.*(Sno2BL-Sno2BF(nx))/LL;
Snh4BFmax=Snh4BF(nx)+2*(nx-(nx-1))*L.*DLSnh4.*(Snh4BL-Snh4BF(nx))/LL;
So2BFmax=So2BF(nx)+2*(nx-(nx-1))*L.*DLSo2.*(So2BL-So2BF(nx))./LL;
faobBC(1)=0.5; fnbBC(1)=0.5;fnspBC(1)=0.5;famxBC(1)=0.5; fhanBC(1)=0.5;
SseBF=[SseBF;SseBFmax];Sno3BF=[Sno3BF;Sno3BFmax];Sno2BF=[Sno2BF;Sno2BFmax];Snh4BF=[Snh4BF;Snh4BFmax];So2BF=[So2BF;So2BFmax];
faob=[faobBC;faob];fnb=[fnbBC;fnb];fnsp=[fnspBC;fnsp];famx=[famxBC;famx];fhan=[fhanBC;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-1
hz=nx-(nx-1);
z=0.012; zeta=z./L;
U=L.*muaver*zeta;
dSseBLdt=Qo.*(Sseo-SseBL)/VL-AL*((DSse/LL)-UL)*(SseBL-SseBF(nx));
dSno3BLdt=Qo.*(Sno3o-Sno3BL)./VL-AL.*((DSno3/LL)-UL).*(Sno3BL-Sno3BF(nx));
dSno2BLdt=Qo.*(Sno2o-Sno2BL)./VL-AL.*((DSno2/LL)-UL).*(Sno2BL-Sno2BF(nx));
dSnh4BLdt=Qo.*(Snh4o-Snh4BL)./VL-AL.*((DSnh4/LL)-UL).*(Snh4BL-Snh4BF(nx));
dSo2BLdt=Qo.*(So2o-So2BL)./VL + KaL.*(8-So2BL) - AL.*((DSo2)-UL).*(So2BL-So2BF(nx));
dXsdt=Qo.*(Xso-Xs)./VL - 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^2) +zeta.*UL(ix).*(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^2.*L^2) +zeta.*UL(ix).*(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)+Sno3BF(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^2.*L^2) +zeta.*UL(ix).*(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+So2BF(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^2.*L^2) +zeta.*UL(ix).*(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^2.*L^2) +zeta.*UL(ix).*(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(ix) -(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(ix)-muaver(ix)).*faob(ix) - (U(ix)-zeta.*UL(ix)).*(faob(ix+1)-faob(ix))./hz;
dfnbdt(ix)=(munb(ix)-muaver(ix)).*fnb(ix) - (U(ix)-zeta.*UL(ix)).*(fnb(ix+1)-fnb(ix))./hz;
dfnspdt(ix)=(munsp(ix)-muaver(ix)).*fnsp(ix) - (U(ix)-zeta.*UL(ix)).*(fnsp(ix+1)-fnsp(ix))./hz;
dfamxdt(ix)=(muamx(ix)-muaver(ix)).*famx(ix) - (U(ix)-zeta.*UL(ix)).*(famx(ix+1)-famx(ix))./hz;
dfhandt(ix)=(muhan(ix)-muaver(ix)).*fhan(ix) - (U(ix)-zeta.*UL(ix)).*(fhan(ix+1)-fhan(ix))./hz;
dLdt=L.*muaver(ix).*zeta+sigma;
end
function fval=BC(t,yBC,nx,muaob,munb,munsp,muamx,muhan,muaver)
fval=zeros(length(yBC),1);
faobBC=yBC(1,:);
fnbBC=yBC(2,:);
fnspBC=yBC(3,:);
famxBC=yBC(4,:);
fhanBC=yBC(5,:);
for i=1
fval=[(muaob(i)-muaver(i))*faobBC(i)
(munb(i)-muaver(i))*fnbBC(i)
(munsp(i)-muaver(i))*fnspBC(i)
(muamx(i)-muaver(i))*famxBC(i)
(muhan(i)-muaver(i))*fhanBC(i)];
end
fval=[faobBC;fnbBC;fnspBC;famxBC;fhanBC];
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 = [0.5;0.5;0.5;0.5;0.5;0.5;1;1;1;1;1;0.4;0.02;0.02;0.3;0.26;0.001];
end
Walter Roberson
on 16 Aug 2023
Edited: Walter Roberson
on 16 Aug 2023
Observe that dSseBLdt is 2000 x 2000.
It is likely that you have accidentally combined row vector calculations with column vector calculations. For example, (1:3)' + [10 20 30] will not result in [11 22 33] and will instead give a 3 x 3 result.
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);
y0 = 0.5*ones(1,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,:).'; y16 = b(16,:).'; y17 = 1*10^-6;
y1=[y0(1); y0(2);y0(3);y0(4); y0(5);y0(6);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);
%main
[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);
Name Size Bytes Class Attributes
AL 1x1 8 double
DLSnh4 1x1 8 double
DLSno2 1x1 8 double
DLSno3 1x1 8 double
DLSo2 1x1 8 double
DLSse 1x1 8 double
DSnh4 1x1 8 double
DSno2 1x1 8 double
DSno3 1x1 8 double
DSo2 1x1 8 double
DSse 1x1 8 double
INBM 1x1 8 double
INXI 1x1 8 double
KH 1x1 8 double
KSsehan 1x1 8 double
KX 1x1 8 double
KaL 1x1 8 double
Knh4amx 1x1 8 double
Knh4aob 1x1 8 double
Kno2amx 1x1 8 double
Kno2han 1x1 8 double
Kno2nb 1x1 8 double
Kno2nsp 1x1 8 double
Kno3amx 1x1 8 double
Kno3aob 1x1 8 double
Kno3han 1x1 8 double
Kno3nb 1x1 8 double
Kno3nsp 1x1 8 double
Ko2amx 1x1 8 double
Ko2aob 1x1 8 double
Ko2han 1x1 8 double
Ko2nb 1x1 8 double
Ko2nsp 1x1 8 double
L 1x1 8 double
LL 1x1 8 double
Qo 1x1 8 double
Snh4BF 2001x1 16008 double
Snh4BFmax 1x1 8 double
Snh4BFmin 1x1 8 double
Snh4BL 1x1 8 double
Snh4o 1x1 8 double
Sno2BF 2001x1 16008 double
Sno2BFmax 1x1 8 double
Sno2BFmin 1x1 8 double
Sno2BL 1x1 8 double
Sno2o 1x1 8 double
Sno3BF 2001x1 16008 double
Sno3BFmax 1x1 8 double
Sno3BFmin 1x1 8 double
Sno3BL 1x1 8 double
Sno3o 1x1 8 double
So2BF 2001x1 16008 double
So2BFmax 1x1 8 double
So2BFmin 1x1 8 double
So2BL 1x1 8 double
So2o 1x1 8 double
SseBF 2001x1 16008 double
SseBFmax 1x1 8 double
SseBFmin 1x1 8 double
SseBL 1x1 8 double
Sseo 1x1 8 double
U 2000x2000 32000000 double
UL 2000x2000 32000000 double
VL 1x1 8 double
VR 1x1 8 double
Xs 1x1 8 double
Xso 1x1 8 double
YNo2han 1x1 8 double
YNo3han 1x1 8 double
Yamx 1x1 8 double
Yaob 1x1 8 double
Yhaer 1x1 8 double
Ynb 1x1 8 double
Ynsp 1x1 8 double
bamx 1x1 8 double
baob 1x1 8 double
bhan 1x1 8 double
bnb 1x1 8 double
bnsp 1x1 8 double
dLdt 1x1 8 double
dSnh4BFdt 2000x1 16000 double
dSnh4BLdt 2000x2000 32000000 double
dSno2BFdt 2000x1 16000 double
dSno2BLdt 2000x2000 32000000 double
dSno3BFdt 2000x1 16000 double
dSno3BLdt 2000x2000 32000000 double
dSo2BFdt 2000x1 16000 double
dSo2BLdt 2000x2000 32000000 double
dSseBFdt 2000x1 16000 double
dSseBLdt 2000x2000 32000000 double
dXsdt 1x1 8 double
dfamxdt 2000x1 16000 double
dfaobdt 2000x1 16000 double
dfhandt 2000x1 16000 double
dfnbdt 2000x1 16000 double
dfnspdt 2000x1 16000 double
dydt 1x1 32 function_handle
fI 1x1 8 double
famx 2001x1 16008 double
famxBC 1x1 8 double
faob 2001x1 16008 double
faobBC 1x1 8 double
fhan 2001x1 16008 double
fhanBC 1x1 8 double
fnb 2001x1 16008 double
fnbBC 1x1 8 double
fnsp 2001x1 16008 double
fnspBC 1x1 8 double
hz 1x1 8 double
ix 1x1 8 double
muamx 2000x2000 32000000 double
muaob 2000x1 16000 double
muaver 2000x2000 32000000 double
muhan 2000x2000 32000000 double
munb 2000x2000 32000000 double
munsp 2000x2000 32000000 double
namx 1x1 8 double
naob 1x1 8 double
nhan 1x1 8 double
nnb 1x1 8 double
nnsp 1x1 8 double
nx 1x1 8 double
sigma 1x1 8 double
t 1x1 8 double
u1 1x1 8 double
u2 1x1 8 double
u3 1x1 8 double
u5 1x1 8 double
u6 1x1 8 double
y 20007x1 160056 double
z 1x1 8 double
zeta 1x1 8 double
Error using vertcat
Dimensions of arrays being concatenated are not consistent.
Dimensions of arrays being concatenated are not consistent.
Error in solution>fun (line 135)
dydt=[dSseBLdt;dSno3BLdt;dSno2BLdt;dSnh4BLdt;dSo2BLdt;dXsdt;dSseBFdt;dSno3BFdt;dSno2BFdt;dSnh4BFdt;dSo2BFdt;dfaobdt;dfnbdt;dfnspdt;dfamxdt;dfhandt;dLdt];
Error in solution>@(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) (line 48)
[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);
Error in odearguments (line 92)
f0 = ode(t0,y0,args{:}); % ODE15I sets args{1} to yp0.
Error in ode15s (line 153)
odearguments(odeIsFuncHandle, odeTreatAsMFile, solver_name, ode, tspan, y0, options, varargin);
%% 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,:);
dydt=@fval;
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);
%BC-2
SseBFmax=SseBF(nx)+2*(nx-(nx-1))*L.*DLSse.*(SseBL-SseBF(nx))./LL;
Sno3BFmax=Sno3BF(nx)+2*(nx-(nx-1))*L.*DLSno3.*(Sno3BL-Sno3BF(nx))./LL;
Sno2BFmax=Sno2BF(nx)+2*(nx-(nx-1))*L.*DLSno2.*(Sno2BL-Sno2BF(nx))/LL;
Snh4BFmax=Snh4BF(nx)+2*(nx-(nx-1))*L.*DLSnh4.*(Snh4BL-Snh4BF(nx))/LL;
So2BFmax=So2BF(nx)+2*(nx-(nx-1))*L.*DLSo2.*(So2BL-So2BF(nx))./LL;
faobBC(1)=0.5; fnbBC(1)=0.5;fnspBC(1)=0.5;famxBC(1)=0.5; fhanBC(1)=0.5;
SseBF=[SseBF;SseBFmax];Sno3BF=[Sno3BF;Sno3BFmax];Sno2BF=[Sno2BF;Sno2BFmax];Snh4BF=[Snh4BF;Snh4BFmax];So2BF=[So2BF;So2BFmax];
faob=[faobBC;faob];fnb=[fnbBC;fnb];fnsp=[fnspBC;fnsp];famx=[famxBC;famx];fhan=[fhanBC;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-1
hz=nx-(nx-1);
z=0.012; zeta=z./L;
U=L.*muaver*zeta;
dSseBLdt=Qo.*(Sseo-SseBL)/VL-AL*((DSse/LL)-UL)*(SseBL-SseBF(nx));
dSno3BLdt=Qo.*(Sno3o-Sno3BL)./VL-AL.*((DSno3/LL)-UL).*(Sno3BL-Sno3BF(nx));
dSno2BLdt=Qo.*(Sno2o-Sno2BL)./VL-AL.*((DSno2/LL)-UL).*(Sno2BL-Sno2BF(nx));
dSnh4BLdt=Qo.*(Snh4o-Snh4BL)./VL-AL.*((DSnh4/LL)-UL).*(Snh4BL-Snh4BF(nx));
dSo2BLdt=Qo.*(So2o-So2BL)./VL + KaL.*(8-So2BL) - AL.*((DSo2)-UL).*(So2BL-So2BF(nx));
dXsdt=Qo.*(Xso-Xs)./VL - 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^2) +zeta.*UL(ix).*(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^2.*L^2) +zeta.*UL(ix).*(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)+Sno3BF(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^2.*L^2) +zeta.*UL(ix).*(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+So2BF(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^2.*L^2) +zeta.*UL(ix).*(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^2.*L^2) +zeta.*UL(ix).*(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(ix) -(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(ix)-muaver(ix)).*faob(ix) - (U(ix)-zeta.*UL(ix)).*(faob(ix+1)-faob(ix))./hz;
dfnbdt(ix)=(munb(ix)-muaver(ix)).*fnb(ix) - (U(ix)-zeta.*UL(ix)).*(fnb(ix+1)-fnb(ix))./hz;
dfnspdt(ix)=(munsp(ix)-muaver(ix)).*fnsp(ix) - (U(ix)-zeta.*UL(ix)).*(fnsp(ix+1)-fnsp(ix))./hz;
dfamxdt(ix)=(muamx(ix)-muaver(ix)).*famx(ix) - (U(ix)-zeta.*UL(ix)).*(famx(ix+1)-famx(ix))./hz;
dfhandt(ix)=(muhan(ix)-muaver(ix)).*fhan(ix) - (U(ix)-zeta.*UL(ix)).*(fhan(ix+1)-fhan(ix))./hz;
dLdt=L.*muaver(ix).*zeta+sigma;
end
function fval=BC(t,yBC,nx,muaob,munb,munsp,muamx,muhan,muaver)
fval=zeros(length(yBC),1);
faobBC=yBC(1,:);
fnbBC=yBC(2,:);
fnspBC=yBC(3,:);
famxBC=yBC(4,:);
fhanBC=yBC(5,:);
for i=1
fval=[(muaob(i)-muaver(i))*faobBC(i)
(munb(i)-muaver(i))*fnbBC(i)
(munsp(i)-muaver(i))*fnspBC(i)
(muamx(i)-muaver(i))*famxBC(i)
(muhan(i)-muaver(i))*fhanBC(i)];
end
fval=[faobBC;fnbBC;fnspBC;famxBC;fhanBC];
end
whos
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 = [0.5;0.5;0.5;0.5;0.5;0.5;1;1;1;1;1;0.4;0.02;0.02;0.3;0.26;0.001];
end
KIPROTICH KOSGEY
on 18 Aug 2023
Hi Walter,
How do you get the sizes of the variables before the program is executed?
If I could manage to do this then I could be checking their attributes as I make corrections.
Torsten
on 18 Aug 2023
KIPROTICH KOSGEY
on 25 Aug 2023
Edited: Torsten
on 25 Aug 2023
Thank you Torsten and Walter.
I have managed to remove all the bugs. But code execution is very slow and takes up a lot of memory. equations are also very stiff.
To save on memory, I have done program execution in bits (timewise). On stiffness, I am re-running the program from the time it fails using the final output from the previous run as input for the new runs. However, program execution is very slow. Could there be a way to speed the execution?
Please see the codes attached.
tstart = 0; tstart1=tSol(end);
The end operator must be used within an array index expression.
tend = 10*24;
nt = 80000;
xstart = 0;
xend = 0.012;
nx = 1000;
t = linspace(tstart,tend,nt); t1 = linspace(tstart1,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);
y0 = 0.5*ones(1,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,:).'; y16 = b(16,:).'; y17 = 1*10^-6;
y1=[y0(1); y0(2);y0(3);y0(4); y0(5);y0(6);y07; y08;y09;y10; y11;y12;y13; y14;y15;y16; y17];
%% execution
options = odeset('RelTol',1e-20,'AbsTol',1e-10,'NonNegative',1);
%main
[tSol, ySol]=ode15s(@(t,y) MBBR(t,y),t,y1);
y2=[ySol(end,1); ySol(end,2);ySol(end,3);ySol(end,4); ySol(end,5);ySol(end,6);ySol(end,7:nx+6)'; ySol(end,nx+7:2*nx+6)';ySol(end,2*nx+7:3*nx+6)';ySol(end,3*nx+7:4*nx+6)'; ySol(end,4*nx+7:5*nx+6)';ySol(end,5*nx+7:6*nx+6)';ySol(end,6*nx+7:7*nx+6)'; ySol(end,7*nx+7:8*nx+6)';ySol(end,8*nx+7:9*nx+6)';ySol(end,9*nx+7:10*nx+6)'; ySol(end,10*nx+7)];
[tSol1, ySol1]=ode15s(@(t1,y) MBBR(t1,y),t1,y2);
writematrix(ySol,'C:\Users\kosgey\Desktop\modelling\MANUSCRIPT\Book1.xls');
%% initial conditions
function u0 = oscic(x)
u0 = [0.5;0.5;0.5;0.5;0.5;0.5;1;1;1;1;1;0.4;0.02;0.02;0.3;0.26;0.001];
end
function dydt=MBBR(t,y)
tstart = 0;
tend = 10*24;
nt = 80000;
xstart = 0;
xend = 0.012;
nx = 1000;
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))];
%% 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; n=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;
KaL=0;
%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,:);
dydt=@fval;
muaob=u1*(So2BF/(Ko2aob+So2BF))*(Snh4BF./(Knh4aob+Snh4BF)); MDaobaer =baob*(So2BF./(Ko2aob+So2BF)); MDaobana =(baob*n*((Ko2aob/(Ko2aob+So2BF))*(Sno2BF+Sno3BF)/(Kno3aob+Sno2BF+Sno3BF)));
munb=u2*(Sno2BF/(Kno2nb+Sno2BF))*(So2BF./(Ko2nb+So2BF)); MDnbaer= bnb*(So2BF./(Ko2nb+So2BF)); MDnbana= (bnb*n*(Ko2nb/(Ko2nb+So2BF))*(Sno2BF+Sno3BF)/(Kno3nb+Sno2BF+Sno3BF));
munsp=u3*(Sno2BF/(Kno2nsp+Sno2BF))*(So2BF./(Ko2nsp+So2BF)); MDnspaer= bnsp*(So2BF./(Ko2nsp+So2BF)); MDnspana= (bnsp*n*(Ko2nsp/(Ko2nsp+So2BF))*(Sno2BF+Sno3BF)/(Kno3nsp+Sno2BF+Sno3BF));
muamx=u5*(Snh4BF./(Knh4amx+Snh4BF)).*(Sno2BF./(Kno2amx+Sno2BF)).*(Ko2amx./(Ko2amx+So2BF)); MDamxaer= bamx*(So2BF./(Ko2amx+So2BF)); MDamxana= (bamx*n*(Ko2amx/(Ko2amx+So2BF))*(Sno2BF+Sno3BF)/(Kno3amx+Sno2BF+Sno3BF));
muhanSno2=u6*(SseBF./(KSsehan+SseBF)).*(Sno2BF./(Kno2han+Sno2BF)).*(Ko2han./(Ko2han+So2BF)); muhanSno3=u6*(SseBF./(KSsehan+SseBF)).*(Sno3BF./(Kno3han+Sno3BF)).*(Ko2han./(Ko2han+So2BF)); muhanSo2=u6*(SseBF./(KSsehan+SseBF)).*(So2BF./(Ko2han+So2BF)); MDhanaer= bhan*(So2BF./(Ko2han+So2BF)); MDhanana = (bhan*n*(Ko2han/(Ko2han+So2BF))*(Sno2BF+Sno3BF)/(Kno3han+Sno2BF+Sno3BF));
muhan=muhanSo2+muhanSno3+muhanSno2;
muaver=faob.*muaob + fnb.*munb + fnsp.*munsp + famx.*muamx + fhan.*muhan;
UL=L*muaver+sigma;
VL=VR-AL*(L + LL);
z=0.012; zeta=z/L;
U=L*muaver*zeta;
m=0.01; %mass/unit_depth (mg/m)
fill_ratio=0.43;
VC=fill_ratio*VR; %volume occupied by carriers
X=5000;
if and (ge(t,0),le(t,5))
Snh4o=916; Qo=0.0452*VL*10^3/(916*24);
elseif and (gt(t,5),le(t,9))
Snh4o=1078; Qo=0.073*VL*10^3/(916*24);
elseif and (gt(t,9), le(t,12))
Snh4o=909;Qo=0.132*VL*10^3/(909*24);
elseif and (gt(t,12), le(t,15))
Snh4o=1078;Qo=0.191*VL*10^3/(1078*24);
elseif and (gt(t,15), le(t,21))
Snh4o=937;Qo=0.185*VL*10^3/(937*24);
elseif and (gt(t,21),le(t,33))
Snh4o=937;Qo=0.2286*VL*10^3/(937*24);
elseif and (gt(t,33),le(t,35))
Snh4o=937;Qo=0.3095*VL*10^3/(937*24);
elseif and (gt(t,35),le(t,37))
Snh4o=937;Qo=0.3158*VL*10^3/(937*24);
elseif and (gt(t,37),le(t,49))
Snh4o=937;Qo=0.348*VL*10^3/(937*24);
elseif and (gt(t,49),le(t,77))
Snh4o=937;Qo=0.448*VL*10^3/(937*24);
elseif and (gt(t,77),le(t,91))
Snh4o=1078; Qo=0.3928*VL*10^3/(1078*24);
elseif and (gt(t,91),le(t,119))
Snh4o=959;Qo=0.4582*VL*10^3/(959*24);
elseif and (gt(t,119),le(t,127))
Snh4o=959;Qo=0.1738*VL*10^3/(959*24);
elseif and (gt(t,127),le(t,134))
Snh4o=959;Qo=0.6246*VL*10^3/(959*24);
elseif and (gt(t,134),le(t,140))
Snh4o=1190;Qo=0.625*VL*10^3/(1190*24);
elseif and (gt(t,140),le(t,148))
Snh4o=1293;Qo=0.3929*VL*10^3/(1293*24);
elseif and (gt(t,148),le(t,155))
Snh4o=1087;Qo=0.8252*VL*10^3/(1087*24);
elseif and (gt(t,155), le(t,161))
Snh4o=1003;Qo=0.47*VL*10^3/(1003*24);
elseif and (gt(t,161),le(t,169))
Snh4o=1078;Qo=0.564*VL*10^3/(1078*24);
else
Snh4o=928; Qo=0.47*VL*10^3/(928*24);
end
while So2BL<0.5
KaL=80;
end
%BC-1
SseBFmin=SseBF(1); Sno3BFmin=Sno3BF(1); Sno2BFmin=Sno2BF(1); Snh4BFmin=Snh4BF(1);
So2BFmin=So2BF(1);
%BC-2
SseBFmax=SseBF(nx-1)+2*(nx-(nx-1))*L*DLSse*(SseBL-SseBF(nx))/LL;
Sno3BFmax=Sno3BF(nx-1)+2*(nx-(nx-1))*L*DLSno3*(Sno3BL-Sno3BF(nx))/LL;
Sno2BFmax=Sno2BF(nx-1)+2*(nx-(nx-1))*L*DLSno2*(Sno2BL-Sno2BF(nx))/LL;
Snh4BFmax=Snh4BF(nx-1)+2*(nx-(nx-1))*L*DLSnh4*(Snh4BL-Snh4BF(nx))/LL;
So2BFmax=So2BF(nx-1)+2*(nx-(nx-1))*L*DLSo2*(So2BL-So2BF(nx))/LL;
faobBC(1)=0.5; fnbBC(1)=0.5;fnspBC(1)=0.5;famxBC(1)=0.5; fhanBC(1)=0.5;
SseBF=[SseBF;SseBFmax];Sno3BF=[Sno3BF;Sno3BFmax];Sno2BF=[Sno2BF;Sno2BFmax];Snh4BF=[Snh4BF;Snh4BFmax];So2BF=[So2BF;So2BFmax];
faob=[faobBC;faob];fnb=[fnbBC;fnb];fnsp=[fnspBC;fnsp];famx=[famxBC;famx];fhan=[fhanBC;fhan];
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);
%% PDEs (7-16) and ODEs (1-6&17)
for ix=2:nx-1
dSseBLdt=Qo*(Sseo-SseBL)/VL-AL*((DSse/LL)-UL(nx))*(SseBL-SseBF(nx));
dSno3BLdt=Qo*(Sno3o-Sno3BL)/VL-AL*((DSno3/LL)-UL(nx))*(Sno3BL-Sno3BF(nx));
dSno2BLdt=Qo*(Sno2o-Sno2BL)/VL-AL*((DSno2/LL)-UL(nx))*(Sno2BL-Sno2BF(nx));
dSnh4BLdt=Qo*(Snh4o-Snh4BL)/VL-AL*((DSnh4/LL)-UL(nx))*(Snh4BL-Snh4BF(nx));
dSo2BLdt=Qo*(So2o-So2BL)/VL + KaL*(8-So2BL) - AL*((DSo2/LL)-UL(nx))*(So2BL-So2BF(nx));
dXsdt=Qo*(Xso-Xs)/VL - KH*((Xs/fhan(nx))/(KX+(Xs/fhan(nx))))*fhan(nx);
dSseBFdt(ix)=DSse*(SseBF(ix+1)-2*SseBF(ix)+(SseBF(ix-1)))/((x(nx+1)-x(nx-1))^2*L^2) +zeta*UL(ix)*(SseBF(ix+1)-SseBF(ix))/L - n*(1/YNo2han)*muhanSno2(ix)*fhan(ix)*X -n*(1/YNo3han)*muhanSno3(ix)*fhan(ix)*X - (1/Yhaer)*muhanSo2(ix)*fhan(ix)*X;
dSno3BFdt(ix)=DSno3*(Sno3BF(ix+1)-2*Sno3BF(ix)+(Sno3BF(ix-1)))/((x(nx+1)-x(nx-1))^2*L^2) +zeta*UL(ix)*(Sno3BF(ix+1)-Sno3BF(ix))/L - n*((1-YNo3han)/(2.86*YNo3han))*muhanSno3(ix)*fhan(ix)*X + (1/1.14)*muamx(ix)*famx(ix)*X + (1/Ynb)*munb(ix)*fnb(ix)*X + (1/Ynsp)*munsp(ix)*fnsp(ix)*X - ((1-fI)/2.86)*MDaobana(ix)*faob(ix)*X - ((1-fI)/2.86)*MDnbana(ix)*fnb(ix) - ((1-fI)/2.86)*MDnspana(ix)*fnsp(ix)*X - ((1-fI)/2.86)*MDamxana(ix)*famx(ix)*X-((1-fI)/2.86)*MDhanana(ix)*fhan(ix)*X;
dSno2BFdt(ix)=DSno2*(Sno2BF(ix+1)-2*Sno2BF(ix)+(Sno2BF(ix-1)))/((x(nx+1)-x(nx-1))^2*L^2) +zeta*UL(ix)*(Sno2BF(ix+1)-Sno2BF(ix))/L - ((1-YNo2han)/(1.71*YNo2han))*muhanSo2(ix)*fhan(ix)*X - ((1/Yamx)+(1/1.14))*muamx(ix)*famx(ix)*X + (1/Yaob)*muaob(ix)*faob(ix)*X - (1/Ynb)*munb(ix)*fnb(ix)*X - (1/Ynsp)*munb(ix)*fnsp(ix)*X;
dSnh4BFdt(ix)=DSnh4*(Snh4BF(ix+1)-2*Snh4BF(ix)+(Snh4BF(ix-1)))/((x(nx+1)-x(nx-1))^2*L^2) +zeta*UL(ix)*(Snh4BF(ix+1)-Snh4BF(ix))/L - (INBM+1/Yaob)*muaob(ix)*faob(ix)*X -(INBM+1/Yamx)*muamx(ix)*famx(ix)*X - INBM*munb(ix)*fnb(ix)*X -INBM*munsp(ix)*fnsp(ix)*X - INBM*n*muhanSno2(ix)*fhan(ix)*X - INBM*n*muhanSno3(ix)*fhan(ix)*X - INBM*muhanSo2(ix)*fhan(ix)*X + (INBM-(fI*INXI))*MDaobana(ix)*faob(ix)*X + (INBM-(fI*INXI))*MDnbana(ix)*fnb(ix)*X +(INBM-(fI*INXI))*MDnspana(ix)*fnsp(ix)*X + (INBM-(fI*INXI))*MDamxana(ix)*famx(ix)*X + (INBM-(fI*INXI))*MDhanana(ix)*fhan(ix)*X +(INBM-(fI*INXI))*MDaobaer(ix)*faob(ix)*X + (INBM-(fI*INXI))*MDnbaer(ix)*fnb(ix)*X +(INBM-(fI*INXI))*MDnspaer(ix)*fnsp(ix)*X + (INBM-(fI*INXI))*MDamxaer(ix)*famx(ix)*X + (INBM-(fI*INXI))*MDhanaer(ix)*fhan(ix)*X;
dSo2BFdt(ix)=DSo2*(So2BF(ix+1)-2*So2BF(ix)+(So2BF(ix-1)))/((x(nx+1)-x(nx-1))^2*L^2) +zeta*UL(ix)*(So2BF(ix+1)-So2BF(ix))/L -((1-Yhaer)/Yhaer)*muhanSo2(ix)*fhan(ix)*X -((3.43-Yaob)/Yaob)*muaob(ix)*faob(ix)*X - ((1.14-Ynb)/Ynb)*munb(ix)*fnb(ix)*X - ((1.14-Ynsp)/Ynsp)*munsp(ix)*fnsp(ix)*X - (1-fI)*MDaobaer(ix)*faob(ix)*X - (1-fI)*MDnbaer(ix)*fnb(ix)*X -(1-fI)*MDnspaer(ix)*fnsp(ix)*X - (1-fI)*MDamxaer(ix)*famx(ix)*X - (1-fI)*MDhanaer(ix)*fhan(ix);
dfaobdt(ix)=(muaob(ix)-muaver(ix))*faob(ix) - (U(ix)-zeta*UL(ix))*(faob(ix+1)-faob(ix))/(x(nx+1)-x(nx-1));
dfnbdt(ix)=(munb(ix)-muaver(ix))*fnb(ix) - (U(ix)-zeta*UL(ix))*(fnb(ix+1)-fnb(ix))/(x(nx+1)-x(nx-1));
dfnspdt(ix)=(munsp(ix)-muaver(ix))*fnsp(ix) - (U(ix)-zeta*UL(ix))*(fnsp(ix+1)-fnsp(ix))/(x(nx+1)-x(nx-1));
dfamxdt(ix)=(muamx(ix)-muaver(ix))*famx(ix) - (U(ix)-zeta*UL(ix))*(famx(ix+1)-famx(ix))/(x(nx+1)-x(nx-1));
dfhandt(ix)=(muhan(ix)-muaver(ix))*fhan(ix) - (U(ix)-zeta*UL(ix))*(fhan(ix+1)-fhan(ix))/(x(nx+1)-x(nx-1));
dLdt=L*muaver(nx)*zeta+sigma;
end
function fval=BC(t,yBC,nx,muaob,munb,munsp,muamx,muhan,muaver)
fval=zeros(length(yBC),1);
faobBC=yBC(1,:);
fnbBC=yBC(2,:);
fnspBC=yBC(3,:);
famxBC=yBC(4,:);
fhanBC=yBC(5,:);
for i=1
fval=[(muaob-muaver)*faobBC(i)
(munb-muaver)*fnbBC(i)
(munsp-muaver)*fnspBC(i)
(muamx-muaver)*famxBC(i)
(muhan-muaver)*fhanBC(i)];
end
fval=[faobBC;fnbBC;fnspBC;famxBC;fhanBC];
end
dydt=[dSseBLdt;dSno3BLdt;dSno2BLdt;dSnh4BLdt;dSo2BLdt;dXsdt;dSseBFdt;dSno3BFdt;dSno2BFdt;dSnh4BFdt;dSo2BFdt;dfaobdt;dfnbdt;dfnspdt;dfamxdt;dfhandt;dLdt];
end
KIPROTICH KOSGEY
on 25 Aug 2023
Edited: Torsten
on 25 Aug 2023
this is the code:
tstart = 0;
tend = 10*24;
nt = 160000;
xstart = 0;
xend = 0.012;
nx = 500;
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);
y0 = 0.5*ones(1,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,:).'; y16 = b(16,:).'; y17 = 1*10^-6;
y1=[y0(1); y0(2);y0(3);y0(4); y0(5);y0(6);y07; y08;y09;y10; y11;y12;y13; y14;y15;y16; y17];
%% execution
options = odeset('RelTol',1e-20,'AbsTol',1e-10,'NonNegative',1);
%main
[tSol, ySol]=ode15s(@(t,y) MBBR(t,y),t,y1);
%writematrix(ySol,'C:\Users\kosgey\Desktop\modelling\MANUSCRIPT\Book1.xls');
%% initial conditions
function u0 = oscic(x)
u0 = [0.5;0.5;0.5;0.5;0.5;0.5;1;1;1;1;1;0.4;0.02;0.02;0.3;0.26;0.001];
end
function dydt=MBBR(t,y)
tstart = 0;
tend = 10*24;
nt = 160000;
xstart = 0;
xend = 0.012;
nx = 500;
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))];
%% 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; n=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;
KaL=0;
%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,:);
dydt=@fval;
muaob=u1*(So2BF/(Ko2aob+So2BF))*(Snh4BF./(Knh4aob+Snh4BF)); MDaobaer =baob*(So2BF./(Ko2aob+So2BF)); MDaobana =(baob*n*((Ko2aob/(Ko2aob+So2BF))*(Sno2BF+Sno3BF)/(Kno3aob+Sno2BF+Sno3BF)));
munb=u2*(Sno2BF/(Kno2nb+Sno2BF))*(So2BF./(Ko2nb+So2BF)); MDnbaer= bnb*(So2BF./(Ko2nb+So2BF)); MDnbana= (bnb*n*(Ko2nb/(Ko2nb+So2BF))*(Sno2BF+Sno3BF)/(Kno3nb+Sno2BF+Sno3BF));
munsp=u3*(Sno2BF/(Kno2nsp+Sno2BF))*(So2BF./(Ko2nsp+So2BF)); MDnspaer= bnsp*(So2BF./(Ko2nsp+So2BF)); MDnspana= (bnsp*n*(Ko2nsp/(Ko2nsp+So2BF))*(Sno2BF+Sno3BF)/(Kno3nsp+Sno2BF+Sno3BF));
muamx=u5*(Snh4BF./(Knh4amx+Snh4BF)).*(Sno2BF./(Kno2amx+Sno2BF)).*(Ko2amx./(Ko2amx+So2BF)); MDamxaer= bamx*(So2BF./(Ko2amx+So2BF)); MDamxana= (bamx*n*(Ko2amx/(Ko2amx+So2BF))*(Sno2BF+Sno3BF)/(Kno3amx+Sno2BF+Sno3BF));
muhanSno2=u6*(SseBF./(KSsehan+SseBF)).*(Sno2BF./(Kno2han+Sno2BF)).*(Ko2han./(Ko2han+So2BF)); muhanSno3=u6*(SseBF./(KSsehan+SseBF)).*(Sno3BF./(Kno3han+Sno3BF)).*(Ko2han./(Ko2han+So2BF)); muhanSo2=u6*(SseBF./(KSsehan+SseBF)).*(So2BF./(Ko2han+So2BF)); MDhanaer= bhan*(So2BF./(Ko2han+So2BF)); MDhanana = (bhan*n*(Ko2han/(Ko2han+So2BF))*(Sno2BF+Sno3BF)/(Kno3han+Sno2BF+Sno3BF));
muhan=muhanSo2+muhanSno3+muhanSno2;
muaver=faob.*muaob + fnb.*munb + fnsp.*munsp + famx.*muamx + fhan.*muhan;
UL=L*muaver+sigma;
VL=VR-AL*(L + LL);
z=0.012; zeta=z/L;
U=L*muaver*zeta;
m=0.01; %mass/unit_depth (mg/m)
fill_ratio=0.43;
VC=fill_ratio*VR; %volume occupied by carriers
X=5000;
if and (ge(t,0),le(t,5))
Snh4o=916; Qo=0.0452*VL*10^3/(916*24);
elseif and (gt(t,5),le(t,9))
Snh4o=1078; Qo=0.073*VL*10^3/(916*24);
elseif and (gt(t,9), le(t,12))
Snh4o=909;Qo=0.132*VL*10^3/(909*24);
elseif and (gt(t,12), le(t,15))
Snh4o=1078;Qo=0.191*VL*10^3/(1078*24);
elseif and (gt(t,15), le(t,21))
Snh4o=937;Qo=0.185*VL*10^3/(937*24);
elseif and (gt(t,21),le(t,33))
Snh4o=937;Qo=0.2286*VL*10^3/(937*24);
elseif and (gt(t,33),le(t,35))
Snh4o=937;Qo=0.3095*VL*10^3/(937*24);
elseif and (gt(t,35),le(t,37))
Snh4o=937;Qo=0.3158*VL*10^3/(937*24);
elseif and (gt(t,37),le(t,49))
Snh4o=937;Qo=0.348*VL*10^3/(937*24);
elseif and (gt(t,49),le(t,77))
Snh4o=937;Qo=0.448*VL*10^3/(937*24);
elseif and (gt(t,77),le(t,91))
Snh4o=1078; Qo=0.3928*VL*10^3/(1078*24);
elseif and (gt(t,91),le(t,119))
Snh4o=959;Qo=0.4582*VL*10^3/(959*24);
elseif and (gt(t,119),le(t,127))
Snh4o=959;Qo=0.1738*VL*10^3/(959*24);
elseif and (gt(t,127),le(t,134))
Snh4o=959;Qo=0.6246*VL*10^3/(959*24);
elseif and (gt(t,134),le(t,140))
Snh4o=1190;Qo=0.625*VL*10^3/(1190*24);
elseif and (gt(t,140),le(t,148))
Snh4o=1293;Qo=0.3929*VL*10^3/(1293*24);
elseif and (gt(t,148),le(t,155))
Snh4o=1087;Qo=0.8252*VL*10^3/(1087*24);
elseif and (gt(t,155), le(t,161))
Snh4o=1003;Qo=0.47*VL*10^3/(1003*24);
elseif and (gt(t,161),le(t,169))
Snh4o=1078;Qo=0.564*VL*10^3/(1078*24);
else
Snh4o=928; Qo=0.47*VL*10^3/(928*24);
end
while So2BL<0.5
KaL=80;
end
%BC-1
SseBFmin=SseBF(1); Sno3BFmin=Sno3BF(1); Sno2BFmin=Sno2BF(1); Snh4BFmin=Snh4BF(1);
So2BFmin=So2BF(1);
%BC-2
SseBFmax=SseBF(nx-1)+2*(nx-(nx-1))*L*DLSse*(SseBL-SseBF(nx))/LL;
Sno3BFmax=Sno3BF(nx-1)+2*(nx-(nx-1))*L*DLSno3*(Sno3BL-Sno3BF(nx))/LL;
Sno2BFmax=Sno2BF(nx-1)+2*(nx-(nx-1))*L*DLSno2*(Sno2BL-Sno2BF(nx))/LL;
Snh4BFmax=Snh4BF(nx-1)+2*(nx-(nx-1))*L*DLSnh4*(Snh4BL-Snh4BF(nx))/LL;
So2BFmax=So2BF(nx-1)+2*(nx-(nx-1))*L*DLSo2*(So2BL-So2BF(nx))/LL;
faobBC(1)=0.5; fnbBC(1)=0.5;fnspBC(1)=0.5;famxBC(1)=0.5; fhanBC(1)=0.5;
SseBF=[SseBF;SseBFmax];Sno3BF=[Sno3BF;Sno3BFmax];Sno2BF=[Sno2BF;Sno2BFmax];Snh4BF=[Snh4BF;Snh4BFmax];So2BF=[So2BF;So2BFmax];
faob=[faobBC;faob];fnb=[fnbBC;fnb];fnsp=[fnspBC;fnsp];famx=[famxBC;famx];fhan=[fhanBC;fhan];
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);
%% PDEs (7-16) and ODEs (1-6&17)
for ix=2:nx-1
dSseBLdt=Qo*(Sseo-SseBL)/VL-AL*((DSse/LL)-UL(nx))*(SseBL-SseBF(nx));
dSno3BLdt=Qo*(Sno3o-Sno3BL)/VL-AL*((DSno3/LL)-UL(nx))*(Sno3BL-Sno3BF(nx));
dSno2BLdt=Qo*(Sno2o-Sno2BL)/VL-AL*((DSno2/LL)-UL(nx))*(Sno2BL-Sno2BF(nx));
dSnh4BLdt=Qo*(Snh4o-Snh4BL)/VL-AL*((DSnh4/LL)-UL(nx))*(Snh4BL-Snh4BF(nx));
dSo2BLdt=Qo*(So2o-So2BL)/VL + KaL*(8-So2BL) - AL*((DSo2/LL)-UL(nx))*(So2BL-So2BF(nx));
dXsdt=Qo*(Xso-Xs)/VL - KH*((Xs/fhan(nx))/(KX+(Xs/fhan(nx))))*fhan(nx);
dSseBFdt(ix)=DSse*(SseBF(ix+1)-2*SseBF(ix)+(SseBF(ix-1)))/((x(nx+1)-x(nx-1))^2*L^2) +zeta*UL(ix)*(SseBF(ix+1)-SseBF(ix))/L - n*(1/YNo2han)*muhanSno2(ix)*fhan(ix)*X -n*(1/YNo3han)*muhanSno3(ix)*fhan(ix)*X - (1/Yhaer)*muhanSo2(ix)*fhan(ix)*X;
dSno3BFdt(ix)=DSno3*(Sno3BF(ix+1)-2*Sno3BF(ix)+(Sno3BF(ix-1)))/((x(nx+1)-x(nx-1))^2*L^2) +zeta*UL(ix)*(Sno3BF(ix+1)-Sno3BF(ix))/L - n*((1-YNo3han)/(2.86*YNo3han))*muhanSno3(ix)*fhan(ix)*X + (1/1.14)*muamx(ix)*famx(ix)*X + (1/Ynb)*munb(ix)*fnb(ix)*X + (1/Ynsp)*munsp(ix)*fnsp(ix)*X - ((1-fI)/2.86)*MDaobana(ix)*faob(ix)*X - ((1-fI)/2.86)*MDnbana(ix)*fnb(ix) - ((1-fI)/2.86)*MDnspana(ix)*fnsp(ix)*X - ((1-fI)/2.86)*MDamxana(ix)*famx(ix)*X-((1-fI)/2.86)*MDhanana(ix)*fhan(ix)*X;
dSno2BFdt(ix)=DSno2*(Sno2BF(ix+1)-2*Sno2BF(ix)+(Sno2BF(ix-1)))/((x(nx+1)-x(nx-1))^2*L^2) +zeta*UL(ix)*(Sno2BF(ix+1)-Sno2BF(ix))/L - ((1-YNo2han)/(1.71*YNo2han))*muhanSo2(ix)*fhan(ix)*X - ((1/Yamx)+(1/1.14))*muamx(ix)*famx(ix)*X + (1/Yaob)*muaob(ix)*faob(ix)*X - (1/Ynb)*munb(ix)*fnb(ix)*X - (1/Ynsp)*munb(ix)*fnsp(ix)*X;
dSnh4BFdt(ix)=DSnh4*(Snh4BF(ix+1)-2*Snh4BF(ix)+(Snh4BF(ix-1)))/((x(nx+1)-x(nx-1))^2*L^2) +zeta*UL(ix)*(Snh4BF(ix+1)-Snh4BF(ix))/L - (INBM+1/Yaob)*muaob(ix)*faob(ix)*X -(INBM+1/Yamx)*muamx(ix)*famx(ix)*X - INBM*munb(ix)*fnb(ix)*X -INBM*munsp(ix)*fnsp(ix)*X - INBM*n*muhanSno2(ix)*fhan(ix)*X - INBM*n*muhanSno3(ix)*fhan(ix)*X - INBM*muhanSo2(ix)*fhan(ix)*X + (INBM-(fI*INXI))*MDaobana(ix)*faob(ix)*X + (INBM-(fI*INXI))*MDnbana(ix)*fnb(ix)*X +(INBM-(fI*INXI))*MDnspana(ix)*fnsp(ix)*X + (INBM-(fI*INXI))*MDamxana(ix)*famx(ix)*X + (INBM-(fI*INXI))*MDhanana(ix)*fhan(ix)*X +(INBM-(fI*INXI))*MDaobaer(ix)*faob(ix)*X + (INBM-(fI*INXI))*MDnbaer(ix)*fnb(ix)*X +(INBM-(fI*INXI))*MDnspaer(ix)*fnsp(ix)*X + (INBM-(fI*INXI))*MDamxaer(ix)*famx(ix)*X + (INBM-(fI*INXI))*MDhanaer(ix)*fhan(ix)*X;
dSo2BFdt(ix)=DSo2*(So2BF(ix+1)-2*So2BF(ix)+(So2BF(ix-1)))/((x(nx+1)-x(nx-1))^2*L^2) +zeta*UL(ix)*(So2BF(ix+1)-So2BF(ix))/L -((1-Yhaer)/Yhaer)*muhanSo2(ix)*fhan(ix)*X -((3.43-Yaob)/Yaob)*muaob(ix)*faob(ix)*X - ((1.14-Ynb)/Ynb)*munb(ix)*fnb(ix)*X - ((1.14-Ynsp)/Ynsp)*munsp(ix)*fnsp(ix)*X - (1-fI)*MDaobaer(ix)*faob(ix)*X - (1-fI)*MDnbaer(ix)*fnb(ix)*X -(1-fI)*MDnspaer(ix)*fnsp(ix)*X - (1-fI)*MDamxaer(ix)*famx(ix)*X - (1-fI)*MDhanaer(ix)*fhan(ix);
dfaobdt(ix)=(muaob(ix)-muaver(ix))*faob(ix) - (U(ix)-zeta*UL(ix))*(faob(ix+1)-faob(ix))/(x(nx+1)-x(nx-1));
dfnbdt(ix)=(munb(ix)-muaver(ix))*fnb(ix) - (U(ix)-zeta*UL(ix))*(fnb(ix+1)-fnb(ix))/(x(nx+1)-x(nx-1));
dfnspdt(ix)=(munsp(ix)-muaver(ix))*fnsp(ix) - (U(ix)-zeta*UL(ix))*(fnsp(ix+1)-fnsp(ix))/(x(nx+1)-x(nx-1));
dfamxdt(ix)=(muamx(ix)-muaver(ix))*famx(ix) - (U(ix)-zeta*UL(ix))*(famx(ix+1)-famx(ix))/(x(nx+1)-x(nx-1));
dfhandt(ix)=(muhan(ix)-muaver(ix))*fhan(ix) - (U(ix)-zeta*UL(ix))*(fhan(ix+1)-fhan(ix))/(x(nx+1)-x(nx-1));
dLdt=L*muaver(nx)*zeta+sigma;
end
function fval=BC(t,yBC,nx,muaob,munb,munsp,muamx,muhan,muaver)
fval=zeros(length(yBC),1);
faobBC=yBC(1,:);
fnbBC=yBC(2,:);
fnspBC=yBC(3,:);
famxBC=yBC(4,:);
fhanBC=yBC(5,:);
for i=1
fval=[(muaob-muaver)*faobBC(i)
(munb-muaver)*fnbBC(i)
(munsp-muaver)*fnspBC(i)
(muamx-muaver)*famxBC(i)
(muhan-muaver)*fhanBC(i)];
end
fval=[faobBC;fnbBC;fnspBC;famxBC;fhanBC];
end
dydt=[dSseBLdt;dSno3BLdt;dSno2BLdt;dSnh4BLdt;dSo2BLdt;dXsdt;dSseBFdt;dSno3BFdt;dSno2BFdt;dSnh4BFdt;dSo2BFdt;dfaobdt;dfnbdt;dfnspdt;dfamxdt;dfhandt;dLdt];
end
KIPROTICH KOSGEY
on 25 Aug 2023
Dear Torsten and Walter,
Kindly, could there be a way of circumventing the error: 'Warning: Failure at t=1.692955e+00. Unable to meet integration tolerances without reducing the step size below the smallest value allowed (3.552714e-15) at time t.
> In ode15s (line 662) In MBBRScript (line 28)' ?
Torsten
on 25 Aug 2023
Edited: Torsten
on 25 Aug 2023
Did you look at the results obtained so far (for times < 1.7 s) ? Do they make sense ?
Sorry, but it's really impossible to debug your code. Split the long expressions in smaller pieces according to the different parts of your equations, insert comments etc.
Don't treat all your equations in one big loop, but each equation separately.
It makes no sense to compute the ODEs for each ix = 2:nx-1 again and again.
dSseBLdt=Qo*(Sseo-SseBL)/VL-AL*((DSse/LL)-UL(nx))*(SseBL-SseBF(nx));
dSno3BLdt=Qo*(Sno3o-Sno3BL)/VL-AL*((DSno3/LL)-UL(nx))*(Sno3BL-Sno3BF(nx));
dSno2BLdt=Qo*(Sno2o-Sno2BL)/VL-AL*((DSno2/LL)-UL(nx))*(Sno2BL-Sno2BF(nx));
dSnh4BLdt=Qo*(Snh4o-Snh4BL)/VL-AL*((DSnh4/LL)-UL(nx))*(Snh4BL-Snh4BF(nx));
dSo2BLdt=Qo*(So2o-So2BL)/VL + KaL*(8-So2BL) - AL*((DSo2/LL)-UL(nx))*(So2BL-So2BF(nx));
dXsdt=Qo*(Xso-Xs)/VL - KH*((Xs/fhan(nx))/(KX+(Xs/fhan(nx))))*fhan(nx);
dLdt=L*muaver(nx)*zeta+sigma;
Further I don't see a call to "BC" in the visible part of your code. Maybe it's hidden somewhere in the endless lines where you compute the derivatives.
KIPROTICH KOSGEY
on 28 Aug 2023
Edited: Torsten
on 28 Aug 2023
Hi Torsten,
I figured out that computing ODEs with SseBF(nx),Sno3BF(nx),Sno2BF(nx),Snh4BF(nx) and So2BF(nx) while ix= 2:nx-1 resulted in a constant values of ODEs since I was using a none-existent values. Therefore I have changed ix=2:nx.
However, there seem to be a problem with the code, as matlab seem to be getting confused since it never finishes computations even after reducing nx and nt to 10 each. I have left MATLAB to run for over 24 hours but still it doesn't finish computing.
I have relooked at the codes and realised that the variables faob,fnb,fnsp,famx,fhan,SseBF,Sno3BF,Sno2BF,Snh4BF,So2BF are (nx+1)x1 but the results of their products (muaob,munb,munsp,muamx,muhan and muaver) are (nx)x1. How is this possible?
I don't know why MATLAB is not throwing an error in the for loop.
SCRPT:
tstart = 0;
tend = 10;
nt = 10;
xstart = 0;
xend = 0.012;
nx = 10;
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);
y0 = 0.5*ones(1,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,:).'; y16 = b(16,:).'; y17 = 1*10^-6;
y1=[y0(1); y0(2);y0(3);y0(4); y0(5);y0(6);y07; y08;y09;y10; y11;y12;y13; y14;y15;y16; y17];
%% execution
options = odeset('RelTol',1e-10,'AbsTol',1e-10,'NonNegative',1);
%main
[tSol, ySol]=ode15s(@(t,y) MBBR(t,y),t,y1);
%writematrix(ySol,'C:\Users\kosgey\Desktop\modelling\MANUSCRIPT\Book1.xls');
%% initial conditions
function u0 = oscic(x)
u0 = [0.5;0.5;0.5;0.5;0.5;0.5;1;1;1;1;1;0.4;0.02;0.02;0.3;0.26;0.001];
end
FUNCTION:
function dydt=MBBR(t,y)
tstart = 0;
tend = 10;
nt = 10;
xstart = 0;
xend = 0.012;
nx = 10;
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))];
%% 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; n=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;
%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,:);
dydt=@fval;
muaob=u1*(So2BF/(Ko2aob+So2BF))*(Snh4BF./(Knh4aob+Snh4BF));
MDaobaer =baob*(So2BF./(Ko2aob+So2BF)); MDaobana =(baob*n*((Ko2aob/(Ko2aob+So2BF))*(Sno2BF+Sno3BF)/(Kno3aob+Sno2BF+Sno3BF)))';
munb=u2*(Sno2BF/(Kno2nb+Sno2BF))*(So2BF./(Ko2nb+So2BF)); MDnbaer= bnb*(So2BF./(Ko2nb+So2BF)); MDnbana= (bnb*n*(Ko2nb/(Ko2nb+So2BF))*(Sno2BF+Sno3BF)/(Kno3nb+Sno2BF+Sno3BF))';
munsp=u3*(Sno2BF/(Kno2nsp+Sno2BF))*(So2BF./(Ko2nsp+So2BF)); MDnspaer= bnsp*(So2BF./(Ko2nsp+So2BF)); MDnspana= (bnsp*n*(Ko2nsp/(Ko2nsp+So2BF))*(Sno2BF+Sno3BF)/(Kno3nsp+Sno2BF+Sno3BF))';
muamx=u5*(Snh4BF./(Knh4amx+Snh4BF)).*(Sno2BF./(Kno2amx+Sno2BF)).*(Ko2amx./(Ko2amx+So2BF)); MDamxaer= bamx*(So2BF./(Ko2amx+So2BF)); MDamxana= (bamx*n*(Ko2amx/(Ko2amx+So2BF))*(Sno2BF+Sno3BF)/(Kno3amx+Sno2BF+Sno3BF))';
muhanSno2=u6*(SseBF./(KSsehan+SseBF)).*(Sno2BF./(Kno2han+Sno2BF)).*(Ko2han./(Ko2han+So2BF));
muhanSno3=u6*(SseBF./(KSsehan+SseBF)).*(Sno3BF./(Kno3han+Sno3BF)).*(Ko2han./(Ko2han+So2BF));
muhanSo2=u6*(SseBF./(KSsehan+SseBF)).*(So2BF./(Ko2han+So2BF)); MDhanaer= bhan*(So2BF./(Ko2han+So2BF)); MDhanana = (bhan*n*(Ko2han/(Ko2han+So2BF))*(Sno2BF+Sno3BF)/(Kno3han+Sno2BF+Sno3BF))';
muhan=muhanSo2+muhanSno3+muhanSno2;
muaver=faob.*muaob + fnb.*munb + fnsp.*munsp + famx.*muamx + fhan.*muhan;
UL=L*muaver+sigma;
VL=VR-AL*(L + LL);
z=0.012; zeta=z/L;
U=L*muaver*zeta;
X=5000; %mg/L
if and (ge(t,0),le(t,5))
Snh4o=916; Qo=0.0452*VL*10^3/(916*24);
elseif and (gt(t,5),le(t,9))
Snh4o=1078; Qo=0.073*VL*10^3/(916*24);
elseif and (gt(t,9), le(t,12))
Snh4o=909;Qo=0.132*VL*10^3/(909*24);
elseif and (gt(t,12), le(t,15))
Snh4o=1078;Qo=0.191*VL*10^3/(1078*24);
elseif and (gt(t,15), le(t,21))
Snh4o=937;Qo=0.185*VL*10^3/(937*24);
elseif and (gt(t,21),le(t,33))
Snh4o=937;Qo=0.2286*VL*10^3/(937*24);
elseif and (gt(t,33),le(t,35))
Snh4o=937;Qo=0.3095*VL*10^3/(937*24);
elseif and (gt(t,35),le(t,37))
Snh4o=937;Qo=0.3158*VL*10^3/(937*24);
elseif and (gt(t,37),le(t,49))
Snh4o=937;Qo=0.348*VL*10^3/(937*24);
elseif and (gt(t,49),le(t,77))
Snh4o=937;Qo=0.448*VL*10^3/(937*24);
elseif and (gt(t,77),le(t,91))
Snh4o=1078; Qo=0.3928*VL*10^3/(1078*24);
elseif and (gt(t,91),le(t,119))
Snh4o=959;Qo=0.4582*VL*10^3/(959*24);
elseif and (gt(t,119),le(t,127))
Snh4o=959;Qo=0.1738*VL*10^3/(959*24);
elseif and (gt(t,127),le(t,134))
Snh4o=959;Qo=0.6246*VL*10^3/(959*24);
elseif and (gt(t,134),le(t,140))
Snh4o=1190;Qo=0.625*VL*10^3/(1190*24);
elseif and (gt(t,140),le(t,148))
Snh4o=1293;Qo=0.3929*VL*10^3/(1293*24);
elseif and (gt(t,148),le(t,155))
Snh4o=1087;Qo=0.8252*VL*10^3/(1087*24);
elseif and (gt(t,155), le(t,161))
Snh4o=1003;Qo=0.47*VL*10^3/(1003*24);
elseif and (gt(t,161),le(t,169))
Snh4o=1078;Qo=0.564*VL*10^3/(1078*24);
else
Snh4o=928; Qo=0.47*VL*10^3/(928*24);
end
KaL=0;
while So2BL<0.5
KaL=80;
end
%Si=SseBF,dSno3BF,dSno2BF, dSnh4BF or dSo2BF
%BC-1 d(SseBF,dSno3BF,dSno2BF, dSnh4BF,dSo2BF)/dx=0;
SseBFmin=SseBF(1); Sno3BFmin=Sno3BF(1); Sno2BFmin=Sno2BF(1); Snh4BFmin=Snh4BF(1);
So2BFmin=So2BF(1);
%SLi=SseBL,dSno3BL,dSno2BL, dSnh4BL,dSo2BL
%BC-2 d(SseBF,dSno3BF,dSno2BF, dSnh4BF,dSo2BF)/dx=L*DLSse*(SLi-Si)/LL;
SseBFmax=SseBF(nx-1)+2*(x(nx)-x(nx-1))*L*DLSse*(SseBL-SseBF(nx))/LL;
Sno3BFmax=Sno3BF(nx-1)+2*(x(nx)-x(nx-1))*L*DLSno3*(Sno3BL-Sno3BF(nx))/LL;
Sno2BFmax=Sno2BF(nx-1)+2*(x(nx)-x(nx-1))*L*DLSno2*(Sno2BL-Sno2BF(nx))/LL;
Snh4BFmax=Snh4BF(nx-1)+2*(x(nx)-x(nx-1))*L*DLSnh4*(Snh4BL-Snh4BF(nx))/LL;
So2BFmax=So2BF(nx-1)+2*(x(nx)-x(nx-1))*L*DLSo2*(So2BL-So2BF(nx))/LL;
faobBC(1)=0.5; fnbBC(1)=0.5;fnspBC(1)=0.5;famxBC(1)=0.5; fhanBC(1)=0.5;
SseBF=[SseBF;SseBFmax];Sno3BF=[Sno3BF;Sno3BFmax];Sno2BF=[Sno2BF;Sno2BFmax];Snh4BF=[Snh4BF;Snh4BFmax];So2BF=[So2BF;So2BFmax];
faob=[faobBC;faob];fnb=[fnbBC;fnb];fnsp=[fnspBC;fnsp];famx=[famxBC;famx];fhan=[fhanBC;fhan];
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);
%% ODEs (1-6&17)
dSseBLdt=Qo*(Sseo-SseBL)/VL-AL*((DSse/LL)-UL(nx))*(SseBL-SseBFmax);
dSno3BLdt=Qo*(Sno3o-Sno3BL)/VL-AL*((DSno3/LL)-UL(nx))*(Sno3BL-Sno3BFmax);
dSno2BLdt=Qo*(Sno2o-Sno2BL)/VL-AL*((DSno2/LL)-UL(nx))*(Sno2BL-Sno2BFmax);
dSnh4BLdt=Qo*(Snh4o-Snh4BL)/VL-AL*((DSnh4/LL)-UL(nx))*(Snh4BL-Snh4BFmax);
dSo2BLdt=Qo*(So2o-So2BL)/VL + KaL*(8-So2BL) - AL*((DSo2/LL)-UL(nx))*(So2BL-So2BFmax);
dXsdt=Qo*(Xso-Xs)/VL - KH*((Xs/fhan(nx))/(KX+(Xs/fhan(nx))))*fhan(nx);
dLdt=L*muaver(nx)*zeta+sigma;
%% PDEs (7-16)
for ix=2:nx
dSseBFdt(ix)=DSse*(SseBF(ix+1)-2*SseBF(ix)+(SseBF(ix-1)))/((x(nx+1)-x(nx-1))^2*L^2) +zeta*UL(ix).*(SseBF(ix+1)-SseBF(ix))/L ...
- n*(1/YNo2han)*muhanSno2(ix).*fhan(ix)*X -n*(1/YNo3han)*muhanSno3(ix).*fhan(ix)*X - (1/Yhaer)*muhanSo2(ix).*fhan(ix)*X; %Sse consumption
dSno3BFdt(ix)=DSno3*(Sno3BF(ix+1)-2*Sno3BF(ix)+(Sno3BF(ix-1)))/((x(nx+1)-x(nx-1))^2*L^2) +zeta*UL(ix).*(Sno3BF(ix+1)-Sno3BF(ix))/L ...
- n*((1-YNo3han)/(2.86*YNo3han))*muhanSno3(ix).*fhan(ix)*X + (1/1.14)*muamx(ix).*famx(ix)*X + (1/Ynb)*munb(ix).*fnb(ix)*X ...
+ (1/Ynsp)*munsp(ix).*fnsp(ix)*X - ((1-fI)/2.86)*MDaobana(ix).*faob(ix)*X - ((1-fI)/2.86)*MDnbana(ix).*fnb(ix) ...
- ((1-fI)/2.86)*MDnspana(ix).*fnsp(ix)*X - ((1-fI)/2.86)*MDamxana(ix).*famx(ix)*X-((1-fI)/2.86)*MDhanana(ix).*fhan(ix)*X; %Sno3 consumption+generation
dSno2BFdt(ix)=DSno2*(Sno2BF(ix+1)-2*Sno2BF(ix)+(Sno2BF(ix-1)))/((x(nx+1)-x(nx-1))^2*L^2) +zeta*UL(ix).*(Sno2BF(ix+1)-Sno2BF(ix))/L ...
- ((1-YNo2han)/(1.71*YNo2han))*muhanSo2(ix).*fhan(ix)*X - ((1/Yamx)+(1/1.14))*muamx(ix).*famx(ix)*X + (1/Yaob)*muaob(ix).*faob(ix)*X ...
- (1/Ynb)*munb(ix).*fnb(ix)*X - (1/Ynsp)*munb(ix).*fnsp(ix)*X; %Sno2 comsumption+generation
dSnh4BFdt(ix)=DSnh4*(Snh4BF(ix+1)-2*Snh4BF(ix)+(Snh4BF(ix-1)))/((x(nx+1)-x(nx-1))^2*L^2) +zeta*UL(ix).*(Snh4BF(ix+1)-Snh4BF(ix))/L ...
- (INBM+1/Yaob)*muaob(ix).*faob(ix)*X -(INBM+1/Yamx)*muamx(ix).*famx(ix)*X - INBM*munb(ix).*fnb(ix)*X -INBM*munsp(ix).*fnsp(ix)*X ...
- INBM*n*muhanSno2(ix).*fhan(ix)*X - INBM*n*muhanSno3(ix).*fhan(ix)*X - INBM*muhanSo2(ix).*fhan(ix)*X + (INBM-(fI*INXI))*MDaobana(ix).*faob(ix)*X ...
+ (INBM-(fI*INXI))*MDnbana(ix).*fnb(ix)*X +(INBM-(fI*INXI))*MDnspana(ix).*fnsp(ix)*X + (INBM-(fI*INXI))*MDamxana(ix).*famx(ix)*X ...
+ (INBM-(fI*INXI))*MDhanana(ix).*fhan(ix)*X +(INBM-(fI*INXI))*MDaobaer(ix).*faob(ix)*X + (INBM-(fI*INXI))*MDnbaer(ix).*fnb(ix)*X ...
+(INBM-(fI*INXI))*MDnspaer(ix).*fnsp(ix)*X + (INBM-(fI*INXI))*MDamxaer(ix).*famx(ix)*X + (INBM-(fI*INXI))*MDhanaer(ix).*fhan(ix)*X;%Snh4 consumption+generation
dSo2BFdt(ix)=DSo2*(So2BF(ix+1)-2*So2BF(ix)+(So2BF(ix-1)))/((x(nx+1)-x(nx-1))^2*L^2) +zeta*UL(ix)*(So2BF(ix+1)-So2BF(ix))/L ...
-((1-Yhaer)/Yhaer)*muhanSo2(ix).*fhan(ix)*X -((3.43-Yaob)/Yaob)*muaob(ix).*faob(ix)*X - ((1.14-Ynb)/Ynb)*munb(ix).*fnb(ix)*X ...
- ((1.14-Ynsp)/Ynsp)*munsp(ix).*fnsp(ix)*X - (1-fI)*MDaobaer(ix).*faob(ix)*X - (1-fI)*MDnbaer(ix).*fnb(ix)*X -(1-fI)*MDnspaer(ix).*fnsp(ix)*X ...
- (1-fI)*MDamxaer(ix).*famx(ix)*X - (1-fI)*MDhanaer(ix).*fhan(ix);%So2 inflow+consumption
dfaobdt(ix)=(muaob(ix)-muaver(ix)).*faob(ix) - (U(ix)-zeta*UL(ix)).*(faob(ix+1)-faob(ix))/((x(nx+1)-x(nx-1))*L);
dfnbdt(ix)=(munb(ix)-muaver(ix)).*fnb(ix) - (U(ix)-zeta*UL(ix)).*(fnb(ix+1)-fnb(ix))/((x(nx+1)-x(nx-1))*L);
dfnspdt(ix)=(munsp(ix)-muaver(ix)).*fnsp(ix) - (U(ix)-zeta*UL(ix)).*(fnsp(ix+1)-fnsp(ix))/((x(nx+1)-x(nx-1))*L);
dfamxdt(ix)=(muamx(ix)-muaver(ix)).*famx(ix) - (U(ix)-zeta*UL(ix)).*(famx(ix+1)-famx(ix))/((x(nx+1)-x(nx-1))*L);
dfhandt(ix)=(muhan(ix)-muaver(ix)).*fhan(ix) - (U(ix)-zeta*UL(ix)).*(fhan(ix+1)-fhan(ix))/((x(nx+1)-x(nx-1))*L);
end
%%left BC for faob,fnb,fnsp,famx,fhan
function fval=BC(t,yBC,muaob,munb,munsp,muamx,muhan,muaver)
fval=zeros(length(yBC),1);
faobBC=yBC(1,:);
fnbBC=yBC(2,:);
fnspBC=yBC(3,:);
famxBC=yBC(4,:);
fhanBC=yBC(5,:);
for i=1
fval=[(muaob(i)-muaver(i))*faobBC(i)
(munb(i)-muaver(i))*fnbBC(i)
(munsp(i)-muaver(i))*fnspBC(i)
(muamx(i)-muaver(i))*famxBC(i)
(muhan(i)-muaver(i))*fhanBC(i)];
end
fval=[faobBC;fnbBC;fnspBC;famxBC;fhanBC];
end
dydt=[dSseBLdt;dSno3BLdt;dSno2BLdt;dSnh4BLdt;dSo2BLdt;dXsdt;dSseBFdt;dSno3BFdt;dSno2BFdt;dSnh4BFdt;dSo2BFdt;dfaobdt;dfnbdt;dfnspdt;dfamxdt;dfhandt;dLdt];
end
Torsten
on 28 Aug 2023
I have relooked at the codes and realised that the variables faob,fnb,fnsp,famx,fhan,SseBF,Sno3BF,Sno2BF,Snh4BF,So2BF are (nx+1)x1 but the results of their products (muaob,munb,munsp,muamx,muhan and muaver) are (nx)x1. How is this possible?
When you read faob,fnb,fnsp,famx,fhan,SseBF,Sno3BF,Sno2BF,Snh4BF,So2BF from the solution vector y, all these arrays become nx x 1 and so are muaob,munb,munsp,muamx,muhan,muaver that you deduce from them. Later on in your code, you most probaby enlarge your solution variable vectors by their boundary value such that they become (nx+1) x 1.
KIPROTICH KOSGEY
on 28 Aug 2023
Thank you Torsten. Your assistance is greatly appreciated, be blessed.
I have taken all your comments into account, though I sometimes come short as I am not a MATLAB expert.
I have moved the codes for computation of muaob,munb,munsp,muamx,muhan and muaver until after the codes for enlargement of faob,fnb,fnsp,famx,fhan,SseBF,Sno3BF,Sno2BF,Snh4BF,So2BF so that all the solution variables become (nx+1)x1. However, this has no effect on the computations since either way MATLAB was self-updating the changes.
As for the endless computations (stretching over 24 hours withouting completing), I have figured out that the while loop caused all these problems:
KaL=0;
while So2BL<0.5
KaL=80;
end
I have then edited these codes as:
KaL=0;
while So2BL<0.5
KaL=80;
break
end
Is this a correct move? Even so, So2BL keeps decreasing yet switching KaL between 0 and 80 is supposed to keep this value within some range.
Also, I am not quite sure on how to incorporate the left BC at x=0 for the pdes d(SseBF,dSno3BF,dSno2BF, dSnh4BF,dSo2BF)/dx=0. The integrals are: SseBFmin=SseBF(1); Sno3BFmin=Sno3BF(1); Sno2BFmin=Sno2BF(1); Snh4BFmin=Snh4BF(1); So2BFmin=So2BF(1);
when I set as SseBF=[SseBFmin;SseBFmax];Sno3BF=[Sno3BFmin;Sno3BFmax];Sno2BF=[Sno2BFmin;Sno2BFmax];Snh4BF=[Snh4BFmin;Snh4BFmax];So2BF=[So2BFmin;So2BFmax]; MATLAB throws an error in the for loop 'index must not exceeed 2' for these variables. So how do i factor in the left BCs?
Sorry to bother this much but I must say that I have learnt a lot from you and Walter throughout this exchange.
Thanks in advance.
Walter Roberson
on 28 Aug 2023
if and (ge(t,0),le(t,5))
Snh4o=916; Qo=0.0452*VL*10^3/(916*24);
elseif and (gt(t,5),le(t,9))
The mathematics of Runge Kutta ODE solvers requires that the code you provide must have continuous first and second derivatives. When that condition is violated, if you are lucky you get an error message about being unable to integrate at a particular time; if you are not lucky then you are not told that your results are invalid.
Your if/elseif looks to be time dependent. You need to recode as a series of calls to ode45(), one for each block of time that the constants such as Snh4o are consistent for, possibly adjusting the boundary conditions (for example if you were injecting chemicals at particular times.)
KIPROTICH KOSGEY
on 28 Aug 2023
Edited: Walter Roberson
on 28 Aug 2023
Thank you Torsten and Walter. Despite my efforts I have failed to achieve my aims.
I have fixed Snh4o and Qo so that none of the inputs vary with time but now is now throwing the error: Warning: Failure at t=1.317414e+00. Unable to meet integration tolerances without reducing the step size below the smallest value allowed (3.552714e-15) at time t.
I could not also code so that KaL is 80 when So2BL is <0.5 and 0 when above this value. I am giving up. Below is my final program, it might be of help to somebody else.
SCRIPT:
tstart = 0;
tend = 10*24;
nt = 30000;
xstart = 0;
xend = 0.012;
nx = 100;
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);
y0 = 0.5*ones(1,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,:).'; y16 = b(16,:).'; y17 = 1*10^-6;
y1=[y0(1); y0(2);y0(3);y0(4); y0(5);y0(6);y07; y08;y09;y10; y11;y12;y13; y14;y15;y16; y17];
%% execution
options = odeset('RelTol',1e-20,'AbsTol',1e-10,'NonNegative',1);
%main
[tSol, ySol]=ode15s(@(t,y) MBBR(t,y),t,y1);
Warning: Failure at t=1.691992e-01. Unable to meet integration tolerances without reducing the step size below the smallest value allowed (4.440892e-16) at time t.
%writematrix(ySol,'C:\Users\kosgey\Desktop\modelling\MANUSCRIPT\Book1.xls');
%% initial conditions
function u0 = oscic(x)
u0 = [0.5;0.5;0.5;0.5;0.5;0.5;1;1;1;1;1;0.4;0.02;0.02;0.3;0.26;0.001];
end
Function:
function dydt=MBBR(t,y)
tstart = 0;
tend = 10;
nt = 30000;
xstart = 0;
xend = 0.012;
nx = 100;
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))];
%% 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; n=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;
%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,:);
dydt=@fval;
VL=VR-AL*(L + LL); %m^3
Snh4o=900; %mg/L
Qo=VL*2;
KaL=80;
%Si=SseBF,dSno3BF,dSno2BF, dSnh4BF or dSo2BF
%BC-1 d(SseBF,dSno3BF,dSno2BF, dSnh4BF,dSo2BF)/dx=0;
SseBFmin=SseBF(1); Sno3BFmin=Sno3BF(1); Sno2BFmin=Sno2BF(1); Snh4BFmin=Snh4BF(1);
So2BFmin=So2BF(1);
%SLi=SseBL,dSno3BL,dSno2BL, dSnh4BL,dSo2BL
%BC-2 d(SseBF,dSno3BF,dSno2BF, dSnh4BF,dSo2BF)/dx=L*DLSse*(SLi-Si)/LL;
SseBFmax=SseBF(nx-1)+2*(x(nx)-x(nx-1))*L*DLSse*(SseBL-SseBF(nx))/LL;
Sno3BFmax=Sno3BF(nx-1)+2*(x(nx)-x(nx-1))*L*DLSno3*(Sno3BL-Sno3BF(nx))/LL;
Sno2BFmax=Sno2BF(nx-1)+2*(x(nx)-x(nx-1))*L*DLSno2*(Sno2BL-Sno2BF(nx))/LL;
Snh4BFmax=Snh4BF(nx-1)+2*(x(nx)-x(nx-1))*L*DLSnh4*(Snh4BL-Snh4BF(nx))/LL;
So2BFmax=So2BF(nx-1)+2*(x(nx)-x(nx-1))*L*DLSo2*(So2BL-So2BF(nx))/LL;
faobBC(1)=0.5; fnbBC(1)=0.5;fnspBC(1)=0.5;famxBC(1)=0.5; fhanBC(1)=0.5;
SseBF=[SseBF;SseBFmax];Sno3BF=[Sno3BF;Sno3BFmax];Sno2BF=[Sno2BF;Sno2BFmax];Snh4BF=[Snh4BF;Snh4BFmax];So2BF=[So2BF;So2BFmax];
faob=[faobBC;faob];fnb=[fnbBC;fnb];fnsp=[fnspBC;fnsp];famx=[famxBC;famx];fhan=[fhanBC;fhan];
muaob=u1*(So2BF./(Ko2aob+So2BF)).*(Snh4BF./(Knh4aob+Snh4BF)); MDaobaer =baob*(So2BF./(Ko2aob+So2BF)); MDaobana =(baob*n*((Ko2aob/(Ko2aob+So2BF))*(Sno2BF+Sno3BF)/(Kno3aob+Sno2BF+Sno3BF)))';
munb=u2*(Sno2BF/(Kno2nb+Sno2BF))*(So2BF./(Ko2nb+So2BF)); MDnbaer= bnb*(So2BF./(Ko2nb+So2BF)); MDnbana= (bnb*n*(Ko2nb/(Ko2nb+So2BF))*(Sno2BF+Sno3BF)/(Kno3nb+Sno2BF+Sno3BF))';
munsp=u3*(Sno2BF/(Kno2nsp+Sno2BF))*(So2BF./(Ko2nsp+So2BF)); MDnspaer= bnsp*(So2BF./(Ko2nsp+So2BF)); MDnspana= (bnsp*n*(Ko2nsp/(Ko2nsp+So2BF))*(Sno2BF+Sno3BF)/(Kno3nsp+Sno2BF+Sno3BF))';
muamx=u5*(Snh4BF./(Knh4amx+Snh4BF)).*(Sno2BF./(Kno2amx+Sno2BF)).*(Ko2amx./(Ko2amx+So2BF)); MDamxaer= bamx*(So2BF./(Ko2amx+So2BF)); MDamxana= (bamx*n*(Ko2amx/(Ko2amx+So2BF))*(Sno2BF+Sno3BF)/(Kno3amx+Sno2BF+Sno3BF))';
muhanSno2=u6*(SseBF./(KSsehan+SseBF)).*(Sno2BF./(Kno2han+Sno2BF)).*(Ko2han./(Ko2han+So2BF));
muhanSno3=u6*(SseBF./(KSsehan+SseBF)).*(Sno3BF./(Kno3han+Sno3BF)).*(Ko2han./(Ko2han+So2BF));
muhanSo2=u6*(SseBF./(KSsehan+SseBF)).*(So2BF./(Ko2han+So2BF)); MDhanaer= bhan*(So2BF./(Ko2han+So2BF)); MDhanana = (bhan*n*(Ko2han/(Ko2han+So2BF))*(Sno2BF+Sno3BF)/(Kno3han+Sno2BF+Sno3BF))';
muhan=muhanSo2+muhanSno3+muhanSno2;
muaver=faob.*muaob + fnb.*munb + fnsp.*munsp + famx.*muamx + fhan.*muhan;
UL=L*muaver+sigma;
z=0.012; zeta=z/L;
U=L*muaver*zeta;
X=5000; %mg/L
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);
%% ODEs (1-6&17)
dSseBLdt=Qo*(Sseo-SseBL)/VL-AL*((DSse/LL)-UL(nx))*(SseBL-SseBFmax);
dSno3BLdt=Qo*(Sno3o-Sno3BL)/VL-AL*((DSno3/LL)-UL(nx))*(Sno3BL-Sno3BFmax);
dSno2BLdt=Qo*(Sno2o-Sno2BL)/VL-AL*((DSno2/LL)-UL(nx))*(Sno2BL-Sno2BFmax);
dSnh4BLdt=Qo*(Snh4o-Snh4BL)/VL-AL*((DSnh4/LL)-UL(nx))*(Snh4BL-Snh4BFmax);
dSo2BLdt=Qo*(So2o-So2BL)/VL + KaL*(8-So2BL) - AL*((DSo2/LL)-UL(nx))*(So2BL-So2BFmax);
dXsdt=Qo*(Xso-Xs)/VL - KH*((Xs/(fhan(nx)*X))/(KX+(Xs/(fhan(nx)*X))))*(fhan(nx)*X);
dLdt=L*muaver(nx)*zeta+sigma;
%% PDEs (7-16)
for ix=2:nx
dSseBFdt(ix)=DSse*(SseBF(ix+1)-2*SseBF(ix)+(SseBF(ix-1)))/((x(nx+1)-x(nx-1))^2*L^2) +zeta*UL(ix).*(SseBF(ix+1)-SseBF(ix))/L ...
- n*(1/YNo2han)*muhanSno2(ix).*fhan(ix)*X -n*(1/YNo3han)*muhanSno3(ix).*fhan(ix)*X - (1/Yhaer)*muhanSo2(ix).*fhan(ix)*X; %Sse consumption
dSno3BFdt(ix)=DSno3*(Sno3BF(ix+1)-2*Sno3BF(ix)+(Sno3BF(ix-1)))/((x(nx+1)-x(nx-1))^2*L^2) +zeta*UL(ix).*(Sno3BF(ix+1)-Sno3BF(ix))/L ...
- n*((1-YNo3han)/(2.86*YNo3han))*muhanSno3(ix).*fhan(ix)*X + (1/1.14)*muamx(ix).*famx(ix)*X + (1/Ynb)*munb(ix).*fnb(ix)*X ...
+ (1/Ynsp)*munsp(ix).*fnsp(ix)*X - ((1-fI)/2.86)*MDaobana(ix).*faob(ix)*X - ((1-fI)/2.86)*MDnbana(ix).*fnb(ix) ...
- ((1-fI)/2.86)*MDnspana(ix).*fnsp(ix)*X - ((1-fI)/2.86)*MDamxana(ix).*famx(ix)*X-((1-fI)/2.86)*MDhanana(ix).*fhan(ix)*X; %Sno3 consumption+generation
dSno2BFdt(ix)=DSno2*(Sno2BF(ix+1)-2*Sno2BF(ix)+(Sno2BF(ix-1)))/((x(nx+1)-x(nx-1))^2*L^2) +zeta*UL(ix).*(Sno2BF(ix+1)-Sno2BF(ix))/L ...
- ((1-YNo2han)/(1.71*YNo2han))*muhanSo2(ix).*fhan(ix)*X - ((1/Yamx)+(1/1.14))*muamx(ix).*famx(ix)*X + (1/Yaob)*muaob(ix).*faob(ix)*X ...
- (1/Ynb)*munb(ix).*fnb(ix)*X - (1/Ynsp)*munb(ix).*fnsp(ix)*X; %Sno2 comsumption+generation
dSnh4BFdt(ix)=DSnh4*(Snh4BF(ix+1)-2*Snh4BF(ix)+(Snh4BF(ix-1)))/((x(nx+1)-x(nx-1))^2*L^2) +zeta*UL(ix).*(Snh4BF(ix+1)-Snh4BF(ix))/L ...
- (INBM+1/Yaob)*muaob(ix).*faob(ix)*X -(INBM+1/Yamx)*muamx(ix).*famx(ix)*X - INBM*munb(ix).*fnb(ix)*X -INBM*munsp(ix).*fnsp(ix)*X ...
- INBM*n*muhanSno2(ix).*fhan(ix)*X - INBM*n*muhanSno3(ix).*fhan(ix)*X - INBM*muhanSo2(ix).*fhan(ix)*X + (INBM-(fI*INXI))*MDaobana(ix).*faob(ix)*X ...
+ (INBM-(fI*INXI))*MDnbana(ix).*fnb(ix)*X +(INBM-(fI*INXI))*MDnspana(ix).*fnsp(ix)*X + (INBM-(fI*INXI))*MDamxana(ix).*famx(ix)*X ...
+ (INBM-(fI*INXI))*MDhanana(ix).*fhan(ix)*X +(INBM-(fI*INXI))*MDaobaer(ix).*faob(ix)*X + (INBM-(fI*INXI))*MDnbaer(ix).*fnb(ix)*X ...
+(INBM-(fI*INXI))*MDnspaer(ix).*fnsp(ix)*X + (INBM-(fI*INXI))*MDamxaer(ix).*famx(ix)*X + (INBM-(fI*INXI))*MDhanaer(ix).*fhan(ix)*X;%Snh4 consumption+generation
dSo2BFdt(ix)=DSo2*(So2BF(ix+1)-2*So2BF(ix)+(So2BF(ix-1)))/((x(nx+1)-x(nx-1))^2*L^2) +zeta*UL(ix)*(So2BF(ix+1)-So2BF(ix))/L ...
-((1-Yhaer)/Yhaer)*muhanSo2(ix).*fhan(ix)*X -((3.43-Yaob)/Yaob)*muaob(ix).*faob(ix)*X - ((1.14-Ynb)/Ynb)*munb(ix).*fnb(ix)*X ...
- ((1.14-Ynsp)/Ynsp)*munsp(ix).*fnsp(ix)*X - (1-fI)*MDaobaer(ix).*faob(ix)*X - (1-fI)*MDnbaer(ix).*fnb(ix)*X -(1-fI)*MDnspaer(ix).*fnsp(ix)*X ...
- (1-fI)*MDamxaer(ix).*famx(ix)*X - (1-fI)*MDhanaer(ix).*fhan(ix);%So2 inflow+consumption
dfaobdt(ix)=(muaob(ix)-muaver(ix)).*faob(ix) - (U(ix)-zeta*UL(ix)).*(faob(ix+1)-faob(ix))/((x(nx+1)-x(nx-1))*L);
dfnbdt(ix)=(munb(ix)-muaver(ix)).*fnb(ix) - (U(ix)-zeta*UL(ix)).*(fnb(ix+1)-fnb(ix))/((x(nx+1)-x(nx-1))*L);
dfnspdt(ix)=(munsp(ix)-muaver(ix)).*fnsp(ix) - (U(ix)-zeta*UL(ix)).*(fnsp(ix+1)-fnsp(ix))/((x(nx+1)-x(nx-1))*L);
dfamxdt(ix)=(muamx(ix)-muaver(ix)).*famx(ix) - (U(ix)-zeta*UL(ix)).*(famx(ix+1)-famx(ix))/((x(nx+1)-x(nx-1))*L);
dfhandt(ix)=(muhan(ix)-muaver(ix)).*fhan(ix) - (U(ix)-zeta*UL(ix)).*(fhan(ix+1)-fhan(ix))/((x(nx+1)-x(nx-1))*L);
end
%%left BC for faob,fnb,fnsp,famx,fhan
function fval=BC(t,yBC,muaob,munb,munsp,muamx,muhan,muaver)
fval=zeros(length(yBC),1);
faobBC=yBC(1,:);
fnbBC=yBC(2,:);
fnspBC=yBC(3,:);
famxBC=yBC(4,:);
fhanBC=yBC(5,:);
for i=1
fval=[(muaob(i)-muaver(i))*faobBC(i)
(munb(i)-muaver(i))*fnbBC(i)
(munsp(i)-muaver(i))*fnspBC(i)
(muamx(i)-muaver(i))*famxBC(i)
(muhan(i)-muaver(i))*fhanBC(i)];
end
fval=[faobBC;fnbBC;fnspBC;famxBC;fhanBC];
end
dydt=[dSseBLdt;dSno3BLdt;dSno2BLdt;dSnh4BLdt;dSo2BLdt;dXsdt;dSseBFdt;dSno3BFdt;dSno2BFdt;dSnh4BFdt;dSo2BFdt;dfaobdt;dfnbdt;dfnspdt;dfamxdt;dfhandt;dLdt];
end
Thanks again and be blessed!
Walter Roberson
on 28 Aug 2023
If you switch to ode23s then you can get up to about 1.32 seconds.
Once you do that, get busy plotting. I suggest plotting 10 at a time, such as
NN = size(ySol,2);
for K = 1 : 10 : NN
idx = K:min(K+9, NN);
figure
plot(tSol, ySol(:,idx));
legend("y_{" + idx + "}");
end
You will see that starting roughly ySol(:,510) that every 100 goes quite steep -- e.g., 519, 619, 719-ish are all steep.
But before that, starting from the second group of 10, ySol(:,11:20) that every plot goes steep around 1.3 for at least some of the outputs.
I am not going to analyze the code to figure out "why" . When it comes to ODEs, sometimes the reason is just "because it does" -- that either the equations are unstable (common!) or else that the equations need extended precision calculations to remain stable (certainly happens, but a lot of the time that loss of precision is suspected as the cause, it turns out that the equations are unstable and high precision calculations just postpones the time thte problem shows up at.)
KIPROTICH KOSGEY
on 29 Aug 2023
I am trying to simulate a biofilm process in which substrates diffuse through some thick layers where microorganisms are residing. This means that none of the variables should be negative. In addition, the key substrate So2 need to range between 0.5 and 1.5 mg/L. To control this value, I have was initially trying to switch KaL between 0 and 80: when So2BL<0.5, KaL=80, when So2BL>1.5, KaL=0, hence the while loop in the below codes. However, MATLAB wzs ignoring this for some reason:
while So2BL<0.5
KaL=80;
break
end
I also included 'NonNegative',1 in options so that none of the solution variables could be negative but then this also didn't help.
I therefore feel that if I could succeeed in setting the two parameters then I could stabilise the calculations.
Walter Roberson
on 29 Aug 2023
Use Event functions to at least detect the case where the values go out of physical range. But you will need to decide whether you should just terminate the system there (that is, it might be telling you that the equations or initial conditions must be wrong) — or if instead you want to adjust the system (for example corresponding to adding chemicals or changing temperature or the like.)
KIPROTICH KOSGEY
on 30 Aug 2023
Edited: Walter Roberson
on 30 Aug 2023
Dear Walter,
Your comments and suggestions are highly appreciated.
I have gone through the codes and didn't find any error relating to the equations. However, I see that faob goes to negative leading to unstable calculations relating to So2BF as a result of -*(-) leading to + when in fact the it should be -.
As suggested, I have introduced an event function in the script but MATLAB doesn't execute it, possibly because of my horrible coding style. I would like to stop the execution when any of the solution variables go to zero.
Sript:
t=[0 180*24];
t1=[0 5*24 9*24 12*24 15*24 21*24 33*24 35*24 37*24 49*24 77*24 91*24 119*24 127*24 134*24 140*24 148*24 155*24 161*24 169*24 180*24];
xstart = 0;
xend = 0.012;
nx=100;
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))];
nx=100;
a=arrayfun(@(i)oscic(x(i)),1:nx,'UniformOutput',false);
b=cell2mat(a);
y0 = 0.5*ones(1,5);
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,:).'; y16 = 1*10^-6;
y1=[y0(1); y0(2);y0(3);y0(4); y0(5);y06;y07; y08;y09;y10; y11;y12;y13; y14;y15;y16];
%% execution
options = odeset('RelTol',1e-20,'AbsTol',1e-10,'events',@eventsfcn);
%main
[tSol, ySol]=ode15s(@(t,y) MBBR(t,y),t,y1,options);
%writematrix(ySol,'C:\Users\kosgey\Desktop\modelling\MANUSCRIPT\Book1.xls');
%% initial conditions
function u0 = oscic(x)
u0 = [0.5;0.5;0.5;0.5;0.5;0.5;1;1;1;1;1;0.4;0.02;0.02;0.3;0.26;0.001];
end
function [value,isterminal,direction]=eventsfcn(t,y)
value=[y(1);y(2);y(3);y(4);y(5);y(6);y(7);y(8);y(9);y(10);y(11);y(12);y(13);y(14);y(15);y(16)];
isterminal = [1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1]; % stop at 0
direction = [-1;-1;-1;-1;-1;-1;-1;-1;-1;-1;-1;-1;-1;-1;-1;-1]; % decreasing
end
Function
function dydt=MBBR(t,y)
nx=100;
xstart = 0;
xend = 0.012;
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))];
%% 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; n=0.5;
KX=1; KH=0.125*exp(0.069*(T-293)); INXI=0.02;
% input parameters
Sseo=300; Sno3o=0.2; 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;
%Define the variables
dydt = zeros(length(y),1);
SseBL=y(1,:);
Sno3BL=y(2,:);
Sno2BL=y(3,:);
Snh4BL=y(4,:);
So2BL=y(5,:);
SseBF=y(6:nx+5,:);
Sno3BF=y(nx+6:2*nx+5,:);
Sno2BF=y(2*nx+6:3*nx+5,:);
Snh4BF=y(3*nx+6:4*nx+5,:);
So2BF=y(4*nx+6:5*nx+5,:);
faob=y(5*nx+6:6*nx+5,:);
fnb=y(6*nx+6:7*nx+5,:);
fnsp=y(7*nx+6:8*nx+5,:);
famx=y(8*nx+6:9*nx+5,:);
fhan=y(9*nx+6:10*nx+5,:);
L=y(10*nx+6,:);
VL=VR-AL*(L + LL); %m^3
Snh4o=900; %mg/L
Qo=VL*2;
KaL=80;
%Si=SseBF,dSno3BF,dSno2BF, dSnh4BF or dSo2BF
%BC-1 d(SseBF,dSno3BF,dSno2BF, dSnh4BF,dSo2BF)/dx=0;
SseBFmin=SseBF(1); Sno3BFmin=Sno3BF(1); Sno2BFmin=Sno2BF(1); Snh4BFmin=Snh4BF(1);
So2BFmin=So2BF(1);
%SLi=SseBL,dSno3BL,dSno2BL, dSnh4BL,dSo2BL
%BC-2 d(SseBF,dSno3BF,dSno2BF, dSnh4BF,dSo2BF)/dx=L*DLSse*(SLi-Si)/LL;
SseBFmax=SseBF(nx-1)+2*(x(nx)-x(nx-1))*L*DLSse*(SseBL-SseBF(nx))/LL;
Sno3BFmax=Sno3BF(nx-1)+2*(x(nx)-x(nx-1))*L*DLSno3*(Sno3BL-Sno3BF(nx))/LL;
Sno2BFmax=Sno2BF(nx-1)+2*(x(nx)-x(nx-1))*L*DLSno2*(Sno2BL-Sno2BF(nx))/LL;
Snh4BFmax=Snh4BF(nx-1)+2*(x(nx)-x(nx-1))*L*DLSnh4*(Snh4BL-Snh4BF(nx))/LL;
So2BFmax=So2BF(nx-1)+2*(x(nx)-x(nx-1))*L*DLSo2*(So2BL-So2BF(nx))/LL;
faobBC(1)=0.5; fnbBC(1)=0.5;fnspBC(1)=0.5;famxBC(1)=0.5; fhanBC(1)=0.5;
SseBF=[SseBF;SseBFmax];Sno3BF=[Sno3BF;Sno3BFmax];Sno2BF=[Sno2BF;Sno2BFmax];Snh4BF=[Snh4BF;Snh4BFmax];So2BF=[So2BF;So2BFmax];
faob=[faobBC;faob];fnb=[fnbBC;fnb];fnsp=[fnspBC;fnsp];famx=[famxBC;famx];fhan=[fhanBC;fhan];
muaob=u1*(So2BF./(Ko2aob+So2BF)).*(Snh4BF./(Knh4aob+Snh4BF)); MDaobaer =baob*(So2BF./(Ko2aob+So2BF)); MDaobana =(baob*n*((Ko2aob/(Ko2aob+So2BF))*(Sno2BF+Sno3BF)/(Kno3aob+Sno2BF+Sno3BF)))';
munb=u2*(Sno2BF/(Kno2nb+Sno2BF))*(So2BF./(Ko2nb+So2BF)); MDnbaer= bnb*(So2BF./(Ko2nb+So2BF)); MDnbana= (bnb*n*(Ko2nb/(Ko2nb+So2BF))*(Sno2BF+Sno3BF)/(Kno3nb+Sno2BF+Sno3BF))';
munsp=u3*(Sno2BF/(Kno2nsp+Sno2BF))*(So2BF./(Ko2nsp+So2BF)); MDnspaer= bnsp*(So2BF./(Ko2nsp+So2BF)); MDnspana= (bnsp*n*(Ko2nsp/(Ko2nsp+So2BF))*(Sno2BF+Sno3BF)/(Kno3nsp+Sno2BF+Sno3BF))';
muamx=u5*(Snh4BF./(Knh4amx+Snh4BF)).*(Sno2BF./(Kno2amx+Sno2BF)).*(Ko2amx./(Ko2amx+So2BF)); MDamxaer= bamx*(So2BF./(Ko2amx+So2BF)); MDamxana= (bamx*n*(Ko2amx/(Ko2amx+So2BF))*(Sno2BF+Sno3BF)/(Kno3amx+Sno2BF+Sno3BF))';
muhanSno2=u6*(SseBF./(KSsehan+SseBF)).*(Sno2BF./(Kno2han+Sno2BF)).*(Ko2han./(Ko2han+So2BF));
muhanSno3=u6*(SseBF./(KSsehan+SseBF)).*(Sno3BF./(Kno3han+Sno3BF)).*(Ko2han./(Ko2han+So2BF));
muhanSo2=u6*(SseBF./(KSsehan+SseBF)).*(So2BF./(Ko2han+So2BF)); MDhanaer= bhan*(So2BF./(Ko2han+So2BF)); MDhanana = (bhan*n*(Ko2han/(Ko2han+So2BF))*(Sno2BF+Sno3BF)/(Kno3han+Sno2BF+Sno3BF))';
muhan=muhanSo2+muhanSno3+muhanSno2;
kinert=0.1; %day^-1
finert=1-(faob+fnb+fnsp+famx+fhan);
muaver=faob.*muaob + fnb.*munb + fnsp.*munsp + famx.*muamx + fhan.*muhan+finert.*kinert*5;
UL=L*muaver+sigma;
z=0.012; zeta=z/L;
U=L*muaver*zeta;
X=5000; %mg/L
if and (ge(t,0),le(t,5))
Snh4o=916; Qo=0.0452*VL*10^3/(916*24);
elseif and (gt(t,5),le(t,9))
Snh4o=1078; Qo=0.073*VL*10^3/(916*24);
elseif and (gt(t,9), le(t,12))
Snh4o=909;Qo=0.132*VL*10^3/(909*24);
elseif and (gt(t,12), le(t,15))
Snh4o=1078;Qo=0.191*VL*10^3/(1078*24);
elseif and (gt(t,15), le(t,21))
Snh4o=937;Qo=0.185*VL*10^3/(937*24);
elseif and (gt(t,21),le(t,33))
Snh4o=937;Qo=0.2286*VL*10^3/(937*24);
elseif and (gt(t,33),le(t,35))
Snh4o=937;Qo=0.3095*VL*10^3/(937*24);
elseif and (gt(t,35),le(t,37))
Snh4o=937;Qo=0.3158*VL*10^3/(937*24);
elseif and (gt(t,37),le(t,49))
Snh4o=937;Qo=0.348*VL*10^3/(937*24);
elseif and (gt(t,49),le(t,77))
Snh4o=937;Qo=0.448*VL*10^3/(937*24);
elseif and (gt(t,77),le(t,91))
Snh4o=1078; Qo=0.3928*VL*10^3/(1078*24);
elseif and (gt(t,91),le(t,119))
Snh4o=959;Qo=0.4582*VL*10^3/(959*24);
elseif and (gt(t,119),le(t,127))
Snh4o=959;Qo=0.1738*VL*10^3/(959*24);
elseif and (gt(t,127),le(t,134))
Snh4o=959;Qo=0.6246*VL*10^3/(959*24);
elseif and (gt(t,134),le(t,140))
Snh4o=1190;Qo=0.625*VL*10^3/(1190*24);
elseif and (gt(t,140),le(t,148))
Snh4o=1293;Qo=0.3929*VL*10^3/(1293*24);
elseif and (gt(t,148),le(t,155))
Snh4o=1087;Qo=0.8252*VL*10^3/(1087*24);
elseif and (gt(t,155), le(t,161))
Snh4o=1003;Qo=0.47*VL*10^3/(1003*24);
elseif and (gt(t,161),le(t,169))
Snh4o=1078;Qo=0.564*VL*10^3/(1078*24);
else
Snh4o=928; Qo=0.47*VL*10^3/(928*24);
end
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);
%% ODEs (1-6&17)
dSseBLdt=Qo*(Sseo-SseBL)/VL-AL*((DSse/LL)-UL(nx))*(SseBL-SseBF(nx));
dSno3BLdt=Qo*(Sno3o-Sno3BL)/VL-AL*((DSno3/LL)-UL(nx))*(Sno3BL-Sno3BF(nx));
dSno2BLdt=Qo*(Sno2o-Sno2BL)/VL-AL*((DSno2/LL)-UL(nx))*(Sno2BL-Sno2BF(nx));
dSnh4BLdt=Qo*(Snh4o-Snh4BL)/VL-AL*((DSnh4/LL)-UL(nx))*(Snh4BL-Snh4BF(nx));
dSo2BLdt=Qo*(So2o-So2BL)/VL + KaL*(8-So2BL) - AL*((DSo2/LL)-UL(nx))*(So2BL-So2BF(nx));
dLdt=L*muaver(nx)*zeta+sigma;
%% PDEs (7-16)
for ix=2:nx
dSseBFdt(ix)=DSse*(SseBF(ix+1)-2*SseBF(ix)+(SseBF(ix-1)))/((x(nx+1)-x(nx-1))^2*L^2) +zeta*UL(ix).*(SseBF(ix+1)-SseBF(ix))/L ...
- n*(1/YNo2han)*muhanSno2(ix).*fhan(ix)*X -n*(1/YNo3han)*muhanSno3(ix).*fhan(ix)*X - (1/Yhaer)*muhanSo2(ix).*fhan(ix)*X; %Sse consumption
dSno3BFdt(ix)=DSno3*(Sno3BF(ix+1)-2*Sno3BF(ix)+(Sno3BF(ix-1)))/((x(nx+1)-x(nx-1))^2*L^2) +zeta*UL(ix).*(Sno3BF(ix+1)-Sno3BF(ix))/L ...
- n*((1-YNo3han)/(2.86*YNo3han))*muhanSno3(ix).*fhan(ix)*X + (1/1.14)*muamx(ix).*famx(ix)*X + (1/Ynb)*munb(ix).*fnb(ix)*X ...
+ (1/Ynsp)*munsp(ix).*fnsp(ix)*X - ((1-fI)/2.86)*MDaobana(ix).*faob(ix)*X - ((1-fI)/2.86)*MDnbana(ix).*fnb(ix) ...
- ((1-fI)/2.86)*MDnspana(ix).*fnsp(ix)*X - ((1-fI)/2.86)*MDamxana(ix).*famx(ix)*X-((1-fI)/2.86)*MDhanana(ix).*fhan(ix)*X; %Sno3 consumption+generation
dSno2BFdt(ix)=DSno2*(Sno2BF(ix+1)-2*Sno2BF(ix)+(Sno2BF(ix-1)))/((x(nx+1)-x(nx-1))^2*L^2) +zeta*UL(ix).*(Sno2BF(ix+1)-Sno2BF(ix))/L ...
- ((1-YNo2han)/(1.71*YNo2han))*muhanSo2(ix).*fhan(ix)*X - ((1/Yamx)+(1/1.14))*muamx(ix).*famx(ix)*X + (1/Yaob)*muaob(ix).*faob(ix)*X ...
- (1/Ynb)*munb(ix).*fnb(ix)*X - (1/Ynsp)*munb(ix).*fnsp(ix)*X; %Sno2 comsumption+generation
dSnh4BFdt(ix)=DSnh4*(Snh4BF(ix+1)-2*Snh4BF(ix)+(Snh4BF(ix-1)))/((x(nx+1)-x(nx-1))^2*L^2) +zeta*UL(ix).*(Snh4BF(ix+1)-Snh4BF(ix))/L ...
- (INBM+1/Yaob)*muaob(ix).*faob(ix)*X -(INBM+1/Yamx)*muamx(ix).*famx(ix)*X - INBM*munb(ix).*fnb(ix)*X -INBM*munsp(ix).*fnsp(ix)*X ...
- INBM*n*muhanSno2(ix).*fhan(ix)*X - INBM*n*muhanSno3(ix).*fhan(ix)*X - INBM*muhanSo2(ix).*fhan(ix)*X + (INBM-(fI*INXI))*MDaobana(ix).*faob(ix)*X ...
+ (INBM-(fI*INXI))*MDnbana(ix).*fnb(ix)*X +(INBM-(fI*INXI))*MDnspana(ix).*fnsp(ix)*X + (INBM-(fI*INXI))*MDamxana(ix).*famx(ix)*X ...
+ (INBM-(fI*INXI))*MDhanana(ix).*fhan(ix)*X +(INBM-(fI*INXI))*MDaobaer(ix).*faob(ix)*X + (INBM-(fI*INXI))*MDnbaer(ix).*fnb(ix)*X ...
+(INBM-(fI*INXI))*MDnspaer(ix).*fnsp(ix)*X + (INBM-(fI*INXI))*MDamxaer(ix).*famx(ix)*X + (INBM-(fI*INXI))*MDhanaer(ix).*fhan(ix)*X;%Snh4 consumption+generation
dSo2BFdt(ix)=DSo2*(So2BF(ix+1)-2*So2BF(ix)+(So2BF(ix-1)))/((x(nx+1)-x(nx-1))^2*L^2) +zeta*UL(ix)*(So2BF(ix+1)-So2BF(ix))/L ...
-((1-Yhaer)/Yhaer)*muhanSo2(ix).*fhan(ix)*X -((3.43-Yaob)/Yaob)*muaob(ix).*faob(ix)*X - ((1.14-Ynb)/Ynb)*munb(ix).*fnb(ix)*X ...
- ((1.14-Ynsp)/Ynsp)*munsp(ix).*fnsp(ix)*X - (1-fI)*MDaobaer(ix).*faob(ix)*X - (1-fI)*MDnbaer(ix).*fnb(ix)*X -(1-fI)*MDnspaer(ix).*fnsp(ix)*X ...
- (1-fI)*MDamxaer(ix).*famx(ix)*X - (1-fI)*MDhanaer(ix).*fhan(ix);%So2 inflow+consumption
dfaobdt(ix)=(muaob(ix)-muaver(ix)).*faob(ix) - (U(ix)-zeta*UL(ix)).*(faob(ix+1)-faob(ix))/((x(nx+1)-x(nx-1))*L);
dfnbdt(ix)=(munb(ix)-muaver(ix)).*fnb(ix) - (U(ix)-zeta*UL(ix)).*(fnb(ix+1)-fnb(ix))/((x(nx+1)-x(nx-1))*L);
dfnspdt(ix)=(munsp(ix)-muaver(ix)).*fnsp(ix) - (U(ix)-zeta*UL(ix)).*(fnsp(ix+1)-fnsp(ix))/((x(nx+1)-x(nx-1))*L);
dfamxdt(ix)=(muamx(ix)-muaver(ix)).*famx(ix) - (U(ix)-zeta*UL(ix)).*(famx(ix+1)-famx(ix))/((x(nx+1)-x(nx-1))*L);
dfhandt(ix)=(muhan(ix)-muaver(ix)).*fhan(ix) - (U(ix)-zeta*UL(ix)).*(fhan(ix+1)-fhan(ix))/((x(nx+1)-x(nx-1))*L);
end
%%left BC for faob,fnb,fnsp,famx,fhan
dydt=@fval;
function fval=BC(t,yBC,muaob,munb,munsp,muamx,muhan,muaver)
fval=zeros(length(yBC),1);
faobBC=yBC(1,:);
fnbBC=yBC(2,:);
fnspBC=yBC(3,:);
famxBC=yBC(4,:);
fhanBC=yBC(5,:);
for i=1
fval=[(muaob(i)-muaver(i))*faobBC(i)
(munb(i)-muaver(i))*fnbBC(i)
(munsp(i)-muaver(i))*fnspBC(i)
(muamx(i)-muaver(i))*famxBC(i)
(muhan(i)-muaver(i))*fhanBC(i)];
end
fval=[faobBC;fnbBC;fnspBC;famxBC;fhanBC];
end
dydt=[dSseBLdt;dSno3BLdt;dSno2BLdt;dSnh4BLdt;dSo2BLdt;dSseBFdt;dSno3BFdt;dSno2BFdt;dSnh4BFdt;dSo2BFdt;dfaobdt;dfnbdt;dfnspdt;dfamxdt;dfhandt;dLdt];
end
Torsten
on 30 Aug 2023
Edited: Torsten
on 30 Aug 2023
y06-y15 all have nx elements, not only one.
Thus
value=[y(1);y(2);y(3);y(4);y(5);y(6);y(7);y(8);y(9);y(10);y(11);y(12);y(13);y(14);y(15);y(16)];
isterminal = [1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1]; % stop at 0
direction = [-1;-1;-1;-1;-1;-1;-1;-1;-1;-1;-1;-1;-1;-1;-1;-1]; % decreasing
is wrong.
But how should the event function help in your case ? The solver will stop integrating almost instantly. If you don't have a strategy on how to continue after the event occurs, you are left with a solution vector after 1 sec or so.
KIPROTICH KOSGEY
on 30 Aug 2023
Edited: Walter Roberson
on 30 Aug 2023
Dear Torsten and Walter.
Thank you so much for giving me hints on possible solutions to the coding mistakes I have been making. I have learnt a lot in the course of this exchange.
I have since edited the script as shown below and coded for the program when So2BL is above or below the given range (since this variable affects all the other variables directly leading to unstable calculations and failures). I have also divided the functions into two to gather for different periods (aerobic (increasing So2BL because KaL=80) and anaerobic (decreasing So2BL because KaL=0)). Although it apperas to be working, I keep getting the error: out of memory. Could there be a way to go around this?
Main functions attached.
Script:
tend=180*24;
t=[0 180*24];
t1=[0 5*24 9*24 12*24 15*24 21*24 33*24 35*24 37*24 49*24 77*24 91*24 119*24 127*24 134*24 140*24 148*24 155*24 161*24 169*24 180*24];
xstart = 0;
xend = 0.012;
nx=100;
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))];
nx=100;
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-20,'AbsTol',1e-10,'events',@eventsfcn_aerobic);
options2 = odeset('RelTol',1e-20,'AbsTol',1e-10,'events',@eventsfcn_anoxic);
%main
[tSol, ySol]=ode15s(@(t,y) MBBR_aerobic(t,y),t,y1,options1);
tplot=tSol;
yplot=ySol;
for i=1:2
y2=[ySol(end,1);ySol(end,2);ySol(end,3);ySol(end,4);ySol(end,5);ySol(end,6:105)';ySol(end,106:205)';ySol(end,206:305)';ySol(end,306:405)';...
ySol(end,406:505)';ySol(end,506:605)';ySol(end,606:705)';ySol(end,706:805)';ySol(end,806:905)';ySol(end,906:1005)';ySol(end,1006)'];
[tSol,ySol]=ode45(@MBBR_anoxic,[tplot(end) tend],y2,options2);
tplot=[tplot;tSol];
yplot=[yplot;ySol];
end
%% initial conditions
function u0 = oscic(x)
u0 = [0.5;0.5;0.5;0.5;0.5;0.5;1;1;1;1;1;0.4;0.02;0.02;0.3;0.26;0.001];
end
%% events aerobic phase
function [value,isterminal,direction]=eventsfcn_aerobic(t,y)
z=y(5);
z1=ceil(1.5-z);
value=z1;
isterminal = 1; % stop at 0
direction = -1; % increasing
end
%% events anoxic phase
function [value,isterminal,direction]=eventsfcn_anoxic(t,y)
z=y(5);
z1=ceil(0.5-z);
value=z1;
isterminal = 1; % stop at 0
direction = 1; % decreasing
end
Torsten
on 30 Aug 2023
Edited: Torsten
on 30 Aug 2023
In which phase of execution of your program do you get the error message ? While using ode15s or ode45 ? If ode45 is called with a two-element vector for tspan as you do, it will collect so many output solutions that you will get a memory overflow.
Why do you use ode45 at all ? It is not suited for your kind of problem.
KIPROTICH KOSGEY
on 30 Aug 2023
I was getting an error in the second phase of execution and I agree with you that ode45 was causing this error! I don't even know why I used ode45 in the first plast. I should have also used yplot instead of ySol in the initialisation of the second phase.
However, after correcting this, I have found that the program is not calling the second function MBBR_anoxic which is supposed to lead to a decrease in So2BL, but rather continues with the first program MBBR_aerobic leading to an exponential increasing in So2BL until the program fails.
Torsten
on 30 Aug 2023
However, after correcting this, I have found that the program is not calling the second function MBBR_anoxic which is supposed to lead to a decrease in So2BL, but rather continues with the first program MBBR_aerobic leading to an exponential increasing in So2BL until the program fails.
Seems you missed to switch a + to a - somewhere in your code for the second phase.
KIPROTICH KOSGEY
on 1 Sep 2023
Edited: Walter Roberson
on 1 Sep 2023
Thank you Torsten, I found the problem leading to the below issue. However, after editting it and modifying the script, MATLAB stops just below y(5) is 1.5 yet I want it to stop when y(5)>1.5.
t=[0 180*24];
xstart = 0;
xend = 0.012;
nx=100;
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))];
nx=100;
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-20,'AbsTol',1e-10,'events',@eventsfcn_aerobic);
options2 = odeset('RelTol',1e-20,'AbsTol',1e-10,'events',@eventsfcn_anoxic);
%main
[tSol, ySol]=ode15s(@(t,y) MBBR_aerobic(t,y),t,y1,options1);
tplot=tSol;
yplot=ySol;
for i=1:2
y2=[yplot(end,1);yplot(end,2);yplot(end,3);yplot(end,4);yplot(end,5);yplot(end,6:105)';yplot(end,106:205)';yplot(end,206:305)';yplot(end,306:405)';...
yplot(end,406:505)';yplot(end,506:605)';yplot(end,606:705)';yplot(end,706:805)';yplot(end,806:905)';yplot(end,906:1005)';yplot(end,1006)];
[tSol,ySol,tev,yev]=ode15s(@MBBR_anoxic,[tplot(end) 5*180],y2,options2);
tplot=[tplot;tSol];
yplot=[yplot;ySol];
end
%% initial conditions
function u0 = oscic(x)
u0 = [0.1;0.1;0.1;0.1;0.1;0.1;0.1;0.1;0.1;0.1;0.1;0.4;0.02;0.02;0.3;0.26;0.001];
end
%% events aerobic phase
function [value,isterminal,direction]=eventsfcn_aerobic(t,y)
z=y(5);
z1=+(z<=1.5);
value=z1;
isterminal = 1; % stop at 0
direction = []; % decreasing
end
%% events anoxic phase
function [value,isterminal,direction]=eventsfcn_anoxic(t,y)
e=y(5);
e1=+(e>0.5);
value=e1;
isterminal = 1; % stop at 0
direction = []; % decreasing
end
Torsten
on 1 Sep 2023
The event function you define must be differentiable with respect to the solution variables. z1 and e1 are not differentiable, but jump functions.
KIPROTICH KOSGEY
on 1 Sep 2023
Edited: Walter Roberson
on 1 Sep 2023
I am learning as I have always done from the time I sought assistance regarding my codes.
with z1 and e1, the program is running and the variable y(5) is changing as expected. However, on the third round, the solution variables are curiously constant when they should have been changing. Could this be related to the issue you have raised?
Basically, I would like to keep y(5) between 0.5 and 1.5. An event as per the examples in MATLAB documentation is recorded when a function crosses zero. So in my case, how can I meet this requirement whereas my variables range above this?
This is my script:
close all; clear; clc;
t1=[0 180*24];
xstart = 0;
xend = 0.012;
nx=100;
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))];
nx=100;
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-20,'AbsTol',1e-10,'events',@eventsfcn_aerobic);
options2 = odeset('RelTol',1e-20,'AbsTol',1e-10,'events',@eventsfcn_anoxic);
%main
[tSol, ySol]=ode15s(@(t,y) MBBR_aerobic(t,y),t1,y1,options1);
tplot=tSol;
yplot=ySol;
for i=1:20
y2=[yplot(end,1);yplot(end,2);yplot(end,3);yplot(end,4);yplot(end,5);yplot(end,6:105)';yplot(end,106:205)';yplot(end,206:305)';yplot(end,306:405)';...
yplot(end,406:505)';yplot(end,506:605)';yplot(end,606:705)';yplot(end,706:805)';yplot(end,806:905)';yplot(end,906:1005)';yplot(end,1006)];
[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:105)';yplot(end,106:205)';yplot(end,206:305)';yplot(end,306:405)';...
yplot(end,406:505)';yplot(end,506:605)';yplot(end,606:705)';yplot(end,706:805)';yplot(end,806:905)';yplot(end,906:1005)';yplot(end,1006)];
[tSol,ySol]=ode15s(@MBBR_aerobic,[tplot(end) 24*180],y3,options1);
tplot=[tplot;tSol];
yplot=[yplot;ySol];
end
%% initial conditions
function u0 = oscic(x)
u0 = [0.1;0.1;0.1;0.1;0.1;0.1;0.1;0.1;0.1;0.1;0.1;0.4;0.02;0.02;0.3;0.26;0.001];
end
%% events aerobic phase
function [value,isterminal,direction]=eventsfcn_aerobic(t,y)
z=y(5);
z1=+(z<=1.5);
value=z1;
isterminal = 1; % stop at 0
direction = -1; % decreasing
end
%% events anoxic phase
function [value,isterminal,direction]=eventsfcn_anoxic(t,y)
e=y(5);
e1=+(e>0.5);
value=e1;
isterminal = 1; % stop at 0
direction = -1; % decreasing
end
Walter Roberson
on 1 Sep 2023
e1=+(e>0.5);
That is not differentiatable. Use
value = [0.5 - e, e - 1.5];
MATLAB will attempt to keep both of those non-positive (so <= 0) .
0.5 - e is less than or equal to 0 when e is 0.5 or greater.
e - 1.5 is less than or equal to 0 when e is no more than 1.5.
Therefore the [] I show will confine e to be in the range [0.5, 1.5]
KIPROTICH KOSGEY
on 2 Sep 2023
Edited: Walter Roberson
on 2 Sep 2023
Thank you Walter and Torsten.
I have edited my script and the program is now running.
However, Since I am trying to call a different function each time an event occurs, is it the correct coding to do the below. I am doing this to mimic process controller that regulates aeration in a biological system that switches aeration on and off depending on the concentration of oxygen y(5) in the bulk solution. The difference between the two functions is that KaL is either 0 or 10 in MBBR-aerobic and MBBR-anoxic:
close all; clear; clc;
t1=[0 180*24];
xstart = 0;
xend = 0.012;
nx=100;
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))];
nx=100;
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
options = odeset('RelTol',1e-20,'AbsTol',1e-10,'events',@eventsfcn);
%main
[tSol, ySol]=ode15s(@(t,y) MBBR_aerobic(t,y),t1,y1,options);
tplot=tSol;
yplot=ySol;
for i=1:200
y2=[yplot(end,1);yplot(end,2);yplot(end,3);yplot(end,4);yplot(end,5);yplot(end,6:105)';yplot(end,106:205)';yplot(end,206:305)';yplot(end,306:405)';...
yplot(end,406:505)';yplot(end,506:605)';yplot(end,606:705)';yplot(end,706:805)';yplot(end,806:905)';yplot(end,906:1005)';yplot(end,1006)];
[tSol,ySol]=ode15s(@MBBR_anoxic,[tplot(end) 24*180],y2,options);
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:105)';yplot(end,106:205)';yplot(end,206:305)';yplot(end,306:405)';...
yplot(end,406:505)';yplot(end,506:605)';yplot(end,606:705)';yplot(end,706:805)';yplot(end,806:905)';yplot(end,906:1005)';yplot(end,1006)];
[tSol, ySol]=ode15s(@(t,y) MBBR_aerobic(t,y),[tplot(end) 24*180],y3,options);
tplot=[tplot;tSol];
yplot=[yplot;ySol];
end
%% initial conditions
function u0 = oscic(x)
u0 = [0.1;0.1;0.1;0.1;0.1;0.1;0.1;0.1;0.1;0.1;0.1;0.4;0.02;0.02;0.3;0.26;0.001];
end
%% events aerobic phase
function [value,isterminal,direction]=eventsfcn(t,y)
z=y(5);
value=[0.5 - z; z - 1.5];
isterminal = [1;1] ; % stop at 0
direction = [1;1]; % increasing
end
KIPROTICH KOSGEY
on 4 Sep 2023
Dear Torsten and Walter,
Thank you so much for the several comments and insights into coding with Matlab.
I have improved my code further by adding an integral for UL and U but its giving me the error:
Error using integral Limits of integration must be double or single scalars.
This are the lines of codes:
UL=integral(@(L) muaver,0,L,'ArrayValued',true)+sigma;
U=integral(@(x) muaver,0,x,'ArrayValued',true);
KIPROTICH KOSGEY
on 4 Sep 2023
Oh sorry, I have realised the mistake I made, I used x instead of xend=0.012 in the limits of integration.
sorry again to have bothered you.
Regards
KIPROTICH KOSGEY
on 4 Sep 2023
Edited: Walter Roberson
on 4 Sep 2023
Hi,
This is my current script. I have set 'NoNegative', 1, but still the y values are going to negative, what could be the problem?
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-20,'AbsTol',1e-10,'events',@eventsfcn_aerobic,'NonNegative',1);
options2 = odeset('RelTol',1e-20,'AbsTol',1e-10,'events',@eventsfcn_anoxic,'NonNegative',1);
%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;1;10;10;10;10;10;1;0.4;0.02;0.02;0.3;0.26;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
Torsten
on 4 Sep 2023
Edited: Torsten
on 4 Sep 2023
If your differential equations are such that the solution becomes negative, you cannot prevent this with the NonNegative option. E.g. the solution of the differential equation y' = -1 with y(0) = 1 is y(t) = 1 - t, and nothing can stop y(t) to become negative for t > 1. It will only be helpful to prevent solution variables to become negative because of numerical issues.
This having set, 'NonNegative',1 only sets the NonNegative option for y(1). To set it for certain or all components of the vector of solution variables, you have to list them all instead of only the 1.
KIPROTICH KOSGEY
on 4 Sep 2023
oh thank you Torsten.
my system consists of odes (y(1:5&16) and pdes (y(6:15)). Would it be correct to set 'NonNegative', [1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1] ?
Or how best to take care of the solution variables for pdes in the mesh?
Torsten
on 4 Sep 2023
odeset('NonNegative',1:6+10*nx)
if you want to keep all your solution components positive.
But as I said: I think your equations are such that the solution becomes negative from your setting - NonNegative cannot help you in this case.
Walter Roberson
on 4 Sep 2023
The 'Nonnegative' option is, as @Torsten indicates, only intended to compensate for numeric roundoff. It does not properly reset boundary conditions to ensure that the integration is valid. It does not have the effect of "bouncing off of 0" (like a ball might bounce on a floor.)
Walter Roberson
on 4 Sep 2023
Is there a newer version of your other .m files since https://www.mathworks.com/matlabcentral/answers/2007287-index-exceeds-the-number-of-array-elements-index-must-not-exceed-1#comment_2866821 ?
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
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
Answers (0)
See Also
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!An Error Occurred
Unable to complete the action because of changes made to the page. Reload the page to see its updated state.
Select a Web Site
Choose a web site to get translated content where available and see local events and offers. Based on your location, we recommend that you select: .
You can also select a web site from the following list
How to Get Best Site Performance
Select the China site (in Chinese or English) for best site performance. Other MathWorks country sites are not optimized for visits from your location.
Americas
- América Latina (Español)
- Canada (English)
- United States (English)
Europe
- Belgium (English)
- Denmark (English)
- Deutschland (Deutsch)
- España (Español)
- Finland (English)
- France (Français)
- Ireland (English)
- Italia (Italiano)
- Luxembourg (English)
- Netherlands (English)
- Norway (English)
- Österreich (Deutsch)
- Portugal (English)
- Sweden (English)
- Switzerland
- United Kingdom (English)
Asia Pacific
- Australia (English)
- India (English)
- New Zealand (English)
- 中国
- 日本Japanese (日本語)
- 한국Korean (한국어)