Details
-
Improvement
-
Status: Closed
-
Major
-
Resolution: Fixed
-
2.0
-
None
-
None
-
Red Hat 5
Description
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