solve a matrix ode45

1 visualizzazione (ultimi 30 giorni)
Rodrigo Blas
Rodrigo Blas il 14 Apr 2020
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);
Having trouble solving this ode45 matrix 3x2.
I split up the equations and the intial conditions into 3 1x2 matrices.
I want to solve each row at the same time.

Risposte (0)

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by