[ https://issues.apache.org/jira/browse/MATH324?page=com.atlassian.jira.plugin.system.issuetabpanels:commenttabpanel&focusedCommentId=12793622#action_12793622
]
Morand Vincent commented on MATH324:

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