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

10 次查看（过去 30 天）
KIPROTICH KOSGEY
KIPROTICH KOSGEY 2020-5-16
评论： KIPROTICH KOSGEY 2020-5-19
Hi,
Kindly, can someone assist to figure out what mistake i made on the attached code that is generating the error below?
'Warning: Rank deficient, rank = 2, tol =
7.141946e-13'
Thanks in advance

请先登录，再回答此问题。

采纳的回答

Walter Roberson
Walter Roberson 2020-5-16
Your function needed a lot of cleanup, and I had to guess about what was desired so you should review it.
With the current code, when you run, at time 0 you will hit a breakpoint in calculating dSno3dt at n = 94, at the new temp2 partial calculation (the lines were too long to debug... copying and pasting is messy in the editor when it exceeds your screen width.)
The reason for the breakpoint is that the calculation will overflow to -infinity . That -infinity will then mix with other calculations leading to an attempt to calculate infinity - infinity which gives nan. Therefore the later fval all return as NaN, leading to an integration failure at time 0.
Please have a careful look at the changes I made. I had to index a number of arrays, and you need to review whether I indexed properly for your calculation purposes. You should also double-check that I got the sign right on adding / subtracting the sub-expressions.
When you work on the code, I recommend that you break the remaining code even further into sub-expressions, to make it easier to debug. Chasing down which operator is causing a problem in a long line of operations is a pain.
By the way:
%% conditions
if le (t,46)
Your functions are discontinuous in time. That is not permitted for the ode* solvers: the mathematics of the Runge-Kutta solvers is only valid for continuous functions. You must terminate the ode15s call at each one of those time boundaries and re-start integration.
function fval=UASBFun6PDE(t,y)
dbstop if naninf
%% Define the variables
Xaob=y(1);
Xnb=y(2);
Xnsp=y(3);
Xcmx=y(4);
Xamx=y(5);
Xhan=y(6);
Xs=y(7);
Sse=y(8);
Sno3=y(9);
Sno2=y(10);
Snh4=y(11);
So2=y(12);
%% parameters
Xo=0; Xe=Xo; Xso=0;
%aob
O1=68/(0.00831*293*303); u1=1.443*exp(O1*15); Ko2aob=0.85*0.594; Knh4aob=0.55*2.4; Yaob=3*0.15; Kno3aob=0.5;
%nb
O2=44/(0.00831*293*303); u2=1.8*0.48*exp(O2*15); Ko2nb=0.98*0.435; Kno2nb=0.6*1.375; Ynb=1.1*0.13;Kno3nb=0.5;
%nsp
O3=44/(0.00831*293*303); u3=1*0.69*exp(O3*15); Ko2nsp=0.95*0.435; Kno2nsp=1.0*0.76; Ynsp=1.0*0.14;Kno3nsp=0.5;
%cmx
O4=44/(0.00831*293*303); u4=0.5*0.72*exp(O4*15); Ko2cmx=0.98*0.43; Kno2cmx=0.15*6.29;Knh4cmx=70*0.011; Ycmx=0.1*0.47;YcmxNo2=2*Ycmx;Kno3cmx=0.5;
%amx
O5=71/(0.00831*293*303); u5=1.2*0.0664*exp(O5*15); Kno2amx=1.2*0.175; Kno3amx=0.5; Knh4amx=2*0.185; Yamx=0.95*0.159; Ko2amx=0.01; Kno3amx=0.5;
%han
u6=7.2*exp(15*0.069); KSsehan=2; Kno2han=0.5; YNo2han=1.8*0.49; Ko2han=0.2;
YNo3han=0.79; Kno3han=0.32;Yhaer=1.8*0.37;
%other constants
baob=0.1296*exp(0.11*15); bnb=0.069*exp(0.11*15);bnsp=0.06*exp(0.11*15);
bcmx=0.06*exp(0.11*15); bamx=0.00312*exp(0.11*15); bhaer=0.192*exp(15*0.11);
bhan=0.192*exp(0.11*15);
INBM=0.0583; fI=0.08; namx=0.5; nhan=0.6;naob=0.5; nnb=0.5;nhaer=0.5;nnsp=0.5;ncmx=0.5;
KX=0.03; KH=3*exp(0.069*15); INXI=0.02;
% input parameters
Sseo=2; Sno3o=0;
%oxygen transfer
KaL=222*(1.02^(-5));
%% conditions
if le (t,46)
Snh4o=50; Sno2o=55; So2G=1;So2o=2; V=0.0038/5; Qo=0.000864; Qe=Qo;
elseif and (gt (t,46),le (t,65))
Snh4o=60.4; Sno2o=25.3; So2G=1;So2o=8.5; V=0.0038/5; Qo=0.000864; Qe=Qo;
elseif and (gt (t,65) , le (t,67))
Snh4o=60.4; Sno2o=25.3; So2G=1; So2o=8.5; V=0.0038/5; Qo=0.000864; Qe=Qo;
elseif and (gt (t,67) , le (t,74))
Snh4o=76.72; Sno2o=45.25; So2G=1;So2o=8.5; V=0.0038/5; Qo=0.002016; Qe=Qo;
elseif and (gt (t,74) , le (t,79))
Snh4o=76.72; Sno2o=45.25; So2G=1; So2o=8.5; V=0.0038/5; Qo=0.002592; Qe=Qo;
elseif and (gt (t,80) , le (t,86))
Snh4o=76.72; Sno2o=45.25; So2o=8.5; V=0.005/5; Qo=0.002592; Qe=Qo; So2G=6.7521; %So2G=14.65-0.41*35+7.99*(10^-3)*(35^2)-7.78*(10^-5)*35^3;
elseif and (gt (t,86) , le (t,87))
Snh4o=76.72; Sno2o=45.25; So2G=0.1;So2o=0.2; V=0.005/5; Qo= 0.002736; Qe=Qo;
elseif and (gt (t,87) , le (t,102))
Snh4o=84.9; Sno2o=50.9; So2G=0.1;So2o=0.2; V=0.005/5; Qo=0.002736; Qe=Qo;
elseif and (gt (t,102) , le (t,108))
Snh4o=76.72; Sno2o=45.25; So2G=0.1;So2o=0.2; V=0.005/5; Qo=0.003168; Qe=Qo;
elseif and (gt (t,108) , le (t,120))
Snh4o=76.72; Sno2o=45.25; So2G=0.1;So2o=0.2; V=0.005/5; Qo=0.004176; Qe=Qo;
elseif and (gt (t,121) , le (t,134))
Snh4o=75.1; Sno2o=62.8; So2G=0.1;So2o=0.2; V=0.005/5; Qo=0.0054; Qe=Qo;
elseif and (gt (t,134) , le (t,137))
Snh4o=61.2; Sno2o=45.7; So2G=0.1;So2o=0.2; V=0.005/5; Qo=0.0054; Qe=Qo;
elseif and (gt (t,137) , le (t,151))
Snh4o=61.2; Sno2o=45.7; So2G=0.1;So2o=0.2; V=0.005/5; Qo=0.006120; Qe=Qo;
elseif and (gt (t,151) , le (t,159))
Snh4o=61.2; Sno2o=45.7; So2G=0.1;So2o=0.2; V=0.005/5; Qo=0.007056; Qe=Qo;
elseif and (gt (t,159) , le (t,176))
Snh4o=68.7; Sno2o=49.7; So2G=0.1;So2o=0.2; V=0.005/5; Qo=0.007920; Qe=Qo;
elseif and (gt (t,176) , le (t,177))
Snh4o=68.7; Sno2o=49.7; So2G=0.1;So2o=0.2; V=0.005/5; Qo=0.006120; Qe=Qo;
elseif and (gt (t,177) , le (t,180))
Snh4o=68.7; Sno2o=49.7; So2G=0.1;So2o=0.2; V=0.005/5; Qo=0.006120; Qe=Qo;
elseif and (gt (t,180) , le (t,197))
Snh4o=50.4; Sno2o=42; So2G=0.1;So2o=0.2; V=0.0038/5; Qo=0.006120; Qe=Qo;
elseif and (gt (t,197) , le (t,208))
Snh4o=59.3; Sno2o=58.2; So2G=0.1;So2o=0.2; V=0.0038/5; Qo=0.006120; Qe=Qo;
elseif and (gt (t,208) , le (t,212))
Snh4o=59.3; Sno2o=58.2; So2G=0.1;So2o=0.2; V=0.0038/5; Qo=0.007920; Qe=Qo;
elseif and (gt (t,212) , le (t,214))
Snh4o=59.3; Sno2o=58.2; So2G=0.1;So2o=0.2; V=0.0038/5; Qo=0.007920; Qe=Qo;
elseif and (gt (t,214) , le (t,221))
Snh4o=65; Sno2o=75; So2G=0.1;So2o=0.2; V=0.0038/5; Qo=0.007920; Qe=Qo;
elseif and (gt (t,221) , le (t,225))
Snh4o=65; Sno2o=75; So2G=0.1;So2o=0.2; V=0.0038/5; Qo=0.007056; Qe=Qo;
elseif and (gt (t,225) , le (t,228))
Snh4o=65; Sno2o=75; So2G=0.1;So2o=0.2; V=0.0038/5; Qo=0.007920; Qe=Qo;
elseif and (gt (t,228) , le (t,237))
Snh4o=72.8; Sno2o=89.6; So2G=0.1;So2o=0.2; V=0.0038/5; Qo=0.007920; Qe=Qo;
elseif and (gt (t,237) , le (t,245))
Snh4o=99.7; Sno2o=116.3; So2G=0.1;So2o=0.2; V=0.0038/5; Qo=0.007920; Qe=Qo;
elseif and (gt (t,245) , lt (t,250))
Snh4o=99.7; Sno2o=116.3; So2G=0.1;So2o=0.2; V=0.0038/5; Qo=0.007920; Qe=Qo;
elseif eq(t,250)
Snh4o=99.7; Sno2o=116.3; So2G=0.1;So2o=0.2; V=0.0038/5; Qo=0.007920; Qe=Qo;
elseif and (gt (t,250) , le (t,254))
Snh4o=112.19; Sno2o=145.38; So2G=0.1;So2o=0.2; V=0.0038/5; Qo=0.01872; Qe=Qo;
elseif and (gt (t,254) , le (t,256))
Snh4o=112.19; Sno2o=145.38; So2G=0.1;So2o=0.2; V=0.0038/5; Qo=0.007920; Qe=Qo;
elseif and (gt (t,256) , le (t,264))
Snh4o=112.19; Sno2o=145.38; So2G=0.1;So2o=0.2; V=0.0038/5; Qo=0.006588; Qe=Qo;
elseif and (gt (t,264) , le (t,278))
Snh4o=112.19; Sno2o=145.38; So2G=0.1;So2o=0.2; V=0.0038/5; Qo=0.007056; Qe=Qo;
elseif and (gt (t,278) , le (t,288))
Snh4o=108.43; Sno2o=140.786; So2G=0.1;So2o=0.2; V=0.0038/5; Qo=0.007056; Qe=Qo;
elseif and (gt (t,288) , le (t,295))
Snh4o=108.43; Sno2o=140.786; So2G=0.1;So2o=0.2; V=0.005/5; Qo=0.007920; Qe=Qo;
elseif and (gt (t,295) , le (t,302))
Snh4o=140.43; Sno2o=196.4; So2G=0.1;So2o=0.2; V=0.005/5; Qo=0.007920; Qe=Qo;
elseif and (gt (t,302) , le (t,311))
Snh4o=160; Sno2o=223.85; So2G=0.1;So2o=0.2; V=0.005/5; Qo=0.007920; Qe=Qo;
elseif and (gt (t,311) , le (t,313))
Snh4o=180.64; Sno2o=252.08; So2G=0.1;So2o=0.2; V=0.005/5; Qo=0.007920; Qe=Qo;
elseif and (gt (t,313) , le (t,322))
Snh4o=180.64; Sno2o=252.08; So2G=0.1;So2o=0.2; V=0.005/5; Qo=0.01; Qe=Qo;
elseif and (gt (t,322) , le (t,327))
Snh4o=200; Sno2o=280; So2G=0.1;So2o=0.2; V=0.005/5; Qo=0.01; Qe=Qo;
elseif and (gt (t,327) , le (t,355))
Snh4o=220; Sno2o=308; So2G=0.1;So2o=0.2; V=0.005/5; Qo=0.01; Qe=Qo;
elseif and (gt (t,355) , le (t,425))
Snh4o=200; Sno2o=264; So2G=0.1;So2o=0.2; V=0.005/5; Qo=0.01; Qe=Qo;
elseif and (gt (t,425) , le (t,434))
Snh4o=110; Sno2o=145; So2G=0.1;So2o=0.2; V=0.005/5; Qo=0.0166667; Qe=Qo;
elseif and (gt (t,434) , le (t,475))
Snh4o=200; Sno2o=264; So2G=0.1;So2o=0.2; V=0.005/5; Qo=0.005; Qe=Qo;
elseif and (gt (t,475) , le (t,491))
Snh4o=200; Sno2o=264; So2G=0.1;So2o=0.2; V=0.005/5; Qo=0.00667; Qe=Qo;
else
Snh4o=200; Sno2o=264; So2G=0.1;So2o=0.2; V=0.005/5; Qo=0.01; Qe=Qo;
end
%% ODE
fval(1)=((0-V*Xaob/250)/V + u1*(So2/(Ko2aob+So2))*(Snh4/(Knh4aob+Snh4))*Xaob - baob*(So2/(Ko2aob+So2))*Xaob - baob*naob*((Ko2aob/(Ko2aob+So2))*(Sno2+Sno3)/(Kno3aob+Sno2+Sno3))*Xaob);
fval(2)=((0-V*Xnb/250)/V + u2*(Sno2/(Kno2nb+Sno2))*(So2/(Ko2nb+So2))*Xnb - bnb*(So2/(Ko2nb+So2))*Xnb - bnb*nnb*((Ko2nb/(Ko2nb+So2))*(Sno2+Sno3)/(Kno3nb+Sno2+Sno3))*Xnb);
fval(3)=((0-V*Xnsp/250)/V + u3*(Sno2/(Kno2nsp+Sno2))*(So2/(Ko2nsp+So2))*Xnsp - bnsp*(So2/(Ko2nsp+So2))*Xnsp - bnsp*nnsp*((Ko2nsp/(Ko2nsp+So2))*(Sno2+Sno3)/(Kno3nsp+Sno2+Sno3))*Xnsp);
fval(4)=((0-V*Xcmx/250)/V + u4*((Sno2/(Kno2cmx+Sno2))*(So2/(Ko2cmx+So2)) + u4*(So2/(Ko2cmx+So2))*(Snh4/(Knh4cmx+Snh4)))*Xcmx - bcmx*(So2/(Ko2cmx+So2))*Xcmx - bcmx*ncmx*((Ko2cmx/(Ko2cmx+So2))*(Sno2+Sno3)/(Kno3cmx+Sno2+Sno3))*Xcmx);
fval(5)=((0-V*Xamx/250)/V + u5*((Snh4/(Knh4amx+Snh4))*(Sno2/(Kno2amx+Sno2))*(Ko2amx/(Ko2amx+So2)))*Xamx-bamx*(So2/(Ko2amx+So2))*Xamx - bamx*namx*((Ko2amx/(Ko2amx+So2))*(Sno2+Sno3)/(Kno3amx+Sno2+Sno3))*Xamx);
fval(6)=((0-V*Xhan/250)/V + (u6*(Sse/(KSsehan+Sse))*(Sno2/(Kno2han+Sno2))*(Ko2han/(Ko2han+So2)) + u6*(Sse/(KSsehan+Sse))*(Sno3/(Kno3han+Sno3))*(Ko2han/(Ko2han+So2))+u6*((Sse/(KSsehan+Sse))*(So2/(Ko2han+So2)))*Xhan)- bhan*(So2/(Ko2han+So2))*Xhan - bhan*nhan*((Ko2han/(Ko2han+So2))*(Sno2+Sno3)/(Kno3han+Sno2+Sno3))*Xhan);
fval(7)=((0-V*Xs/250)/V - KH*((Xs/Xhan)/(KX+(Xs/Xhan)))*Xhan + (1-fI)*(baob*Xaob+bnb*Xnb+bnsp*Xnsp+bcmx*Xcmx+bamx*Xamx+bhan*Xhan));
%% PDE
L=0.31;
N=535/1;
J=5;
hz=L/J;
for j=1:J+1; Sse(1,j)=0.0005; end
for n=1:N+1; Sse(n,1)=Sseo; end
for j=1:J+1; Sno3(1,j)=0.0005; end
for n=1:N+1; Sno3(n,1)=Sno3o; end
for j=1:J+1; Sno2(1,j)=55; end
for n=1:N+1; Sno2(n,1)=Sno2o; end
for j=1:J+1; Snh4(1,j)=50; end
for n=1:N+1; Snh4(n,1)=Snh4o; end
for j=1:J+1; So2(1,j)=0.5; end
for n=1:N+1; So2(n,1)=So2o; end
dSsedt=zeros(N+1,J+1);
dSno3dt=zeros(N+1,J+1);
dSno2dt=zeros(N+1,J+1);
dSnh4dt=zeros(N+1,J+1);
dSo2dt=zeros(N+1,J+1);
Dma=3;
U=1;
for n=1:N
for j=2:J
temp1 = -U*(dSsedt(n+1,j)-dSsedt(n,j));
temp2 = Dma*(dSsedt(n,j+1)-2*dSsedt(n,j)+(dSsedt(n,j-1)))/hz^2;
temp3 = (u6*nhan*(1/YNo2han)*(Sse(n,J)./(KSsehan+Sse(n,j)))*(Sno2(n,j)./(Kno2han+Sno2(n,j)))*(Ko2han/(Ko2han+So2(n,j))))*Xhan;
temp4 = (u6*nhan*(1/YNo3han)*(Sse(n,j)/(KSsehan+Sse(n,j)))*(Sno3(n,j)/(Kno3han+Sno3(n,j)))*(Ko2han/(Ko2han+So2(n,j))))*Xhan;
dSsedt(n+1,j)= temp1 + temp2 - temp3 - temp4 - ((u6/Yhaer)*(Sse(n,j)/(KSsehan+Sse(n,j)))*(So2(n,j)/(Ko2han+So2(n,j)))*Xhan) + KH*(Xs/(Xhan))/(KX+(Xs/(Xhan)))*(Xhan)+0.0*(Xaob+Xnb+Xnsp+Xcmx+Xamx+Xhan);
temp1 = -U*(dSno3dt(n+1,j)-dSno3dt(n,j));
temp2 = Dma*(dSno3dt(n,j+1)-2*dSno3dt(n,j)+dSno3dt(n,j-1))/hz^2;
temp3 = (u6*nhan*((1-YNo3han)/(2.86*YNo3han))*(Sse(n,j)/(KSsehan+Sse(n,j)))*(Sno3(n,j)/(Kno3han+Sno3(n,j)))*(Ko2han/(Ko2han+So2(n,j))))*Xhan;
temp4 = (u5/1.14)*(Snh4(n,j)/(Knh4amx+Snh4(n,j)))*(Sno2(n,j)/(Kno2amx+Sno2(n,j)))*(Ko2amx/(Ko2amx+So2(n,j)))*Xamx;
dSno3dt(n+1,j) = temp1 + temp2 - temp3 + temp4 + (u2/Ynb*(Sno2(n,j)/(Kno2nb+Sno2(n,j)))*(So2(n,j)./(Ko2nb+So2(n,j)))*Xnb) + ((u3/Ynsp)*(Sno2(n,j)/(Kno2nsp+Sno2(n,j)))*(So2(n,j)./(Ko2nsp+So2(n,j)))*Xnsp) + ((u4/YcmxNo2)*(Sno2(n,j)/(Kno2cmx+Sno2(n,j)))*(So2(n,j)/(Ko2cmx+So2(n,j)))*Xcmx)+((u4/Ycmx)*(Snh4(n,j)/(Knh4cmx+Snh4(n,j)))*(So2(n,j)./(Ko2cmx+So2(n,j)))*Xcmx) - ((1-fI)/2.86)*baob*naob*(Ko2aob/(Ko2aob+So2(n,j)))*(Sno2(n,j)+Sno3(n,j))/(Kno3aob+Sno2(n,j)+Sno3(n,j))*Xaob - ((1-fI)/2.86)*bnb*nnb*(Ko2nb/(Ko2nb+So2(n,j)))*(Sno2(n,j)+Sno3(n,j))/(Kno3nb+Sno2(n,j)+Sno3(n,j))*Xnb - ((1-fI)/2.86)*bnsp*nnsp*(Ko2nsp/(Ko2nsp+So2(n,j)))*(Sno2(n,j)+Sno3(n,j))/(Kno3nsp+Sno2(n,j)+Sno3(n,j))*Xnsp - ((1-fI)/2.86)*bcmx*ncmx*(Ko2cmx./(Ko2cmx+So2(n,j)))*(Sno2(n,j)+Sno3(n,j))/(Kno3cmx+Sno2(n,j)+Sno3(n,j))*Xcmx - ((1-fI)/2.86)*bamx*namx*(Ko2amx/(Ko2amx+So2(n,j)))*(Sno2(n,j)+Sno3(n,j))/(Kno3amx+Sno2(n,j)+Sno3(n,j))*Xamx-((1-fI)/2.86)*bhan*nhan*(Ko2han/(Ko2han+So2(n,j)))*(Sno2(n,j)+Sno3(n,j))/(Kno3han+Sno2(n,j)+Sno3(n,j))*Xhan;
temp1 = -U*(dSno2dt(n+1,j)-dSno2dt(n,j));
temp2 = Dma*(dSno2dt(n,j+1)-2*dSno2dt(n,j)+dSno2dt(n,j-1))/hz^2;
temp3 = (((1-YNo2han)/(1.71*YNo2han))*u6*nhan*(Sse(n,j)/(KSsehan+Sse(n,j)))*(Sno2(n,j)/(Kno2han+Sno2(n,j)))*(Ko2han/(Ko2han+So2(n,j)))*Xhan);
temp4 = (((u5/Yamx)+u5/1.14)*(Snh4(n,j)/(Knh4amx+Snh4(n,j)))*(Sno2(n,j)/(Kno2amx+Sno2(n,j)))*(Ko2amx/(Ko2amx+So2(n,j))))*Xamx;
dSno2dt(n+1,j) = temp1 + temp2 - temp3 - temp4 + ((u1/Yaob)*(So2(n,j)/(Ko2aob+So2(n,j)))*(Snh4(n,j)/(Knh4aob+Snh4(n,j))))*Xaob - ((u2/Ynb)*(Sno2(n,j)/(Kno2nb+Sno2(n,j)))*(So2(n,j)/(Ko2nb+So2(n,j)))*Xnb) - ((u3/Ynsp)*(Sno2(n,j)/(Kno2nb+Sno2(n,j)))*(So2(n,j)/(Ko2nb+So2(n,j))))*Xnsp - ((u4/YcmxNo2)*(Sno2(n,j)/(Kno2nb+Sno2(n,j)))*(So2(n,j)/(Ko2nb+So2(n,j))))*Xcmx;
temp1 = -U*(dSnh4dt(n+1,j)-dSnh4dt(n,j));
temp2 = Dma*(dSnh4dt(n,j+1)-2*dSnh4dt(n,j)+dSnh4dt(n,j-1))/hz^2;
temp3 = (u1*(INBM+1/Yaob)*(So2(n,j)/(Ko2aob+So2(n,j)))*(Snh4(n,j)/(Knh4aob+Snh4(n,j))))*Xaob;
temp4 = (u4*(INBM+1/Ycmx)*(So2(n,j)/(Ko2cmx+So2(n,j)))*(Snh4(n,j)/(Knh4cmx+Snh4(n,j))))*Xcmx;
dSnh4dt(n+1,j) = temp1 + temp2 - temp3 - temp4 -(u5*(INBM+1/Yamx)*(Snh4(n,j)/(Knh4amx+Snh4(n,j)))*(Sno2(n,j)/(Kno2amx+Sno2(n,j)))*(Ko2amx/(Ko2amx+So2(n,j))))*Xamx - (u2*INBM*(Sno2(n,j)/(Kno2nb+Sno2(n,j)))*(So2(n,j)/(Ko2nb+So2(n,j))))*Xnb -(u3*INBM*(Sno2(n,j)/(Kno2nsp+Sno2(n,j)))*(So2(n,j)/(Ko2nsp+So2(n,j))))*Xnsp - INBM*u6*nhan*((Sse(n,j)/(KSsehan+Sse(n,j)))*(Sno2(n,j)/(Kno2han+Sno2(n,j)))*(Ko2han/(Ko2han+So2(n,j))))*Xhan - u6*INBM*nhan*((Sse(n,j)/(KSsehan+Sse(n,j)))*(Sno3(n,j)/(Kno3han+Sno3(n,j)))*(Ko2han/(Ko2han+So2(n,j))))*Xhan - u6*INBM*((Sse(n,j)/(KSsehan+Sse(n,j)))*(So2(n,j)/(Ko2han+So2(n,j))))*Xhan + (INBM-(fI*INXI))*baob*naob*((Ko2aob/(Ko2aob+So2(n,j)))*(Sno2(n,j)+Sno3(n,j))/(Kno3aob+Sno2(n,j)+Sno3(n,j)))*Xaob + (INBM-(fI*INXI))*bnb*nnb*((Ko2nb/(Ko2nb+So2(n,j)))*(Sno2(n,j)+Sno3(n,j))/(Kno3nsp+Sno2(n,j)+Sno3(n,j)))*Xnb +(INBM-(fI*INXI))*bnsp*nnsp*((Ko2nsp/(Ko2nsp+So2(n,j)))*(Sno2(n,j)+Sno3(n,j))/(Kno3nsp+Sno2(n,j)+Sno3(n,j)))*Xnsp + (INBM-(fI*INXI))*(bcmx*ncmx*(Ko2cmx/(Ko2cmx+So2(n,j)))*(Sno2(n,j)+Sno3(n,j))/(Kno3cmx+Sno2(n,j)+Sno3(n,j)))*Xcmx +(INBM-(fI*INXI))*bamx*namx*((Ko2amx/(Ko2amx+So2(n,j)))*(Sno2(n,j)+Sno3(n,j))/(Kno3amx+Sno2(n,j)+Sno3(n,j)))*Xamx + (INBM-(fI*INXI))*bhan*nhan*((Ko2han/(Ko2han+So2(n,j)))*(Sno2(n,j)+Sno3(n,j))/(Kno3han+Sno2(n,j)+Sno3(n,j)))*Xhan +(INBM-(fI*INXI))*baob*(So2(n,j)/(Ko2aob+So2(n,j)))*Xaob + (INBM-(fI*INXI))*bnb*(So2(n,j)/(Ko2nb+So2(n,j)))*Xnb +(INBM-(fI*INXI))*bnsp*(So2(n,j)/(Ko2nsp+So2(n,j)))*Xnsp + (INBM-(fI*INXI))*bcmx*(So2(n,j)/(Ko2cmx+So2(n,j)))*Xcmx + (INBM-(fI*INXI))*bamx*(So2(n,j)/(Ko2amx+So2(n,j)))*Xamx + (INBM-(fI*INXI))*bhan*(So2(n,j)/(Ko2han+So2(n,j)))*Xhan;
temp1 = -U*(dSo2dt(n+1,j)-dSo2dt(n,j));
temp2 = Dma*(dSo2dt(n,j+1)-2*dSo2dt(n,j)+dSo2dt(n,j-1))/hz^2;
temp3 = (KaL*(So2G-So2(n,j)))-(((1-Yhaer)/Yhaer)*u6*(Sse(n,j)/(KSsehan+Sse(n,j)))*(So2(n,j)/(Ko2han+So2(n,j))))*Xhan;
temp4 = (((3.43-Yaob)/Yaob)*u1*(So2(n,j)/(Ko2aob+So2(n,j)))*(Snh4(n,j)/(Knh4aob+Snh4(n,j))))*Xaob;
dSo2dt(n+1,j) = temp1 + temp2 + temp3 - temp4 - (((1.14-Ynb)/Ynb)*u2*(Sno2(n,j)/(Kno2nb+Sno2(n,j)))*(So2(n,j)/(Ko2nb+So2(n,j))))*Xnb - (((1.14-Ynsp)/Ynsp)*u3*(Sno2(n,j)/(Kno2nsp+Sno2(n,j)))*(So2(n,j)/(Ko2nsp+So2(n,j))))*Xnsp - (((4.57-Ycmx)/Ycmx)*u4*(So2(n,j)/(Ko2cmx+So2(n,j)))*(Snh4(n,j)/(Knh4cmx+Snh4(n,j))))*Xcmx -(((1.14-YcmxNo2)/YcmxNo2)*u4*(Sno2(n,j)/(Kno2cmx+Sno2(n,j)))*(So2(n,j)/(Ko2cmx+So2(n,j))))*Xcmx- (1-fI)*baob*(So2(n,j)/(Ko2aob+So2(n,j)))*Xaob - (1-fI)*bnb*(So2(n,j)/(Ko2nb+So2(n,j)))*Xnb -(1-fI)*bnsp*(So2(n,j)/(Ko2nsp+So2(n,j)))*Xnsp - (1-fI)*bcmx*(So2(n,j)/(Ko2cmx+So2(n,j)))*Xcmx - (1-fI)*bamx*(So2(n,j)/(Ko2amx+So2(n,j)))*Xamx - (1-fI)*bhan*(So2(n,j)/(Ko2han+So2(n,j)))*Xhan;
end
end
fval(8)=dSsedt(n+1,j);
fval(9)=dSno3dt(n+1,j);
fval(10)=dSno2dt(n+1,j);
fval(11)=dSnh4dt(n+1,j);
fval(12)=dSo2dt(n+1,j);
fval = fval(:);
dbclear if naninf
end
  5 个评论
显示 3更早的评论
Walter Roberson
Walter Roberson 2020-5-19
Do not use events for this.
y0=[0.5*18.16;0.5*0.52;0.5*25.64;0.5*3.6;0.5*172.6;0.5*179.48;0.0005;2;0.0005;55;50;0.5];
tspans = [0, 46, 65, 67, 74, 80, 86, 87, 102, 108, 121, 134, 137, 151, 159, 176, 177, 180, 197, 208, 212, 214, 221, 225, 228, 237, 245, 250, 250.1, 254, 256, 264, 278, 288, 295, 302, 311, 313, 322, 327, 355, 425, 434, 475, 491, 535];
nspans = length(tspans) - 1;
tout = cell(nspans,1);
yout = cell(nspans,1);
hold off
for idx = 1 : nspans
% Solve until the first terminal event.
[t,y] = ode15s(@(t,y) MBBRFun5(t,y), linspace(tspans(idx),tspans(idx+1),25), y0);
% Accumulate output. This could be passed out as output arguments.
if idx == 1
tout{idx} = t;
yout{idx} = y;
else
tout{idx} = t(2:end);
yout{idx} = y(2:end,:);
end
for col = 1 : 12
ax = subplot(3,4,col);
plot(ax, t, y(:,col));
title(ax, sprintf('y(:,%d)', col) );
hold(ax, 'on');
end
drawnow
%set new initial conditions
y0 = y(end,:);
end
all_t = vertcat(tout{:});
all_y = vertcat(yout{:});
The slowest part is the plotting. If you were to postpone that to the end, then it would be faster.
KIPROTICH KOSGEY
KIPROTICH KOSGEY 2020-5-19
Thank you soo much Walter for your efforts and time.
I realise that this codes are giving exactly the same results as i had before stopping at the discontinuities. However, my codes were difficult to execute because of high stiffness so i kept altering the tolerances. My previous script is attached.
As for the coupled PDEs and ODEs, i bet i will have to dig a lot deep in future.
I cannot thank you enough!

请先登录，再进行评论。

更多回答（0 个）

类别

在 Help Center 和 File Exchange 中查找有关 Partial Differential Equation Toolbox 的更多信息

标签

Community Treasure Hunt

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

Start Hunting!

Translated by