Solving simultaneous 3 ODE : I can't use the results in the equations themselves.

Hi everyone.

I want to solve 3 ODE equations simultaneously and I need to use the values calculated in the ODE itself. My variable is z, and dependent values are X, T, P. In the function file, I want to use the resulting X, T, P values but I obtained NaN values only. Also, if I write dep = dep' at the end of function file in order to convert it into column vector, I received an error that I should convert it to column vector.

Here is my main program :

clc , clear all

global nu a b c d e f h A Ea R eps rocat D Mw dcat mugas f0

%%data
nu = [-1 -0.5 0 1]; % stoichiometric matrix
R = 8.314; % [J.mol-1.K]
T0 = 703.15; %feed temperature [K]
F0 = 1.2676e3; % [mol.s-1] calculated by hand 
P0 = 101e3; % feed pressure [Pa]
y0 = [0.095 0.115 0.79 0 ]; % initial feed composition
f0 = F0.*y0; % inlet feed concentration of components [mol.s-1]
robed = 5e2; %bed/bulk density [kg.m-3]
eps = 0.40; %porosity
rocat = robed/(1-eps); %catalyst density [kg.m-3]
dcat = 8e-3; % catalyst diameter [m]
mugas = 36.6*10^-6; %[Pa.s] 
Mw = [64.066 16.000 14.007 80.006]*10^-3; % molecular weight [kg.mol-1]
X0 = 0;

a = [21.43029 30.03235 19.50583 24.02503];
b = [74.35094 8.772972 19.88705 119.4607];
c = [-57.75217 -3.988133 -8.598535 -94.38686];
d = [16.35534 0.788313 1.369784 26.96237];
e = [0.086731 -0.741599 0.527601 -0.117517];
f = [-305.7688 -11.32468 -4.935202 -407.8526];
h = [-296.8422 0 0 -395.7654];
a = a';
b = b';
c = c';
d = d';
e = e';
f = f';
h = h';
A =  5.3428;
Ea = 7.9670e+04; %[J.mol-1.K-1]

%%initial condition    
dep0(1) = X0;
dep0(2) = T0;
dep0(3) = P0;

%%ODE solver
L = 5; % [m]
D = 0.5; % [m]
wcat = (3.14/4.*D.^2).*L./robed; % catalyst weight [kg]
zrange = [0 L];
[z, dep] = ode15s(@newodesolv, zrange, dep0);
X = dep(1);
T = dep(2);
P = dep(3);

fi = f0 - nu.*f0(1)*X;

fprintf('%12.4e %12.4e %12.4e %12.4e\n', wcat, X, T, P)
fprintf('%12.4e %12.4e %12.4e %12.4e\n', fi)

------------------------------------------------------------------------------------------------

And this is my function file.

------------------------------------------------------------------------------------------------

function dep = newodesolv(z,dep)

%dep(1) = X
%dep(2) = T
%dep(3) = P

global nu a b c d e f h A Ea R eps rocat D Mw dcat mugas f0
X = dep(1);
T = dep(2);
P = dep(3);
fi = f0 - nu.*f0(1)*dep(1);
y = fi/(fi(1)+fi(2)+fi(3)+fi(4));
Mwmix = (Mw*y'); % molecular weight of mixture
rogas = (P * Mwmix)/(R*T); %[kg.m-3]
t = T/1000; %[K]
cp = a+b.*t+c.*(t^2)+d.*(t^3)+e./(t^2); % [J.mol-1.K-1] 
Formation_ent = (a.*t+(b.*t^2)/2+(c.*t^3)/3+(d.*t^4)/4-(e./t)+ f - h)*1000; % [J.mol-1]
delta_H_rxn = (nu*Formation_ent)*10^3; % [J.mol-1]
k = A*exp(-Ea/(R*T))*1000/101325; %[mol.kg,cat-1.Pa-1]
Kp = 7.97*10^-5*exp(12100./T); % [MPa^-0.5]
ficp = fi*cp; %[J.K-1.s-1]
fiMw = fi*Mw'; % [kg.s-1] 
r = k.*((P*y(2)*y(1)^(0.5))/(y(4)^(0.5))-((y(4)^(1.5))/((y(1)^(1.5))*(Kp^2)*10^6))); % [mol.kg,cat-1.s-1]
dep(1) = r.*(1-eps).*rocat*3.14/4*(D.^2)./f0(1);
dep(2) = -(r.*delta_H_rxn.*(1-eps).*rocat.*(3.14/4)*(D.^2))/(ficp);
dep(3) = -(fiMw)./((rogas).*(3.14/4).*(D.^2).*dcat)*...
    ((1-eps)/(eps^3)).*((150*(1-eps).*mugas./dcat) + (1.75.*(fiMw)./(3.14/4*(D.^2))));
dep;

Could you please help me? Thank you.

Answer 1

[Star Strider]

First, remove the semicolons from your dep(1)...dep(3) statements. It could be that one or more of them are returning matrices. If so, you need to figure out why and how to fix them so they will return elements of column vectors.

The NaNs are usually the result of (0/0) calculations. One NaN in the wrong place can ruin all your calculations. You need to find out what is generating them and fix it.

I don’t know what you’re doing, or what you can or can’t change to get the result you want. Only you can determine effectively what’s wrong with your code and how to fix it.

I suggest you set zrange = [0 0.1 0.2]; until you get your code running successfully.

[annisette]

Thank you for your answer. After I removed the semicolumns, I don't get the vector columns error annymore. But still I have NaN values. I think the reason of (0/0) calculation is, T, X and P values are the result of ODEs and they should be used before the ODE solution and also inside the ODEs. I don't know how to fix that part and need help on it.

[Star Strider]

My pleasure!

Removing the semicolons simply reports the output of those calculations to the Command Window. Are the values for dep(1) and the rest vectors or matrices? They should each be individual elements (scalars) and not vectors at each step.

For X, T, and P, the first step is to find the elements of those that are either zero, Inf or NaN, since those may also cause the problem. The easiest way is to use the find function for that.