24小时热门版块排行榜    

CyRhmU.jpeg
查看: 882  |  回复: 3

郑美琴琴

金虫 (著名写手)

[求助] ??? Attempted to access C(2); index out of bounds because numel(C)=1.已有2人参与

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));
回复此楼

» 猜你喜欢

» 本主题相关价值贴推荐,对您同样有帮助:

已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

材料廖

木虫 (正式写手)

【答案】应助回帖

★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★
感谢参与,应助指数 +1
郑美琴琴: 金币+20, ★★★★★最佳答案 2015-06-14 15:24:47
C(it)=(C(it-1)+3*m*kf*dt/(R*pp*V)*Cs(it))/(1+3*m*kf*dt/(R*pp*V));
C=C(it);
第二条语句有问题,把c变成了一个数,而不是向量了
2楼2015-06-14 09:44:33
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

zhuguiqiu

禁虫 (文坛精英)

感谢参与,应助指数 +1
本帖内容被屏蔽

3楼2015-06-14 09:50:58
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

郑美琴琴

金虫 (著名写手)

引用回帖:
2楼: Originally posted by 材料廖 at 2015-06-14 09:44:33
C(it)=(C(it-1)+3*m*kf*dt/(R*pp*V)*Cs(it))/(1+3*m*kf*dt/(R*pp*V));
C=C(it);
第二条语句有问题,把c变成了一个数,而不是向量了

大神,太谢谢你了!感激不尽!
4楼2015-06-14 15:25:05
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖
相关版块跳转 我要订阅楼主 郑美琴琴 的主题更新
信息提示
请填处理意见