24小时热门版块排行榜

返回列表

郑美琴琴

金虫 (著名写手)

应助: 3 (幼儿园)
金币: 1607.8
散金: 84
红花: 1
帖子: 1002
在线: 147.4小时
虫号: 1765560
注册: 2012-04-19
性别: MM
专业: 化工系统工程

[求助] matlab程序出错

大神帮帮忙，感激不尽！！！

M文件
function yt=ped_a(t,y)
global nr nz dr dz drs dzs...
      r  z Ds DL c q...
      cs v e R  kf K...
      c0 p n ncall
   %1D to 2D矩阵
   for i=1:nz
      for j=1:nr
         ij=(i-1)*nr+j;
         c(i)=y(i);
         q(i,j)=y(ij+nz);
      end
   end
   %r,z网格化
   for i=1:nz
      for j=1:nr
         %2/r*qr
         if(j==1) %r=0
               qr(i,j)=2.0*(q(i,j+1)-q(i,j))/drs;
         elseif(j==nr) %r=R
               qr(i,j)=2/R*kf/(p*Ds)*(c(i)-cs);
         else
               qr(i,j)=2/r(j)*(q(i,j+1)-q(i,j-1))/(2*dr);
         end
         %qrr
         if(j==1)
               qrr(i,j)=2.0*(q(i,j+1)-q(i,j))/drs;
         elseif(j==nr)
               qf(i,j)=q(i,j-1)+2.0*dr*kf/(p*Ds)*(c(i)-cs);
               qrr(i,j)=(qf(i,j)-2.0*q(i,j)+q(i,j-1))/drs;
         else
               qrr(i,j)=(q(i,j+1)-2.0*q(i,j)+q(i,j-1))/drs;
         end
         %cz
         if(i==1) %z=0
               cz(i)=(c(i)-c0)/dz;
         elseif(i==nz) %z=L
               cz(i)=0;
         else
               cz(i)=(c(i)-c(i-1))/dz;
         end
         %czz
         czz(i)=(c(i+1)-2.0*c(i)+c(i-1))/dzs;
         %PDEs
         qt(i,j)=Ds*(qrr(i,j)+qr(i,j));
         ct(i)=DL*czz(i)-v*cz(i)-3*(1-e)/(e*R)*kf*(c(i)-cs(i,j));
         q(i,j)=K*cs(i,j)^(1/n);
      end
   end
   %2D to 1D
   for i=1:nz
      for j=1:nr
         ij=(i-1)*nr+j;
         yt(i)=ct(i);
         yt(ij+nz)=qt(i,j);
      end
   end
   %转置和计数
   yt=yt';
   ncall=ncall+1;

Command window：
%全部区域
>> global nr nz dr dz drs dzs...
      r  z Ds DL c q...
      cs v e R  kf K...
      p n ncall
>> %模型参数
>> c0=250;
>> L=0.15;
>> R=0.001;
>> p=436.8;
>> v=5.8*10^(-5);
>> e=0.363;
>> DL=1.66*10^(-7);
>> Ds=2.4281*10^(-12);
>> K=0.14;
>> n=1;
>> %z方向网格
>> nz=20;
>> dz=L/nz;
>> for i=1:nz
z(i)=i*dz;
end
>> dzs=dz^2;
>> %r方向网格
>> nr=7;
>> dr=R/(nr-1);
>> for j=1:nr
r(j)=(j-1)*dr;
end
>> drs=dr^2;
%ODE集成
>> tf=200;
>> tout=[0.0:50:tf];
>> nout=5;
>> ncall=0;
>> %初始条件
>> for i=1:nz
for j=1:nr
if(i==1)
c(i)=c0
else
c(i)=0
end
q(i,j)=0
y0(i)=c(i);
y0((i-1)*nr+j+nz)=q(i,j)
end
end
%ODE集成
>> reltol=1.0e-04; abstol=1.0e-04;
options=odeset('RelTol',reltol,'AbsTol',abstol);
[t,y]=ode15s(@ped_a,tout,y0,options);

模拟结果：
???  In an assignment  A(I) = B, the number of elements in B and
I must be the same.

Error in ==> ped_a at 36
               cz(i)=(c(i)-c0)/dz;

Error in ==> funfun\private\odearguments at 110
f0 = feval(ode,t0,y0,args{:}); % ODE15I sets args{1} to yp0.

Error in ==> ode15s at 227
[neq, tspan, ntspan, next, t0, tfinal, tdir, y0, f0, odeArgs, ...

回复此楼