MATLAB: Error code drying of solid

solutions please

when running the program it reflects the error to me

code

%% PARAMETROS DE ENTRADA AL MODELO
a=0.2755; b=-0.34; D0=0.008 ;r0=D0/2; as=3/r0;
%Cálculo del flujo másico de aire
l1=0.53;
l2=0.51;
A=l1*l2;
L1=100/100;
V=A*L1;
v=4; %Velocidad del aire
ma=0.028;
Ga=(ma/A)*3600;
N=4;
L=0.5;
x=L/N;
h=0.5;
tfinal=24*4.5;
Mi=0.2409
v=1;
Sp=0.685*0.51;
%Coeficientes de transferencia de calor por conducción
lamb_madera=0.13;
lamb_vidrio=0.05;
lamb_yeso=0.35;
em=2.5/100;
ev=6/100;
ey=1.5/100;
hcp=(lamb_madera/em)+(lamb_vidrio/ev)+(lamb_yeso/ey);
hvv=(5.67+3.86*v);
hvfpi=4.1;
%Coeficiente de Fourier
a0=29.86 
a11=-2.997 
b1 =-6.595 
a22 =0.2317 
b2 =1.949 
a33 =-0.2656 
b3 =-0.4193 
a4 =-0.3443 
b4 =-0.3929 
a5 =0.2462 
b5 =0.3422 
a6 =-0.1787 
b6 =-0.4957 
a7 =-0.3067 
b7 =0.2473 
w =0.1488
%Temperatura de ambiente
a0a=25.91 
a1a=-2.422 
b1a =-3.467 
a2a =0.4955 
b2a =1.07 
a3a =-0.4853 
b3a =0.4878 
a4a =0.1001 
b4a =0.07658 
a5a =0.2295 
b5a =-0.04174 
a6a =-0.0225 
b6a =0.02728 
a7a =-0.07581 
b7a =-0.1319 
a8a =-0.1274 
b8a =0.1078 
wa =0.1298
%% DETERMINACION DEL NUMERO DE NODOS PARA LA SIMULACION
N=L/x;
N=ceil(N);
X=tfinal/h;
X=ceil(X);
%% CONDICIONES INICIALES 
M0=56;
M0=100*M0/(100-M0);
for n=2:N
 M(n,13)=M0;
 Tg(n,13)=20;
end
% load('Tf0.mat');
for i=14:X
 H(1,i)=0.01;
end
for ia=14:(X)
 
T(1,ia)=a0+a11*cos(ia*w)+b1*sin(ia*w)+a22*cos(2*ia*w)+b2*sin(2*ia*w)+a33*cos(3*ia*w)+b3*sin(3*ia*w)+a4*cos(4*ia*w)+b4*sin(4*ia*w)+a5*cos(5*ia*w)+b5*sin(5*ia*w)+a6*cos(6*ia*w)+b6*sin(6*ia*w)+a7*cos(7*ia*w)+b7*sin(7*ia*w);
Ta(1,ia)=a0a+a1a*cos(ia*wa)+b1a*sin(ia*wa)+a2a*cos(2*ia*wa)+b2a*sin(2*ia*wa)+a3a*cos(3*ia*wa)+b3a*sin(3*ia*wa)+a4a*cos(4*ia*wa)+b4a*sin(4*ia*wa)+a5a*cos(5*ia*wa)+b5a*sin(5*ia*wa)+a6a*cos(6*ia*wa)+b6a*sin(6*ia*wa)+a7a*cos(7*ia*wa)+b7a*sin(7*ia*wa)+a8a*cos(8*ia*wa)+b8a*sin(8*ia*wa);
 Tsky(1,ia)=Ta(1,ia)-6;
end
%% CALCULO DE LOS PARAMETROS INTERMEDIOS y H, M,T,Tg
for j=14:X
 for n=2:N
     
 end
 %Calculo de M(n,j)
 RH(n-1,j)=1.626*H(n-1,j)*(10^5)/((1+1.608*H(n-1,j))*6892.4*exp(54.633-12301.69/(1.8*T(n-1,j)+492)-5.169*log(1.8*T(n-1,j)+492)));
 RH(n-1,j)=abs(RH(n-1,j));
 Me(n-1,j)=100*(0.01087324+0.23758842*RH(n-1,j)-0.63848052*(RH(n-1,j)^2)+0.56197989*(RH(n-1,j)^3))*exp((0.02965424-0.33506858*RH(n-1,j)+0.95779665*(RH(n-1,j)^2)-1.27182954*(RH(n-1,j)^3)+0.573953345*(RH(n-1,j)^4))*(T(n-1,j)-57.2913401));
 Me(n-1,j)=abs(Me(n-1,j));
 K(n-1,j)=(pi*pi/(9))*as*as*3600*8.64*(10^-4)*exp(-39.94*1000/(8.314*(Tg(n,j-1)+273)));
 K(n-1,j)=abs(K(n-1,j));
 
 k1(n-1,j)=-K(n-1,j)*(M(n,j-1)-Me(n-1,j));
 k2(n-1,j)=-K(n-1,j)*((M(n,j-1)+0.05*k1(n-1,j))-Me(n-1,j));
 k3(n-1,j)=-K(n-1,j)*((M(n,j-1)+0.05*k2(n-1,j))-Me(n-1,j));
 k4(n-1,j)=-K(n-1,j)*((M(n,j-1)+0.1*k3(n-1,j))-Me(n-1,j));
 M(n,j)=M(n,j-1)+0.0167*(k1(n-1,j)+2*k2(n-1,j)+2*k3(n-1,j)+k4(n-1,j));
 dM(n,j)=0.0167*(k1(n-1,j)+2*k2(n-1,j)+2*k3(n-1,j)+k4(n-1,j));
 %Calculo de H(n,j)
 rhop(n,j-1)=726.4-4.135*-4.135*(M(n,j-1)/(100+M(n,j-1)));
 rhop(n,j-1)=abs(rhop(n,j-1));
 
 H(n,j)=H(n-1,j)-0.01*(rhop(n,j-1)*x/(Ga*2*h))*(M(n,j)-M(n,j-1));
 %Calculo de T(n,j)
 Ca(n-1,j)=1003.4+0.178*T(n-1,j);
 Ca(n-1,j)=abs(Ca(n-1,j));
 
 Cv(n-1,j)=1859+0.236*T(n-1,j);
 Cv(n-1,j)=abs(Cv(n-1,j));
 
 ua(n-1,j)=(1.691*(10^-5)+4.987*(10^-8)*T(n-1,j)-3.187*(10^-11)*(T(n-1,j)^2)+1.319*(10^-14)*(T(n-1,j)^3));
 ua(n-1,j)=abs(ua(n-1,j));
 
 ha(n-1,j)=172.2*(Ga^0.5)*as;
 ha(n-1,j)=abs(ha(n-1,j));
 hrc(1,j)=0.93*(5.670373*(10^-8))*(T(n-1,j)+Tsky(1,j))*((T(n-1,j)^2)+(Tsky(1,j)^2));
 Tap(1,j)=(Ta(1,j)*hvv+Tsky(1,j)*hrc)/(hvv+hrc);
 Tpi(1,j)=((hvfpi*(hcp+hrc(1,j)+hvv)*T(n-1,j))+(hcp*(hrc(1,j)+hvv)*Tap(1,j)))/((hcp*(hrc(1,j)+hvv))+(hvfpi*(hcp+hrc(1,j)+hvv)));
 
 per=x*3600*(hvfpi*Sp*(T(n-1,j)-Tpi(1,j)))/(Ga*Ca(n-1,j)+Ga*H(n,j)*Cv(n-1,j));
 T(n,j)=T(n-1,j)-(ha(n-1,j)*x/(Ga*Ca(n-1,j)+Ga*H(n,j)*Cv(n-1,j)))*(T(n-1,j)-Tg(n,j-1));
 %Calculo de Tg(n,j)
 Cp(n,j-1)=1755+23.45*0.01*M(n,j-1);
 Cp(n,j-1)=abs(Cp(n,j-1));
 
 hfg(n,j-1)=2501.6+(1840-4180)*(Tg(n,j-1)+273);
 hfg(n,j-1)=abs(hfg(n,j-1));
 
 Cw=4200;
 
 Tg(n,j)=Tg(n,j-1)+(10/6)*(h*ha(n-1,j)/(rhop(n,j-1)*Cp(n,j-1)+rhop(n,j-1)*Cw*0.01*M(n,j-1)))*(T(n-1,j)-Tg(n,j-1))-(10/6)*Ga*h*(hfg(n,j-1)+(T(n,j)-Tg(n,j-1))*Cv(n-1,j))*(H(n,j)-H(n-1,j))/(x*(rhop(n,j-1)*Cp(n,j-1)+rhop(n,j-1)*Cw*0.01*M(n,j-1)));
 LV(n,j)=hfg(n,j-1)+(T(n,j)-Tg(n,j-1))*Cv(n-1,j);
end
t=6:h:(tfinal-0.5);
t1=6:h:(tfinal-1);
M4=M(4,:);
M4=M4(13:1:X);
M3=M(3,:);
M3=M3(13:1:X);
M2=M(2,:);
M2=M2(13:1:X);
plot(t,M4,'b','linewidth',3)
hold on
plot(t,M3,'r','linewidth',3)
plot(t,M2,'k','linewidth',3)
legend('Bandeja 3','Bandeja 2','Bandeja 1')
xlabel('tiempo(hrs del dia)')
ylabel('Humedad del cacao(%db)')
grid on
figure(2)
dM1=dM(2,:);
dM1=dM1(2:1:X);
divM=M/M0;
divM1=divM(2,:);
divM1=divM1(1:1:X-1);
plot(divM1(13:1:X-1),-dM1(13:1:X-1)/0.5,'k','linewidth',3)
hold on
dM2=dM(3,:);
dM2=dM2(2:1:X);
divM2=M/M0;
divM2=divM2(3,:);
divM2=divM2(1:1:X-1);
plot(divM2(13:1:X-1),-dM2(13:1:X-1)/0.5,'r','linewidth',3)
dM3=dM(4,:);
dM3=dM3(2:1:X);
divM3=M/M0;
divM3=divM3(4,:);
divM3=divM3(1:1:X-1);
plot(divM3(13:1:X-1),-dM3(13:1:X-1)/0.5,'b','linewidth',3)
legend('Bandeja 1','Bandeja 2','Bandeja 3')
xlabel('M/M0 (db)')
ylabel('-dM/dt (kg de agua/kg de cacao *h^-1)')
grid on
figure(3)
dM11=dM(2,:);
dM11=dM11(14:1:X);
dM22=dM(3,:);
dM22=dM22(14:1:X);
dM33=dM(4,:);
dM33=dM33(14:1:X);
plot(t1,-dM11/0.5,'k','linewidth',3)
hold on
plot(t1,-dM22/0.5,'r','linewidth',3)
plot(t1,-dM33/0.5,'b','linewidth',3)
legend('Bandeja 1','Bandeja 2','Bandeja 3')
xlabel('tiempo(hrs del dia)')
ylabel('-dM/dt (kg de agua/kg de cacao *h^-1)')
grid on
figure(4)
plot(H(2,:),'.k')
hold on
plot(H(3,:),'.r')
plot(H(4,:),'.b')
legend('layer 1','layer 2','layer 3','layer 4')
xlabel('tiempo(hrs)')
ylabel('Humedad especifica del aire(kg/kg)')
grid on
figure(5)
Tent=T(1,:);
Tent=Tent(14:1:X)
T1=T(2,:);
T1=T1(14:1:X);
T2=T(3,:);
T2=T2(14:1:X);
T3=T(4,:);
T3=T3(14:1:X);
plot(t1,Tent,'g','linewidth',2)
hold on
plot(t1,T1,'k','linewidth',2)
plot(t1,T2,'r','linewidth',2)
plot(t1,T3,'b','linewidth',2)
legend('Entrada a la camara','Bandeja 1','Bandeja 2','Bandeja 3')
xlabel('tiempo(hrs)')
ylabel('Temperatura del aire desecante (C)')
grid on
figure(6)
Tg1=Tg(2,:);
Tg1=Tg1(13:1:X)
Tg2=Tg(3,:);
Tg2=Tg2(13:1:X)
Tg3=Tg(4,:);
Tg3=Tg3(13:1:X)
plot(t,Tg1,'k','linewidth',3)
hold on
plot(t,Tg2,'r','linewidth',3)
plot(t,Tg3,'b','linewidth',3)
legend('Bandeja 1','Bandeja 2','Bandeja 3')
xlabel('tiempo(hrs)')
ylabel('Temperatura del cacao (C)')
grid on
figure(7)
H1=H(2,:);
H1=H1(14:1:X);
H2=H(3,:);
H2=H2(14:1:X);
H3=H(4,:);
H3=H3(14:1:X);
plot(t1,H1,'k','linewidth',3)
hold on
plot(t1,H2,'r','linewidth',3)
plot(t1,H3,'b','linewidth',3)
legend('Bandeja 1','Bandeja 2','Bandeja 3')
xlabel('tiempo(hrs)')
ylabel('Humedad especifica del aire desecante')
grid on
axis([6 75 0.005 0.035])
error:
Index in position 1 exceeds array bounds (must not exceed 1).
Error in secado (line 94)
 RH(n-1,j)=1.626*H(n-1,j)*(10^5)/((1+1.608*H(n-1,j))*6892.4*exp(54.633-12301.69/(1.8*T(n-1,j)+492)-5.169*log(1.8*T(n-1,j)+492))); 
 
 solutions please :(

Best Answer

 for n=2:N
     
 end

You have nothing inside that loop. The loop will execute, doing nothing. Afterwards, the variable n will be left at the last value that was assigned to it, namely the value of N. So on the next line you try to access an array indexing with n-1 but n is N so that is index N-1 . However you have not assigned to that row.

Related Solutions

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

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

MATLAB: Someone correct me please! How can i plot with this codes to get the following graph

Here's a start, but what is N and x? And you didn't call plot(), so what, in this alphabet soup of a program, is being plotted along the x and y axes?

N = 100;  % Just a guess!!!
U = zeros(N*2+2);
A = zeros(N*2+2);
B = zeros(N*2+2);
C = zeros(N*2+2);
D = zeros(N*2+2);
Y = zeros(N*2+2);
onesN=1.0/N;
Re=5000;
EE=1.0;
t=0;
dt=EE*Re*N^2;
AA = -0.5 * EE;
BB = 1.0 + EE;
KKEND=12;
KKMOD=1;
Y(1)=0;
U(1)=0.0;
U(x)=1.0;
N=20
for j=2:N
  U(j)=0.0;
  A(1) = 1.0;
  B(1) = 1.0;
  C(1) = 1.0;
  D(1) = 1.0;
end
for j=2:N
  Y(j)=Y(j-1)+N;
  A(j)=AA;
  if j==N
    A(j)=0.0;
    D(j)=BB;
    B(j)=AA;
  end
  if j==2
    B(j)=0.0;
    C(j)=1.0-EE*U(j)+0.5*EE*U(j+1)+U(j-1);
  end
  if j==N
    C(j)=C(j)-AA*U(x);
  end
end
for j=3:N
  D(j)=D(j)-B(j)*A(j-1)/D(j-1);
  C(j)=C(j)-C(j-1)*B(j)/D(j-1);
  for k=2:N
    M=N-(k-2);
    U(M)=(C(M)-A(M)*U(M+1))/D(M);
    Y(1)=0.0;
    Y(x)=Y(N)+N;
    t=t+dt;
  end
end

Best Answer

Related Solutions

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

MATLAB: Someone correct me please! How can i plot with this codes to get the following graph

Related Question