Commons Math
  1. Commons Math
  2. MATH-995

Adaptive division of segments in Quadrature Legendre-Gauss

    Details

    • Type: Bug Bug
    • Status: Closed
    • Priority: Major Major
    • Resolution: Fixed
    • Affects Version/s: None
    • Fix Version/s: 3.3
    • Labels:
      None

      Description

      I think the existing Legendre-Gauss object fails for certain integrals. An example of failure and a solution that divides segments based on error is provided. Please let me know if I'm not using the Legendre-Gauss object correctly.

      1. patch-test
        4 kB
        Ajo Fod
      2. patch-code
        6 kB
        Ajo Fod
      3. patch-code
        3 kB
        Ajo Fod
      4. patch-code
        3 kB
        Ajo Fod
      5. gaussian__sigma_1000.png
        12 kB
        Gilles
      6. gaussian__sigma_1000_zoom.png
        16 kB
        Gilles
      7. gaussian__sigma_1.png
        18 kB
        Gilles

        Issue Links

          Activity

          Hide
          Ajo Fod added a comment -

          The attached files show how to use AdaptiveQuadrature and how the existing method fails. I've wrapped the IterativeLegendreGaussIntegrator in an InfiniteIntegral object to show a specific instance of a failure of the class. The AdaptiveQuadrature object is more efficient at solving problems (in function evaluation counts) because it selectively increases resolution where the error is high.

          This problem is not limited to infinite integrals because the underlying IterativeLegendreGaussIntegrator is integrating in the region [-1,1].

          The attached solution uses 1st and 2nd order polynomials, but it can be generalized to a higher order polynomial solutions.

          Show
          Ajo Fod added a comment - The attached files show how to use AdaptiveQuadrature and how the existing method fails. I've wrapped the IterativeLegendreGaussIntegrator in an InfiniteIntegral object to show a specific instance of a failure of the class. The AdaptiveQuadrature object is more efficient at solving problems (in function evaluation counts) because it selectively increases resolution where the error is high. This problem is not limited to infinite integrals because the underlying IterativeLegendreGaussIntegrator is integrating in the region [-1,1] . The attached solution uses 1st and 2nd order polynomials, but it can be generalized to a higher order polynomial solutions.
          Hide
          Gilles added a comment -

          Your "InfiniteIntegral" class instantiates an "IterativeLegendreGaussIntegrator" with only 3 sample points.
          Could you test whether the problem persists when increasing the number of points?

          Show
          Gilles added a comment - Your "InfiniteIntegral" class instantiates an "IterativeLegendreGaussIntegrator" with only 3 sample points. Could you test whether the problem persists when increasing the number of points?
          Hide
          Ajo Fod added a comment - - edited

          These patches show the failure of the LGQ and the success of the Adaptive method. The adaptive method uses a fixed resolution quadrature (3) ... the concept itself is more general and work with higher orders. I've left system printline code to show failure/success ... feel free to remove/edit the code if you chose to accept it.

          Show
          Ajo Fod added a comment - - edited These patches show the failure of the LGQ and the success of the Adaptive method. The adaptive method uses a fixed resolution quadrature (3) ... the concept itself is more general and work with higher orders. I've left system printline code to show failure/success ... feel free to remove/edit the code if you chose to accept it.
          Hide
          Gilles added a comment -

          In revision 1497904, I've added a unit test.
          One step at a time: can we first figure out what exactly is going wrong?

          Show
          Gilles added a comment - In revision 1497904, I've added a unit test. One step at a time: can we first figure out what exactly is going wrong?
          Hide
          Ajo Fod added a comment -

          Look at the digits in accuracy per evaluation for LGQ vs AQ

          Show
          Ajo Fod added a comment - Look at the digits in accuracy per evaluation for LGQ vs AQ
          Hide
          Ajo Fod added a comment -

          An improvement in accuracy as suggested by Konstantin. This avoids the redundant computations in of the middle function values.

          Show
          Ajo Fod added a comment - An improvement in accuracy as suggested by Konstantin. This avoids the redundant computations in of the middle function values.
          Hide
          Gilles added a comment -

          I attached plots that hint to why the "IterativeGaussLegendreIntegrator" is not adapted to this particular change of variable.
          It appears to be a limitation of the algorithm but not a bug in the implementation.

          Show
          Gilles added a comment - I attached plots that hint to why the "IterativeGaussLegendreIntegrator" is not adapted to this particular change of variable. It appears to be a limitation of the algorithm but not a bug in the implementation.
          Hide
          Gilles added a comment -

          The ML contains a lengthy discussion about this issue, that is drifting away from the original report about a bug in CM.

          A warning has been added (revision 1499765) in the class's Javadoc to avoid the kind of usage that hits the limitation of this algorithm.

          Adaptive strategies are nice features to have in CM; but such a request belongs to another report, and how to best implement them in CM should be discussed in a specific thread on the "dev" ML.

          Show
          Gilles added a comment - The ML contains a lengthy discussion about this issue, that is drifting away from the original report about a bug in CM. A warning has been added (revision 1499765) in the class's Javadoc to avoid the kind of usage that hits the limitation of this algorithm. Adaptive strategies are nice features to have in CM; but such a request belongs to another report, and how to best implement them in CM should be discussed in a specific thread on the "dev" ML.
          Hide
          Ajo Fod added a comment -

          I like the plots ... does explain what is happening.

          However, if the IterativeGaussLegendreIntegrator is retained in the codebase, I'd suggest at the very least reporting on the non-convergence of the integration process or throwing an exception when the integration does'nt converge. Otherwise, any user of the code is likely to have to deal with expensive debugging session ... not good for confidence.

          Show
          Ajo Fod added a comment - I like the plots ... does explain what is happening. However, if the IterativeGaussLegendreIntegrator is retained in the codebase, I'd suggest at the very least reporting on the non-convergence of the integration process or throwing an exception when the integration does'nt converge. Otherwise, any user of the code is likely to have to deal with expensive debugging session ... not good for confidence.
          Hide
          Gilles added a comment -

          [...] if the IterativeGaussLegendreIntegrator is retained in the codebase [...]

          It is not the only instance where CM contains an algorithm that has shortcomings or drawbacks. Not every condition that leads to (numerical) problems can be easily characterized, short of seeing that the result is incorrect.
          [I recall a discussion about the "Regula falsi" root solver where the algorithm was stuck in an infinite loop; in that case it was possible to cheaply test for the condition and throw an exception.]

          The Javadoc now draws attention that the algorithm is not 100% fool-proof.

          Show
          Gilles added a comment - [...] if the IterativeGaussLegendreIntegrator is retained in the codebase [...] It is not the only instance where CM contains an algorithm that has shortcomings or drawbacks. Not every condition that leads to (numerical) problems can be easily characterized, short of seeing that the result is incorrect. [I recall a discussion about the "Regula falsi" root solver where the algorithm was stuck in an infinite loop; in that case it was possible to cheaply test for the condition and throw an exception.] The Javadoc now draws attention that the algorithm is not 100% fool-proof.
          Hide
          Ajo Fod added a comment -

          Well, I suppose if there is nothing one can do, you might as well mark this issue closed.

          Show
          Ajo Fod added a comment - Well, I suppose if there is nothing one can do, you might as well mark this issue closed.
          Hide
          Gilles added a comment -

          When some issue has been dealt with, it is first marked as "resolved"; it is "closed" post release.

          Show
          Gilles added a comment - When some issue has been dealt with, it is first marked as "resolved"; it is "closed" post release.
          Hide
          Gilles added a comment -

          Warning added (in r1499765) in the code documentation (cf. comment above).

          Show
          Gilles added a comment - Warning added (in r1499765) in the code documentation (cf. comment above).
          Hide
          Luc Maisonobe added a comment -

          Closing all resolved issue now available in released 3.3 version.

          Show
          Luc Maisonobe added a comment - Closing all resolved issue now available in released 3.3 version.

            People

            • Assignee:
              Unassigned
              Reporter:
              Ajo Fod
            • Votes:
              0 Vote for this issue
              Watchers:
              3 Start watching this issue

              Dates

              • Created:
                Updated:
                Resolved:

                Development