# ode23tb and non-linear curve fitting

2 views (last 30 days)
Tim on 29 Apr 2020
Answered: Tim on 6 May 2020
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 on 30 Apr 2020
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
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.

Tim on 6 May 2020
Dear All,
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.

R2016b

### Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!