I have two arrays $f(z)$ and $z$ both indexed by k and I want to solve $\frac{df}{dz}=\mu(1-f)$ to find $\mu(z)$
What would be the best numerical method to solve this equation?
$\begingroup$
$\endgroup$
-
1$\begingroup$ this is a ODE parameter identification problem: given the two data arrays $(z_k,f_k)$ (indexed by $k$), find the coefficient $\mu$ that best fits the data, where best is defined in a statistically precise way. My first approach to this would be to write it as an ODE-constrained optimization problem with a properly regularized cost function, then solve its discrete equivalent using a numerical method. The slides I linked to should get you started if you decide to go this route. $\endgroup$ – GoHokies Aug 29 '17 at 14:24
-
1$\begingroup$ alternatively, MATLAB can help with this as well $\endgroup$ – GoHokies Aug 29 '17 at 14:24
-
1$\begingroup$ or, since you have mentioned FORTRAN, there is the PDEFIT package that seems to do the trick: paper and software $\endgroup$ – GoHokies Aug 29 '17 at 14:36
-
$\begingroup$ I'm not sure about that recommendation in the slides. For nonlinear ODEs, solving parameter estimation problems with Newton's method or other local optimizers is usually a bad idea and will not find the local minima (though this case may be simple enough). Here I wrote up a general overview of how the algorithms tend to work though, and you can stick any optimizer and diffeq solver of your choice into that general framework. Of course, Julia has this all setup, but MATLAB has AMIGO2 as well $\endgroup$ – Chris Rackauckas Aug 29 '17 at 18:54
-
$\begingroup$ @ChrisRackauckas usually a bad idea - that's a bold statement. Newton will usually converge to a local minimum. The "quality" of that minimum will depend on your prior, on the data quality ($p(\rm{data}|\mu)$) and the exact choice of local optimizer, among other factors. $\endgroup$ – GoHokies Aug 29 '17 at 20:00
$\begingroup$
$\endgroup$
What you have is a linear ODE-constrained parameter identification problem: given the two data arrays $(z_k,f_k)$ (indexed by $k$) and the ODE model (the constraint), find the parameter $\mu$ that best fits the data, where the notion of best is defined in a statistically precise way. A very simple example would be a least-squares type regularized functional
$$ \arg \min_\mu \sum_k \left( f_k - f(z_k) \right)^2 + \| \mu \|^2 \; {\rm subject \; to\;} f_z = \mu (1-f).$$
One approach to this would be to write it as an ODE-constrained optimization problem with a properly regularized cost function, then solve its discrete equivalent using a numerical method and off-the shelf optimizer, as exemplified in these slides (1). Note that the optimization algorithm used there is a local Newton-type procedure. This means, for most practical problems, that you will need a reasonably good initial guess for the Newton procedure to converge to a "good-quality" minimum point, i.e., to a $\bar\mu$ not "too far" (in some norm) from the unknown "true" parameter $\mu^{\rm true}$.
If local minima are a major concern, and you do not have a good "guess" at $\mu^{\rm true}$ to start your Newton routine from, then global optimization is an alternative (caveat emptor: these procedures are more computationally expensive than a local routine). Or you can bring out the best of the two approaches into a combined global/local optimization strategy.
If you don't have the time to implement all this, then there are other options to consider:
- The MATLAB system identification toolbox can help (Google
idgrey,greyest). The workflow for your problem is described in detail here. - There are some FORTRAN libraries that do this. For example PDEFIT and/or EASY-FIT. I'd suggest you to read the manual to get acquainted with the features (and limitations) of each.
- Other programming languages appear to have at least some support for parameter identification: R, python with numpy/scipy.
