Thank you for your detailed answer.

JIRA issue MATH-1225 was created.

> 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:

[hidden email]
> For additional commands, e-mail:

[hidden email]
>

>