Warning: Matrix is singular, close to singular or badly scaled.

Hi,

I am writing a progrem with Matlab for the solution of a system of differential equations.

The program is:

% Orto-xylene oxidation in an adiabatic PFR.
% ethylbenzene -> Anhydrated-xylene + H2
% Water as an inert in the system
global y0 F0 deltaHR0 P0 ni n 
y0(1) = 0.0769; % ethylbenzene
y0(2) = 0.; % styrene
y0(3) = 0.; % hydrogen
y0(4) = 0.923; % water
Feed  = 1105.; % F01 + FH2O Feed rate [mol/s]
F0 = Feed .* y0; % Inlet feed rate of components [mol/s]
Tin = 925.;
n = 4.;        % Number of reacting compound
deltaHR0 = 124850.; % heat of reaction at standard conditions[J/mol]
P0 = 240000.; % pressure [Pa]
W = 25400.; % reactor weight [kg]
ni = [-1 1 1 0]; % Stoichiometric matrix
% Initial conditions (feed and T0)
u0 = F0;
u0(n+1) = Tin;
wspan = [0 W];
[w_adiab,u] = ode15s(@dfuns,wspan,u0);
conv_adiab = 1 - u(:,1)/F0(1);
T = u(:,n+1);

%function file dfuns 
function f = dfuns (w,U)
global y0 F0 deltaHR0 P0 ni n 
y = U(1:n)/sum(U(1:n)); % molar fractions
T = U(n+1);
k = 3.46* 10e8.* exp(-10980./T); %Needs to be in Pascal 
K = 8.2.* 10e11.* exp(-15200./T); %Needs to be in Pascal 
rm  = k.* (y(1)* P0- (y(2).* y(3).*(P0^2))./ K);
f(1:n) = -ni(1:n).* rm;
Cp1 = 37.778 + 87.940e-3.* T + 113.436e-5.* T^2; % ethylbenzene 
Cp4 = 36.540 - 34.827e-3.*T + 11.681e-5.*T^2; % water 
deltaCp = 0.;
deltaHR = deltaHR0;
f(n+1) = (1./((F0(1)+F0(4)).* (y(1).* Cp1 + y(4).* Cp4 +y(1).* deltaCp)*...
    (1 - U(1)/F0(1)))).*-deltaHR.* rm;
f = f'; % transformation to column vector

I get a warning stating:

Warning: Matrix is singular, close to singular or badly scaled. Results may be inaccurate. RCOND = NaN. 
and then, after many other warnings:

Warning: Failure at t=0.000000e+00.  Unable to meet integration tolerances without reducing the step size below the
smallest value allowed (7.905050e-323) at time t.

I think there might be a problem with the dimensions of some of the vectors I defined but i struggle to see which one and why.

Can someone please help me? Thank you very much

[Torsten]

f(n+1) becomes -Infinity in the first call to dfuns since (1-U(1)/F0(1))=0.

Best wishes

Torsten.

Answer 1

[Walter Roberson]

You probably have a division by 0. To track it down, give the command

dbstop if naninf

and run the program. It should stop at the point where inf or nan was produced.

Note: instead of using global, it is advised that you parameterize your calls

[Serena Casanova]

Dear Walter,

thank you for answering.

Unfortunately i'm not so much into how matlab approaches this kindo of problem so I don't know how to interpret this analysis produced by the command you suggested to make:

function varargout = dataviewerhelper(whichcall, varargin)
%DATAVIEWERHELPER Helper functions for Workspace, Variable Editor, and other tools

%   Copyright 2008 The MathWorks, Inc.

switch whichcall
    case 'upconvertIntegralType',
        varargout = {upconvertIntegralType(varargin{:})};
    case 'isUnsignedIntegralType',
        varargout = {isUnsignedIntegralType(varargin{:})};
    otherwise
        error(message('MATLAB:dataviewerhelper:unknownOption'));
end

%********************************************************************
function converted = upconvertIntegralType(value)
converted = value;
if ~isfloat(value)

       if isa(value, 'uint8')
          converted = uint16(value);
      elseif isa(value, 'uint16')
          converted = uint32(value);
      elseif isa(value, 'uint32')
          converted = uint64(value);
      elseif isa(value, 'uint64')
          converted = uint64(value);
      elseif isa(value, 'int8')
          converted = int16(value);
      elseif isa(value, 'int16')
          converted = int32(value);
      elseif isa(value, 'int32')
          converted = int64(value);
      elseif isa(value,'int64') 
      %Necessary for referring to the parent class when passing through the
      %JAVA interface
          converted = int64(value);
      end

end

%********************************************************************
function unsigned = isUnsignedIntegralType(value)
unsigned = false;
if ~isfloat(value)
    unsigned = isa(value, 'uint8') || isa(value, 'uint16') || ...
        isa(value, 'uint32') || isa(value, 'uint64');
end

The line with the arrow finally is:

    if ~isfloat(value)

[Walter Roberson]

Umm, I think I would

dbcont

at that point as it appears to have to do with the workspace browser.

[Serena Casanova]

If it was the case of dividing by zero, where do you think it could be?

In here:

   y = U(1:n)/sum(U(1:n)); % molar fractions

the sum of U(1:n) is always positive, even for the initial conditions.

Here:

rm  = k.* (y(1)* P0- (y(2).* y(3).*(P0^2))./ K); 
K is varies exponentially with the temperature- no zero

And:

    f(n+1) = (1./((F0(1)+F0(4)).* (y(1).* Cp1 + y(4).* Cp4 +y(1).* deltaCp)*...
        (1 - U(1)/F0(1)))).*-deltaHR.* rm;

Again the denominator is positive because at least

y(4).*Cp4

is always positive.

I'm lost :-S

[Walter Roberson]

I suggest tossing in an OutputFcn option, probably together with OutputSel; search down in the documentation near "odeplot".

For example if the outputs turned out to be essentially identical then as the function attempted to establish a gradient at time 0, you could end up with a 0 determinant (infinite rcond)

[Serena Casanova]

Thank you very much,

In this case, where do you think the problem with my code could be? I understand what you mean but I don't see how to solve it in my specific case.

[Walter Roberson]

I do not know. First you need to verify whether the problem is happening in every run or if it is happening when it goes to integrate the results of several runs, and if it is the result of several runs then you need to compare the values to get an idea of how close the outputs are.