[
https://issues.apache.org/jira/browse/MATH-324?page=com.atlassian.jira.plugin.system.issuetabpanels:all-tabpanel
]
Luc Maisonobe resolved MATH-324.
--------------------------------
Resolution: Fixed
Fixed in subversion repository as of r893281.
Thanks for the bug report, the analysis and the fix.
> 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.
