commons-issues mailing list archives

Site index · List index
Message view « Date » · « Thread »
Top « Date » · « Thread »
From "Phil Steitz (JIRA)" <j...@apache.org>
Subject [jira] Closed: (MATH-324) Error in the the convergence control for GraggBulirschStoerIntegrator (GBS)
Date Sat, 03 Apr 2010 20:39:29 GMT

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


Mime
View raw message