Unable to meet integration tolerances

I have a system of 13 stiff ODEs (DAEs) with an initial value for 13 variables. I have tried the attached code with different tolerance levels but I still get the same warning and I am struggling to fix this warning. This problem is similar to the Robertson's ODE system. I read some suggestions on this forum where they say the integral needs to be split. I am not sure if that is necessary for my system. I have played around with the time interval of integration and even changing the constants 'k'. None of these seem to work and I look forward to some productive feedback! Thank you.

Here's the warning:

Warning: Failure at t=1.905321e-19. Unable to meet integration tolerances without reducing the

step size below the smallest value allowed (3.851860e-34) at time t.

> In ode15s (line 730)

In mkm_solver (line 18)

The ODEs of the system are:

Conservation eqn:

My code is:

function hb1dae
M = [1 0 0 0 0 0 0 0 0 0 0 0 0
   0 1 0 0 0 0 0 0 0 0 0 0 0 
   0 0 1 0 0 0 0 0 0 0 0 0 0
   0 0 0 1 0 0 0 0 0 0 0 0 0 
   0 0 0 0 1 0 0 0 0 0 0 0 0 
   0 0 0 0 0 1 0 0 0 0 0 0 0
   0 0 0 0 0 0 1 0 0 0 0 0 0
   0 0 0 0 0 0 0 1 0 0 0 0 0
   0 0 0 0 0 0 0 0 1 0 0 0 0 
   0 0 0 0 0 0 0 0 0 1 0 0 0 
   0 0 0 0 0 0 0 0 0 0 1 0 0
   0 0 0 0 0 0 0 0 0 0 0 1 0
   0 0 0 0 0 0 0 0 0 0 0 0 0];
tspan = [0 10];
y0 = [1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0];
options = odeset('Mass',M,'RelTol',1e-10,'AbsTol',[1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10], ...
   'Vectorized','on');
[t,y] = ode15s(@f,tspan,y0,options);
y(:,2) = y(:,2);
figure;
semilogx(t,y);
legend('*','A','L','D','E','G','I','J','K','CO','OH','H','H2O');
ylabel('1e4 * y(:,2)');
title('Robertson DAE problem with a Conservation Law, solved by ODE15S');
xlabel('This is equivalent to the stiff ODEs coded in HB1ODE.');
end
function out = f(t,y)
 %Forward rate constants
  kf0 = 4.25e+18;
  kf1 = 1.09e+19;
  %kf2 = 2.59e+15;
  kf3 = 1.58e+21;
  kf4 = 2.87e+16;
  %kf5 = 5.29e+21;
  %kf6 = 4.74e+13;
  kf7 = 4.35e+17;
  kf8 = 2.79e+14;
  kf9 = 8.38e+13;
  kf10 = 3.79e+19;
  kf11 = 8.34e+12;
  kf12 = 1.00e+13;
  kf13 = 2.24e+52;
  kf14 = 2.51e+23;
  kf15 = 3.58e+17;
  kf16 = 1.00e+23;          
  kf17 = 1.00e+23; 
 %Reverse rate constants
  kb0 = 5.8e+13;  
  kb1 = 1.41e+18;
  %kb2 = 3.38e+12;
  kb3 = 1.66e+15;
  kb4 = 3.41e+13;
  %kb5 = 7.14e+16;
  %kb6 = 1.47e+09;
  kb7 = 1.25e+06;
  kb8 = 8.32e+06;
  kb9 = 1.05e+01;
  kb10 = 2.89e+28;
  kb11 = 1.55e+21;
  kb12 = 3.52e+13;
  kb13 = 6.72e+18;       
  kb14 = 3.26e+15;
  kb15 = 6.96e+18;
  kb16 = 1.00e+13;        
  kb17 = 1.00e+13;      
  pa = 1.0;               %partial pressure of A
  pk = 0.0;
  pCO = 0.0;
  pH2O = 0.0;
%These are the ODEs:
%yA
    y1 = (kf0*pa*y(1,:).*y(1,:))-(kb0*y(2,:))-(kf1*y(2,:).*y(1,:))+(kb1*y(3,:).*y(11,:));
%yL
    y2 = (kf1*y(2,:).*y(1,:))-(kb1*y(3,:).*y(11,:))-(kf3*y(3,:))+(kb3*y(8,:).*y(10,:))-(kf4*y(3,:).*y(1,:).*y(1,:))+(kb4*y(4,:).*y(12,:));
%yD
    y3 = (kf4*y(3,:).*y(1,:).*y(1,:))-(kb4*y(4,:).*y(12,:))-(kf7*y(4,:))+(kb7*y(7,:).*y(10,:))-(kf8*y(4,:).*y(1,:))+(kb8*y(5,:).*y(12,:));
%yE
    y4 = (kf8*y(4,:).*y(1,:))-(kb8*y(5,:).*y(12,:))-(kf9*y(5,:))+(kb9*y(6,:).*y(10,:).*y(1,:));
%yG
    y5 = (kf9*y(5,:))-(kb9*y(6,:).*y(10,:).*y(1,:))-(kf10*y(6,:).*y(12,:))+(kb10*y(7,:));
%yI
    y6 = (kf7*y(4,:))-(kb7*y(7,:).*y(10,:))+(kf10*y(6,:).*y(12,:))-(kb10*y(7,:))-(kf11*y(7,:).*y(12,:))+(kb11*y(8,:).*y(1,:).*y(1,:));
%yJ
    y7 = (kf3*y(3,:))-(kb3*y(8,:).*y(10,:))+(kf11*y(7,:).*y(12,:))-(kb11*y(8,:).*y(1,:).*y(1,:))-(kf12*y(8,:).*y(12,:))+(kb12*y(9,:).*y(1,:));
%yK
    y8 = (kf12*y(8,:).*y(12,:))-(kb12*y(9,:).*y(1,:))-(kf13*y(9,:))+(kb13*pk*y(1,:));
%yCO
    y9 = (kf3*y(3,:))-(kb3*y(8,:).*y(10,:))+(kf7*y(4,:))-(kb7*y(7,:).*y(10,:))+(kf9*y(5,:))-(kb9*y(6,:).*y(10,:).*y(1,:))-(kf15*y(10,:))+(kb15*pCO*y(1,:));
%yOH
    y10 = (kf1*y(2,:).*y(1,:))-(kb1*y(3,:).*y(11,:))-(kf16*y(11,:).*y(12,:))+(kb16*y(13,:).*y(1,:));
%yH
    y11 = (kf4*y(3,:).*y(1,:).*y(1,:))-(kb4*y(4,:).*y(12,:))+(kf8*y(4,:).*y(1,:))-(kb8*y(5,:).*y(12,:))-(kf10*y(6,:).*y(12,:))+(kb10*y(7,:))-(kf11*y(7,:).*y(12,:))+(kb11*y(8,:).*y(1,:).*y(1,:))-(kf12*y(8,:).*y(12,:))+(kb12*y(9,:).*y(1,:))-(kf16*y(11,:).*y(12,:))+(kb16*y(13,:).*y(1,:));
%yH2O
    y12 = (kf16*y(11,:).*y(12,:))-(kb16*y(13,:).*y(1,:))-(kf17*y(13,:))+(kb17*pH2O*y(1,:));
%yf
    %y13 = (kf9*y(5,:))-(kb9*y(6,:).*y(10,:).*y(1,:))+(kf11*y(7,:).*y(12,:))-(kb11*y(8,:).*y(1,:).*y(1,:))+(kf12*y(8,:).*y(12,:))-(kb12*y(9,:).*y(1,:))+(kf13*y(9,:))-(kb13*pk*y(1,:))+(kf15*y(10,:))-(kb15*pCO*y(1,:))+(kf16*y(11,:).*y(12,:))-(kb16*y(13,:).*y(1,:))+(kf17*y(13,:))-(kb17*pH2O*y(1,:))-(kf0*pa*y(1,:).*y(1,:))+(kb0*y(2,:))-(kf1*y(2,:).*y(1,:))+(kb1*y(3,:).*y(11,:))-(kf4*y(3,:).*y(1,:).*y(1,:))+(kb4*y(4,:).*y(12,:))-(kf8*y(4,:).*y(1,:))+(kb8*y(5,:).*y(12,:));
%conservation equation below or the above y13 can be considered. The below
%eqn makes the mass matrix singular! 
    y13 = y(1,:) + y(2,:) + y(3,:) + y(4,:) + y(5,:) + y(6,:) + y(7,:) + y(8,:) + y(9,:) + y(10,:) + y(11,:) + y(12,:) + y(13,:)- 1;
out = [y1; y2; y3; y4; y5; y6; y7; y8; y9; y10; y11; y12; y13];
end

[Bjorn Gustavsson]

It might help if you use the code-markup icon (the left-mose icon in above "code" in the editing frame) to make your code more readable, if possible it might also be easier for people to get what your DAEs model if you can present the mathematical equations (the capital Sigma allows you to insert latex-formated equations (reasonably exhaustive latex-support)), then we would not read and interpret the mathematical meaning from your code.

[Keerthan Raghavendra Rao]

Awesome! Thank you. I have updated the post with eqns and code.

[Keerthan Raghavendra Rao]

After more investigation, it looks like some of the y values become negative to offset other y values>1 to satisfy the conservation equation (y13). I am not sure why it would work for a three-equation problem like the Robertson ODE and not for this system. Is there a work-around for this?

Answer 1

[Bjorn Gustavsson]

You might get enough control of ode15s by telling it to keep the solutions nonegative, try setting the "NonNegative" field to a vector ones(12,1) - but then you'll have to abandon the conservation equation (which should be trivially satisfied? If your reaction-rates are correctly implemented and your losses of some theta matches a source of one of the other they should cancel out in the sum of the components, right?).

[Keerthan Raghavendra Rao]

Thank you, I tried the NonNegative option but ode15s ignores it! This option doesn't seem to be available for stiff solvers. Any other suggestions?

[Bjorn Gustavsson]

In the help to odeset it claimed that ode15s would listen to that - but it couldn't be combined with a mass-matrix. So to use that you might have to add in the ode for y13 and remove the mass-conservation equation, outside of that I have no other suggestions right now...