[
https://issues.apache.org/jira/browse/MATH-324?page=com.atlassian.jira.plugin.system.issuetabpanels:comment-tabpanel&focusedCommentId=12793622#action_12793622
]
Morand Vincent commented on MATH-324:
-------------------------------------
Hi Luc,
you have perfectly understood the error, targetIter is the index we have
predicted in the previous step and k is the current index, whereas Hairer uses
only "k" to explain his theory.
Before sending my bug report i ran some test (those provided by apache) and
also found a larger number of function evaluations. I think it is explained by
the fact that the asymptotic development of the error is just an estimation of
the chance of convergence in k+1. Sometimes it gives the good information,
sometimes not, depending on the function itself.
Whatever, applying the fix allows us to follow the theory described by Hairer.
Thanks for the efficiency of your answer
Vincent
> Error in the the convergence control for GraggBulirschStoerIntegrator (GBS)
> ---------------------------------------------------------------------------
>
> Key: MATH-324
> URL: https://issues.apache.org/jira/browse/MATH-324
> Project: Commons Math
> Issue Type: Improvement
> Affects Versions: 2.0
> Environment: Red Hat 5
> Reporter: Morand Vincent
>
> By reading the source code and making a compareason with the theory described
> by Hairer in his book (Solving Ordinary Differential Equations 1 : Nonstiff
> Problems) I have found an error in the convergence control.
> EFFECTS :
> The mistake is unvisible from a user's point of view but makes the
> integration slower than it should be. Some steps are rejected whereas they
> could have offer convergence. The theory discribed by Hairer isn't correctly
> used.
> LOCATION :
> The problem comes from the line 693 of GraggBulirschStoerIntegrator.java :
> "final double ratio = ((double) sequence [k] * sequence[k+1]) / (sequence[0]
> * sequence[0]);"
> This line should be replaced by :
> final double ratio = ((double) sequence [k+1] * sequence[k+2]) / (sequence[0]
> * sequence[0]);
> or (which is the same) : final double ratio = ((double) sequence [targetIter]
> * sequence[targetIter + 1]) / (sequence[0] * sequence[0]);
> EXPLANATION :
> The corresponding chapter in Hairer's book is page 233 : Order and Step Size
> Control (for GBS method) :
> Following the theory, we compute k-1 modified midpoint integration to obtain
> the k-1lines of the extrapolation tableau where k is the predicted index in
> which there should be convergence. (In the java code source the variable is
> 'targetIter')
> If there is convergence (errk-1 < 1) we accept the step and continue the
> integration.
> If not we use the asymptotic evolution of error to evaluate if there is
> convergence at least in line k+1 (so two integration later)
> This stage is wrongly done in the java code source by the line 693.
> Following the theory the estimation of convergence in line k+1 (when you are
> in line k-1) is
> errk-1 > (nk+1 * nk / (n1 * n1) )²
> So you shall use the number of substeps (nk and nk+1) that haven't been used
> yet. you shall use the number of substep of the next two integration.
> In the java code source :
> * the predicted index for convergence is targetIter
> * the current integration number is k
> * nk, number of substep for the integration is sequence[ ]
> When we are in line 693 :
> Here k = targetIter - 1 (case -1 of the switch)
> To evaluate the convergence at least in targetIter + 1 is is done :
> "final double ratio = ((double) sequence [k] * sequence[k+1]) /
> (sequence[0] * sequence[0]);"
> It seems like a nonsense to use sequence[k] to evaluate convergence in next
> lines whereas sequence[k] has already been used in the last call to
> tryStep(...,k,..)
> we should do :
> final double ratio = ((double) sequence [targetIter] *
> sequence[targetIter + 1]) / ( sequence[0] * sequence[0]);
> Or (because targetIter = k+1) final double ratio = ((double) sequence [k+1]
> *sequence[k+2]) /(sequence[0] * sequence[0]);
> ie to use the number os substeps for the next integrations, as described in
> the Hairer's book.
> COMMENT :
> 1) in the java code for case = 0 we have k = targetIter and the estimation
> of the convergence at least in targetIter + 1 is done by
> final double ratio = ((double) sequence[k+1]) / sequence[0];
> Since k = targetIter the instruction is correct.
> 2) The compareason with the fortran code delivered by Hairer (easily found on
> his website page) confirms the error in the java code.
> Sorry for the length of this message, looking forward to hearing from you soon
> Vincent Morand
--
This message is automatically generated by JIRA.
-
You can reply to this email to add a comment to the issue online.
