| 查看: 1596 | 回复: 12 | ||
| 当前只显示满足指定条件的回帖，点击这里查看本话题的所有回帖 | ||
郑美琴琴金虫 (著名写手)
|
[求助]
求助，程序中的错误时什么原因 已有1人参与
|
|
|
这是Command window的程序： global nr nz dr dz drs dzs... r z Dc Dt ca Tk... cae Tke h k E R... rk0 v rho Cp Tw dH... ncall %模型参数 ca0=0; cae=0.01; Tk0=305.0; Tke=305.0; Tw=355.0; r0=2.0; z1=100.0; v=1; Dc=0.1; Dt=0.1; k=0.01; h=0.01; rho=1.0; Cp=0.5; rk0=1.5e+09; dH=-10000.0; E=15000.0; R=1.987; %x轴网格 nz=20; dz=z1/nz; for i=1:nz; z(i)=i*dz; end %半径网格化 nr=7; dr=r0/(nr-1); for j=1:nr; r(j)=(j-1)*dr; end drs=dr^2; %自变量 tf=200.0; tout=[0:50:tf]'; nout=5; ncall=0; %初始条件 for i=1:nz; for j=1:nr ca(i,j)=ca0; Tk(i,j)=Tk0; y0((i-1)*nr+j)=ca(i,j); y0((i-1)*nr+j*nz*nr)=Tk(i,j); end end %ODE集成 reltol=1.0e-04;abstol=1.0e-04; options=odeset('RelTol',reltol,'AbsTol',abstol); [t,y]=ode15s(@pde_13,tout,y0,options); 这是function pde_的程序： function yt=pde_13(t,y) global nr nz dr dz drs dzs... r z Dc Dt ca Tk... cae Tke h k E R... rk0 v rho Cp Tw dH... ncall for i=1:nz for j=1:nr ij=(i-1)*nr+j; ca(i,j)=y(ij); Tk(i,j)=y(ij+nr*nz); end end for i=1:nz for j=1:nr if(j==1) car(i,j)=2*(ca(i,j+1)-ca(i,j))/drs; Tkr(i,j)=2*(Tk(i,j+1)-Tk(i,j))/drs; elseif(j==nr) car(i,j)=0.0; Tkr(i,j)=(1/r(j))*(h/k)*(Tw-Tk(i,j)); else car(i,j)=(1/r(j))*(ca(i,j+1)-ca(i,j-1))/(2*dr); Tkr(i,j)=(1/r(j))*(Tk(i,j+1)-Tk(i,j-1))/(2*dr); end if(j==1) carr(i,j)=2*(ca(i,j+1)-ca(i,j))/drs; Tkrr(i,j)=2*(Tk(i,j+1)-Tk(i,j))/drs; elseif(j==nr) carr(i,j)=2*(ca(i,j-1)-ca(i,j))/drs; Tkf(i,j)=Tk(i,j-1)+2*dr*h/k*(Tw-Tk(i,j)); Tkrr(i,j)=(Tkf-2.0*Tk(i,j)+Tk(i,j-1))/drs; else carr(i,j)=(ca(i,j+1)-2.0*ca(i,j)+ca(i,j-1))/drs; Tkrr(i,j)=(Tk(i,j+1)-2.0*Tk(i,j)+Tk(i,j-1))/drs; end if(i==1) caz(i,j)=(ca(i,j)-cae)/dz; Tkz(i,j)=(Tk(i,j)-Tke)/dz; else caz(i,j)=(ca(i,j)-ca(i-1,j))/dz; Tkz(i,j)=(Tk(i,j)-Tk(i-1,j))/dz; end rk=rk0*exp(-E/(R*Tk(i,j)))*ca(i,j)^2; cat(i,j)=Dc*(carr(i,j)+car(i,j))-v*caz(i,j)-rk; Tkt(i,j)=Dt*(Tkrr(i,j)+Tkr(i,j))-v*Tkz(i,j)-dH/(rho*Cp)*rk; end end for i=1:nz for j=1:nr ij=(i-1)*nr+j; yt(ij)=cat(i,j); yt(ij+nr*nz)=Tkt(i,h); end end yt=yt'; ncall=ncall+1; 运行结果： Warning: Divide by zero. > In pde_13 at 44 In funfun\private\odearguments at 110 In ode15s at 227 Warning: Divide by zero. > In pde_13 at 44 In funfun\private\odearguments at 110 In ode15s at 227 Warning: Divide by zero. > In pde_13 at 44 In funfun\private\odearguments at 110 In ode15s at 227 Warning: Divide by zero. > In pde_13 at 44 In funfun\private\odearguments at 110 In ode15s at 227 Warning: Divide by zero. > In pde_13 at 44 In funfun\private\odearguments at 110 In ode15s at 227 Warning: Divide by zero. > In pde_13 at 44 In funfun\private\odearguments at 110 In ode15s at 227 ??? Subscripted assignment dimension mismatch. Error in ==> pde_13 at 32 Tkrr(i,j)=(Tkf-2.0*Tk(i,j)+Tk(i,j-1))/drs; 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, ... 各位大神，帮帮忙。感激不尽！ |
回复此楼» 猜你喜欢
极限拉扯，如何抉择
已经有15人回复
-
申请科研项目纯靠实力能成功吗？
已经有8人回复
-
中科院海洋所招收海洋材料相关领域博士研究生及研究助理
已经有4人回复
-
40+中年人失业后，还能找到好工作吗？
已经有3人回复
-
二审4个多月了还是Awaiting Reviewer Scores
已经有4人回复
-
OpenAI又攻克数学难题了，离AI做有机合成还有多远？
已经有13人回复
-
课题组会帮你通过考核么？
已经有7人回复
-
入职后第一次管课题组的专利缴费，对账才发现官费比别人多掏了一倍
已经有6人回复
-
售一区SCI文章T0P，我:8O.551.O54,科目全,可十急
已经有9人回复
-
售SCI一区文章，我:8O5.5.1.O5.4,科目全,可伽急
已经有10人回复
郑美琴琴
金虫 (著名写手)
- 应助: 3 (幼儿园)
- 金币: 1607.8
- 散金: 84
- 红花: 1
- 帖子: 1002
- 在线: 147.4小时
- 虫号: 1765560
- 注册: 2012-04-19
- 性别: MM
- 专业: 化工系统工程
11楼2015-04-26 19:58:35
★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ...
感谢参与，应助指数 +1
郑美琴琴: 金币+100, ★★★★★最佳答案 2015-04-19 16:14:21
感谢参与，应助指数 +1
郑美琴琴: 金币+100, ★★★★★最佳答案 2015-04-19 16:14:21
|
本帖内容被屏蔽 |
2楼2015-04-19 11:23:19
郑美琴琴
金虫 (著名写手)
- 应助: 3 (幼儿园)
- 金币: 1607.8
- 散金: 84
- 红花: 1
- 帖子: 1002
- 在线: 147.4小时
- 虫号: 1765560
- 注册: 2012-04-19
- 性别: MM
- 专业: 化工系统工程
3楼2015-04-19 16:15:14
|
本帖内容被屏蔽 |
4楼2015-04-19 16:39:45
200