[jira] [Created] (MATH-797) Single step integrators

classic Classic list List threaded Threaded
26 messages Options
12
Reply | Threaded
Open this post in threaded view
|

[jira] [Created] (MATH-797) Single step integrators

ASF GitHub Bot (Jira)
Gilles created MATH-797:
---------------------------

             Summary: Single step integrators
                 Key: MATH-797
                 URL: https://issues.apache.org/jira/browse/MATH-797
             Project: Commons Math
          Issue Type: Wish
    Affects Versions: 3.0
            Reporter: Gilles
            Assignee: Gilles
            Priority: Trivial
             Fix For: 3.1


CM assumes that the user wants to integrate a complex function on a large interval, so the large interval has to be subdivided into many subintervals. CM does the partition, and performs convergence checks, using an iterative approach.

However, if the function is smooth enough, no subdivision of the integration interval is required. Those use-cases could benefit from the efficiency gain of not performing a convergence check.

The proposal is to provide a new interface "UnivariateSingleStepIntegrator":
{code}
interface SingleIntervalIntegrator {
    /**
     * Method for implementing a single interval integration.
     * There is no convergence checks because it is not iterative.
     *
     * @param f Function to integrate.
     * @param lower Lower bound of the interval over which to integrate.
     * @param upper Upper bound of the interval over which to integrate.
     * @return the integrated value.
     */
    double integrate(UnivariateFunction f,
                     double lower,
                     double upper);
}
{code}

In effect, the implementation of the above "integrate" method of a new "LegendreGaussIntegratorSingleStepIntegrator" would the equivalent of "stage(1)" in the current "LegendreGaussIntegrator".


--
This message is automatically generated by JIRA.
If you think it was sent incorrectly, please contact your JIRA administrators: https://issues.apache.org/jira/secure/ContactAdministrators!default.jspa
For more information on JIRA, see: http://www.atlassian.com/software/jira

       
Reply | Threaded
Open this post in threaded view
|

[jira] [Commented] (MATH-797) Single step integrators

ASF GitHub Bot (Jira)

    [ https://issues.apache.org/jira/browse/MATH-797?page=com.atlassian.jira.plugin.system.issuetabpanels:comment-tabpanel&focusedCommentId=13287998#comment-13287998 ]

Sébastien Brisard commented on MATH-797:
----------------------------------------

I don't really follow. Are you proposing to duplicate all existing integrators (one single step, one adaptative implementation)? The immediate answer would be to decide that all integrators are single step, and define one iterative integrator, which would be provided with a single step integrator at the time of construction. But that does not work (see below).

I think that we can't have one single policy for all integrators
* trapezoid, Simpson and Romberg integrators are iterative in nature, and IIRC, there are some clever optimisations which prevent the separation single-step/iterative: the next iteration has to have access to internal calculations done in the previous iteration.
* Gauss-like integrators are in my view single step integrators. Indeed, I think that people use Gauss integration schemes in order to have absolute control on the number of operations which are carried out. This goes opposite to the current implementation of these integrators in CM3.

So, to sum up, I think that if we were to create a second interface, maybe we could think about creating a GaussIntegrator, dedicated to Gaussian integration schemes. These schemes would not be iterative, like what is proposed in this ticket. All other existing integrators would remain unchanged.

By the way, the maximum number of integration points in the Gauss-Legendre scheme is limiting. I have some reliable (hi accuracy) code to compute the points and weights at any order. Also, I have an implementation for other Gauss-like schemes, which I'm happy to contribute if we can agree on a common interface with Gaussian integration schemes (which I find very useful and efficient).
               

> Single step integrators
> -----------------------
>
>                 Key: MATH-797
>                 URL: https://issues.apache.org/jira/browse/MATH-797
>             Project: Commons Math
>          Issue Type: Wish
>    Affects Versions: 3.0
>            Reporter: Gilles
>            Assignee: Gilles
>            Priority: Trivial
>             Fix For: 3.1
>
>
> CM assumes that the user wants to integrate a complex function on a large interval, so the large interval has to be subdivided into many subintervals. CM does the partition, and performs convergence checks, using an iterative approach.
> However, if the function is smooth enough, no subdivision of the integration interval is required. Those use-cases could benefit from the efficiency gain of not performing a convergence check.
> The proposal is to provide a new interface "UnivariateSingleStepIntegrator":
> {code}
> interface SingleIntervalIntegrator {
>     /**
>      * Method for implementing a single interval integration.
>      * There is no convergence checks because it is not iterative.
>      *
>      * @param f Function to integrate.
>      * @param lower Lower bound of the interval over which to integrate.
>      * @param upper Upper bound of the interval over which to integrate.
>      * @return the integrated value.
>      */
>     double integrate(UnivariateFunction f,
>                      double lower,
>                      double upper);
> }
> {code}
> In effect, the implementation of the above "integrate" method of a new "LegendreGaussIntegratorSingleStepIntegrator" would the equivalent of "stage(1)" in the current "LegendreGaussIntegrator".

--
This message is automatically generated by JIRA.
If you think it was sent incorrectly, please contact your JIRA administrators: https://issues.apache.org/jira/secure/ContactAdministrators!default.jspa
For more information on JIRA, see: http://www.atlassian.com/software/jira


Reply | Threaded
Open this post in threaded view
|

[jira] [Commented] (MATH-797) Single step integrators

ASF GitHub Bot (Jira)
In reply to this post by ASF GitHub Bot (Jira)

    [ https://issues.apache.org/jira/browse/MATH-797?page=com.atlassian.jira.plugin.system.issuetabpanels:comment-tabpanel&focusedCommentId=13288028#comment-13288028 ]

Gilles commented on MATH-797:
-----------------------------

Actually, I filed this issue on behalf of another team within the project I work for.
They indeed are currently using a stripped-down (single step) version of the "LegendreGaussIntegrator".
Their problem is mainly an efficiency one because for their use-case, they know in advance that a single iteration is enough.

So, your proposal is quite fine (i.e. focus only on  Gauss-like integrators and make them more flexible). Thanks!

               

> Single step integrators
> -----------------------
>
>                 Key: MATH-797
>                 URL: https://issues.apache.org/jira/browse/MATH-797
>             Project: Commons Math
>          Issue Type: Wish
>    Affects Versions: 3.0
>            Reporter: Gilles
>            Assignee: Gilles
>            Priority: Trivial
>             Fix For: 3.1
>
>
> CM assumes that the user wants to integrate a complex function on a large interval, so the large interval has to be subdivided into many subintervals. CM does the partition, and performs convergence checks, using an iterative approach.
> However, if the function is smooth enough, no subdivision of the integration interval is required. Those use-cases could benefit from the efficiency gain of not performing a convergence check.
> The proposal is to provide a new interface "UnivariateSingleStepIntegrator":
> {code}
> interface SingleIntervalIntegrator {
>     /**
>      * Method for implementing a single interval integration.
>      * There is no convergence checks because it is not iterative.
>      *
>      * @param f Function to integrate.
>      * @param lower Lower bound of the interval over which to integrate.
>      * @param upper Upper bound of the interval over which to integrate.
>      * @return the integrated value.
>      */
>     double integrate(UnivariateFunction f,
>                      double lower,
>                      double upper);
> }
> {code}
> In effect, the implementation of the above "integrate" method of a new "LegendreGaussIntegratorSingleStepIntegrator" would the equivalent of "stage(1)" in the current "LegendreGaussIntegrator".

--
This message is automatically generated by JIRA.
If you think it was sent incorrectly, please contact your JIRA administrators: https://issues.apache.org/jira/secure/ContactAdministrators!default.jspa
For more information on JIRA, see: http://www.atlassian.com/software/jira

       
Reply | Threaded
Open this post in threaded view
|

[jira] [Commented] (MATH-797) Single step integrators

ASF GitHub Bot (Jira)
In reply to this post by ASF GitHub Bot (Jira)

    [ https://issues.apache.org/jira/browse/MATH-797?page=com.atlassian.jira.plugin.system.issuetabpanels:comment-tabpanel&focusedCommentId=13288125#comment-13288125 ]

Luc Maisonobe commented on MATH-797:
------------------------------------

Big +1 from me.
We also hit some problems with iterations with the current integrators at work, and had to recode some stripped down versions with fixed steps. Having deterministic integration without convergence checks for some simple cases using Gauss-Legendre would be fine.

We also dicussed some time ago about extending the order of Gauss-Legendre in another issue (I don't remember which one) but did not implement it at that time. It would be nice to do it now.

               

> Single step integrators
> -----------------------
>
>                 Key: MATH-797
>                 URL: https://issues.apache.org/jira/browse/MATH-797
>             Project: Commons Math
>          Issue Type: Wish
>    Affects Versions: 3.0
>            Reporter: Gilles
>            Assignee: Gilles
>            Priority: Trivial
>             Fix For: 3.1
>
>
> CM assumes that the user wants to integrate a complex function on a large interval, so the large interval has to be subdivided into many subintervals. CM does the partition, and performs convergence checks, using an iterative approach.
> However, if the function is smooth enough, no subdivision of the integration interval is required. Those use-cases could benefit from the efficiency gain of not performing a convergence check.
> The proposal is to provide a new interface "UnivariateSingleStepIntegrator":
> {code}
> interface SingleIntervalIntegrator {
>     /**
>      * Method for implementing a single interval integration.
>      * There is no convergence checks because it is not iterative.
>      *
>      * @param f Function to integrate.
>      * @param lower Lower bound of the interval over which to integrate.
>      * @param upper Upper bound of the interval over which to integrate.
>      * @return the integrated value.
>      */
>     double integrate(UnivariateFunction f,
>                      double lower,
>                      double upper);
> }
> {code}
> In effect, the implementation of the above "integrate" method of a new "LegendreGaussIntegratorSingleStepIntegrator" would the equivalent of "stage(1)" in the current "LegendreGaussIntegrator".

--
This message is automatically generated by JIRA.
If you think it was sent incorrectly, please contact your JIRA administrators: https://issues.apache.org/jira/secure/ContactAdministrators!default.jspa
For more information on JIRA, see: http://www.atlassian.com/software/jira

       
Reply | Threaded
Open this post in threaded view
|

[jira] [Comment Edited] (MATH-797) Single step integrators

ASF GitHub Bot (Jira)
In reply to this post by ASF GitHub Bot (Jira)

    [ https://issues.apache.org/jira/browse/MATH-797?page=com.atlassian.jira.plugin.system.issuetabpanels:comment-tabpanel&focusedCommentId=13288163#comment-13288163 ]

Sébastien Brisard edited comment on MATH-797 at 6/3/12 1:58 PM:
----------------------------------------------------------------

I have a jQuadrature project sitting in my workspace for almost two years now, which I'm very happy to refactor and contribute if we settle on a proper interface for Gauss integrators. I have currently implemented: Gauss-Legendre, Gauss-Chebyshev (first and second kinds), Gauss-Hermite. Integration points and weights are computed and cached if needed, so the order is not limited.

In NR for example, the integration points were computed with built-in Newton iterations. I chose bracketing instead. It may be a bit slower, but it's not called very often (only when new Gaussian rules are needed), and accuracy of the computed roots is then garanteed (to within 1ulp). It's more difficult to insure the accuracy of the weights (when you compute them on-the-fly), and I wrote a high accuracy computation of integration points and weights for the Gauss-Legendre rule, using {{BigDecimal}} (at that time, I knew but did not use Commons-Math, shame on me!).

If any of this seems interesting to you, it's right here !

One point worth noting: with Gauss integration schemes, I always wondered whether calling the integration method {{integrate(f, a, b)}} was a good idea. Indeed, generic integration schemes do not integrate the provided univariate function {{f( x )}}, but the weighted function {{f( x ) * w( x )}}, which might be confusing. Of course, there is no such ambiguity with Gauss-Legendre...
               
      was (Author: celestin):
    I have a jQuadrature project sitting in my workspace for almost two years now, which I'm very happy to refactor and contribute if we settle on a proper interface for Gauss integrators. I have currently implemented: Gauss-Legendre, Gauss-Chebyshev (first and second kinds), Gauss-Hermite. Integration points and weights are computed and cached if needed, so the order is not limited.

In NR for example, the integration points were computed with built-in Newton iterations. I chose bracketing instead. It may be a bit slower, but it's not called very often (only when new Gaussian rules are needed), and accuracy of the computed roots is then garanteed (to within 1ulp). It's more difficult to insure the accuracy of the weights (when you compute them on-the-fly), and I wrote a high accuracy computation of integration points and weights for the Gauss-Legendre rule, using {{BigDecimal}} (at that time, I knew but did not use Commons-Math, shame on me!).

If any of this seems interesting to you, it's right here !

One point worth noting: with Gauss integration schemes, I always wondered whether calling the integration method {{integrate(f, a, b)}} was a good idea. Indeed, generic integration schemes do not integrate the provided univariate function {{f(x)}}, but the weighted function {{f(x) * w(x)}}, which might be confusing. Of course, there is no such ambiguity with Gauss-Legendre...
                 

> Single step integrators
> -----------------------
>
>                 Key: MATH-797
>                 URL: https://issues.apache.org/jira/browse/MATH-797
>             Project: Commons Math
>          Issue Type: Wish
>    Affects Versions: 3.0
>            Reporter: Gilles
>            Assignee: Gilles
>            Priority: Trivial
>             Fix For: 3.1
>
>
> CM assumes that the user wants to integrate a complex function on a large interval, so the large interval has to be subdivided into many subintervals. CM does the partition, and performs convergence checks, using an iterative approach.
> However, if the function is smooth enough, no subdivision of the integration interval is required. Those use-cases could benefit from the efficiency gain of not performing a convergence check.
> The proposal is to provide a new interface "UnivariateSingleStepIntegrator":
> {code}
> interface SingleIntervalIntegrator {
>     /**
>      * Method for implementing a single interval integration.
>      * There is no convergence checks because it is not iterative.
>      *
>      * @param f Function to integrate.
>      * @param lower Lower bound of the interval over which to integrate.
>      * @param upper Upper bound of the interval over which to integrate.
>      * @return the integrated value.
>      */
>     double integrate(UnivariateFunction f,
>                      double lower,
>                      double upper);
> }
> {code}
> In effect, the implementation of the above "integrate" method of a new "LegendreGaussIntegratorSingleStepIntegrator" would the equivalent of "stage(1)" in the current "LegendreGaussIntegrator".

--
This message is automatically generated by JIRA.
If you think it was sent incorrectly, please contact your JIRA administrators: https://issues.apache.org/jira/secure/ContactAdministrators!default.jspa
For more information on JIRA, see: http://www.atlassian.com/software/jira


Reply | Threaded
Open this post in threaded view
|

[jira] [Commented] (MATH-797) Single step integrators

ASF GitHub Bot (Jira)
In reply to this post by ASF GitHub Bot (Jira)

    [ https://issues.apache.org/jira/browse/MATH-797?page=com.atlassian.jira.plugin.system.issuetabpanels:comment-tabpanel&focusedCommentId=13288163#comment-13288163 ]

Sébastien Brisard commented on MATH-797:
----------------------------------------

I have a jQuadrature project sitting in my workspace for almost two years now, which I'm very happy to refactor and contribute if we settle on a proper interface for Gauss integrators. I have currently implemented: Gauss-Legendre, Gauss-Chebyshev (first and second kinds), Gauss-Hermite. Integration points and weights are computed and cached if needed, so the order is not limited.

In NR for example, the integration points were computed with built-in Newton iterations. I chose bracketing instead. It may be a bit slower, but it's not called very often (only when new Gaussian rules are needed), and accuracy of the computed roots is then garanteed (to within 1ulp). It's more difficult to insure the accuracy of the weights (when you compute them on-the-fly), and I wrote a high accuracy computation of integration points and weights for the Gauss-Legendre rule, using {{BigDecimal}} (at that time, I knew but did not use Commons-Math, shame on me!).

If any of this seems interesting to you, it's right here !

One point worth noting: with Gauss integration schemes, I always wondered whether calling the integration method {{integrate(f, a, b)}} was a good idea. Indeed, generic integration schemes do not integrate the provided univariate function {{f(x)}}, but the weighted function {{f(x) * w(x)}}, which might be confusing. Of course, there is no such ambiguity with Gauss-Legendre...
               

> Single step integrators
> -----------------------
>
>                 Key: MATH-797
>                 URL: https://issues.apache.org/jira/browse/MATH-797
>             Project: Commons Math
>          Issue Type: Wish
>    Affects Versions: 3.0
>            Reporter: Gilles
>            Assignee: Gilles
>            Priority: Trivial
>             Fix For: 3.1
>
>
> CM assumes that the user wants to integrate a complex function on a large interval, so the large interval has to be subdivided into many subintervals. CM does the partition, and performs convergence checks, using an iterative approach.
> However, if the function is smooth enough, no subdivision of the integration interval is required. Those use-cases could benefit from the efficiency gain of not performing a convergence check.
> The proposal is to provide a new interface "UnivariateSingleStepIntegrator":
> {code}
> interface SingleIntervalIntegrator {
>     /**
>      * Method for implementing a single interval integration.
>      * There is no convergence checks because it is not iterative.
>      *
>      * @param f Function to integrate.
>      * @param lower Lower bound of the interval over which to integrate.
>      * @param upper Upper bound of the interval over which to integrate.
>      * @return the integrated value.
>      */
>     double integrate(UnivariateFunction f,
>                      double lower,
>                      double upper);
> }
> {code}
> In effect, the implementation of the above "integrate" method of a new "LegendreGaussIntegratorSingleStepIntegrator" would the equivalent of "stage(1)" in the current "LegendreGaussIntegrator".

--
This message is automatically generated by JIRA.
If you think it was sent incorrectly, please contact your JIRA administrators: https://issues.apache.org/jira/secure/ContactAdministrators!default.jspa
For more information on JIRA, see: http://www.atlassian.com/software/jira


Reply | Threaded
Open this post in threaded view
|

[jira] [Comment Edited] (MATH-797) Single step integrators

ASF GitHub Bot (Jira)
In reply to this post by ASF GitHub Bot (Jira)

    [ https://issues.apache.org/jira/browse/MATH-797?page=com.atlassian.jira.plugin.system.issuetabpanels:comment-tabpanel&focusedCommentId=13288163#comment-13288163 ]

Sébastien Brisard edited comment on MATH-797 at 6/3/12 2:04 PM:
----------------------------------------------------------------

I have a jQuadrature project sitting in my workspace for almost two years now, which I'm very happy to refactor and contribute if we settle on a proper interface for Gauss integrators. I have currently implemented: Gauss-Legendre, Gauss-Chebyshev (first and second kinds), Gauss-Hermite. Integration points and weights are computed and cached if needed, so the order is not limited.

In NR for example, the integration points were computed with built-in Newton iterations. I chose bracketing instead. It may be a bit slower, but it's not called very often (only when new Gaussian rules are needed), and accuracy of the computed roots is then garanteed (to within 1ulp). It's more difficult to insure the accuracy of the weights (when you compute them on-the-fly), and I wrote a high accuracy computation of integration points and weights for the Gauss-Legendre rule, using {{BigDecimal}} (at that time, I knew but did not use Commons-Math, shame on me!).

If any of this seems interesting to you, it's right here !

One point worth noting: with Gauss integration schemes, I always wondered whether naming the integration method {{integrate(f, a, b)}} was a good idea. Indeed, generic integration schemes do not integrate the provided univariate function {{f( x )}}, but the weighted function {{f( x ) * w( x )}}, which might be confusing. Of course, there is no such ambiguity with Gauss-Legendre...
               
      was (Author: celestin):
    I have a jQuadrature project sitting in my workspace for almost two years now, which I'm very happy to refactor and contribute if we settle on a proper interface for Gauss integrators. I have currently implemented: Gauss-Legendre, Gauss-Chebyshev (first and second kinds), Gauss-Hermite. Integration points and weights are computed and cached if needed, so the order is not limited.

In NR for example, the integration points were computed with built-in Newton iterations. I chose bracketing instead. It may be a bit slower, but it's not called very often (only when new Gaussian rules are needed), and accuracy of the computed roots is then garanteed (to within 1ulp). It's more difficult to insure the accuracy of the weights (when you compute them on-the-fly), and I wrote a high accuracy computation of integration points and weights for the Gauss-Legendre rule, using {{BigDecimal}} (at that time, I knew but did not use Commons-Math, shame on me!).

If any of this seems interesting to you, it's right here !

One point worth noting: with Gauss integration schemes, I always wondered whether calling the integration method {{integrate(f, a, b)}} was a good idea. Indeed, generic integration schemes do not integrate the provided univariate function {{f( x )}}, but the weighted function {{f( x ) * w( x )}}, which might be confusing. Of course, there is no such ambiguity with Gauss-Legendre...
                 

> Single step integrators
> -----------------------
>
>                 Key: MATH-797
>                 URL: https://issues.apache.org/jira/browse/MATH-797
>             Project: Commons Math
>          Issue Type: Wish
>    Affects Versions: 3.0
>            Reporter: Gilles
>            Assignee: Gilles
>            Priority: Trivial
>             Fix For: 3.1
>
>
> CM assumes that the user wants to integrate a complex function on a large interval, so the large interval has to be subdivided into many subintervals. CM does the partition, and performs convergence checks, using an iterative approach.
> However, if the function is smooth enough, no subdivision of the integration interval is required. Those use-cases could benefit from the efficiency gain of not performing a convergence check.
> The proposal is to provide a new interface "UnivariateSingleStepIntegrator":
> {code}
> interface SingleIntervalIntegrator {
>     /**
>      * Method for implementing a single interval integration.
>      * There is no convergence checks because it is not iterative.
>      *
>      * @param f Function to integrate.
>      * @param lower Lower bound of the interval over which to integrate.
>      * @param upper Upper bound of the interval over which to integrate.
>      * @return the integrated value.
>      */
>     double integrate(UnivariateFunction f,
>                      double lower,
>                      double upper);
> }
> {code}
> In effect, the implementation of the above "integrate" method of a new "LegendreGaussIntegratorSingleStepIntegrator" would the equivalent of "stage(1)" in the current "LegendreGaussIntegrator".

--
This message is automatically generated by JIRA.
If you think it was sent incorrectly, please contact your JIRA administrators: https://issues.apache.org/jira/secure/ContactAdministrators!default.jspa
For more information on JIRA, see: http://www.atlassian.com/software/jira


Reply | Threaded
Open this post in threaded view
|

[jira] [Commented] (MATH-797) Single step integrators

ASF GitHub Bot (Jira)
In reply to this post by ASF GitHub Bot (Jira)

    [ https://issues.apache.org/jira/browse/MATH-797?page=com.atlassian.jira.plugin.system.issuetabpanels:comment-tabpanel&focusedCommentId=13288195#comment-13288195 ]

Luc Maisonobe commented on MATH-797:
------------------------------------

Sounds good.

If you intend to contribute a bunch of code that is already existing and not developed within the Apache Commons project, then we would need some paperwork before including it. All IP holders of this code should sign a Software Grant that allow the foundation to distribute this under the Apache V2 license. See [http://www.apache.org/licenses/#grants].
               

> Single step integrators
> -----------------------
>
>                 Key: MATH-797
>                 URL: https://issues.apache.org/jira/browse/MATH-797
>             Project: Commons Math
>          Issue Type: Wish
>    Affects Versions: 3.0
>            Reporter: Gilles
>            Assignee: Gilles
>            Priority: Trivial
>             Fix For: 3.1
>
>
> CM assumes that the user wants to integrate a complex function on a large interval, so the large interval has to be subdivided into many subintervals. CM does the partition, and performs convergence checks, using an iterative approach.
> However, if the function is smooth enough, no subdivision of the integration interval is required. Those use-cases could benefit from the efficiency gain of not performing a convergence check.
> The proposal is to provide a new interface "UnivariateSingleStepIntegrator":
> {code}
> interface SingleIntervalIntegrator {
>     /**
>      * Method for implementing a single interval integration.
>      * There is no convergence checks because it is not iterative.
>      *
>      * @param f Function to integrate.
>      * @param lower Lower bound of the interval over which to integrate.
>      * @param upper Upper bound of the interval over which to integrate.
>      * @return the integrated value.
>      */
>     double integrate(UnivariateFunction f,
>                      double lower,
>                      double upper);
> }
> {code}
> In effect, the implementation of the above "integrate" method of a new "LegendreGaussIntegratorSingleStepIntegrator" would the equivalent of "stage(1)" in the current "LegendreGaussIntegrator".

--
This message is automatically generated by JIRA.
If you think it was sent incorrectly, please contact your JIRA administrators: https://issues.apache.org/jira/secure/ContactAdministrators!default.jspa
For more information on JIRA, see: http://www.atlassian.com/software/jira

       
Reply | Threaded
Open this post in threaded view
|

[jira] [Commented] (MATH-797) Single step integrators

ASF GitHub Bot (Jira)
In reply to this post by ASF GitHub Bot (Jira)

    [ https://issues.apache.org/jira/browse/MATH-797?page=com.atlassian.jira.plugin.system.issuetabpanels:comment-tabpanel&focusedCommentId=13288210#comment-13288210 ]

Sébastien Brisard commented on MATH-797:
----------------------------------------

No, no. That's personal code I've never managed to really finish...
               

> Single step integrators
> -----------------------
>
>                 Key: MATH-797
>                 URL: https://issues.apache.org/jira/browse/MATH-797
>             Project: Commons Math
>          Issue Type: Wish
>    Affects Versions: 3.0
>            Reporter: Gilles
>            Assignee: Gilles
>            Priority: Trivial
>             Fix For: 3.1
>
>
> CM assumes that the user wants to integrate a complex function on a large interval, so the large interval has to be subdivided into many subintervals. CM does the partition, and performs convergence checks, using an iterative approach.
> However, if the function is smooth enough, no subdivision of the integration interval is required. Those use-cases could benefit from the efficiency gain of not performing a convergence check.
> The proposal is to provide a new interface "UnivariateSingleStepIntegrator":
> {code}
> interface SingleIntervalIntegrator {
>     /**
>      * Method for implementing a single interval integration.
>      * There is no convergence checks because it is not iterative.
>      *
>      * @param f Function to integrate.
>      * @param lower Lower bound of the interval over which to integrate.
>      * @param upper Upper bound of the interval over which to integrate.
>      * @return the integrated value.
>      */
>     double integrate(UnivariateFunction f,
>                      double lower,
>                      double upper);
> }
> {code}
> In effect, the implementation of the above "integrate" method of a new "LegendreGaussIntegratorSingleStepIntegrator" would the equivalent of "stage(1)" in the current "LegendreGaussIntegrator".

--
This message is automatically generated by JIRA.
If you think it was sent incorrectly, please contact your JIRA administrators: https://issues.apache.org/jira/secure/ContactAdministrators!default.jspa
For more information on JIRA, see: http://www.atlassian.com/software/jira


Reply | Threaded
Open this post in threaded view
|

[jira] [Commented] (MATH-797) Single step integrators

ASF GitHub Bot (Jira)
In reply to this post by ASF GitHub Bot (Jira)

    [ https://issues.apache.org/jira/browse/MATH-797?page=com.atlassian.jira.plugin.system.issuetabpanels:comment-tabpanel&focusedCommentId=13401373#comment-13401373 ]

Gilles commented on MATH-797:
-----------------------------

Sébastien,

Do you intend to work on this (in which case we'll change the assignee), or shall I implement the "single-step Legendre-Gauss" part (using the code provided by the team I referred to above)?

               

> Single step integrators
> -----------------------
>
>                 Key: MATH-797
>                 URL: https://issues.apache.org/jira/browse/MATH-797
>             Project: Commons Math
>          Issue Type: Wish
>    Affects Versions: 3.0
>            Reporter: Gilles
>            Assignee: Gilles
>            Priority: Trivial
>             Fix For: 3.1
>
>
> CM assumes that the user wants to integrate a complex function on a large interval, so the large interval has to be subdivided into many subintervals. CM does the partition, and performs convergence checks, using an iterative approach.
> However, if the function is smooth enough, no subdivision of the integration interval is required. Those use-cases could benefit from the efficiency gain of not performing a convergence check.
> The proposal is to provide a new interface "UnivariateSingleStepIntegrator":
> {code}
> interface SingleIntervalIntegrator {
>     /**
>      * Method for implementing a single interval integration.
>      * There is no convergence checks because it is not iterative.
>      *
>      * @param f Function to integrate.
>      * @param lower Lower bound of the interval over which to integrate.
>      * @param upper Upper bound of the interval over which to integrate.
>      * @return the integrated value.
>      */
>     double integrate(UnivariateFunction f,
>                      double lower,
>                      double upper);
> }
> {code}
> In effect, the implementation of the above "integrate" method of a new "LegendreGaussIntegratorSingleStepIntegrator" would the equivalent of "stage(1)" in the current "LegendreGaussIntegrator".

--
This message is automatically generated by JIRA.
If you think it was sent incorrectly, please contact your JIRA administrators: https://issues.apache.org/jira/secure/ContactAdministrators!default.jspa
For more information on JIRA, see: http://www.atlassian.com/software/jira


Reply | Threaded
Open this post in threaded view
|

[jira] [Commented] (MATH-797) Single step integrators

ASF GitHub Bot (Jira)
In reply to this post by ASF GitHub Bot (Jira)

    [ https://issues.apache.org/jira/browse/MATH-797?page=com.atlassian.jira.plugin.system.issuetabpanels:comment-tabpanel&focusedCommentId=13401388#comment-13401388 ]

Sébastien Brisard commented on MATH-797:
----------------------------------------

I would be happy to get involved in this, but certainly do not want to tread on anyone's toes...
I guess it all really depends on how urgently you need this feature (I'm quite busy with testing the vector implementations, right now, plus summer vacations are approaching...).

My understanding of the above conclusion was that it is better to create a new hierarchy of classes for Gauss-like integrators. Personally, I'd rather we take some time to discuss the interface common to these integrators. Implementation could then be quite fast, as I could easily refactor what I've already written, and provide various implementations (see above), including tests; also we could refactor what has already been done on your side.

So I would say: OK to get involved (but not tomorrow!), but let's discuss first !

May be to fire up the discussion, I could include the interfaces I had designed a few years back. It's just to initiate the discussion, though... So here goes. I started by defining a Gauss integration rule
{code:java}
public interface GaussianRule {
        /**
         * Copies the abscissae of all integration points in the specified array,
         * starting at the specified location. This method will generate an
         * {@link ArrayIndexOutOfBoundsException} if the specified array is not
         * large enough.
         *
         * @param x
         *            the array in which the nodes are to be copied
         * @param pos
         *            the index of the first cell to be filled
         */
        void copyNodes(final double[] x, final int pos);

        /**
         * Copies the weights of all integration points in the specified array,
         * starting at the specified location. This method will generate an
         * {@link ArrayIndexOutOfBoundsException} if the specified array is not
         * large enough.
         *
         * @param w
         *            the array in which the weights are to be copied
         * @param pos
         *            the index of the first cell to be filled
         */
        void copyWeights(final double[] w, final int pos);

        /**
         * Returns the abscissa of the <code>i</code>-th integration point.
         *
         * @param i
         *            the index of the node
         * @return the abscissa of the node
         */
        double getNode(final int i);

        /**
         * Returns the number of integration points of this quadrature rule.
         *
         * @return the number of integration points
         */
        int getNumNodes();

        /**
         * Returns the value of the weight of the <code>i</code>-th integration
         * point.
         *
         * @param i
         *            the index of the node
         * @return the weight of the node
         */
        double getWeight(final int i);

        /**
         * Returns the weighted integral of the specified function between the
         * natural bounds of this quadrature rule.
         *
         * @param f
         *            the univariate real function to be integrated
         * @return the estimate of the integral
         */
        double integrate(final UnivariateRealFunction f);

        /**
         * Returns the maximum degree of polynomials for which integration with this
         * quadrature rule is exact.
         *
         * @return the maximum degree
         */
        int getMaximumDegreeForExactness();

        /**
         * Returns the lower-bound of the integration interval for which this rule
         * is designed. Appropriate change of the integration variable may allow
         * integration on a different interval.
         *
         * @return the lower bound of the integration interval
         */
        double getRangeOfIntegrationLowerBound();

        /**
         * Returns the upper-bound of the integration interval for which this rule
         * is designed. Appropriate change of the integration variable may allow
         * integration on a different interval.
         *
         * @return the upper bound of the integration interval
         */
        double getRangeOfIntegrationUpperBound();
}
{code}

Then, each integration rule (Gauss-Legendre, Gauss-Chebyshev, etc...) is created with a factory
{code:java}
public interface GaussianRuleFactory {
    /**
     * Returns a new <code>n</code>-point Gaussian quadrature rule. If no rule
     * has been defined for the specified number of integration points, this
     * method returns <code>null</code>.
     *
     * @param n
     *            the number of integration points
     * @return the quadrature rule
     */
    GaussianRule create(final int n);
}

{code}

Example for Gauss-Legendre ({{@mathml}} tags are taglets I've written to include MathML code into Javadocs, please do not mind them).
{code:java}
package org.jquadrature.gauss;

import org.changemyname.WeightedPoint1d;

import static org.celestin.constants.DoubleConstants.ONE_HALF;

//@formatter:off
/**
 * <p>
 * This factory creates Gauss-type quadrature rule using Legendre polynomials.
 * Such a quadrature rule allows the exact calculation of integrals
 * {@mathML doc-files/gauss-legendre-01.mml} where {@inlineMathML <mi>f</mi>}
 * is a polynomial of degree lower than or equal to
 * {@mathML doc-files/gauss-legendre-02.mml}, and {@inlineMathML n} is the
 * number of integration points. The lower- and upper-bounds of the natural
 * interval of integration are {@inlineMathML <mo>-</mo><mn>1</mn>} and
 * {@inlineMathML <mn>1</mn>} in this implementation.
 * </p>
 * <p>
 * First, the roots of the Legendre polynomial are computed (using bisection)
 * with half-ulp accuracy. Legendre polynomials are evaluated using the
 * recurrence relation (
 * <a href="{@docRoot}/doc-files/references.xhtml#ABRA1964">
 * Abramowitz and Stegun, 1964
 * </a>, 22.7.10) {@mathML doc-files/gauss-legendre-03.mml}
 * </p>
 * <p>
 * Relations 25.4.29 and 22.8.5 from
 * <a href="{@docRoot}/doc-files/references.xhtml#ABRA1964">
 * Abramowitz and Stegun (1964)</a> are then combined to compute the weights
 * {@mathML doc-files/gauss-legendre-04.mml} where {@inlineMathML
 * <msub><mi>x</mi><mi>i</mi></msub>} is the {@inlineMathML
 * <mi>i</mi>}-th node, and {@inlineMathML <msub><mi>w</mi><mi>i</mi></msub>}
 * the corresponding weight. The accuracy of this computation is somewhat worse,
 * that is why all quadrature rules up to order 20 are precomputed using
 * <a href="doc-files/gauss-legendre.mac">this script</a> for the
 * <a href="http://maxima.sourceforge.net"/>Maxima</a> computer algebra system.
 * </p>
 * <p>
 * The root-finding method is based on the property of interlacing of roots of
 * orthogonal polynomials. Therefore, in order to compute the nodes and weights
 * for the {@inlineMathML <mi>n</mi>}-point quadrature rule, the roots of
 * <em>all</em> {@mathML doc-files/gauss-legendre-05.mml} Legendre
 * polynomials must first be computed. These roots are saved for further use.
 * </p>
 *
 * @author S&eacute;bastien Brisard
 */
//@formatter:on
public class GaussLegendreFactory extends DefaultGaussianRuleFactory {
        /**
         * The unique instance of this class (singleton pattern).
         */
        private static GaussianRuleFactory factory;

        /**
         * The last rule computed.
         */
        private WeightedPoint1d[] oldRule;

        /**
         * Creates a new instance of this class.
         */
        private GaussLegendreFactory() {
                super();

                // 1-point quadrature rule
                oldRule = new WeightedPoint1d[1];
                oldRule[0] = new WeightedPoint1d(0.0E0, 2.0E0);
                addRule(oldRule);

                // 2-point quadrature rule
                oldRule = new WeightedPoint1d[2];
                oldRule[0] = new WeightedPoint1d(-5.7735026918962576451E-1, 1.0E0);
                oldRule[1] = new WeightedPoint1d(5.7735026918962576451E-1, 1.0E0);
                addRule(oldRule);

                // 3-point quadrature rule
                oldRule = new WeightedPoint1d[3];
                oldRule[0] = new WeightedPoint1d(-7.7459666924148337703E-1,
                                5.5555555555555555556E-1);
                oldRule[1] = new WeightedPoint1d(0.0E0, 8.8888888888888888889E-1);
                oldRule[2] = new WeightedPoint1d(7.7459666924148337703E-1,
                                5.5555555555555555556E-1);
                addRule(oldRule);

                // 4-point quadrature rule
                oldRule = new WeightedPoint1d[4];
                oldRule[0] = new WeightedPoint1d(-8.6113631159405257523E-1,
                                3.4785484513745385737E-1);
                oldRule[1] = new WeightedPoint1d(-3.399810435848562648E-1,
                                6.5214515486254614263E-1);
                oldRule[2] = new WeightedPoint1d(3.399810435848562648E-1,
                                6.5214515486254614263E-1);
                oldRule[3] = new WeightedPoint1d(8.6113631159405257523E-1,
                                3.4785484513745385737E-1);
                addRule(oldRule);

                // 5-point quadrature rule
                oldRule = new WeightedPoint1d[5];
                oldRule[0] = new WeightedPoint1d(-9.061798459386639928E-1,
                                2.3692688505618908751E-1);
                oldRule[1] = new WeightedPoint1d(-5.3846931010568309104E-1,
                                4.7862867049936646804E-1);
                oldRule[2] = new WeightedPoint1d(0.0E0, 5.6888888888888888889E-1);
                oldRule[3] = new WeightedPoint1d(5.3846931010568309104E-1,
                                4.7862867049936646804E-1);
                oldRule[4] = new WeightedPoint1d(9.061798459386639928E-1,
                                2.3692688505618908751E-1);
                addRule(oldRule);

                // 6-point quadrature rule
                oldRule = new WeightedPoint1d[6];
                oldRule[0] = new WeightedPoint1d(-9.3246951420315202782E-1,
                                1.7132449237917034503E-1);
                oldRule[1] = new WeightedPoint1d(-6.6120938646626451366E-1,
                                3.6076157304813860757E-1);
                oldRule[2] = new WeightedPoint1d(-2.3861918608319690863E-1,
                                4.6791393457269104739E-1);
                oldRule[3] = new WeightedPoint1d(2.3861918608319690863E-1,
                                4.6791393457269104739E-1);
                oldRule[4] = new WeightedPoint1d(6.6120938646626451366E-1,
                                3.6076157304813860757E-1);
                oldRule[5] = new WeightedPoint1d(9.3246951420315202782E-1,
                                1.7132449237917034503E-1);
                addRule(oldRule);

                // 7-point quadrature rule
                oldRule = new WeightedPoint1d[7];
                oldRule[0] = new WeightedPoint1d(-9.4910791234275852452E-1,
                                1.2948496616886969328E-1);
                oldRule[1] = new WeightedPoint1d(-7.4153118559939443987E-1,
                                2.797053914892766679E-1);
                oldRule[2] = new WeightedPoint1d(-4.0584515137739716691E-1,
                                3.8183005050511894495E-1);
                oldRule[3] = new WeightedPoint1d(0.0E0, 4.1795918367346938776E-1);
                oldRule[4] = new WeightedPoint1d(4.0584515137739716691E-1,
                                3.8183005050511894495E-1);
                oldRule[5] = new WeightedPoint1d(7.4153118559939443987E-1,
                                2.797053914892766679E-1);
                oldRule[6] = new WeightedPoint1d(9.4910791234275852452E-1,
                                1.2948496616886969328E-1);
                addRule(oldRule);

                // 8-point quadrature rule
                oldRule = new WeightedPoint1d[8];
                oldRule[0] = new WeightedPoint1d(-9.6028985649753623169E-1,
                                1.0122853629037625915E-1);
                oldRule[1] = new WeightedPoint1d(-7.9666647741362673959E-1,
                                2.2238103445337447055E-1);
                oldRule[2] = new WeightedPoint1d(-5.2553240991632898582E-1,
                                3.1370664587788728734E-1);
                oldRule[3] = new WeightedPoint1d(-1.8343464249564980494E-1,
                                3.6268378337836198297E-1);
                oldRule[4] = new WeightedPoint1d(1.8343464249564980494E-1,
                                3.6268378337836198297E-1);
                oldRule[5] = new WeightedPoint1d(5.2553240991632898582E-1,
                                3.1370664587788728734E-1);
                oldRule[6] = new WeightedPoint1d(7.9666647741362673959E-1,
                                2.2238103445337447055E-1);
                oldRule[7] = new WeightedPoint1d(9.6028985649753623169E-1,
                                1.0122853629037625915E-1);
                addRule(oldRule);

                // 9-point quadrature rule
                oldRule = new WeightedPoint1d[9];
                oldRule[0] = new WeightedPoint1d(-9.6816023950762608983E-1,
                                8.1274388361574411976E-2);
                oldRule[1] = new WeightedPoint1d(-8.360311073266357943E-1,
                                1.8064816069485740406E-1);
                oldRule[2] = new WeightedPoint1d(-6.1337143270059039731E-1,
                                2.6061069640293546232E-1);
                oldRule[3] = new WeightedPoint1d(-3.2425342340380892904E-1,
                                3.1234707704000284007E-1);
                oldRule[4] = new WeightedPoint1d(0.0E0, 3.3023935500125976317E-1);
                oldRule[5] = new WeightedPoint1d(3.2425342340380892904E-1,
                                3.1234707704000284007E-1);
                oldRule[6] = new WeightedPoint1d(6.1337143270059039731E-1,
                                2.6061069640293546232E-1);
                oldRule[7] = new WeightedPoint1d(8.360311073266357943E-1,
                                1.8064816069485740406E-1);
                oldRule[8] = new WeightedPoint1d(9.6816023950762608983E-1,
                                8.1274388361574411976E-2);
                addRule(oldRule);

                // 10-point quadrature rule
                oldRule = new WeightedPoint1d[10];
                oldRule[0] = new WeightedPoint1d(-9.7390652851717172008E-1,
                                6.6671344308688137597E-2);
                oldRule[1] = new WeightedPoint1d(-8.6506336668898451074E-1,
                                1.4945134915058059314E-1);
                oldRule[2] = new WeightedPoint1d(-6.7940956829902440623E-1,
                                2.19086362515982044E-1);
                oldRule[3] = new WeightedPoint1d(-4.333953941292471908E-1,
                                2.6926671930999635509E-1);
                oldRule[4] = new WeightedPoint1d(-1.4887433898163121088E-1,
                                2.9552422471475287018E-1);
                oldRule[5] = new WeightedPoint1d(1.4887433898163121088E-1,
                                2.9552422471475287018E-1);
                oldRule[6] = new WeightedPoint1d(4.333953941292471908E-1,
                                2.6926671930999635509E-1);
                oldRule[7] = new WeightedPoint1d(6.7940956829902440623E-1,
                                2.19086362515982044E-1);
                oldRule[8] = new WeightedPoint1d(8.6506336668898451074E-1,
                                1.4945134915058059314E-1);
                oldRule[9] = new WeightedPoint1d(9.7390652851717172008E-1,
                                6.6671344308688137597E-2);
                addRule(oldRule);

                // 11-point quadrature rule
                oldRule = new WeightedPoint1d[11];
                oldRule[0] = new WeightedPoint1d(-9.7822865814605699281E-1,
                                5.5668567116173666479E-2);
                oldRule[1] = new WeightedPoint1d(-8.8706259976809529907E-1,
                                1.2558036946490462464E-1);
                oldRule[2] = new WeightedPoint1d(-7.3015200557404932409E-1,
                                1.8629021092773425143E-1);
                oldRule[3] = new WeightedPoint1d(-5.1909612920681181592E-1,
                                2.3319376459199047992E-1);
                oldRule[4] = new WeightedPoint1d(-2.6954315595234497233E-1,
                                2.6280454451024666218E-1);
                oldRule[5] = new WeightedPoint1d(0.0E0, 2.7292508677790063072E-1);
                oldRule[6] = new WeightedPoint1d(2.6954315595234497233E-1,
                                2.6280454451024666218E-1);
                oldRule[7] = new WeightedPoint1d(5.1909612920681181592E-1,
                                2.3319376459199047992E-1);
                oldRule[8] = new WeightedPoint1d(7.3015200557404932409E-1,
                                1.8629021092773425143E-1);
                oldRule[9] = new WeightedPoint1d(8.8706259976809529907E-1,
                                1.2558036946490462464E-1);
                oldRule[10] = new WeightedPoint1d(9.7822865814605699281E-1,
                                5.5668567116173666479E-2);
                addRule(oldRule);

                // 12-point quadrature rule
                oldRule = new WeightedPoint1d[12];
                oldRule[0] = new WeightedPoint1d(-9.8156063424671925069E-1,
                                4.71753363865118272E-2);
                oldRule[1] = new WeightedPoint1d(-9.0411725637047485668E-1,
                                1.0693932599531843096E-1);
                oldRule[2] = new WeightedPoint1d(-7.6990267419430468704E-1,
                                1.6007832854334622634E-1);
                oldRule[3] = new WeightedPoint1d(-5.873179542866174473E-1,
                                2.0316742672306592175E-1);
                oldRule[4] = new WeightedPoint1d(-3.6783149899818019376E-1,
                                2.3349253653835480876E-1);
                oldRule[5] = new WeightedPoint1d(-1.2523340851146891547E-1,
                                2.49147045813402785E-1);
                oldRule[6] = new WeightedPoint1d(1.2523340851146891547E-1,
                                2.49147045813402785E-1);
                oldRule[7] = new WeightedPoint1d(3.6783149899818019376E-1,
                                2.3349253653835480876E-1);
                oldRule[8] = new WeightedPoint1d(5.873179542866174473E-1,
                                2.0316742672306592175E-1);
                oldRule[9] = new WeightedPoint1d(7.6990267419430468704E-1,
                                1.6007832854334622634E-1);
                oldRule[10] = new WeightedPoint1d(9.0411725637047485668E-1,
                                1.0693932599531843096E-1);
                oldRule[11] = new WeightedPoint1d(9.8156063424671925069E-1,
                                4.71753363865118272E-2);
                addRule(oldRule);

                // 13-point quadrature rule
                oldRule = new WeightedPoint1d[13];
                oldRule[0] = new WeightedPoint1d(-9.8418305471858814948E-1,
                                4.0484004765315879512E-2);
                oldRule[1] = new WeightedPoint1d(-9.1759839922297796521E-1,
                                9.2121499837728447916E-2);
                oldRule[2] = new WeightedPoint1d(-8.015780907333099128E-1,
                                1.3887351021978723846E-1);
                oldRule[3] = new WeightedPoint1d(-6.4234933944034022064E-1,
                                1.7814598076194573828E-1);
                oldRule[4] = new WeightedPoint1d(-4.4849275103644685288E-1,
                                2.0781604753688850231E-1);
                oldRule[5] = new WeightedPoint1d(-2.3045831595513479406E-1,
                                2.2628318026289723841E-1);
                oldRule[6] = new WeightedPoint1d(0.0E0, 2.3255155323087391019E-1);
                oldRule[7] = new WeightedPoint1d(2.3045831595513479406E-1,
                                2.2628318026289723841E-1);
                oldRule[8] = new WeightedPoint1d(4.4849275103644685288E-1,
                                2.0781604753688850231E-1);
                oldRule[9] = new WeightedPoint1d(6.4234933944034022064E-1,
                                1.7814598076194573828E-1);
                oldRule[10] = new WeightedPoint1d(8.015780907333099128E-1,
                                1.3887351021978723846E-1);
                oldRule[11] = new WeightedPoint1d(9.1759839922297796521E-1,
                                9.2121499837728447916E-2);
                oldRule[12] = new WeightedPoint1d(9.8418305471858814948E-1,
                                4.0484004765315879512E-2);
                addRule(oldRule);

                // 14-point quadrature rule
                oldRule = new WeightedPoint1d[14];
                oldRule[0] = new WeightedPoint1d(-9.8628380869681233884E-1,
                                3.5119460331751863032E-2);
                oldRule[1] = new WeightedPoint1d(-9.2843488366357351733E-1,
                                8.0158087159760209808E-2);
                oldRule[2] = new WeightedPoint1d(-8.2720131506976499319E-1,
                                1.2151857068790318469E-1);
                oldRule[3] = new WeightedPoint1d(-6.8729290481168547015E-1,
                                1.5720316715819353457E-1);
                oldRule[4] = new WeightedPoint1d(-5.1524863635815409197E-1,
                                1.8553839747793781374E-1);
                oldRule[5] = new WeightedPoint1d(-3.1911236892788976044E-1,
                                2.0519846372129560397E-1);
                oldRule[6] = new WeightedPoint1d(-1.0805494870734366207E-1,
                                2.152638534631577902E-1);
                oldRule[7] = new WeightedPoint1d(1.0805494870734366207E-1,
                                2.152638534631577902E-1);
                oldRule[8] = new WeightedPoint1d(3.1911236892788976044E-1,
                                2.0519846372129560397E-1);
                oldRule[9] = new WeightedPoint1d(5.1524863635815409197E-1,
                                1.8553839747793781374E-1);
                oldRule[10] = new WeightedPoint1d(6.8729290481168547015E-1,
                                1.5720316715819353457E-1);
                oldRule[11] = new WeightedPoint1d(8.2720131506976499319E-1,
                                1.2151857068790318469E-1);
                oldRule[12] = new WeightedPoint1d(9.2843488366357351733E-1,
                                8.0158087159760209808E-2);
                oldRule[13] = new WeightedPoint1d(9.8628380869681233884E-1,
                                3.5119460331751863032E-2);
                addRule(oldRule);

                // 15-point quadrature rule
                oldRule = new WeightedPoint1d[15];
                oldRule[0] = new WeightedPoint1d(-9.8799251802048542849E-1,
                                3.0753241996117268357E-2);
                oldRule[1] = new WeightedPoint1d(-9.3727339240070590431E-1,
                                7.0366047488108124707E-2);
                oldRule[2] = new WeightedPoint1d(-8.482065834104272162E-1,
                                1.0715922046717193501E-1);
                oldRule[3] = new WeightedPoint1d(-7.2441773136017004742E-1,
                                1.3957067792615431445E-1);
                oldRule[4] = new WeightedPoint1d(-5.7097217260853884754E-1,
                                1.6626920581699393355E-1);
                oldRule[5] = new WeightedPoint1d(-3.941513470775633699E-1,
                                1.8616100001556221103E-1);
                oldRule[6] = new WeightedPoint1d(-2.011940939974345223E-1,
                                1.9843148532711157646E-1);
                oldRule[7] = new WeightedPoint1d(0.0E0, 2.0257824192556127288E-1);
                oldRule[8] = new WeightedPoint1d(2.011940939974345223E-1,
                                1.9843148532711157646E-1);
                oldRule[9] = new WeightedPoint1d(3.941513470775633699E-1,
                                1.8616100001556221103E-1);
                oldRule[10] = new WeightedPoint1d(5.7097217260853884754E-1,
                                1.6626920581699393355E-1);
                oldRule[11] = new WeightedPoint1d(7.2441773136017004742E-1,
                                1.3957067792615431445E-1);
                oldRule[12] = new WeightedPoint1d(8.482065834104272162E-1,
                                1.0715922046717193501E-1);
                oldRule[13] = new WeightedPoint1d(9.3727339240070590431E-1,
                                7.0366047488108124707E-2);
                oldRule[14] = new WeightedPoint1d(9.8799251802048542849E-1,
                                3.0753241996117268357E-2);
                addRule(oldRule);

                // 16-point quadrature rule
                oldRule = new WeightedPoint1d[16];
                oldRule[0] = new WeightedPoint1d(-9.894009349916499326E-1,
                                2.7152459411754094847E-2);
                oldRule[1] = new WeightedPoint1d(-9.4457502307323257608E-1,
                                6.225352393864789286E-2);
                oldRule[2] = new WeightedPoint1d(-8.6563120238783174388E-1,
                                9.5158511682492784808E-2);
                oldRule[3] = new WeightedPoint1d(-7.5540440835500303389E-1,
                                1.2462897125553387205E-1);
                oldRule[4] = new WeightedPoint1d(-6.1787624440264374845E-1,
                                1.4959598881657673208E-1);
                oldRule[5] = new WeightedPoint1d(-4.5801677765722738634E-1,
                                1.6915651939500253819E-1);
                oldRule[6] = new WeightedPoint1d(-2.8160355077925891323E-1,
                                1.8260341504492358887E-1);
                oldRule[7] = new WeightedPoint1d(-9.5012509837637440184E-2,
                                1.8945061045506849629E-1);
                oldRule[8] = new WeightedPoint1d(9.5012509837637440184E-2,
                                1.8945061045506849629E-1);
                oldRule[9] = new WeightedPoint1d(2.8160355077925891323E-1,
                                1.8260341504492358887E-1);
                oldRule[10] = new WeightedPoint1d(4.5801677765722738634E-1,
                                1.6915651939500253819E-1);
                oldRule[11] = new WeightedPoint1d(6.1787624440264374845E-1,
                                1.4959598881657673208E-1);
                oldRule[12] = new WeightedPoint1d(7.5540440835500303389E-1,
                                1.2462897125553387205E-1);
                oldRule[13] = new WeightedPoint1d(8.6563120238783174388E-1,
                                9.5158511682492784808E-2);
                oldRule[14] = new WeightedPoint1d(9.4457502307323257608E-1,
                                6.225352393864789286E-2);
                oldRule[15] = new WeightedPoint1d(9.894009349916499326E-1,
                                2.7152459411754094847E-2);
                addRule(oldRule);

                // 17-point quadrature rule
                oldRule = new WeightedPoint1d[17];
                oldRule[0] = new WeightedPoint1d(-9.9057547531441733567E-1,
                                2.4148302868547931965E-2);
                oldRule[1] = new WeightedPoint1d(-9.5067552176876776122E-1,
                                5.5459529373987201132E-2);
                oldRule[2] = new WeightedPoint1d(-8.8023915372698590213E-1,
                                8.5036148317179180882E-2);
                oldRule[3] = new WeightedPoint1d(-7.8151400389680140692E-1,
                                1.118838471934039711E-1);
                oldRule[4] = new WeightedPoint1d(-6.5767115921669076585E-1,
                                1.3513636846852547329E-1);
                oldRule[5] = new WeightedPoint1d(-5.1269053708647696789E-1,
                                1.5404576107681028808E-1);
                oldRule[6] = new WeightedPoint1d(-3.512317634538763153E-1,
                                1.6800410215645004451E-1);
                oldRule[7] = new WeightedPoint1d(-1.7848418149584785585E-1,
                                1.7656270536699264633E-1);
                oldRule[8] = new WeightedPoint1d(0.0E0, 1.7944647035620652546E-1);
                oldRule[9] = new WeightedPoint1d(1.7848418149584785585E-1,
                                1.7656270536699264633E-1);
                oldRule[10] = new WeightedPoint1d(3.512317634538763153E-1,
                                1.6800410215645004451E-1);
                oldRule[11] = new WeightedPoint1d(5.1269053708647696789E-1,
                                1.5404576107681028808E-1);
                oldRule[12] = new WeightedPoint1d(6.5767115921669076585E-1,
                                1.3513636846852547329E-1);
                oldRule[13] = new WeightedPoint1d(7.8151400389680140692E-1,
                                1.118838471934039711E-1);
                oldRule[14] = new WeightedPoint1d(8.8023915372698590213E-1,
                                8.5036148317179180882E-2);
                oldRule[15] = new WeightedPoint1d(9.5067552176876776122E-1,
                                5.5459529373987201132E-2);
                oldRule[16] = new WeightedPoint1d(9.9057547531441733567E-1,
                                2.4148302868547931965E-2);
                addRule(oldRule);

                // 18-point quadrature rule
                oldRule = new WeightedPoint1d[18];
                oldRule[0] = new WeightedPoint1d(-9.9156516842093094673E-1,
                                2.1616013526483310319E-2);
                oldRule[1] = new WeightedPoint1d(-9.5582394957139775518E-1,
                                4.9714548894969796457E-2);
                oldRule[2] = new WeightedPoint1d(-8.926024664975557392E-1,
                                7.6425730254889056531E-2);
                oldRule[3] = new WeightedPoint1d(-8.0370495897252311568E-1,
                                1.0094204410628716556E-1);
                oldRule[4] = new WeightedPoint1d(-6.9168704306035320787E-1,
                                1.2255520671147846019E-1);
                oldRule[5] = new WeightedPoint1d(-5.5977083107394753461E-1,
                                1.406429146706506512E-1);
                oldRule[6] = new WeightedPoint1d(-4.1175116146284264604E-1,
                                1.5468467512626524493E-1);
                oldRule[7] = new WeightedPoint1d(-2.5188622569150550959E-1,
                                1.6427648374583272299E-1);
                oldRule[8] = new WeightedPoint1d(-8.4775013041735301242E-2,
                                1.6914238296314359184E-1);
                oldRule[9] = new WeightedPoint1d(8.4775013041735301242E-2,
                                1.6914238296314359184E-1);
                oldRule[10] = new WeightedPoint1d(2.5188622569150550959E-1,
                                1.6427648374583272299E-1);
                oldRule[11] = new WeightedPoint1d(4.1175116146284264604E-1,
                                1.5468467512626524493E-1);
                oldRule[12] = new WeightedPoint1d(5.5977083107394753461E-1,
                                1.406429146706506512E-1);
                oldRule[13] = new WeightedPoint1d(6.9168704306035320787E-1,
                                1.2255520671147846019E-1);
                oldRule[14] = new WeightedPoint1d(8.0370495897252311568E-1,
                                1.0094204410628716556E-1);
                oldRule[15] = new WeightedPoint1d(8.926024664975557392E-1,
                                7.6425730254889056531E-2);
                oldRule[16] = new WeightedPoint1d(9.5582394957139775518E-1,
                                4.9714548894969796457E-2);
                oldRule[17] = new WeightedPoint1d(9.9156516842093094673E-1,
                                2.1616013526483310319E-2);
                addRule(oldRule);

                // 19-point quadrature rule
                oldRule = new WeightedPoint1d[19];
                oldRule[0] = new WeightedPoint1d(-9.9240684384358440319E-1,
                                1.9461788229726477028E-2);
                oldRule[1] = new WeightedPoint1d(-9.6020815213483003085E-1,
                                4.4814226765699600332E-2);
                oldRule[2] = new WeightedPoint1d(-9.0315590361481790165E-1,
                                6.9044542737641226579E-2);
                oldRule[3] = new WeightedPoint1d(-8.2271465653714282498E-1,
                                9.1490021622449999465E-2);
                oldRule[4] = new WeightedPoint1d(-7.2096617733522937862E-1,
                                1.1156664554733399472E-1);
                oldRule[5] = new WeightedPoint1d(-6.0054530466168102347E-1,
                                1.2875396253933622768E-1);
                oldRule[6] = new WeightedPoint1d(-4.6457074137596094572E-1,
                                1.4260670217360661178E-1);
                oldRule[7] = new WeightedPoint1d(-3.1656409996362983199E-1,
                                1.5276604206585966678E-1);
                oldRule[8] = new WeightedPoint1d(-1.6035864564022537587E-1,
                                1.5896884339395434765E-1);
                oldRule[9] = new WeightedPoint1d(0.0E0, 1.6105444984878369598E-1);
                oldRule[10] = new WeightedPoint1d(1.6035864564022537587E-1,
                                1.5896884339395434765E-1);
                oldRule[11] = new WeightedPoint1d(3.1656409996362983199E-1,
                                1.5276604206585966678E-1);
                oldRule[12] = new WeightedPoint1d(4.6457074137596094572E-1,
                                1.4260670217360661178E-1);
                oldRule[13] = new WeightedPoint1d(6.0054530466168102347E-1,
                                1.2875396253933622768E-1);
                oldRule[14] = new WeightedPoint1d(7.2096617733522937862E-1,
                                1.1156664554733399472E-1);
                oldRule[15] = new WeightedPoint1d(8.2271465653714282498E-1,
                                9.1490021622449999465E-2);
                oldRule[16] = new WeightedPoint1d(9.0315590361481790165E-1,
                                6.9044542737641226579E-2);
                oldRule[17] = new WeightedPoint1d(9.6020815213483003085E-1,
                                4.4814226765699600332E-2);
                oldRule[18] = new WeightedPoint1d(9.9240684384358440319E-1,
                                1.9461788229726477028E-2);
                addRule(oldRule);

                // 20-point quadrature rule
                oldRule = new WeightedPoint1d[20];
                oldRule[0] = new WeightedPoint1d(-9.9312859918509492479E-1,
                                1.7614007139152118312E-2);
                oldRule[1] = new WeightedPoint1d(-9.6397192727791379127E-1,
                                4.0601429800386941329E-2);
                oldRule[2] = new WeightedPoint1d(-9.1223442825132590587E-1,
                                6.267204833410906357E-2);
                oldRule[3] = new WeightedPoint1d(-8.391169718222188234E-1,
                                8.3276741576704748724E-2);
                oldRule[4] = new WeightedPoint1d(-7.4633190646015079262E-1,
                                1.0193011981724043504E-1);
                oldRule[5] = new WeightedPoint1d(-6.3605368072651502545E-1,
                                1.1819453196151841731E-1);
                oldRule[6] = new WeightedPoint1d(-5.10867001950827098E-1,
                                1.316886384491766269E-1);
                oldRule[7] = new WeightedPoint1d(-3.7370608871541956067E-1,
                                1.4209610931838205133E-1);
                oldRule[8] = new WeightedPoint1d(-2.2778585114164507808E-1,
                                1.4917298647260374679E-1);
                oldRule[9] = new WeightedPoint1d(-7.6526521133497333755E-2,
                                1.527533871307258507E-1);
                oldRule[10] = new WeightedPoint1d(7.6526521133497333755E-2,
                                1.527533871307258507E-1);
                oldRule[11] = new WeightedPoint1d(2.2778585114164507808E-1,
                                1.4917298647260374679E-1);
                oldRule[12] = new WeightedPoint1d(3.7370608871541956067E-1,
                                1.4209610931838205133E-1);
                oldRule[13] = new WeightedPoint1d(5.10867001950827098E-1,
                                1.316886384491766269E-1);
                oldRule[14] = new WeightedPoint1d(6.3605368072651502545E-1,
                                1.1819453196151841731E-1);
                oldRule[15] = new WeightedPoint1d(7.4633190646015079262E-1,
                                1.0193011981724043504E-1);
                oldRule[16] = new WeightedPoint1d(8.391169718222188234E-1,
                                8.3276741576704748724E-2);
                oldRule[17] = new WeightedPoint1d(9.1223442825132590587E-1,
                                6.267204833410906357E-2);
                oldRule[18] = new WeightedPoint1d(9.6397192727791379127E-1,
                                4.0601429800386941329E-2);
                oldRule[19] = new WeightedPoint1d(9.9312859918509492479E-1,
                                1.7614007139152118312E-2);
                addRule(oldRule);
        }

        /**
         * Returns an instance of this class.
         *
         * @return an instance
         */
        public static GaussianRuleFactory getInstance() {
                if (factory == null) {
                        factory = new GaussLegendreFactory();
                }
                return factory;
        }

        @Override
        public GaussianRule create(final int n) {
                while (oldRule.length < n) {
                        incrementOrder();
                }
                return super.create(n);
        }

        /**
         * Computes the nodes and weights of the next Gaussian quadrature rule.
         */
        private void incrementOrder() {
                // Degree of the last polynomial computed
                final int n = oldRule.length;
                final WeightedPoint1d[] newRule = new WeightedPoint1d[n + 1];
                // Lower-bound of the interval, and corresponding values of P[j-1],
                // P[j], P[j+1].
                double a, pma, pa, ppa;
                // Upper-bound of the interval, and corresponding values of P[j-1],
                // P[j], P[j+1].
                double b, pmb, pb, ppb;
                // Middle of the interval, and corresponding values of P[j-1], P[j],
                // and P[j+1].
                double x, pm, p, pp;
                final int imax = (n + 1) / 2;
                // Find i-th root of P[n+1] by bracketing.
                for (int i = 0; i < imax; i++) {
                        if (i == 0) {
                                a = -1.;
                        } else {
                                a = oldRule[i - 1].getX();
                        }
                        if (i == n) {
                                b = 1.;
                        } else {
                                b = oldRule[i].getX();
                        }
                        pma = 1.;
                        pa = a;
                        pmb = 1.;
                        pb = b;
                        for (int j = 1; j <= n; j++) {
                                // Compute P[j+1](a) and P[j+1](b)
                                ppa = ((2 * j + 1) * a * pa - j * pma) / (j + 1);
                                ppb = ((2 * j + 1) * b * pb - j * pmb) / (j + 1);
                                pma = pa;
                                pa = ppa;
                                pmb = pb;
                                pb = ppb;
                        }
                        // Now pa = P[n+1](a), and pma = P[n](a) (same holds for b).
                        x = ONE_HALF * (a + b);
                        pm = 1.;
                        p = x;
                        boolean done = false;
                        while (!done) {
                                done = b - a <= Math.ulp(x);
                                pm = 1.;
                                p = x;
                                for (int j = 1; j <= n; j++) {
                                        // Compute P[j+1](c)
                                        pp = ((2 * j + 1) * x * p - j * pm) / (j + 1);
                                        pm = p;
                                        p = pp;
                                }
                                // Now pc = P[n+1](c) and pmc = P[n](c).
                                if (!done) {
                                        if (pa * p < 0.) {
                                                b = x;
                                                pmb = pm;
                                                pb = p;
                                        } else {
                                                a = x;
                                                pma = pm;
                                                pa = p;
                                        }
                                        x = ONE_HALF * (a + b);
                                }
                        }
                        double dummy = (n + 1) * (pm - x * p);
                        double w = 2. * (1. - x * x) / (dummy * dummy);
                        newRule[i] = new WeightedPoint1d(x, w);
                        newRule[n - i] = new WeightedPoint1d(-x, w);
                }
                // If (n+1) is odd, 0 is a root.
                if (n % 2 == 0) {
                        pm = 1.;
                        for (int j = 1; j < n; j += 2) {
                                pm = -j * pm / (j + 1.);
                        }
                        double dummy = (n + 1) * pm;
                        double w = 2. / (dummy * dummy);
                        newRule[imax] = new WeightedPoint1d(0., w);
                }
                addRule(newRule);
                oldRule = newRule;
        }
}
{code}

It should be noted that Gauss-Legendre rules of any order can be computed and stored. Also, I chose the midpoint rule instead of Newton in order to guarantee the accuracy of the result.

I would be happy to read your thoughts on all that.
Sébastien
               

> Single step integrators
> -----------------------
>
>                 Key: MATH-797
>                 URL: https://issues.apache.org/jira/browse/MATH-797
>             Project: Commons Math
>          Issue Type: Wish
>    Affects Versions: 3.0
>            Reporter: Gilles
>            Assignee: Gilles
>            Priority: Trivial
>             Fix For: 3.1
>
>
> CM assumes that the user wants to integrate a complex function on a large interval, so the large interval has to be subdivided into many subintervals. CM does the partition, and performs convergence checks, using an iterative approach.
> However, if the function is smooth enough, no subdivision of the integration interval is required. Those use-cases could benefit from the efficiency gain of not performing a convergence check.
> The proposal is to provide a new interface "UnivariateSingleStepIntegrator":
> {code}
> interface SingleIntervalIntegrator {
>     /**
>      * Method for implementing a single interval integration.
>      * There is no convergence checks because it is not iterative.
>      *
>      * @param f Function to integrate.
>      * @param lower Lower bound of the interval over which to integrate.
>      * @param upper Upper bound of the interval over which to integrate.
>      * @return the integrated value.
>      */
>     double integrate(UnivariateFunction f,
>                      double lower,
>                      double upper);
> }
> {code}
> In effect, the implementation of the above "integrate" method of a new "LegendreGaussIntegratorSingleStepIntegrator" would the equivalent of "stage(1)" in the current "LegendreGaussIntegrator".

--
This message is automatically generated by JIRA.
If you think it was sent incorrectly, please contact your JIRA administrators: https://issues.apache.org/jira/secure/ContactAdministrators!default.jspa
For more information on JIRA, see: http://www.atlassian.com/software/jira


Reply | Threaded
Open this post in threaded view
|

[jira] [Commented] (MATH-797) Single step integrators

ASF GitHub Bot (Jira)
In reply to this post by ASF GitHub Bot (Jira)

    [ https://issues.apache.org/jira/browse/MATH-797?page=com.atlassian.jira.plugin.system.issuetabpanels:comment-tabpanel&focusedCommentId=13401433#comment-13401433 ]

Gilles commented on MATH-797:
-----------------------------

I'm not familiar _at all_ with those algorithms, so I wouldn't be able to comment much.

However, I'm indeed eager to propose something quite soon (because time for code changes starts to run out for some teams in the project I work for). If not, they will keep their "fork" from the current CM code (and we run the risk of implementing something that nobody will actually test).

Thus, what I would propose is to create a code that performs what is requested, using a stripped-down version of your code (i.e. no interface, minimum number of accessors, ...), maybe hiding everything in a "SingleStepGaussLegendreIntegrator" class. This will enable to quickly assess correctness (e.g. comparing with the current CM code) and performance (the reason for the reporter's request).
Then we can gradually add more flexibility/refinements ("public" factories, etc.) as we agree on their API.

How does that sound?

               

> Single step integrators
> -----------------------
>
>                 Key: MATH-797
>                 URL: https://issues.apache.org/jira/browse/MATH-797
>             Project: Commons Math
>          Issue Type: Wish
>    Affects Versions: 3.0
>            Reporter: Gilles
>            Assignee: Gilles
>            Priority: Trivial
>             Fix For: 3.1
>
>
> CM assumes that the user wants to integrate a complex function on a large interval, so the large interval has to be subdivided into many subintervals. CM does the partition, and performs convergence checks, using an iterative approach.
> However, if the function is smooth enough, no subdivision of the integration interval is required. Those use-cases could benefit from the efficiency gain of not performing a convergence check.
> The proposal is to provide a new interface "UnivariateSingleStepIntegrator":
> {code}
> interface SingleIntervalIntegrator {
>     /**
>      * Method for implementing a single interval integration.
>      * There is no convergence checks because it is not iterative.
>      *
>      * @param f Function to integrate.
>      * @param lower Lower bound of the interval over which to integrate.
>      * @param upper Upper bound of the interval over which to integrate.
>      * @return the integrated value.
>      */
>     double integrate(UnivariateFunction f,
>                      double lower,
>                      double upper);
> }
> {code}
> In effect, the implementation of the above "integrate" method of a new "LegendreGaussIntegratorSingleStepIntegrator" would the equivalent of "stage(1)" in the current "LegendreGaussIntegrator".

--
This message is automatically generated by JIRA.
If you think it was sent incorrectly, please contact your JIRA administrators: https://issues.apache.org/jira/secure/ContactAdministrators!default.jspa
For more information on JIRA, see: http://www.atlassian.com/software/jira

       
Reply | Threaded
Open this post in threaded view
|

[jira] [Commented] (MATH-797) Single step integrators

ASF GitHub Bot (Jira)
In reply to this post by ASF GitHub Bot (Jira)

    [ https://issues.apache.org/jira/browse/MATH-797?page=com.atlassian.jira.plugin.system.issuetabpanels:comment-tabpanel&focusedCommentId=13401444#comment-13401444 ]

Sébastien Brisard commented on MATH-797:
----------------------------------------

Fine by me. Only, I would not call it a {{SingleStepLegendreGaussIntegrator}}, as Gauss integrators _are_ single-step by nature (of course, you can use refined Gauss-like integrators, like Gauss-Kronrod, but that's another story). It's a shame that the name {{LegendreGaussIntegrator}} should already be taken, by a class that is not that useful...

What we could do is create {{SingleStepLegendreGaussIntegrator}} for the time being in 3.1. Then, I think the best we could do is get rid in 4.0 of the existing {{LegendreGaussIntegrator}}, rename {{SingleStepLegendreGaussIntegrator}} and include it in a hierarchy of Gauss-like integrators.

What do you think?

Sébastien
               

> Single step integrators
> -----------------------
>
>                 Key: MATH-797
>                 URL: https://issues.apache.org/jira/browse/MATH-797
>             Project: Commons Math
>          Issue Type: Wish
>    Affects Versions: 3.0
>            Reporter: Gilles
>            Assignee: Gilles
>            Priority: Trivial
>             Fix For: 3.1
>
>
> CM assumes that the user wants to integrate a complex function on a large interval, so the large interval has to be subdivided into many subintervals. CM does the partition, and performs convergence checks, using an iterative approach.
> However, if the function is smooth enough, no subdivision of the integration interval is required. Those use-cases could benefit from the efficiency gain of not performing a convergence check.
> The proposal is to provide a new interface "UnivariateSingleStepIntegrator":
> {code}
> interface SingleIntervalIntegrator {
>     /**
>      * Method for implementing a single interval integration.
>      * There is no convergence checks because it is not iterative.
>      *
>      * @param f Function to integrate.
>      * @param lower Lower bound of the interval over which to integrate.
>      * @param upper Upper bound of the interval over which to integrate.
>      * @return the integrated value.
>      */
>     double integrate(UnivariateFunction f,
>                      double lower,
>                      double upper);
> }
> {code}
> In effect, the implementation of the above "integrate" method of a new "LegendreGaussIntegratorSingleStepIntegrator" would the equivalent of "stage(1)" in the current "LegendreGaussIntegrator".

--
This message is automatically generated by JIRA.
If you think it was sent incorrectly, please contact your JIRA administrators: https://issues.apache.org/jira/secure/ContactAdministrators!default.jspa
For more information on JIRA, see: http://www.atlassian.com/software/jira


Reply | Threaded
Open this post in threaded view
|

[jira] [Commented] (MATH-797) Single step integrators

ASF GitHub Bot (Jira)
In reply to this post by ASF GitHub Bot (Jira)

    [ https://issues.apache.org/jira/browse/MATH-797?page=com.atlassian.jira.plugin.system.issuetabpanels:comment-tabpanel&focusedCommentId=13401805#comment-13401805 ]

Gilles commented on MATH-797:
-----------------------------

bq. [...] I would not call it a {{SingleStepLegendreGaussIntegrator}}, as Gauss integrators are single-step  [...]

Maybe just {{GaussIntegrator}} then?

               

> Single step integrators
> -----------------------
>
>                 Key: MATH-797
>                 URL: https://issues.apache.org/jira/browse/MATH-797
>             Project: Commons Math
>          Issue Type: Wish
>    Affects Versions: 3.0
>            Reporter: Gilles
>            Assignee: Gilles
>            Priority: Trivial
>             Fix For: 3.1
>
>
> CM assumes that the user wants to integrate a complex function on a large interval, so the large interval has to be subdivided into many subintervals. CM does the partition, and performs convergence checks, using an iterative approach.
> However, if the function is smooth enough, no subdivision of the integration interval is required. Those use-cases could benefit from the efficiency gain of not performing a convergence check.
> The proposal is to provide a new interface "UnivariateSingleStepIntegrator":
> {code}
> interface SingleIntervalIntegrator {
>     /**
>      * Method for implementing a single interval integration.
>      * There is no convergence checks because it is not iterative.
>      *
>      * @param f Function to integrate.
>      * @param lower Lower bound of the interval over which to integrate.
>      * @param upper Upper bound of the interval over which to integrate.
>      * @return the integrated value.
>      */
>     double integrate(UnivariateFunction f,
>                      double lower,
>                      double upper);
> }
> {code}
> In effect, the implementation of the above "integrate" method of a new "LegendreGaussIntegratorSingleStepIntegrator" would the equivalent of "stage(1)" in the current "LegendreGaussIntegrator".

--
This message is automatically generated by JIRA.
If you think it was sent incorrectly, please contact your JIRA administrators: https://issues.apache.org/jira/secure/ContactAdministrators!default.jspa
For more information on JIRA, see: http://www.atlassian.com/software/jira

       
Reply | Threaded
Open this post in threaded view
|

[jira] [Commented] (MATH-797) Single step integrators

ASF GitHub Bot (Jira)
In reply to this post by ASF GitHub Bot (Jira)

    [ https://issues.apache.org/jira/browse/MATH-797?page=com.atlassian.jira.plugin.system.issuetabpanels:comment-tabpanel&focusedCommentId=13401938#comment-13401938 ]

Sébastien Brisard commented on MATH-797:
----------------------------------------

Not really, no. Gauss integrator is the generic name for a bunch of integration schemes which read
{noformat}
 b                n - 1
/                 ====
[                 \
I  f(t) w(t) dt =  >    w  f(t )
]                 /      i    i
/                 ====
a                i = 0
{noformat}
where {{w(t)}} is a weighting function. The Gauss points t[i] and weights w[i] are chosen so that the above formula is *exact* for {{f(t)}} polynomial, {{deg(f) <= 2 * n - 1}}. Various choices of {{w(t)}} lead to various Gauss integration rules
* {{w(t) = 1}} is the so-called Gauss-Legendre rule,
* {{w(t) = (1 - x^2)^(-1 / 2)}}: Gauss-Chebyshev,
* etc...

So, we basically need a mother class to handle Gauss points and Gauss weights, as well as the above formula. Then, to implement a new Gauss integration rule, we would only need to derive it from this mother class.
               

> Single step integrators
> -----------------------
>
>                 Key: MATH-797
>                 URL: https://issues.apache.org/jira/browse/MATH-797
>             Project: Commons Math
>          Issue Type: Wish
>    Affects Versions: 3.0
>            Reporter: Gilles
>            Assignee: Gilles
>            Priority: Trivial
>             Fix For: 3.1
>
>
> CM assumes that the user wants to integrate a complex function on a large interval, so the large interval has to be subdivided into many subintervals. CM does the partition, and performs convergence checks, using an iterative approach.
> However, if the function is smooth enough, no subdivision of the integration interval is required. Those use-cases could benefit from the efficiency gain of not performing a convergence check.
> The proposal is to provide a new interface "UnivariateSingleStepIntegrator":
> {code}
> interface SingleIntervalIntegrator {
>     /**
>      * Method for implementing a single interval integration.
>      * There is no convergence checks because it is not iterative.
>      *
>      * @param f Function to integrate.
>      * @param lower Lower bound of the interval over which to integrate.
>      * @param upper Upper bound of the interval over which to integrate.
>      * @return the integrated value.
>      */
>     double integrate(UnivariateFunction f,
>                      double lower,
>                      double upper);
> }
> {code}
> In effect, the implementation of the above "integrate" method of a new "LegendreGaussIntegratorSingleStepIntegrator" would the equivalent of "stage(1)" in the current "LegendreGaussIntegrator".

--
This message is automatically generated by JIRA.
If you think it was sent incorrectly, please contact your JIRA administrators: https://issues.apache.org/jira/secure/ContactAdministrators!default.jspa
For more information on JIRA, see: http://www.atlassian.com/software/jira


Reply | Threaded
Open this post in threaded view
|

[jira] [Commented] (MATH-797) Single step integrators

ASF GitHub Bot (Jira)
In reply to this post by ASF GitHub Bot (Jira)

    [ https://issues.apache.org/jira/browse/MATH-797?page=com.atlassian.jira.plugin.system.issuetabpanels:comment-tabpanel&focusedCommentId=13402095#comment-13402095 ]

Gilles commented on MATH-797:
-----------------------------

Sorry, I don't get the "Not really, no."

When you say:
bq. [...] we basically need a mother class to handle Gauss points and Gauss weights, [...]
I can only suggest {{GaussIntegrator}} as the name for this class.

*IIUC*, the constructor would be something like
{code}
public GaussIntegrator(double[] points, double[] weights) { /* ... */ }
{code}

And there would be a factory:
{code}
public class GaussIntegratorFactory {
  public static GaussIntegrator legendre(int order) { /* ... */ }
  public static GaussIntegrator chebyshev(int order) { /* ... */ }
  // etc.
}
{code}

At first sight, this is all what should be in the public API. Or am I mistaken?

               

> Single step integrators
> -----------------------
>
>                 Key: MATH-797
>                 URL: https://issues.apache.org/jira/browse/MATH-797
>             Project: Commons Math
>          Issue Type: Wish
>    Affects Versions: 3.0
>            Reporter: Gilles
>            Assignee: Gilles
>            Priority: Trivial
>             Fix For: 3.1
>
>
> CM assumes that the user wants to integrate a complex function on a large interval, so the large interval has to be subdivided into many subintervals. CM does the partition, and performs convergence checks, using an iterative approach.
> However, if the function is smooth enough, no subdivision of the integration interval is required. Those use-cases could benefit from the efficiency gain of not performing a convergence check.
> The proposal is to provide a new interface "UnivariateSingleStepIntegrator":
> {code}
> interface SingleIntervalIntegrator {
>     /**
>      * Method for implementing a single interval integration.
>      * There is no convergence checks because it is not iterative.
>      *
>      * @param f Function to integrate.
>      * @param lower Lower bound of the interval over which to integrate.
>      * @param upper Upper bound of the interval over which to integrate.
>      * @return the integrated value.
>      */
>     double integrate(UnivariateFunction f,
>                      double lower,
>                      double upper);
> }
> {code}
> In effect, the implementation of the above "integrate" method of a new "LegendreGaussIntegratorSingleStepIntegrator" would the equivalent of "stage(1)" in the current "LegendreGaussIntegrator".

--
This message is automatically generated by JIRA.
If you think it was sent incorrectly, please contact your JIRA administrators: https://issues.apache.org/jira/secure/ContactAdministrators!default.jspa
For more information on JIRA, see: http://www.atlassian.com/software/jira

       
Reply | Threaded
Open this post in threaded view
|

[jira] [Commented] (MATH-797) Single step integrators

ASF GitHub Bot (Jira)
In reply to this post by ASF GitHub Bot (Jira)

    [ https://issues.apache.org/jira/browse/MATH-797?page=com.atlassian.jira.plugin.system.issuetabpanels:comment-tabpanel&focusedCommentId=13402148#comment-13402148 ]

Sébastien Brisard commented on MATH-797:
----------------------------------------

Yes, that's much better than what I had in mind. I thought of something like {{GaussIntegrator}} as a mother class, with various impls: {{GaussLegendreIntegrator extends GaussIntegrator}}, {{GaussChebyshevIntegrator extends GaussIntegrator}}, and so on. That's why I was saying that {{GaussIntegrator}} was a bad name for {{GaussLegendreIntegrator}}.

But I like your way of seeing things very much. Plus it would be fairly easy to start with a basic implementation, only legendre, limited order, and build on this later on.

I think that {{GaussIntegrator}} should not implement {{UnivariateIntegrator}}, because the implicit weighting function would otherwise be misleading. What do you think?
               

> Single step integrators
> -----------------------
>
>                 Key: MATH-797
>                 URL: https://issues.apache.org/jira/browse/MATH-797
>             Project: Commons Math
>          Issue Type: Wish
>    Affects Versions: 3.0
>            Reporter: Gilles
>            Assignee: Gilles
>            Priority: Trivial
>             Fix For: 3.1
>
>
> CM assumes that the user wants to integrate a complex function on a large interval, so the large interval has to be subdivided into many subintervals. CM does the partition, and performs convergence checks, using an iterative approach.
> However, if the function is smooth enough, no subdivision of the integration interval is required. Those use-cases could benefit from the efficiency gain of not performing a convergence check.
> The proposal is to provide a new interface "UnivariateSingleStepIntegrator":
> {code}
> interface SingleIntervalIntegrator {
>     /**
>      * Method for implementing a single interval integration.
>      * There is no convergence checks because it is not iterative.
>      *
>      * @param f Function to integrate.
>      * @param lower Lower bound of the interval over which to integrate.
>      * @param upper Upper bound of the interval over which to integrate.
>      * @return the integrated value.
>      */
>     double integrate(UnivariateFunction f,
>                      double lower,
>                      double upper);
> }
> {code}
> In effect, the implementation of the above "integrate" method of a new "LegendreGaussIntegratorSingleStepIntegrator" would the equivalent of "stage(1)" in the current "LegendreGaussIntegrator".

--
This message is automatically generated by JIRA.
If you think it was sent incorrectly, please contact your JIRA administrators: https://issues.apache.org/jira/secure/ContactAdministrators!default.jspa
For more information on JIRA, see: http://www.atlassian.com/software/jira


Reply | Threaded
Open this post in threaded view
|

[jira] [Commented] (MATH-797) Single step integrators

ASF GitHub Bot (Jira)
In reply to this post by ASF GitHub Bot (Jira)

    [ https://issues.apache.org/jira/browse/MATH-797?page=com.atlassian.jira.plugin.system.issuetabpanels:comment-tabpanel&focusedCommentId=13402643#comment-13402643 ]

Gilles commented on MATH-797:
-----------------------------

bq. I think that GaussIntegrator should not implement UnivariateIntegrator, because the implicit weighting function would otherwise be misleading.

If I refer to your nice ASCII equation, it indeed seems that {{UnivariateIntegrator}} and {{GaussIntegrator}} are different concept that happen to intersect at the "Legendre point". :-)
Anyways, it will be easier to start by not trying to make the concepts depend on each other or otherwise overlap. Maybe "adapters" can be added later on.

So, are we set to go?

               

> Single step integrators
> -----------------------
>
>                 Key: MATH-797
>                 URL: https://issues.apache.org/jira/browse/MATH-797
>             Project: Commons Math
>          Issue Type: Wish
>    Affects Versions: 3.0
>            Reporter: Gilles
>            Assignee: Gilles
>            Priority: Trivial
>             Fix For: 3.1
>
>
> CM assumes that the user wants to integrate a complex function on a large interval, so the large interval has to be subdivided into many subintervals. CM does the partition, and performs convergence checks, using an iterative approach.
> However, if the function is smooth enough, no subdivision of the integration interval is required. Those use-cases could benefit from the efficiency gain of not performing a convergence check.
> The proposal is to provide a new interface "UnivariateSingleStepIntegrator":
> {code}
> interface SingleIntervalIntegrator {
>     /**
>      * Method for implementing a single interval integration.
>      * There is no convergence checks because it is not iterative.
>      *
>      * @param f Function to integrate.
>      * @param lower Lower bound of the interval over which to integrate.
>      * @param upper Upper bound of the interval over which to integrate.
>      * @return the integrated value.
>      */
>     double integrate(UnivariateFunction f,
>                      double lower,
>                      double upper);
> }
> {code}
> In effect, the implementation of the above "integrate" method of a new "LegendreGaussIntegratorSingleStepIntegrator" would the equivalent of "stage(1)" in the current "LegendreGaussIntegrator".

--
This message is automatically generated by JIRA.
If you think it was sent incorrectly, please contact your JIRA administrators: https://issues.apache.org/jira/secure/ContactAdministrators!default.jspa
For more information on JIRA, see: http://www.atlassian.com/software/jira

       
Reply | Threaded
Open this post in threaded view
|

[jira] [Commented] (MATH-797) Single step integrators

ASF GitHub Bot (Jira)
In reply to this post by ASF GitHub Bot (Jira)

    [ https://issues.apache.org/jira/browse/MATH-797?page=com.atlassian.jira.plugin.system.issuetabpanels:comment-tabpanel&focusedCommentId=13402854#comment-13402854 ]

Sébastien Brisard commented on MATH-797:
----------------------------------------

{quote}
If I refer to your nice ASCII equation
{quote}
[Maxima|http://maxima.sourceforge.net/] is a very faithful friend of mine...

I now remember why I implemented {{GaussLegendreFactory}}, {{GaussChebyshevFactory}}, and the likes, instead of one big factory qs you suggested. The reason for this is that n-point rules are computed recursively: you need first to compute the (n-1), (n-2), ... 1-point rules. I therefore cached all previously computed rules, which can make one single factory messy, because you would need one cache for each rule. However
* I am not sure that caching is indeed necessary: computing a new rule will presumably not occur too often in particular applications, and we can probably allow for the CPU-cost,
* we could also think of a factory method which would take the (n-1) point rule and return the n-point rule.
* the two approaches (one big factory vs. multiple factories) is anyway tractable: we can still define multiple factories, and have one {{GaussIntegratorFactory}} which holds singleton references to each of these specific factories. I have to say I find your idea of one big single factory very neat.

{quote}
So, are we set to go?
{quote}
I think so. I guess you will start with implementing the stripped down G-L integrator, where all points and weights have been precomputed (not computed on-the-fly), so you are not likely to have the above problem.
               

> Single step integrators
> -----------------------
>
>                 Key: MATH-797
>                 URL: https://issues.apache.org/jira/browse/MATH-797
>             Project: Commons Math
>          Issue Type: Wish
>    Affects Versions: 3.0
>            Reporter: Gilles
>            Assignee: Gilles
>            Priority: Trivial
>             Fix For: 3.1
>
>
> CM assumes that the user wants to integrate a complex function on a large interval, so the large interval has to be subdivided into many subintervals. CM does the partition, and performs convergence checks, using an iterative approach.
> However, if the function is smooth enough, no subdivision of the integration interval is required. Those use-cases could benefit from the efficiency gain of not performing a convergence check.
> The proposal is to provide a new interface "UnivariateSingleStepIntegrator":
> {code}
> interface SingleIntervalIntegrator {
>     /**
>      * Method for implementing a single interval integration.
>      * There is no convergence checks because it is not iterative.
>      *
>      * @param f Function to integrate.
>      * @param lower Lower bound of the interval over which to integrate.
>      * @param upper Upper bound of the interval over which to integrate.
>      * @return the integrated value.
>      */
>     double integrate(UnivariateFunction f,
>                      double lower,
>                      double upper);
> }
> {code}
> In effect, the implementation of the above "integrate" method of a new "LegendreGaussIntegratorSingleStepIntegrator" would the equivalent of "stage(1)" in the current "LegendreGaussIntegrator".

--
This message is automatically generated by JIRA.
If you think it was sent incorrectly, please contact your JIRA administrators: https://issues.apache.org/jira/secure/ContactAdministrators!default.jspa
For more information on JIRA, see: http://www.atlassian.com/software/jira


Reply | Threaded
Open this post in threaded view
|

[jira] [Commented] (MATH-797) Single step integrators

ASF GitHub Bot (Jira)
In reply to this post by ASF GitHub Bot (Jira)

    [ https://issues.apache.org/jira/browse/MATH-797?page=com.atlassian.jira.plugin.system.issuetabpanels:comment-tabpanel&focusedCommentId=13403231#comment-13403231 ]

Gilles commented on MATH-797:
-----------------------------

I'd rather include the algorithm for computing the points and weights. So, if you don't mind (and the paperwork is all done), could you please post the code which you wanted to donate?

               

> Single step integrators
> -----------------------
>
>                 Key: MATH-797
>                 URL: https://issues.apache.org/jira/browse/MATH-797
>             Project: Commons Math
>          Issue Type: Wish
>    Affects Versions: 3.0
>            Reporter: Gilles
>            Assignee: Gilles
>            Priority: Trivial
>             Fix For: 3.1
>
>
> CM assumes that the user wants to integrate a complex function on a large interval, so the large interval has to be subdivided into many subintervals. CM does the partition, and performs convergence checks, using an iterative approach.
> However, if the function is smooth enough, no subdivision of the integration interval is required. Those use-cases could benefit from the efficiency gain of not performing a convergence check.
> The proposal is to provide a new interface "UnivariateSingleStepIntegrator":
> {code}
> interface SingleIntervalIntegrator {
>     /**
>      * Method for implementing a single interval integration.
>      * There is no convergence checks because it is not iterative.
>      *
>      * @param f Function to integrate.
>      * @param lower Lower bound of the interval over which to integrate.
>      * @param upper Upper bound of the interval over which to integrate.
>      * @return the integrated value.
>      */
>     double integrate(UnivariateFunction f,
>                      double lower,
>                      double upper);
> }
> {code}
> In effect, the implementation of the above "integrate" method of a new "LegendreGaussIntegratorSingleStepIntegrator" would the equivalent of "stage(1)" in the current "LegendreGaussIntegrator".

--
This message is automatically generated by JIRA.
If you think it was sent incorrectly, please contact your JIRA administrators: https://issues.apache.org/jira/secure/ContactAdministrators!default.jspa
For more information on JIRA, see: http://www.atlassian.com/software/jira

       
12