[ https://issues.apache.org/jira/browse/MATH324?page=com.atlassian.jira.plugin.system.issuetabpanels:alltabpanel
]
Luc Maisonobe resolved MATH324.

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: MATH324
> URL: https://issues.apache.org/jira/browse/MATH324
> 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 k1 modified midpoint integration to obtain the k1lines
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 (errk1 < 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
k1) is
> errk1 > (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.
