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.
