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

warning: imaginary parts of complex x and/or y arguments ignoredwile e. coyote

I believe this error may be coming from line 92, where I apply the transformation to KG to get KH. When I remove the semicolon and don't suppress output for line 91, KG always results in a real number, but when I look at line 92, KH jumps around and results in complex numbers.
Line 92 isn't really that complex, it's just taking KG and multiplying by the integral form of the van't Hoff equation. I have a suspicion that maybe the temperature term (x(4)) is causing some sort of error.
Any insight into what's going on would be greatly appreciated.
The code:
function [] = THT2()
clc;
clear all;
close all;
guess = 1000; % Initial guess of T(z=0).  in K
% Define constants
P = 300; % atm
Qf = 0.05713; % m^3 /s
XAf = 0.015;
XNf = 0.985/4;
XHf = 0.985*3/4;
Ac = 1; % m^2
L = 12; % m
Taf = 323; % K

gamma = 0.5; % 1/m
beta = -2.342; % (m^2*s*K) /mol
deltaG = -4.25*10^3; % cal/mol

deltaH = -1.2*10^4; % cal/mol
kminus10 = 7.794*10^11; % atm/s
EoverR = 2*10^4; % K
Rcal = 1.987; % cal /(K*mol)
Ratm = 8.205 * 10^-5; % m^3 atm /(K*mol)
NTf = (Qf*P)/(Ratm*guess);     % Total feed molar flux mol/s, N = QP/RT
NAf = NTf*XAf;  % feed molar flux of ammonia, mol/s
NHf = NTf*XHf;  % feed molar flux of hydrogen, mol/s
NNf = NTf*XNf;  % feed molar flux of nitrogen, mol/s
% Length vector
v0 = 0; % meters, top of column, entering the reactor
vf = 12; % meters, bottom of column, extiting reactor core
tspan = linspace(v0, vf, 7001);
initcons = [NNf, NHf, NAf, guess, guess]; % initial conditions entering reactor core (z=0)
% 1=N, 2=H, 3=A, 4=T, 5=Ta
params = [P Qf Ac Taf gamma beta deltaG deltaH kminus10 EoverR L Rcal Ratm NTf NAf NHf NNf guess];
[v,X] = ode15s(@(v,x)ODEfun(v,x,params),tspan,initcons);
subplot(2,1,1);
plot(v,X(:,4),'r', v, X(:,5),'b')
xlabel('Length (m)');
ylabel('Temperature (K)');
title('Temperature versus Length');
legend('Annulus', 'Reactor core');
subplot(2,1,2);
plot(v, X(:,3),'r', v, X(:,1),'b', v, X(:,2),'g')
xlabel('Length (m)');
ylabel('Flux (mol/s)');
title('Concentration vs Length');
legend('Ammonia','Nitrogen','Hydrogen');
end
function dX = ODEfun (v,x, params);
P       = params(1);
Qf      = params(2);
Ac      = params(3);
Taf     = params(4);
gamma   = params(5);
beta    = params(6);
deltaG  = params(7);
deltaH  = params(8);
kminus10= params(9);
EoverR  = params(10);
L       = params(11);
Rcal    = params(12);
Ratm    = params(13);
NTf     = params(14);
NAf     = params(15);
NHf     = params(16);
NNf     = params(17);
guess   = params(18);
%NN2 = x(1),  Molar flowrate of nitrogen
%NH2 = x(2),  Molar flowrate of hydrogen
%NA = x(3),   Molar flowrate of ammonia
%T = x(4),   Temperature in the reactor core
%Ta = x(5),  Temperature in the annulus
krev = (kminus10 * exp(-EoverR/x(4)));  % atm/s, solve for k reverse
KG = (exp(-deltaG/(Rcal*guess)));  % equilibrium constant
KH = KG*exp(-(deltaH/Rcal)*((1/x(4))-(1/guess))); % Van't Hoff equation
Q = (Ratm*x(4)/P)*(x(1)+x(2)+x(3)); % Volumetric Flowrate
PN = Ratm*x(4)*x(1)/Q; % Partial Pressure of Nitrogen


PH = Ratm*x(4)*x(2)/Q; % Partial Pressure of Nitrogen
PA = Ratm*x(4)*x(3)/Q; % Partial Pressure of Nitrogen
KHsquared = KH^2;
PHroot = PH^(3/2);
% Reaction rate equation, in mol / (s*m^3)
r = (krev/Ratm*x(4)) *((KHsquared*PN*PHroot/PA) - (PA/PHroot));
dX = zeros(5,1);
dX(1) = NNf -.5*(x(3)-NAf);
dX(2) = NHf - 1.5*(x(3)-NAf);
dX(3) = 2*Ac*r;  % Ammonia flux through reactor
dX(4) = -beta*r + gamma*(x(5) - x(4));  % Temperature profile through reactor core
dX(5) = gamma * (x(5) - x(4)); % Temperature profile through annulus = gamma * (Ta - T)
end

Best Answer

The problem with the code is that eventually (at about 10.47) it tries to project to a location with some negative x values. That leads to negative partial pressures, which leads to complex solutions.
You can add an odeset() options structure with 'NonNegative', 1:5 to try to force the entries to be non-negative. Unfortunately, the handling of the NonNegative option in the ode* routines does not prevent projected negative values from being passed in: it only clamps negative outputs.
I worked around this by adding
x = max(0, x);
at the top of your ode routine, which gets rid of the imaginaries, but at potential cost to the accuracy of the solution.
Related Question