MATLAB: Ode23tb and non-linear curve fitting

Question

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.

Answer 1

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.