Index exceeds the number of array elements. Index must not exceed 1.
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
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
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
Torsten
on 11 Aug 2023
Loop from ix = 2:nx-1 and set the boundary conditions separately.
KIPROTICH KOSGEY
on 11 Aug 2023
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
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
KIPROTICH KOSGEY
on 16 Aug 2023
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
KIPROTICH KOSGEY
on 16 Aug 2023
Edited: Walter Roberson
on 16 Aug 2023
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);
%% 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
KIPROTICH KOSGEY
on 25 Aug 2023
Edited: Torsten
on 25 Aug 2023
Torsten
on 25 Aug 2023
I can't make the two .m files work.
KIPROTICH KOSGEY
on 25 Aug 2023
Edited: Torsten
on 25 Aug 2023
KIPROTICH KOSGEY
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
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
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
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
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
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
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
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
KIPROTICH KOSGEY
on 1 Sep 2023
fahad isah adam
on 1 Sep 2023
Moved: John D'Errico
on 1 Sep 2023
Kudos to you, Sir.
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
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
KIPROTICH KOSGEY
on 4 Sep 2023
KIPROTICH KOSGEY
on 4 Sep 2023
KIPROTICH KOSGEY
on 4 Sep 2023
Edited: Walter Roberson
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
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
Answers (0)
Categories
Find more on Eigenvalue Problems in Help Center and File Exchange
on 10 Aug 2023
on 5 Sep 2023
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
