帮看下这个程序 Input argument "y" is undefined.

lxx244lxx

function PDEs2DS_324

clear 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 = [543.3        574.5203462        607.6242169        640.1415255        671.1218255        700.2984569        727.6492903        753.2367257        777.1515236        799.4921566        820.3568002 0        0.055093867        0.113920607        0.171932185        0.227264584        0.279360997        0.328151221        0.373737811        0.416283021        0.455966554        0.492968961 0        0.052761405        0.106914317        0.15911271        0.208568428        0.255208536        0.299132063        0.340477585        0.379389327        0.416007726        0.450466717];  % y=[T1 T2 T3 T4 x1 x2 x3]
tspan=0:1:3000;
[t,y]=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 = [dTdR dx1dt dx2dt]';

% ------------------------------------------------------------------
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);
% ------------------------------------------------------------------

ChaChing

回复 1 # lxx244lxx 的帖子

1.程序是给齐了, 但报错没给齐, 大牛们会懒得出手的
??? Input argument "y" is undefined.
Error in ==> Euqations at 4
T = y(1:N);
...
Error in ==> zza at 26
[t,y]=ode45(@Euqations,tspan,y0);

2.参数不是这样传的
3.看看
4F, 常见的程序出错问题整理 (eight)
http://forum.vibunion.com/forum/thread-46001-1-1.html
关于求解变参数微分方程
http://forum.vibunion.com/forum/ ... id=44972&page=1

lxx244lxx

回复 2 # ChaChing 的帖子

谢谢了，这个问题我已经解决了。我还想问下，关于多元非线性微分方程组的解法，例如有30元。。什么方法比较好。。。

ChaChing

回复 3 # lxx244lxx 的帖子

这方面不熟, 待真正高手路过!

或搜搜老帖, 或明显些发帖询问!

ChaChing

lxx244lxx 发表于 2011-3-27 11:57

                               
登录/注册后可看大图

回复 2 # ChaChing 的帖子

谢谢了，这个问题我已经解决了。我还想问下，关于多元非线性微分方程组的解法， ...

恭贺解决了! 建议与大家分享成果, 共同进步!

西西想幸福

回复 1 # lxx244lxx 的帖子

我也想要解一个动力学的微分方程组，跟你的差不多！！请问你是怎么解决这个问题的？