I have a RKF45 numerical integrator that simulates polymerization of proteins using CUDA.  It does so by tracking the populations of discrete length polymers, e.g. monomers, dimers, trimers, etc. all the way up to 8192-mers using the equations of the type that follow.$$\frac{dc_r}{dt}= k_nc_r^{n_c}+2k_a(c_{r-1}-c_r) + k_m(2\sum^{max}_{s=r+2}c_s-(r-1)c_r) $$where r represents the length of the polymer and c its concentration at that given time.  The k's are the parameters that I'd like to minimize an objective function with respect to.  These equations are integrated over about 20,000-30,000 timesteps.  My goal is to start fitting the model to experimental data.  The first type of fit I'm trying is fitting a value derived from these polymer populations at arbitrary timesteps.  For instance, I want to calculate the average length of polymers at five different times in the simulation/integration and compare them to externally supplied data.  The average length is derived for each timestep from the concentrations like this : $$L(t) = \frac{M(t)}{N(t)}$$$$M(t) = \sum^{max}_{s=n_c}s*c_s(t)$$$$N(t) = \sum^{max}_{s=n_c}c_s(t)$$The objective function I'm trying to minimize thusly looks like this : $$y(t)=\sqrt{\sum_{fit data}(L_{calculated}-L_{experimental})^2}$$I've already implemented a solution to do this and run it through the following norm and minimize said norm/objective function using Scipy's Nelder-Mead optimize method.  It seems to work, but if it's possible to make it perform better, I'd like to make it happen.  I want to try and utilize more robust and swifter converging methods like L-BFGS-B, especially because they also allow me to impose bounds on the parameters, that need to remain positive.Now, the Runge-Kutta algorithm, for the unfamiliar, is here : http://en.wikipedia.org/wiki/Runge%E2%80%93Kutta_methodsDue to the recursive nature of the Runge Kutta calculations such as $$k_1=f(y_n)$$$$k_2=f(y_n+a_{21}k_1)$$$$k_3=f(y_n+a_{31}k_1+a_{32}k_2)$$I'm having issues figuring out how to calculate the Jacobian, which I believe would be this.$$J(t) = [\frac{dy}{dk_a} \frac{dy}{dk_m} \frac{dy}{dk_n}]$$So, if I framed this question properly, is it possible at some point during the RKF calculation to construct a Jacobian to pass out to Scipy along with the objective function evaluation?  I'd like to get an optimized optimization script up and running ASAP.
Is there a relatively simple way to extract the Jacobian from a Runge-Kutta 4/5 integrator?
optimization;runge kutta;jacobian
So, if I framed this question properly, is it possible at some point during the RKF calculation to construct a Jacobian to pass out to Scipy along with the objective function evaluation?The Runge-Kutta-Fehlberg 4(5) method you describe is an explicit numerical method for solving an ODE. As such, it does not require the solution of a nonlinear algebraic equation; Jacobian matrices are not required for that algorithm.Is it possible to construct a Jacobian? Sure; the RKF4(5) method has nothing to do with it.What you want are parameter sensitivities, namely $\partial{c_{r}}/\partial{k_{a}}$, $\partial{c_{r}}/\partial{k_{m}}$, etc., so that you can calculate the derivatives you need.These can be calculated numerically in a few ways, to my knowledge; all of them involve solving a set of auxiliary differential equations:the adjoint method; see Steven G. Johnson's notes, and this paper by Petzold, et al., among others by Petzold. Adjoint methods work best when you have fewer state variables than parameters. (That is, the dimension of your $c_{r}$ is smaller than the number of parameters you have.)the simultaneous direct method, where you solve the sensitivity ODE system and your original state ODE system at the same time as a larger coupled ODE systemthe staggered direct method (see M. Caracotsios, W.E. Stewart, Sensitivity analysis of initial value problems with mixed ODEs and algebraic equations, Comput. Chem. Engrg. 9 (4) (1985) 359365)the simultaneous corrector method, where you solve for your state variables and sensitivities at the same time (see T. Maly, L.R. Petzold, Numerical methods and software for sensitivity analysis of differentialalgebraicsystems, Appl. Numer. Math. 20 (1997) 5779)the staggered corrector method (by Barton and Feehery; disclaimer: Barton was one of my advisors); see Feehery's thesis, Chapter 4, also see W.F. Feehery, J.E. Tolsma, P.I. Barton, Efficient sensitivity analysis of large-scale differentialalgebraic systems, Appl. Numer. Math. 25 (1997) 4154In chemistry and chemical engineering, there are a few other types of sensitivity analysis available, but I don't think of them will give you the Jacobian matrix information you need for parameter fitting.One final warning: ODE-constrained optimization problems are notoriously tricky, and generally nonconvex. It is possible to apply a standard nonlinear programming solver to these problems and obtain local minima that are globally suboptimal. See the thesis of Joseph K. Scott for examples. (Disclaimer: Dr. Scott was a former colleague.)