MATLAB: ODE45 returns NaN values

nanode45

I am having trouble identifying why ODE45 returns or values for the following model I'm trying to compute.
I have experimented with a variety of cases which I will list below, but I still have no idea why it will return or values. I have also presented the relevant functions and scripts I have been running below.
RHS is the function for which I'm applying ODE45. RHS is being computed over a grid (the unit square). I don't believe the RHS is part of the problem since there are no undefined values (such as dividing by 0 or ). If you believe so, please let me know why.
function uv = RHS(t,y)
% defines the right hand side of the Eulerian velocity field obtained
% from the analytical solution of the Stommel equation.
% Relevant equations and assumptions can be found from the paper:
% NUMERICAL SIMULATIONS OF PASSIVE TRACERS DISPERSION IN THE SEA, page 36
day = 11.6;
a = 1e7;
bta = 2*10^-11;
r = 1/(day*86400);
scalefac = 0.05;
es = r/(a*bta);
uv = [-(1-exp(-y(1)./es)-y(1))*pi^2.*cos(pi.*y(2))*scalefac;...
    (exp(-y(1)./es)./es -1)*pi*sin(pi.*y(2))*scalefac];
(The following are the cases that I've tried, to hopefully see the cause of this problem.)
Case 1: For some initial conditons, ODE45 computes some iterations but ends up with values.
tspan = [0 3];
y0 = [0.125; 0.265];
    [t,y] = ode45(@RHS, tspan, y0);
    figure(1)
    plot(y(:,1),y(:,2),'-o');
    title('Trajectory of particle');
    xlabel('x'); ylabel('y');
    axis([0 1 0 1]);
The result is a solution y of size where its 14th entry suddenly becomes . continues till the end (53rd entry). Interestingly, the 17th entry presents .
This result seems to be 'fixed' when I apply option tolerances as shown below. (For an arbitrary tolerance choice).
tspan = [0 3];
y0 = [0.125; 0.265];
opts = odeset('RelTol',1e-8,'AbsTol',1e-8);
    [t,y] = ode45(@RHS, tspan, y0, opts);
    figure(1)
    plot(y(:,1),y(:,2),'-o');
    title('Trajectory of particle');
    xlabel('x'); ylabel('y');
    axis([0 1 0 1]);
Case 2: For some initial conditions, the solution when applying option tolerances produces values.
tspan = [0 3];
y0 = [0.005; 0.005];
opts = odeset('RelTol',1e-8,'AbsTol',1e-8);
    [t,y] = ode45(@RHS, tspan, y0, opts);
    figure(1)
    plot(y(:,1),y(:,2),'-o');
    title('Trajectory of particle');
    xlabel('x'); ylabel('y');
    legend('t = [0 3]');
    axis([0 1 0 1]);
This time it goes straight to values by the 2nd iteration. The result y is size filled with values except the initial condition. And of course, when running this without option tolerances, the trajectories are computed just fine. This has become a problem because I need to compute many trajectories and I can't just go specifying certain conditions for many different points. There is a more general problem that I'm not seeing.
Case 3: The above two cases seem to be 'fixed' when I ignore ODE45's default timestepping, and instead manually state a time step, as in the following, I let increase from 0 by steps of . We try Case 2 with this:
tspan = [0:0.01:3];
y0 = [0.005; 0.005];
opts = odeset('RelTol',1e-8,'AbsTol',1e-8);
    [t,y] = ode45(@RHS, tspan, y0, opts);
    figure(1)
    plot(y(:,1),y(:,2),'-o');
    title('Trajectory of particle');
    xlabel('x'); ylabel('y');
    legend('t = [0 3]');
    axis([0 1 0 1]);
It computes the trajectory just fine, and the same results occur with Case 1. With or without option tolerances. It seems to be that the adaptive timestep used by ODE45 may be the problem, but I eventually found another point that began to cause problems.
Case 4: Some initial positions produce NaN values when specifying a timestep.
tspan = [0:0.01:3];
y0 = [0.7635; 0.0935];
opts = odeset('RelTol',1e-8,'AbsTol',1e-8);
    [t,y] = ode45(@RHS, tspan, y0, opts);
    figure(1)
    plot(y(:,1),y(:,2),'-o');
    title('Trajectory of particle');
    xlabel('x'); ylabel('y');
    legend('t = [0 3]');
    axis([0 1 0 1]);
The result y of size seems to be computed just fine until the 143rd entry, where it begins to present values until the end. Where the last, 151st entry is . This occurs with or without option tolerances and with or without explicit timesteps.
It seems that the above 4 cases rule out each of the possibilities to do with timestepping, options and that there is definitely something that I'm overlooking. Any help is appreciated.

Best Answer

ode45 does not simply compute the function at "every" point and do something like trapazoidal rule. Instead, it starts from a particular set of conditions, and does careful probing around with different times and different initial conditions, hoping the calculated error tolerences come out low enough that it considers that it has been successful in projecting over a range of values. If the error it calculates comes out small enough, it accepts the projected time and boundary conditions and works forward from there. If, though, the error it calculates comes out too large, then it declares the step failed, and goes back to where it was and takes a smaller step. As steps are successful, it tends to grow the step size, and as steps fail, it shrinks the step size. In this way, it can find rapid changes in slope and be more careful in those cases; however, this algorithm can also fail to notice a sufficiently narrow peak if the step size happens to be large enough at the time and the value on the other side of the peak is similar enough to what it is expecting.
The locations it samples at to decide whether to accept or reject are according to mixing coefficients carefully calculated from theory -- there is no randomness or arbitrary choices involved other than the initial step size. It uses the theoretically best coefficients, according to certain assumptions about continuity. There is no room for Mathworks to "tune" the coefficients for choosing where to project to; these are the best locations.
However, since it is doing projections and cannot know ahead of time whether any one projection is reasonable (since the ode function is a "black box" to it, not something it can examine the mathematics of), on projections that are going to end up failing, it can end up asking about the curve values at combinations of boundary conditions that are in locations steeply off of the real curve -- locations that you wouldn't think to ask about if you were taking smaller steps, but which are reasonable to ask about when you realize that the steps are designed to process the curve as rapidly as is compatible with being able to keep reasonable error bounds and this is just probing hoping this direction might be a good one. Projection is successful most of the time, but sometimes you can be close enough to boundary conditions that the function value has to suddenly be large.
Your function happens to have an exp(-y(1)./es), and if y is large enough and negative, then the -y(1)/es can be large and positive and the exp() can go to infinity, leading to an infinite uv result. Because exp() is non-linear, it does not necessarily take much curve projection to go from "hiking along smoothly" to "fell off of the cliff!" if you take too ambitious of a step.
When you generate the inf, then it gets cycled into forming the next set of input coefficients, typically resulting in inf in several input coefficients, and it then being common to end up subtracting or dividing infinities, and both of those lead to nan. And once you have a nan, everything goes out the window, "nan poisoning".
Now, in such cases where you overflow, if somehow you stopped short of infinity, but with a sufficiently large finite value, then the calculation of error would proceed, and MATLAB would see the error was too big and would fail the step, and reduce the stepsize and try again, and the step that would have produced the infinity would be thrown away as irrelevant. The infinity is not a problem because it is infinite: it is a problem because it typically leads to multiple infinities that lead to calculations between infinities that give nan, and nan is a problem. If the function produced an infinity at just the right point or if ode45() somehow noticed it in time, then the step would just be rejected and the algorithm would be back up and try again as designed.
So, what you can do is detect the infinity yourself and return a large enough value to "poison" the calculation to be reasonably sure that the error generated would be too large and that the code will back up. In the case of your function, you can do that by adding the lines:
if any(~isfinite(uv))
    uv = [1e10;1e10] .* sign(uv) ;
end
This keeps the result finite but large enough to be over error tolerances (usually.) This line of code is not fool-proof though; in other equations it might have to be adjusted for different error values.
So.. detect the infinity, send out a large value to create a situation where the error will be large enough that ode45 will back up, but without having sent out an actual infinity that can lead to infinity mixes that lead to nan poisoning.
Related Question