PDE
in [Matlab]
|
From: Luca Zanotti Fragonara on 27 Feb 2010 15:12 Hello, it is very hard to understand what's going on in your code if we don't see the code. It seems that one of your variables is not well defined, and that's why you get your error.
From: Nga Do on 2 Mar 2010 02:42
Torsten Hennig <Torsten.Hennig(a)umsicht.fhg.de> wrote in message <1616844385.279521.1267270324231.JavaMail.root(a)gallium.mathforum.org>... > > On Feb 26, 7:03 pm, Torsten Hennig > > <Torsten.Hen...(a)umsicht.fhg.de> > > wrote: > > > > Good afternoon, > > > > > > > I am using PDEs to solve the mathematic problem, > > heat > > > > and mass transfer. However, i met this mistake, > > how > > > > can i overcome it? can you help me? > > > > Warning: Failure at t=0.000000e+000. Unable to > > meet > > > > integration > > > > tolerances without reducing the step size below > > the > > > > smallest value > > > > allowed (7.905050e-323) at time t. > > > > > In ode15s at 753 > > > > In pdepe at 320 > > > > > > > Warning: Time integration has failed. > > > > > > > and > > > > > > > ??? Attempted to access sR(1); index out of > > bounds > > > > because > > > > numel(sR)=0. > > > > > > > i am waiting for you reply. > > > > > > .. and we are waiting for your code. > > > > > > > Regards. > > > > > > Best wishes > > > Torsten.- Hide quoted text - > > > > > > - Show quoted text - > > > > oh, i just see that > > > > The hardest one to debug is if you have used / > > instead of ./ somewhere > > in a function. If you integrate the function, this > > may cause quad to > > give a very odd error message: > > ??? Attempted to access y(7); index out of bounds > > because numel(y)=1. > > > > How can i overcome that mistake, i do not think i > > used / instead ./ > > With my answer I wanted to say: > We cannot answer your questions if we do not see > exactly what you coded. > > Best wishes > Torsten. This is the part which i met the errors, when i change it with other equation, it runs smoothly. function [c,f,s] = pdex_apde(x,t,u,DuDx) global .... CH4t=pCH4_0*(1-u(2)); H2t=(pH2_0+2*pCH4_0*u(4)); COt=(pCO_0+pCH4_0*u(3)); CO2t=pCH4_0*(u(2)-u(3)); O2t=(pO2_0+pCH4_0*(2*u(4)-2*u(2)+.5*u(3))+.5*pCO_0+.5*pH2_0); H2Ot=2*pCH4_0*(u(2)-u(4)); N2t=pN2_0; n_sumt=CH4t+H2t+COt+CO2t+O2t+H2Ot+N2t; dent=1/n_sumt*P_0_sum; pCH4=CH4t*dent; pH2=H2t*dent; pCO=COt*dent; pCO2=CO2t*dent; pO2=O2t*dent; pH2O=H2Ot*dent; pN2=N2t*dent; M_t=(pCH4*16+2*pH2+28*pCO+44*pCO2+32*pO2+18*pH2O+28*pN2)/P_0_sum; rhogT=rhog; v_t=v; omega_1=cpg*rhog/(epxg*cpg*rhogT+epxs*rhos*cps); Pe=(epxg*cpg*rhogT+epxs*rhos*cps)*L*v_t/lamdac; omega_2=1/((epxg*cpg*rhogT+epxs*rhos*cps)*v_t); Pe_2=2200;%moi them vao Jun 24 06 k01=0.0013554;%973 k02=5.922*10^8;% k03=6.028*10^-4;% %k04=.55*2.83*10^14;%2.83 10^14 beta1=(-deltaH1)*L*gama0/((epxg*cpg*rhogT+epxs*rhos*cps)*v_t*T0); beta2=(-deltaH2)*L*gama0/((epxg*cpg*rhogT+epxs*rhos*cps)*v_t*T0); beta3=(-deltaH3)*L*gama0/((epxg*cpg*rhogT+epxs*rhos*cps)*v_t*T0); %beta4=(-deltaH4)*L*gama0/((epxg*cpg*rhogT+epxs*rhos*cps)*v_t*T0); phi=L/v_t; k1=k01*exp(-gama1/((u(1)/gama0)+1)); k2=k02*exp(-gama2/(u(1)/gama0+1)); k3=k03*exp(-gama3/(u(1)/gama0+1)); K_CH4_ox=19.6116*exp((-21.731/(R*T0))*(1/(u(1)/gama0+1)-T0/973)); K_O2_ox=45.04*exp((-55.513/(R*T0))*(1/(u(1)/gama0+1)-T0/973)); K_p2=(1.198*10^17)/((103)^2)*exp(-gamap2/(u(1)/gama0+1)); K_p3=(1.767*10^-2)/((103)^2)*exp(-gamap3/(u(1)/gama0+1)); K_H=5.68*10^-10*exp(-gamah2/(u(1)/gama0+1)); K_CO=5.127*10^-13*exp(-gamaco/(u(1)/gama0+1)); K_H2O=9.521*exp(-gamah2o/(u(1)/gama0+1)); T_eta=u(1)*T0/(E0/R/T0)+T0; eta_1= polyval(eta_1_fun,T_eta); eta_2= polyval(eta_2_fun,T_eta); eta_4= polyval(eta_4_fun,T_eta); A=(1+(K_CH4_ox*pCH4)+sqrt(K_O2_ox*pO2)); B=1+K_CO*pCO+K_H*(pH2)^0.5+K_H2O*(pH2O/pH2); r1=((k1*K_CH4_ox*pCH4)*sqrt(K_O2_ox*pO2))./A^2; r2=k2*(pCH4*(pH2O^0.5)/(pH2^(1.25))*(1-(pCO*(pH2^3)/(K_p2*pCH4*pH2O))))./B^2; r3=k3*(pCO*(pH2O^0.5)/(pH2^0.5)*(1-(pCO2+pH2/(K_p3*pCO*pH2O))))./B^2; %--------- c = [1;1;1;1]; if ((x>=a)&&(x<=b)) f = [DuDx(1)/Pe;DuDx(2)/Pe_2;DuDx(3)/Pe_2;DuDx(4)/Pe_2]; (1)/gama0+1)-Tw))+... s1 =-omega_1*DuDx(1)+omega_2*(-alpha*area*gama0*L/T0*(T0*(u(1)/gama0+1)-Tw))+... (epxs*rhos*((beta1*r1)+(beta2*r2)+(beta3*r3))); s2=-(1./epxs)*DuDx(2)+ M_CH4*epxs*rhos/(epxg*epxs)*phi*(r1+r2);%CH4 s3=-(1./epxs)*DuDx(3)+ M_CO*epxs*rhos/(epxg*epxs)*phi*(r2-r3);%CO s4=-(1./epxs)*DuDx(4)+ M_H2*epxs*rhos/(epxg*epxs)*phi*(3*r2+r3);%H2 s =[s1;s2;s3;s4]; else f = [DuDx(1)/Pe;0;0;0]; s1=omega_2*(-alpha*area*gama0*L/T0*(T0*(u(1)/gama0+1)-Tw))-omega_1*DuDx(1); s =[s1;0;0;0]; end length(s); % ------------------------------------------------------------------------- function u0 = pdex_aic(x) global utc stepsize xtotal a b stepsize stepsize_fine x1=xtotal;%(0:stepsize:1) p_tem=[1 0]; if x<=a pos_x=x/stepsize+1; else if x<=b pos_x=a/stepsize+(x-a)/stepsize_fine +1; else pos_x=a/stepsize+(b-a)/stepsize_fine+(x-b)/stepsize+1; end end pos_x=floor(pos_x); u0=[utc(pos_x);0;0;0]; function [pl,ql,pr,qr] = pdex_abc(xl,ul,xr,ur,t) global omega_1 pl = [-omega_1*ul(1);ul(2);ul(3);ul(4)];%-omega_1*ul(1) ql = [1;0;0;0]; pr=[0;0;0;0];%cho ddieu kien bien dt/dz=0 qr=[1;1;1;1]; Best wishes Quynh Nga |