solve a matrix ode45

1 visualización (últimos 30 días)
Rodrigo Blas
Rodrigo Blas el 14 de Abr. de 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.

Respuestas (0)

Categorías

Más información sobre Ordinary Differential Equations en Help Center y File Exchange.

Etiquetas

Community Treasure Hunt

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

Start Hunting!

Translated by