You are viewing a plain text version of this content. The canonical link for it is here.
Posted to issues@commons.apache.org by "Morand Vincent (JIRA)" <ji...@apache.org> on 2009/12/21 16:40:18 UTC
[jira] Created: (MATH-324) Error in the the convergence control for
GraggBulirschStoerIntegrator (GBS)
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.
[jira] Updated: (MATH-324) Error in the the convergence control for
GraggBulirschStoerIntegrator (GBS)
Posted by "Phil Steitz (JIRA)" <ji...@apache.org>.
[ https://issues.apache.org/jira/browse/MATH-324?page=com.atlassian.jira.plugin.system.issuetabpanels:all-tabpanel ]
Phil Steitz updated MATH-324:
-----------------------------
Fix Version/s: 2.1
> 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
> Fix For: 2.1
>
>
> 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.
[jira] Commented: (MATH-324) Error in the the convergence control
for GraggBulirschStoerIntegrator (GBS)
Posted by "Luc Maisonobe (JIRA)" <ji...@apache.org>.
[ https://issues.apache.org/jira/browse/MATH-324?page=com.atlassian.jira.plugin.system.issuetabpanels:comment-tabpanel&focusedCommentId=12793339#action_12793339 ]
Luc Maisonobe commented on MATH-324:
------------------------------------
Hi Vincent,
I'm trying to follow your reasoning. From what you say, it appears that what we call targetIter is what is called k in Hairer (with the classical index shift due to array beginning at 0 and not 1), and what we call k is only a current index when Hairer only specifically identity k-1, k and k+1. Am I right ? If yes, then I agree with you that the current implementation is wrong.
What puzzled me first was that applying the fix lead to a larger number of functions evaluations in the few test that were affected by the change. After more investigation, I think this is due to the fact the order is not reduced when the step is accepted so the number of evaluations per step is kept higher. In one case, I get a reduced number of steps but a larger number of calls, and at the end a slightly increased accuracy.
If you can confirm the understanding of our k and targetIter variables with respect to the theory, I'll commit 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.
[jira] Closed: (MATH-324) Error in the the convergence control for
GraggBulirschStoerIntegrator (GBS)
Posted by "Phil Steitz (JIRA)" <ji...@apache.org>.
[ https://issues.apache.org/jira/browse/MATH-324?page=com.atlassian.jira.plugin.system.issuetabpanels:all-tabpanel ]
Phil Steitz closed MATH-324.
----------------------------
> 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
> Fix For: 2.1
>
>
> 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.
[jira] Commented: (MATH-324) Error in the the convergence control
for GraggBulirschStoerIntegrator (GBS)
Posted by "Morand Vincent (JIRA)" <ji...@apache.org>.
[ 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.
[jira] Resolved: (MATH-324) Error in the the convergence control
for GraggBulirschStoerIntegrator (GBS)
Posted by "Luc Maisonobe (JIRA)" <ji...@apache.org>.
[ 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.