ode23tb and non-linear curve fitting

Dear All,

I think that I found an interesting bug-like thing while I was using ode23tb (MATLAB 2016b).

According to the Documentation :

[t,y] = ode23tb(odefun,tspan,y0,options) %where t is the same as tspan,if tspan has more than 2 elements

My tspan has 1185 elements.

I use the ode23tb in a function,which I want to minimize in a non-linear curve fitting procedure. I use fminsearchcon to find the minimum. There is no problem to find a nice-looking minimum.

The problem first occured when I called the function jacobianest in order to calculate the the confidence intervals. This function calls multiple times the fit function and therefore ode23tb, but with different (fit) parameters because the jacobianest calculates the partial derivatives. Near the end of the jacobianest routine, ode23tb returns with t that has only 759 elements while the tspan has still 1185 elements. It seems that a special configuration of parameters for ode23tb reduces the dimension of the returning values. Of course, later this results in an error because of dimension mismatch. This is very weird and I could not find the problematic part in my code.

When I got this dimension loss of t, I saved the paramteres from the jacobianest , and I just plugged them into ode23tb without any fitting routine and it gave the same error.

The order of the parameters were correct but one obtained a negative sign, which is not realistic, though mathematically it should not be an issue.

Do you have any idea what is happening here?

Just for clarity: I have 5 fit parameters, which define the solution of the ode23tb, none of them is the initial parameter. I have three equations, the ode23tb returns three columns. Then, another function makes a single column out of these three. This is the vector that I want to fit to my data. I guess the issue is not with the fitting routine but with ode23tb.

[Ameer Hamza]

I guess the shortening of vector 't' is accompanied by a warning in the command window stating that the ode23tb is unable to meet integration tolerance or something similar. My guess is that your system of ODE is unstable at that point. MATLAB ode functions return partial solutions if the response of ODE becomes unstable.

[Walter Roberson]

That would be my suspicion too.

[Tim]

Thank you for the comments. Indeed, I got the warning:Warning: Failure at t=4.164915e-10. Unable to meet integration tolerances without reducing the step size below the smallest value allowed (1.479675e-24) at time t.

I overlooked this, due to the other errors that I get.

I reduced the step size in tspan drastically (down to 10 points) and I still get the same failure. Could you please comment on this?

My original problem is still the confidence interval calculation, and it seems to me that jacobianest function drives the ode23tb to the territory where the ODE is unstable. Do you know other methods where I can get confidence intervals in such a complicated case?

[Walter Roberson]

The step size in your tspan is not the problem; ode23tb calculates where-ever it feels like calculating to find locations that it feels comfortable are within error tolerance of the previous location, and it uses those where-ever-it-feels-like points to predict the values at the tspan values you provide. You are not using a fixed step solver, so the solver does not bother solving at the tspan locations, just close enough to them that it is happy with the error estimates.

What would make a difference is reducing your maximum t below 4.164915e-10 to avoid the discontinuity. And if that is not possible, then you have a problem.

[Ameer Hamza]

Step-size in tspan has no effect on the internal working of the ode functions. These functions use their internal algorithm to calculate a solution, and then use interpolation to return the solution vector at points in tspan.

You may need to analyze your system of ODE to find out where they can become unstable and then set some bounds on the parameter passed to jacobianest such that they remain inside the stability region of ode function.

[Are Mjaavatten]

I think the comments by Walter and Ameer are very relevant. When I encounter the "unable to meet integration limits" error in solving an ODE system it generally means that I hit a singularity. This may be a real singularity in my problem, but usually it means that I have made some error. Since the error comes very early in your case, I would first run odefun(0,y0) and see if all derivatives seem reasonable. If they seem OK, I would follow Walter's suggestion and use tspan = [0, 4.16e-10], or shorter if this crashes. I would then examine the time series. One or more components may start to diverge, which may give you a hint as to what is happening.

[Walter Roberson]

A difficulty that I have encountered several times recently with ode45, is a situation where ode45 is projecting points (there is a mathematically optimal way of projecting for each different ode* routine) to calculate error, and it happens to project to a location where the return value from odefun includes nan (typically because the previous point produced inf). When that happens, the ode*() routine fails and cannot continue. The ode*() functions can survive +/- inf in themselves, but not nan -- but if you generate an inf then it is highly likely that a nan will soon follow.

The expection of the ode* routines is that if they project too far, that the return values will be sufficiently different from the current values that the calculated error is large, which the ode* routines would detect and reject the location as a candidate location, shrinking back to a less ambitious step, Over-projecting is normal for the ode*() routines, and as long as the odefun does not return nan recovery is fine. But the logic does not permit recovery when a nan is encountered.

The work-around for this is to put into your code tests for the calculated derivatives being nan or inf, and if that is found, emitting something large but finite that will "surely" lead to the ode*() routine deciding that the error is too much and that the step should be rejected.

[Tim]

Thank you to all for the thoughtful replies. I have learnt a lot.

The problem with Walter`s and Are`s solution with the shorter tspan is that I loose the purpose of this calculation because the fit happens in the full range of the original tspan, so I need to know what is happening in the full range.

Regarding Ameer`s idea: I found at which iteration of jacobianest the system is unstable. I wrote a condition in the jacobianest function saying that if there is a shortening of t, then the difference which gives a second order estimate, is not evaluated, instead it goes to the next index of that loop. With this ``dirty``` solution I managed to run the jacobianest function fully and finally get the confidence intervals. Though, I doubt that it is the correct way to proceed.

Regarding Walter`s last idea: do you mean to identify the Jacobian of the ODE system, check when is inf or nan and then to modify somehow the ode23tb solver itself?

[Ameer Hamza]

Tim, The idea suggested by Walter is to filter out Inf and NaN values inside your odefun, so that the ode23tb never get to see those. Something like this

function dydt = odefun(t,y)
    % your current code which calculates dydt
    
    dydt(isinf(dydt)) = 1e50; % 1e50 or a very large number as compared to the normal 
                              % value of your function
    dydt(isnan(dydt)) = 1e50; % same for nan.
end

Add these lines at the end of your odefun. It will replace all the infinities in dydt with a very large number. This technique will signal the ode23tb that it is approaching infinity if it projects too far from the current point while prevent the messing up of its internal states because of some unexpected NaN value caused by infinities.

Answer 1

[Tim]

Dear All,

Thanks for all your comments. They were really helpful.

I. It turned out that the problem is not with infinities or nan values. They never occured during the calculation.

However, I always got this message: Failure at t=4.164915e-10. Unable to meet integration tolerances without reducing the step size below the smallest value allowed (1.479675e-24) at time t

before the shortening of vector 't'. My time values were in the order of 1e-12 seconds. I thought this is small value so it would better if I reformulate my ODE system such that it uses picoseconds instead of 1e-12 seconds. After the implementation of this change I have not got any failure.

But I got very long running time. Even after 6 hours the ode solver (I tried all) did not even calculated the half of the points I needed.

As you said before, the ode solver decides the stepsize during the calculation, and there I found a problem: line 542-543 in ode15s :

% Predict a solution at t+h.

tnew = t + h;

The h value went below 1e-12 while the solver calculated a point between -3 and -4 ps. But I do not need better resolution than 1 ps unit.

It seems that some sets of parameters (I do not mean here the initial paramteres of the ODE system) gives the solver very hard time to make reasonable stepsize. This part I still do not understand.

II. The solution: using a solver with the fixed stepsize. I chose ode5 and it worked remarkably well. Not only for solving the ODE system but also using within the jacobianest function.

I performed the fitting and I got a consistent story in the end. Thanks again very much, I have really learnt a lot.