Change step size in MATLAB for ODE45

Hi, is it possible to avoid this:

    Warning: Failure at t=-3.000000e+00.  Unable to meet integration tolerances without reducing the step size below the smallest value
    allowed (7.105427e-15) at time t. 
    > In ode45 (line 308)
      In S_plot (line 6)

I would need to go far lower, type 10^-104

[Torsten]

Depends on the ODE you are trying to solve.

What do you mean by "I would need to go far lower, type 10^-104" ?

[Sergio Manzetti]

The ODE I try to solve is this:

syms h x
fun = @(x,y)[y(2);y(3);(y(1)*(2.56E-20-x^2))/1.1728E-102];
y0 = [1 0 1]; %y, y', y''
xspan = [-3 3];
[X,Y] = ode45(fun,xspan,y0);
plot(X,Y(:,1)),,% X,Y(:,2),X,Y(:,3),X,Y(:,1).* (2 + X .^ 2));

by far lower I mean the warning message:

Warning: Failure at t=-3.000000e+00.  Unable to meet integration tolerances without reducing the step size below the smallest value
allowed (7.105427e-15) at time t. 
> In ode45 (line 308)
  In S_plot (line 6)

"allowed (7.105427e-15) at time t. "

It must most likely go beyond that limit.

Answer 1

[Jan]

If the step size controller of ODE45 reaches 7e-15, the integration will take many years of processing time: Remember, that a day has less than 1e5 seconds only and even if ODE45 would get 1 million iterations per second, the number of steps is still huge.

Such a tiny step size is a secure indicator of either a discontinuity, a pole or a stiff equation. Reducing the step size is not a valid option, because this increases the run time (see above) and the accumulated rounding errors due to the massive number of steps. You have to examine the function mathematically instead. Maybe a stiff solver is better or avoiding the pole.

[Sergio Manzetti]

Thanks Jan Simon, I will try some of these

https://se.mathworks.com/help/matlab/math/solve-stiff-odes.html

[Steven Lord]

I don't know where this ODE came from but I would recommend double-checking that you have correctly copied it from your source or have correctly derived it from your problem. Even a stiffer solver probably isn't going to help you. Look at the ODEs you're trying to solve.

fun = @(x,y)[y(2);y(3);(y(1)*(2.56E-20-x^2))/1.1728E-102];

Another way of writing that third component is:

(y(1)*(2.56e-20-x^2))*(8.5266e101)

If y(1) is even slightly different from zero or x is slightly different from sqrt(2.56e-20), this component is going to take off like the Starship Enterprise going to warp speed. For your actual problem, the initial conditions are y(1) = 1 and x = -3 and this component turns out to be -7.6739e102.

That value of this component will, at later time steps, feed back into the expression for the derivative of the other components. After only a couple time steps the values of these components would exceed what can be represented in double precision and you'd start getting Inf and NaN values in the components.

[Sergio Manzetti]

Thanks indeed there was an innecessary use of high exponentials. I have changed the ODE now and it appears normal. However, I can't really fathom (being new to numerical analysis) why we talk about a time component in the error message, when the first dimension is x, and the second is y. Where is time in all this?

Thanks!

PS: Nice with Starship Enterprise! I take Picard would be impressed with this function!

All the best

[Steven Lord]

The warning message talks about t rather than time or x because of the way the help text and documentation describes the input. From the documentation for ode45:

"[t,y] = ode45(odefun,tspan,y0), where tspan = [t0 tf], integrates the system of differential equations y'=f(t,y) from t0 to tf with initial conditions y0."

and

"The function dydt = odefun(t,y), for a scalar t and a column vector y, must return a column vector dydt of data type single or double that corresponds to f(t,y)."

While it may be possible in a lot of cases to actually figure out the name the function the user passed into ode45 internally uses for its first input argument, there are also cases where it wouldn't be possible. One example where it wouldn't be possible is if the function accepts varargin. The error message can't be always in sync with the code the user wrote and passed into ode45, but it can be in sync with the terminology used in the help text and documentation.