Graph of ODE comes out wrong

%%data
R=62.364;            %%torr L/gmol K
T=473;               %%K
vo=100/3600;         %%L/s
pnh3=.046;          %%torr
po2=.068;           %%torr
%%intial conditions {gmol/s}
icnh3=(pnh3*vo)/(R*T);
ico2=(po2*vo)/(R*T);

ic=[icnh3,ico2];
zspan=[0 1];

[z,f]=ode45(@Ljallo_f_q1,zspan,ic);
plot(z,f(:,1),z,f(:,2));
xlabel('Length');
ylabel('Flowrate')
legend('Fnh3','Fo2');

function diffyq1=Ljallo_f_q1(~,X)
x=X(1); y=X(2);
%%Data
R1=1.987 ;           %%cal/gmol*K
R=62.364;            %%torr L/gmol K
T=473;               %%K
vo=100/3600;         %%L/s
vo2=-5/2; vnh3=-1;   %%coeff
%%Perimeter cm
PL=40*10^-4;PW=300*10^-4;  
Per=2*PL+2*PW;
%%Concentration and Pressure
cnh3=y/vo; pnh3=cnh3*R*T;   
co2=y/vo; po2=co2*R*T;
%%Rate Components 
A=3.4*10^-4*exp(21700/(R1*T))*pnh3*po2^(1/2);
B=1+8*10^-2*exp(4400/(R1*T))*po2^(1/2);
C=1+1.6*10^-3*exp(25500/(R1*T))*pnh3;
%%Rate 
rnh3=A/(B*C);
%%Diffy eq
diffyq1(1)=Per*vnh3*rnh3;
diffyq1(2)=Per*vo2*rnh3;
diffyq1=diffyq1';
end

This graph is wrong and is telling me that the intial conditions start at 0 when I plot it. This shouldnt be the case. Not sure where I'm wrong

Answer 1

[Walter Roberson]

LJallo_f_q1 does not depend upon X(1), only upon X(2). If you happen to have the symbolic toolbox, you can find the expression for the results

If you assume that X(2) is non-negative, then the first element comes out (to one decimal place) with the form

-(2.7e+14*X2^(3/2))/((1.0e+15*X2 + 1.0)*(8.9e+3*X2^(1/2) + 1.0))

with non-negative X2, the numerator 2.7e+14*X2^(3/2) is positive and real-valued, and so is the denominator (1.0e+15*X2 + 1.0)*(8.9e+3*X2^(1/2) + 1.0) because each of the terms in the denominator is positive times X2 plus positive constant and that can never be less than the positive constant for non-negative X2.

So you have -positive/positive and that is -positive and you have to conclude that for any non-negative initial X2 that the function will continue to produce negative values that will decrease the boundary condition value, until eventually X(2) is negative.

Once you do not assume that X(2) is non-negative, then the first elemetn comes out (to one decimal place) with the form

 -(1.2e+38*X2*conj(X2^(1/2)))/((5.4e+23*X2 + 5.2e+8)*(7.8e+18*conj(X2^(1/2)) + 8.8e+14))

Notice the X2^(1/2) . WIth negative X2 that is going to generate a complex result, and taking the complex conjugate of that is still going to be complex.

Putting these two together, any positive initial conditions are eventually going to lead to complex-valued results. The only question is whether your initial values are large enough that within your time span

If my calculations are right, then for that ic(1) value, your ic(2) would need to be no larger than about 6e-15 for you not to encounter complex values by time 1, but it is about 10^7 times larger than that.

Your plot does not show the plot starting from 0. Notice the scale of the plot on the y axes, that the marks are 1E-5 apart. Your plot starts from about 5E-8 which on the scale of the plot is not distinguishable from 0 unless you zoom in a fair way.

If you notice in the command window, after your plot you have

Warning: Imaginary parts of complex X and/or Y arguments ignored.

That is because your y is complex-valued.

By the way, please remember in your line of code

 diffyq1=diffyq1';

that you are taking the conjugate transpose, not the regular transpose. That will not affect whether the results are complex or not, but it will affect the complex values returned. The plot of the imaginary component is considerably more interesting if you use the regular transpose .' instead of the conjugate transpose '

[Rodrigo Blas]

%%data
R=62.364;            %%torr L/gmol K
T=473;               %%K
vo=100/3600;         %%L/s
pnh3=.046;          %%torr
po2=.068;           %%torr
%%intial conditions {gmol/s}
icnh3=(pnh3*vo)/(R*T);
ico2=(po2*vo)/(R*T);

ic=[icnh3,ico2];
zspan=[0 1];

[z,f]=ode45(@Ljallo_f_q1,zspan,ic);
plot(z,f(:,1),z,f(:,2));
xlabel('Length');
ylabel('Flowrate')
legend('Fnh3','Fo2');

function diffyq1=Ljallo_f_q1(~,X)
fnh3=X(1); fo2=X(2);
%%Data
R1=1.987 ;           %%cal/gmol*K
T=473;               %%K
vo2=-5/2; vnh3=-1;   %%coeff
P=0.115;            %%torr
%%Perimeter cm
PL=.004;PW=.03;  
Per=2*PL+2*PW;
%%ratio,P,Flowrate
ft=fnh3+fo2;
pnh3=(fnh3/ft)*P; 
po2=(fo2/ft)*P;
%%Rate Components 
A=3.4*10^-4*exp(21700/(R1*T))*pnh3*po2^(1/2);
B=1+8*10^-2*exp(4400/(R1*T))*po2^(1/2);
C=1+1.6*10^-3*exp(25500/(R1*T))*pnh3;
%%Rate 
rnh3=A/(B*C);
%%Diffy eq
diffyq1(1)=Per*vnh3*rnh3;
diffyq1(2)=Per*vo2*rnh3;
diffyq1=diffyq1';
end

I have adjusted the code. However, I stil get the same kind of graph.

[Walter Roberson]

Again, your values start from about 1e-8, which just is not a visible difference from 0 compared to the plot scale of 1e-5. If you zoom the plot enough you will see that it does not start from 0.

Negative flow rates do seem unlikely, but that is what you compute. We would need the original equations in mathematical form to look more closely.