| 查看: 294 | 回复: 0 | ||
郑美琴琴金虫 (著名写手)
|
[求助]
程序出错
|
|
function qt=HSDMJ(t,q) global nr dr drs ... r z Ds C Cs q... v R kf K... C0 p n ncall qr qrr pp m V it y % r网格化 for j=1:nr % 2/r*qr if(j==1) % r=0 qr(j)=2.0*(q(j+1)-q(j))/drs; elseif(j==nr) % r=R qr(j)=2/R*kf/(pp*Ds)*(C-Cs); else qr(j)=2/r(j)*(q(j+1)-q(j-1))/(2*dr); end % qrr if(j==1) qrr(j)=2.0*(q(j+1)-q(j))/drs; elseif(j==nr) qf(j)=q(j-1)+2.0*dr*kf/(pp*Ds)*(C-Cs); qrr(j)=(qf(j)-2.0*q(j)+q(j-1))/drs; else qrr(j)=(q(j+1)-2.0*q(j)+q(j-1))/drs; end % PDEs qt(j)=Ds*(qrr(j)+qr(j)); end % 转置和计数 qt=qt'; ncall=ncall+1; clc; clear all; close all; %% 定义全局变量 global nr dr drs ... r z Ds C q... v R kf K... C0 p n ncall qr qrr F pp V x y m it y Cs %% 模型参数 C0=300; % mg/L F=30; % mL/min m=23; % g d=0.02; % m R=0.0015; % m p=1100; %kg/m^3 pp=718.6; %kg/m^3 Ds=60*4.187e-12; % m^2/min K=32.948; n=2.822; kf=60*2.495e-5; % m/min V=0.025; %L %% r方向网格 % r的范围为0-R nr=7; r=linspace(0, R, nr); dr=R/(nr-1); for j=1:nr r(j)=(j-1)*dr; end drs=dr^2; tf=180; tout=0:5:tf; nout=37; ncall=0; dt=tf/(nout-1); Cs(1)=0; C(1)=C0; %% 初始条件 for j=1:nr q(j)=0; q0(j)=q(j); end C=zeros(nout,1); for it=1:nout if(it==1) C(it)=C0; else C(it)=(C(it-1)+3*m*kf*dt/(R*pp*V)*Cs(it))/(1+3*m*kf*dt/(R*pp*V)); end C=C(it); Cs=Cs(it); %% ode参数设置及方程求解 reltol=1.0e-04; abstol=1.0e-04; options=odeset('RelTol',reltol,'AbsTol',abstol); [t,q]=ode15s(@HSDMJ,tout,q0,options); Cs(it+1)=(q(it,nr)/K)^n; end 报错： ??? Attempted to access C(2); index out of bounds because numel(C)=1. Error in ==> HSDMJ_1 at 47 C(it)=(C(it-1)+3*m*kf*dt/(R*pp*V)*Cs(it))/(1+3*m*kf*dt/(R*pp*V)); |
回复此楼» 猜你喜欢
2019年青年基金涵评意见，大家看看几个A，几个B?
已经有4人回复
-
应该是93bebmhtak前后十一个字符比较关键
已经有29人回复
-
奇怪，两个人的filecode固定段从头到尾一模一样
已经有9人回复
-
是这周出结果还是下周出结果？
已经有3人回复
-
小木虫上这么多卖论文的，真有人买论文么？感觉没必要啊
已经有9人回复
-
有时候，自然基金真的不能太认真 （我的申报经验）
已经有4人回复
-
帮忙看看fileCode
已经有13人回复
-
为什么网上很多人说本周 12号出结果
已经有7人回复
-
读博的好处
已经有3人回复
-
什么时候出结果，有咨询渠道？？？
已经有3人回复
100