ODE problem : Unable to meet integration tolerances

Hello,

I'm trying to solve a differential non linear equation with very high values (each constant is like 10^-20 or 10^20 or even more). And i got problems of tolerances. I've tried that :

options=odeset('RelTol',1e-30,'AbsTol',1e-30)
xvalues=30:200;
[x,y]=ode45('fun',xvalues,1, options)

Where the function is, as i told you, with very big values (for example, a term with exp(x)...)

""Warning: RelTol has been increased to
2.22045e-14. 

Warning: Failure at t=3.000000e+01.  Unable
to meet integration tolerances without
reducing the step size below the smallest
value allowed (1.136868e-13) at time t. ""

If you have any idea, help me !!!

Thank you

Fred

[Walter Roberson]

Could I ask you to show your ODE in symbolic form? I have not done much work with transforming back from ode45 calls into formulae and I suspect I got it wrong when I did.

Could you confirm your boundary condition is that the derivative evaluated at the initial coordinate, xvalues(1), 35, should be 1? The documentation is not explicit about the fact that the initial (t,y) will have t be the first time and y be the initial condition, but my tests confirm it in practice. (An alternative reading would be that the initial values should be for an implicit time of 0, which would seem less practical, but I wanted to be sure I got the theory right.)

[Fred Dugli]

I haven't written it another way than this one I can try to make it more clear (sadly it looks i can't make LaTeX on this forum :( ) :

dy/dx= -C1*x^(-5/2)/sqrt(C2+C3/x) * ( y^2 - (1-y)*x^(3/2)*e(-x)/C4 )

With C1, C2, C3, C4 constants. I hope it can help you

Boundary condition is that y=1 for x small (for example y(x=35)=0 ; in my problem, I study x from 30 or 35, as u can see with my code, y don't change between 30 and 35)

But I think it worked well since the result really looks like what I have to find (the equation was already solved before, but I want to test influence of some parameters). From this example I really think the first alternative you are speaking about must be true.

If you want more informations, don't hesitate to ask me, i'll to my best to answer you.

Answer 1

[Star Strider]

I would increase the tolerance (to perhaps 1E-8 or more) rather than decrease it, but that may not solve your problem.

The warning is that at t=30, your ODE is most likely encountering a singularity. First, plot it to see what it is doing, then start it again at some value greater than t=30 (perhaps t=35) and see what it does, then experiment further to get the result you want. You may have to do a piecewise integration. That may not work in the end (you may not be able to integrate it beyond t=30 regardless), but unless you experiment with it, you will never know.

[Fred Dugli]

Thank you very much for your answer. I tried to change from 1E-5 to 1E-30 for tolerance, without success. I also tried to start from t=35 and t=40 but still same problem.

I found an approximate of the same function by neglecting some terms and this time i haven't got this problem but answer is Xe =

1.000000000000000                     
              NaN +               NaNi
              NaN +               NaNi
              NaN +               NaNi
              NaN +               NaNi
              NaN +               NaNi
              NaN +               NaNi
              NaN +               NaNi
              NaN +               NaNi
              ................. etc etc

Here is the code of the approximate function :

    function y = Approx( y , x )
    alpha2=7.1472*10^(-14);
    H=(4.5198*10^(-11)*sqrt(0.15+2.4047*1/x))*1/(x^(3/2));
    n= 3.3206e+08 * 1/(x^3);
    S=2.19216e-16*exp(x)*1/(x^(3/2));
    y=-alpha2*n*(y^2-(1-y)/S)*1/(x*H);
    end

And my main is :

options=odeset('RelTol',1e-40,'AbsTol',1e-40);
xvalues=35:160;
[x,y]=ode45('Approx',xvalues,1, options)

I have also changed initial point and tolerance, but still same problem. It must be due to the extremely big figures As you can see in the Approx function...

But I have no idea how what I could try ...

[Star Strider]

I vectorised your ODE function (necessary because you are likely doing array rather than matrix operations):

    function y = Approx( y , x )
    alpha2 = 7.1472E-14;
    H = 4.5198E-11*sqrt(0.15+2.4047./x)./(x.^(3/2));
    n= 3.3206e+08 ./ (x.^3);
    S=2.19216e-16*exp(x)./(x.^(3/2));
    y=-alpha2*n.*(y.^2-(1-y)./S)./(x.*H);
    end

and:

options=odeset('RelTol',1e-8,'AbsTol',1e-8);
xvalues=[0 160];
[x,y]=ode45(@Approx,xvalues,1, options);
%
figure(1)
plot(x, y)
grid

and found that it only appears to exist on x from 0 to 1. I cannot get it to integrate and produce other than NaN values outside those limits. I haven’t analysed it, but the extreme values of your constants likely are the reason.

[Fred Dugli]

Thank you for your help.

I know the solution must be more or less this: y going down from 1 (almost) to 0 (or at least less than 10e-2) when x is going from 35 to 160. I'll check my values, to be sure my constants are the good ones.

[Star Strider]

My pleasure.

Usually, NaN values are the result of attempting a ‘0/0’ or ‘Inf/Inf’ operation. It’s not obvious to me by looking at your function how it might be generating them outside the (0,1) interval.

[Walter Roberson]

Star Strider, note you ran the ode from time 0 rather than from time 35, but you did not change the initial condition from being 1.

[Star Strider]

Neither did Fred.

Answer 2

[Fred Dugli]

I don't really know why, i made a few changes (like inverting x and y, changing tolerance). And now it works !!!!!!!!!!!!!!!!

I get the nice curve i expected. Don't ask me why it now works ...

In case people go on this post and have a similar problem (hope it could help them). Here is the code that worked for me :

function y2 = Approx2( x , y )
    alpha2 = 7.1472E-14;
    H = 4.5198E-11*sqrt(0.15+2.4047./x)./(x.^(3/2));
    n= 3.3206e+08 ./ (x.^3);
    S=2.19216e-16*exp(x)./(x.^(3/2));
    y2=-alpha2*n.*(y.^2-(1-y)./S)./(x.*H);
    end

And the main is :

options=odeset('RelTol',1e-12,'AbsTol',1e-12);
xvalues=[35, 160];
[x,y]=ode45(@Approx2,xvalues,1, options);
y
%
figure(1)
plot(x, y)
grid

Thank you all for your help, and time for reading/answering me (and sorry for my mistakes in english !), if I find more details about why it now works, i'll let you know

Answer 3

[Walter Roberson]

You can never get a relative tolerance of 1E-30 with double precision floating point numbers. The absolute least relative tolerance that could make sense is eps(1) = 2.22044604925031e-16 as that represents the relative difference between adjacent double precision floating point numbers.

An absolute tolerance of 1E-30 might make sense in some situations, but only when the values being computed are in the range of roughly 1E-18 to 1E-24.

If you have large values such as 1E+20 in your calculation, to get down below 1E+4 (eps(1) times the large value) you need to be dividing by large values or you need to be hoping for extreme cancellation of values.

I think you need to be re-thinking your tolerances.

[Fred Dugli]

Thank you for your answer, I knew tolerance was too small since it's less than eps, but I just wanted to be sure to have the minimum tolerance possible. It's not possible to change the value of eps ?? I assume no ... even if that surprises me, it's not rare for practical calculus to have values under eps.

And i think it would be hard to divide anything since for example my S variable takes values from 10e3 to 10e50 ...