solve a matrix ode45

Question

Rodrigo Blas 2020-4-14

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/517719-solve-a-matrix-ode45

pnh3=.046;      %%torr
po2=.068;       %%torr
R1=1.987 ;       %%cal/gmol*K
R=62.364;       %%torr L/gmol K
vo=100/3600;    %%L/s
P=.115;     %%torr
Per=40*300; %%micro meters 
zf=1;        %%cm
vnh3=1;
vo2=5/2;
zspan=[0 1];    %%cm 
T=zeros(3,1);
T(1)=473;
T(2)=673;
T(3)=800;
fio=Ljallo_f2_q1(T,R,vo,pnh3,po2);
ic1=fio(1,:);
ic2=fio(2,:);
ic3=fio(3,:);
syms y
DFDT=Ljallo_f1_q1(T,R1,Per,vnh3,vo2,y,vo,R);
diffy1=DFDT(1,:);
diffy2=DFDT(2,:);
diffy3=DFDT(3,:);
f1=@(z,y) diffy1;
f2=@(z,y) diffy2;
f3=@(z,y) diffy3;
[z,y]=ode45(@(z,y) f1,zspan,ic1);
[z,y]=ode45(@(z,y) f2,zspan,ic2);
[z,y]=ode45(@(z,y) f3,zspan,ic3);
function dfdt=Ljallo_f1_q1(T,R1,Per,vnh3,vo2,y,vo,R)
A=zeros(3,1);
B=zeros(3,1);
C=zeros(3,1);
rn=zeros(3,1);
dfdt=zeros(3,2);
for i=1:3
A(i)=3.4*10^-8*exp(21700./(R1.*T(i)))*(y(1)./vo*R.*T(i)).*(y(2)./vo*R.*T(i)).^(.5);
B(i)=1+8*10^-2*exp(4400/(R1.*T(i)))*(y(2)/vo*R*T(i))^(.5);
C(i)=1+1.6*10^-3*exp(25500/(R1.*T(i)))*(y(1)/vo*T(i));
rn(i)=A(i)/(B(i)*C(i));
dfdt(i,1)=rn(i).*Per*vnh3;
dfdt(i,2)=rn(i).*Per*vo2;
end 
function iFlows=Ljallo_f2_q1(T,R,vo,pnh3,po2)
cnh3=zeros(3,1);
finh3o=zeros(3,1);
co2=zeros(3,2);
fio2o=zeros(3,2);
iFlows=zeros(3,2);
for i=1:3
cnh3(i,1)=pnh3/(R.*T(i));
finh3o(i,1)=cnh3(i,1).*vo;
iFlows(i,1)=finh3o(i,1);
co2(i,2)=po2/(R.*T(i));
fio2o(i,2)=co2(i,2).*vo;    
iFlows(i,2)=fio2o(i,2);
end
>> Ljallo_s_q1
Index exceeds the number of array elements (1).
Error in sym/subsref (line 890)
            R_tilde = builtin('subsref',L_tilde,Idx);
Error in Ljallo_f1_q1 (line 8)
A(i)=3.4*10^-8*exp(21700./(R1.*T(i)))*(y(1)./vo*R.*T(i)).*(y(2)./vo*R.*T(i)).^(.5);
Error in Ljallo_s_q1 (line 23)
DFDT=Ljallo_f1_q1(T,R1,Per,vnh3,vo2,y,vo,R);
 