Le 14/05/2015 16:40, Bernard GODARD a écrit : > Dear all, Hi Bernard,
> > The user guide on ordinary differential equations > http://commons.apache.org/proper/commons-math/userguide/ode.html > is very useful to understand how to use ODEWithJacobians and > ParameterizedODE. > but not up to date. You are right, we need to update this. Thanks for pointing the problem. Could you open a JIRA issue at <https://issues.apache.org/jira/browse/MATH> so we don't forget about it? > >>From the latest documentation, it does not seem clear to me how to use > ParameterizedODE and ParameterJacobianProvider in math 3.5. > > Moreover ParameterJacobianProvider-> computeParameterJacobian looks like > it can only handle one parameter at a time and is computed independently > from the main dynamic equation. For the use case I am interested in > (integration of the orbit of a satellite with complex dynamic model of > several hundred parameters) I prefer to compute together the force model > and its partial derivatives with respect to all parameters together > because these computations share a lot of intermediate results. Otherwise I > have to either recompute or use a cache mechanism with a test on input > parameter change (which also may be invalidated by the ODE integrator > depending on the order of the evaluation calls). This design then seems > inefficient to me but there might be a good reason for it. I would > appreciate if someone could explain this choice. > > I would rather have a method like: > > void computeDerivativesWithJacobian(double t, double[] y, > parameterList,double[] yDot,double[][] Jacobian ) > with parameterList an input specifying the ordered list of parameters for > which to compute the Jacobian. The rationale for splitting the computation in several parts (mainly the yDot from the Jacobian) was twofold. First, it allowed design of a simpler class when Jacobians are not desired, which could then be extended by subclasses that add Jacobians support. The second reason was that it allowed the Apache Commons Math library to provide by itself a generic implementation using finite differences for computing the Jacobians when the user would still provide only the base differential equation (this is a case of what Hairer, Nørsett and Wanner called internal differentiation). The rationale for splitting the Jacobians itself with one column at a time for one parameter is also twofold. First it frees the user (who provides this Jacobian) from worrying about putting the right parameter in the right column. The matrix is built automatically with proper columns assignment even if for some runs all parameters were selected and in other runs only a subset of the known parameters are used for differentiation, the other parameters being considered fixed. The second reason is that the providers for the various columns (i.e. for the various parameters) can be different objects. You can even have one different class for each parameter if you want. This allows for example to handle ODE for which some parameters are computed analytically and some other parameters are computed by finite differences. In such a case, user would call JacobianMatrices.addParameterJacobianProvider for the parameters for which he knows how to differentiate, and he would rather call JacobianMatrices.setParameterStep for the parameters for which he doesn't know how to differentiate. The use case for these two aspect was also orbit determination: user can decide at run time which parameters are estimated and which are fixed, hence increasing or decreasing the number of columns in the Jacobian, and some force model may be complex enough that no analytical expression is available, hence needing finite differences. The current design can be used in your case and benefit from already computed values. This is not clear in the documentation, though and should be better explained. In fact, the various calls to the computeParameterJacobian that are used to retrieve all the columns one at a time are all done in a small loop and without any change to the state, i.e. they correspond to the same step. The scheduling is as follows: for each evaluation (i.e. for one value of t and y[]): - one call to the computeMainStateJacobian method (which belongs to either the MainStateJacobianProvider interface or the MainStateJacobianProvider interface, depending on which JacobianMatrices constructor you use) - for each selected parameter, iterating in the exact same order the parameters were specified in JacobianMatrices constructor: * one call to computeParameterJacobian for the current parameter, using the proper provider (which may be a finite differences provider for some parameters and an analytical provider for other parameters) - end for end for With the calls sequence above, you can for example perform all the computation early (say when computeMainStateJacobian is called) and store the Jacobians columns in an internal array, and then only copy the column back when the computeParameterJacobian method is called in the internal loop. Of course, this implies the implementations of MainStateJacobianProvider and ParameterJacobianProvider know each other. This is even easier if they correspond to the same object, but it is not required. So in your case, if you already know how to compute the full Jacobian by yourself, I would suggest to put everything in a top level class that would implement: - either MainStateJacobianProvider or MainStateJacobianProvider - AND ParameterJacobianProvider This could be done as follows (note that the computeParameterJacobian is complete here, it can really be restricted to two lines of code): MyClass implements MainStateJacobianProvider, ParameterJacobianProvider { private double [][] jacobian; private Map<String, int> nameToIndex; public MyClass() { // set up the nameToIndex map } public void computeDerivatives(double t, double[] y, double[] yDot) { // compute yDot // compute and store Jacobian // (beware to store it in transposed form ...) } void computeParameterJacobian(double t, double[] y, double[] yDot, String paramName, double[] dFdP) { // this method is known to be called after computeDerivatives // we simply pick up the column that was already computed before // and that was stored transposed final double[] precomputed = jacobian[nameToIndex.get(name)]; System.arraycopy(precomputed, 0, dFdP, 0, precomputed.length); } } Hope this helps, Luc > > > Thank you, > > Kind regards, > > > Bernard Godard > --------------------------------------------------------------------- To unsubscribe, e-mail: user-unsubscr...@commons.apache.org For additional commands, e-mail: user-h...@commons.apache.org
