Unable to meet integration tolerances
5 Ansichten (letzte 30 Tage)
Ältere Kommentare anzeigen
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
3 Kommentare
Antworten (1)
Bjorn Gustavsson
am 11 Aug. 2021
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?).
2 Kommentare
Bjorn Gustavsson
am 11 Aug. 2021
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...
Siehe auch
Kategorien
Mehr zu Ordinary Differential Equations finden Sie in Help Center und File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!