帮看下这个程序 Input argument "y" is undefined.
function PDEs2DS_324clear all
clc
global r Fco Fco2 F H1 H2 Cpm N u L kg1a pb pa kg2a
% Parameters
r = 1; % W/(m K)
Fco = 14; % kg/(m^2 h)
Fco2 = 4; % kJ/(kg K)
F = 200; % J/mol
H1=206100; % kg/m^3
H2=164700; % radius of reactor, m
Cpm =31.02577;% u0*c0 obtained from u0*c0*A=0.069 kmol/h
N=11;% CP', J/(kg K)
u=1.4385; % kg/s
L=3;
kg1a=6.76;
pb=1175;
pa=2;
kg2a=5.83;
y0 = ;% y=
tspan=0:1:3000;
=ode45(@Euqations,tspan,y0);
T = y(:,1:N);
x1 = y(:,N+1:2*N);
x2 = y(:,2*N+1:3*N);
% ------------------------------------------------------------------
function dydx = Euqations(x1,x2,y)
global r Fco Fco2 F H1 H2 Cpm N u L
detaz=L/(N-1);
T = y(1:N);
x1 = y(N+1:2*N);
x2 = y(2*N+1:3*N);
rco= ReactionRate1(T(1:N),x1);
rco2=ReactionRate2(T(1:N),x2);
dTdt(1) = u*((T(1)-T(2))/detaz+(pi*r*r/F)*(rco(1)*H1+rco2(1)*H2)/Cpm);
dx1dt(1) =u*((x1(1)-x1(2))/detaz+(pi*r*r/Fco)*rco(1));
dx2dt(1) =u*((x2(1)-x2(2))/detaz+(pi*r*r/Fco2)*rco2(1));
dTdt(N) = u*((T(N-1)-T(N))/detaz+(pi*r*r/F)*(rco(N)*H1+rco2(N)*H2)/Cpm);
dx1dt(N) = u*((x1(N-1)-x1(N))/detaz+(pi*r*r/Fco)*rco(N));
dx2dt(N) = u*((x1(N-1)-x1(N))/detaz+(pi*r*r/Fco2)*rco(N));
for i = 2:N-1
dTdt(i) = u*((T(i-1)-T(i+1))/2/detaz+(pi*r*r/F)*(rco(i)*H1+rco2(i)*H2)/Cpm);
dx1dt(i) =u*((x1(i-1)-x1(i+1))/2/detaz+(pi*r*r/Fco)*rco(i));
dx2dt(i) =u*((x1(i-1)-x1(i+1))/2/detaz+(pi*r*r/Fco)*rco(i));
end
dydx = ';
% ------------------------------------------------------------------
function f1 = ReactionRate1(T,x,T0) % 计算反应速度
global kg1a pb pa
CA0=pa*100000/8.314/T0*0.07;
f1=(kg1a*CA0*(1-x)*pb*533488.642.*exp(-76856.4/8.314./T).*(CA0*(1-x)).^0.673833)/20./(kg1a*CA0*(1-x)+pb*533529*exp(-76857/8.314./T).*(CA0*(1-x)).^0.674);
function f2 = ReactionRate2(T,x,T0) % 计算反应速度
global kg2a pb pa
CB0=pa*100000/8.314/T0*0.02;
f2=(kg2a*CB0*(1-x)*pb*296452638.*exp(-105858/8.314./T).*(CB0*(1-x)).^0.008)/20./(kg2a*CB0*(1-x)+pb*296452638*exp(-105858/8.314./T).*(CB0*(1-x)).^0.008);
% ------------------------------------------------------------------
回复 1 # lxx244lxx 的帖子
1.程序是给齐了, 但报错没给齐, 大牛们会懒得出手的
??? Input argument "y" is undefined.
Error in ==> Euqations at 4
T = y(1:N);
...
Error in ==> zza at 26
=ode45(@Euqations,tspan,y0);
2.参数不是这样传的
3.看看
4F, 常见的程序出错问题整理 (eight)
http://forum.vibunion.com/forum/thread-46001-1-1.html
关于求解变参数微分方程
http://forum.vibunion.com/forum/viewthread.php?tid=44972&page=1
回复 2 # ChaChing 的帖子
谢谢了，这个问题我已经解决了。我还想问下，关于多元非线性微分方程组的解法，例如有30元。。什么方法比较好。。。 本帖最后由 ChaChing 于 2011-3-27 12:33 编辑
回复 3 # lxx244lxx 的帖子
这方面不熟, 待真正高手路过!
或搜搜老帖, 或明显些发帖询问!
lxx244lxx 发表于 2011-3-27 11:57 static/image/common/back.gif
回复 2 # ChaChing 的帖子
谢谢了，这个问题我已经解决了。我还想问下，关于多元非线性微分方程组的解法， ...
恭贺解决了! 建议与大家分享成果, 共同进步! 回复 1 # lxx244lxx 的帖子
我也想要解一个动力学的微分方程组，跟你的差不多！！请问你是怎么解决这个问题的？
页:
[1]