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

lxx244lxx

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

qibbxxt

你的方程定义的有问题，x1和x2可不可以放在一个数组中？