[jira] Closed: (MATH-324) Error in the the convergence control for GraggBulirschStoerIntegrator (GBS)

["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.