commons-commits mailing list archives

Site index · List index
Message view « Date » · « Thread »
Top « Date » · « Thread »
From l..@apache.org
Subject svn commit: r735876 - in /commons/proper/math/trunk/src: java/org/apache/commons/math/ java/org/apache/commons/math/analysis/integration/ mantissa/src/org/spaceroots/mantissa/quadrature/scalar/ mantissa/tests-src/org/spaceroots/mantissa/quadrature/scal...
Date Mon, 19 Jan 2009 23:40:15 GMT
Author: luc
Date: Mon Jan 19 15:40:14 2009
New Revision: 735876

URL: http://svn.apache.org/viewvc?rev=735876&view=rev
Log:
added a Legendre-Gauss integrator

Added:
    commons/proper/math/trunk/src/java/org/apache/commons/math/analysis/integration/LegendreGaussIntegrator.java
  (contents, props changed)
      - copied, changed from r735530, commons/proper/math/trunk/src/mantissa/src/org/spaceroots/mantissa/quadrature/scalar/GaussLegendreIntegrator.java
    commons/proper/math/trunk/src/test/org/apache/commons/math/analysis/integration/LegendreGaussIntegratorTest.java
  (contents, props changed)
      - copied, changed from r735530, commons/proper/math/trunk/src/mantissa/tests-src/org/spaceroots/mantissa/quadrature/scalar/GaussLegendreIntegratorTest.java
Removed:
    commons/proper/math/trunk/src/mantissa/src/org/spaceroots/mantissa/quadrature/scalar/GaussLegendreIntegrator.java
    commons/proper/math/trunk/src/mantissa/tests-src/org/spaceroots/mantissa/quadrature/scalar/GaussLegendreIntegratorTest.java
Modified:
    commons/proper/math/trunk/src/java/org/apache/commons/math/MessagesResources_fr.java
    commons/proper/math/trunk/src/site/xdoc/changes.xml
    commons/proper/math/trunk/src/site/xdoc/userguide/analysis.xml

Modified: commons/proper/math/trunk/src/java/org/apache/commons/math/MessagesResources_fr.java
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/java/org/apache/commons/math/MessagesResources_fr.java?rev=735876&r1=735875&r2=735876&view=diff
==============================================================================
--- commons/proper/math/trunk/src/java/org/apache/commons/math/MessagesResources_fr.java (original)
+++ commons/proper/math/trunk/src/java/org/apache/commons/math/MessagesResources_fr.java Mon
Jan 19 15:40:14 2009
@@ -300,6 +300,12 @@
    { "invalid iteration limits: min={0}, max={1}",
      "limites d''it\u00e9rations invalides : min = {0}, max = {1}" },
 
+   // org.apache.commons.math.analysis.integration.LegendreGaussIntegrator
+   { "{0} points Legendre-Gauss integrator not supported," +
+     " number of points must be in the {1}-{2} range",
+     "int\u00e9grateur de Legendre-Gauss non support\u00e9 en {0} points, " +
+     "le nombre de points doit \u00eatre entre {1} et {2}" },
+
    // org.apache.commons.math.fraction.Fraction
    { "zero denominator in fraction {0}/{1}",
      "d\u00e9nominateur null dans le nombre rationnel {0}/{1}" },

Copied: commons/proper/math/trunk/src/java/org/apache/commons/math/analysis/integration/LegendreGaussIntegrator.java
(from r735530, commons/proper/math/trunk/src/mantissa/src/org/spaceroots/mantissa/quadrature/scalar/GaussLegendreIntegrator.java)
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/java/org/apache/commons/math/analysis/integration/LegendreGaussIntegrator.java?p2=commons/proper/math/trunk/src/java/org/apache/commons/math/analysis/integration/LegendreGaussIntegrator.java&p1=commons/proper/math/trunk/src/mantissa/src/org/spaceroots/mantissa/quadrature/scalar/GaussLegendreIntegrator.java&r1=735530&r2=735876&rev=735876&view=diff
==============================================================================
--- commons/proper/math/trunk/src/mantissa/src/org/spaceroots/mantissa/quadrature/scalar/GaussLegendreIntegrator.java
(original)
+++ commons/proper/math/trunk/src/java/org/apache/commons/math/analysis/integration/LegendreGaussIntegrator.java
Mon Jan 19 15:40:14 2009
@@ -1,151 +1,238 @@
-// Licensed to the Apache Software Foundation (ASF) under one
-// or more contributor license agreements.  See the NOTICE file
-// distributed with this work for additional information
-// regarding copyright ownership.  The ASF licenses this file
-// to you under the Apache License, Version 2.0 (the
-// "License"); you may not use this file except in compliance
-// with the License.  You may obtain a copy of the License at
-// 
-//   http://www.apache.org/licenses/LICENSE-2.0
-// 
-// Unless required by applicable law or agreed to in writing,
-// software distributed under the License is distributed on an
-// "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY
-// KIND, either express or implied.  See the License for the
-// specific language governing permissions and limitations
-// under the License.
-
-package org.spaceroots.mantissa.quadrature.scalar;
-
-import org.spaceroots.mantissa.functions.scalar.ComputableFunction;
-import org.spaceroots.mantissa.functions.FunctionException;
-
-/** This class implements a Gauss-Legendre integrator.
+/*
+ * Licensed to the Apache Software Foundation (ASF) under one or more
+ * contributor license agreements.  See the NOTICE file distributed with
+ * this work for additional information regarding copyright ownership.
+ * The ASF licenses this file to You under the Apache License, Version 2.0
+ * (the "License"); you may not use this file except in compliance with
+ * the License.  You may obtain a copy of the License at
+ *
+ *      http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+package org.apache.commons.math.analysis.integration;
 
- * <p>Gauss-Legendre integrators are efficient integrators that can
+import org.apache.commons.math.ConvergenceException;
+import org.apache.commons.math.FunctionEvaluationException;
+import org.apache.commons.math.MathRuntimeException;
+import org.apache.commons.math.MaxIterationsExceededException;
+import org.apache.commons.math.analysis.UnivariateRealFunction;
+
+/**
+ * Implements the <a href="http://mathworld.wolfram.com/Legendre-GaussQuadrature.html">
+ * Legendre-Gauss</a> quadrature formula.
+ * <p>
+ * Legendre-Gauss integrators are efficient integrators that can
  * accurately integrate functions with few functions evaluations. A
- * Gauss-Legendre integrator using an n-points quadrature formula can
- * integrate exactly 2n-1 degree polynoms.</p>
+ * Legendre-Gauss integrator using an n-points quadrature formula can
+ * integrate exactly 2n-1 degree polynomialss.
+ * </p>
+ * <p>
+ * These integrators evaluate the function on n carefully chosen
+ * abscissas in each step interval (mapped to the canonical [-1  1] interval).
+ * The evaluation abscissas are not evenly spaced and none of them are
+ * at the interval endpoints. This implies the function integrated can be
+ * undefined at integration interval endpoints.
+ * </p>
+ * <p>
+ * The evaluation abscissas x<sub>i</sub> are the roots of the degree n
+ * Legendre polynomial. The weights a<sub>i</sub> of the quadrature formula
+ * integrals from -1 to +1 &int; Li<sup>2</sup> where Li (x) =
+ * &prod; (x-x<sub>k</sub>)/(x<sub>i</sub>-x<sub>k</sub>)
for k != i.
+ * </p>
+ * <p>
+ * @version $Revision$ $Date$
+ * @since 1.2
+ */
 
- * <p>These integrators evaluate the function on n carefully chosen
- * points in each step interval. These points are not evenly
- * spaced. The function is <emph>never</emph> evaluated at the
- * boundary points, which means it can be undefined at these
- * points.</p>
+public class LegendreGaussIntegrator extends UnivariateRealIntegratorImpl {
 
- * @version $Id$
- * @author L. Maisonobe
+    /** Serializable version identifier. */
+    private static final long serialVersionUID = -331962723352824098L;
 
- */
+    /** Abscissas for the 2 points method. */
+    private static final double[] ABSCISSAS_2 = {
+        -1.0 / Math.sqrt(3.0),
+         1.0 / Math.sqrt(3.0)
+    };
+
+    /** Weights for the 2 points method. */
+    private static final double[] WEIGHTS_2 = {
+        1.0,
+        1.0
+    };
+
+    /** Abscissas for the 3 points method. */
+    private static final double[] ABSCISSAS_3 = {
+        -Math.sqrt(0.6),
+         0.0,
+         Math.sqrt(0.6)
+    };
+
+    /** Weights for the 3 points method. */
+    private static final double[] WEIGHTS_3 = {
+        5.0 / 9.0,
+        8.0 / 9.0,
+        5.0 / 9.0
+    };
+
+    /** Abscissas for the 4 points method. */
+    private static final double[] ABSCISSAS_4 = {
+        -Math.sqrt((15.0 + 2.0 * Math.sqrt(30.0)) / 35.0),
+        -Math.sqrt((15.0 - 2.0 * Math.sqrt(30.0)) / 35.0),
+         Math.sqrt((15.0 - 2.0 * Math.sqrt(30.0)) / 35.0),
+         Math.sqrt((15.0 + 2.0 * Math.sqrt(30.0)) / 35.0)
+    };
+
+    /** Weights for the 4 points method. */
+    private static final double[] WEIGHTS_4 = {
+        (90.0 - 5.0 * Math.sqrt(30.0)) / 180.0,
+        (90.0 + 5.0 * Math.sqrt(30.0)) / 180.0,
+        (90.0 + 5.0 * Math.sqrt(30.0)) / 180.0,
+        (90.0 - 5.0 * Math.sqrt(30.0)) / 180.0
+    };
+
+    /** Abscissas for the 5 points method. */
+    private static final double[] ABSCISSAS_5 = {
+        -Math.sqrt((35.0 + 2.0 * Math.sqrt(70.0)) / 63.0),
+        -Math.sqrt((35.0 - 2.0 * Math.sqrt(70.0)) / 63.0),
+         0.0,
+         Math.sqrt((35.0 - 2.0 * Math.sqrt(70.0)) / 63.0),
+         Math.sqrt((35.0 + 2.0 * Math.sqrt(70.0)) / 63.0)
+    };
+
+    /** Weights for the 5 points method. */
+    private static final double[] WEIGHTS_5 = {
+        (322.0 - 13.0 * Math.sqrt(70.0)) / 900.0,
+        (322.0 + 13.0 * Math.sqrt(70.0)) / 900.0,
+        128.0 / 225.0,
+        (322.0 + 13.0 * Math.sqrt(70.0)) / 900.0,
+        (322.0 - 13.0 * Math.sqrt(70.0)) / 900.0
+    };
+
+    /** Abscissas for the current method. */
+    private final double[] abscissas;
+
+    /** Weights for the current method. */
+    private final double[] weights;
+
+    /** Build a Legendre-Gauss integrator.
+     * @param n number of points desired (must be between 2 and 5 inclusive)
+     * @param defaultMaximalIterationCount maximum number of iterations
+     * @exception IllegalArgumentException if the number of points is not
+     * in the supported range
+     */
+    public LegendreGaussIntegrator(final int n, final int defaultMaximalIterationCount)
+        throws IllegalArgumentException {
+        super(defaultMaximalIterationCount);
+        switch(n) {
+        case 2 :
+            abscissas = ABSCISSAS_2;
+            weights   = WEIGHTS_2;
+            break;
+        case 3 :
+            abscissas = ABSCISSAS_3;
+            weights   = WEIGHTS_3;
+            break;
+        case 4 :
+            abscissas = ABSCISSAS_4;
+            weights   = WEIGHTS_4;
+            break;
+        case 5 :
+            abscissas = ABSCISSAS_5;
+            weights   = WEIGHTS_5;
+            break;
+        default :
+            throw MathRuntimeException.createIllegalArgumentException(
+                    "{0} points Legendre-Gauss integrator not supported, " +
+                    "number of points must be in the {1}-{2} range",
+                    new Object[] { n, 2, 5 });
+        }
 
-public class GaussLegendreIntegrator
-  implements ComputableFunctionIntegrator {
-  /** Build a Gauss-Legendre integrator.
-
-   * <p>A Gauss-Legendre integrator is a formula like:
-   * <pre>
-   *    int (f) from -1 to +1 = Sum (ai * f(xi))
-   * </pre>
-   * </p>
-   *
-   * <p>The coefficients of the formula are computed as follow:
-   * <pre>
-   *   let n be the desired number of points
-   *   the xi are the roots of the degree n Legendre polynomial
-   *   the ai are the integrals int (Li^2) from -1 to +1
-   *   where Li (x) = Prod (x-xk)/(xi-xk) for k != i
-   * </pre>
-   * </p>
-   *
-   * <p>A formula in n points can integrate exactly polynoms of degree
-   * up to 2n-1.</p>
-   *
-   * @param minPoints minimal number of points desired
-   * @param rawStep raw integration step (the precise step will be
-   * adjusted in order to have an integer number of steps in the
-   * integration range).
-   * */
-  public GaussLegendreIntegrator(int minPoints, double rawStep) {
-    if (minPoints <= 2) {
-      weightedRoots = new double[][] {
-        { 1.0, -1.0 / Math.sqrt(3.0) },
-        { 1.0,  1.0 / Math.sqrt(3.0) }
-      };
-    } else if (minPoints <= 3) {
-      weightedRoots = new double[][] {
-        { 5.0 / 9.0, -Math.sqrt(0.6) },
-        { 8.0 / 9.0,            0.0  },
-        { 5.0 / 9.0,  Math.sqrt(0.6) }
-      };
-    } else if (minPoints <= 4) {
-      weightedRoots = new double[][] {
-        { (90.0 - 5.0 * Math.sqrt(30.0)) / 180.0,
-             -Math.sqrt((15.0 + 2.0 * Math.sqrt(30.0)) / 35.0) },
-        { (90.0 + 5.0 * Math.sqrt(30.0)) / 180.0,
-             -Math.sqrt((15.0 - 2.0 * Math.sqrt(30.0)) / 35.0) },
-        { (90.0 + 5.0 * Math.sqrt(30.0)) / 180.0,
-              Math.sqrt((15.0 - 2.0 * Math.sqrt(30.0)) / 35.0) },
-        { (90.0 - 5.0 * Math.sqrt(30.0)) / 180.0,
-              Math.sqrt((15.0 + 2.0 * Math.sqrt(30.0)) / 35.0) }
-      };
-    } else {
-      weightedRoots = new double[][] {
-        { (322.0 - 13.0 * Math.sqrt(70.0)) / 900.0,
-             -Math.sqrt((35.0 + 2.0 * Math.sqrt(70.0)) / 63.0) },
-        { (322.0 + 13.0 * Math.sqrt(70.0)) / 900.0,
-             -Math.sqrt((35.0 - 2.0 * Math.sqrt(70.0)) / 63.0) },
-        { 128.0 / 225.0,
-              0.0 },
-        { (322.0 + 13.0 * Math.sqrt(70.0)) / 900.0,
-              Math.sqrt((35.0 - 2.0 * Math.sqrt(70.0)) / 63.0) },
-        { (322.0 - 13.0 * Math.sqrt(70.0)) / 900.0,
-              Math.sqrt((35.0 + 2.0 * Math.sqrt(70.0)) / 63.0) }
-      };
     }
 
-    this.rawStep = rawStep;
+    /** {@inheritDoc} */
+    @Deprecated
+    public double integrate(final double min, final double max)
+        throws ConvergenceException,  FunctionEvaluationException, IllegalArgumentException
{
+        return integrate(f, min, max);
+    }
 
-  }
+    /** {@inheritDoc} */
+    public double integrate(final UnivariateRealFunction f,
+            final double min, final double max)
+        throws ConvergenceException,  FunctionEvaluationException, IllegalArgumentException
{
+        
+        clearResult();
+        verifyInterval(min, max);
+        verifyIterationCount();
+
+        // compute first estimate with a single step
+        double oldt = stage(f, min, max, 1);
+
+        int n = 2;
+        for (int i = 0; i < maximalIterationCount; ++i) {
+
+            // improve integral with a larger number of steps
+            final double t = stage(f, min, max, n);
+
+            // estimate error
+            final double delta = Math.abs(t - oldt);
+            final double limit =
+                Math.max(absoluteAccuracy,
+                         relativeAccuracy * (Math.abs(oldt) + Math.abs(t)) * 0.5);
+
+            // check convergence
+            if ((i + 1 >= minimalIterationCount) && (delta <= limit)) {
+                setResult(t, i);
+                return result;
+            }
+
+            // prepare next iteration
+            double ratio = Math.min(4, Math.pow(delta / limit, 0.5 / abscissas.length));
+            n = Math.max((int) (ratio * n), n + 1);
+            oldt = t;
 
-  /** Get the number of functions evaluation per step.
-   * @return number of function evaluation per step
-   */
-  public int getEvaluationsPerStep() {
-    return weightedRoots.length;
-  }
-
-  public double integrate(ComputableFunction f, double a, double b)
-    throws FunctionException {
-
-    // swap the bounds if they are not in ascending order
-    if (b < a) {
-      double tmp = b;
-      b          = a;
-      a          = tmp;
-    }
+        }
 
-    // adjust the step according to the bounds
-    long   n     = Math.round(0.5 + (b - a) / rawStep);
-    double step  = (b - a) / n;
-
-    // integrate over all elementary steps
-    double halfStep = step / 2.0;
-    double midPoint = a + halfStep;
-    double sum = 0.0;
-    for (long i = 0; i < n; ++i) {
-      for (int j = 0; j < weightedRoots.length; ++j) {
-        sum += weightedRoots[j][0]
-          * f.valueAt(midPoint + halfStep * weightedRoots[j][1]);
-      }
-      midPoint += step;
-    }
+        throw new MaxIterationsExceededException(maximalIterationCount);
 
-    return halfStep * sum;
+    }
 
-  }
+    /**
+     * Compute the n-th stage integral.
+     * @param f the integrand function
+     * @param min the lower bound for the interval
+     * @param max the upper bound for the interval
+     * @param n number of steps
+     * @return the value of n-th stage integral
+     * @throws FunctionEvaluationException if an error occurs evaluating the
+     * function
+     */
+    private double stage(final UnivariateRealFunction f,
+                         final double min, final double max, final int n)
+        throws FunctionEvaluationException {
+
+        // set up the step for the current stage
+        final double step     = (max - min) / n;
+        final double halfStep = step / 2.0;
+
+        // integrate over all elementary steps
+        double midPoint = min + halfStep;
+        double sum = 0.0;
+        for (int i = 0; i < n; ++i) {
+            for (int j = 0; j < abscissas.length; ++j) {
+                sum += weights[j] * f.value(midPoint + halfStep * abscissas[j]);
+            }
+            midPoint += step;
+        }
 
-  double[][] weightedRoots;
+        return halfStep * sum;
 
-  double rawStep;
+    }
 
 }

Propchange: commons/proper/math/trunk/src/java/org/apache/commons/math/analysis/integration/LegendreGaussIntegrator.java
------------------------------------------------------------------------------
    svn:eol-style = native

Propchange: commons/proper/math/trunk/src/java/org/apache/commons/math/analysis/integration/LegendreGaussIntegrator.java
------------------------------------------------------------------------------
    svn:keywords = Author Date Id Revision

Propchange: commons/proper/math/trunk/src/java/org/apache/commons/math/analysis/integration/LegendreGaussIntegrator.java
------------------------------------------------------------------------------
    svn:mergeinfo = 

Modified: commons/proper/math/trunk/src/site/xdoc/changes.xml
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/site/xdoc/changes.xml?rev=735876&r1=735875&r2=735876&view=diff
==============================================================================
--- commons/proper/math/trunk/src/site/xdoc/changes.xml (original)
+++ commons/proper/math/trunk/src/site/xdoc/changes.xml Mon Jan 19 15:40:14 2009
@@ -39,6 +39,9 @@
   </properties>
   <body>
     <release version="2.0" date="TBD" description="TBD">
+      <action dev="luc" type="add" >
+        Added a Legendre-Gauss integrator.
+      </action>
       <action dev="psteitz" type="fix" issue="MATH-240" due-to="Christian Semrau">
         Fixed error in factorial computation for 17 <= n <= 20.
       </action>

Modified: commons/proper/math/trunk/src/site/xdoc/userguide/analysis.xml
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/site/xdoc/userguide/analysis.xml?rev=735876&r1=735875&r2=735876&view=diff
==============================================================================
--- commons/proper/math/trunk/src/site/xdoc/userguide/analysis.xml (original)
+++ commons/proper/math/trunk/src/site/xdoc/userguide/analysis.xml Mon Jan 19 15:40:14 2009
@@ -312,6 +312,8 @@
           Simpson's</a></li>
           <li><a href="../apidocs/org/apache/commons/math/analysis/integration/TrapezoidIntegrator.html">
           trapezoid method</a></li>
+          <li><a href="../apidocs/org/apache/commons/math/analysis/integration/LegendreGaussIntegrator.html">
+          Legendre-Gauss method</a></li>
           </ul>      
         </p>
       </subsection>

Copied: commons/proper/math/trunk/src/test/org/apache/commons/math/analysis/integration/LegendreGaussIntegratorTest.java
(from r735530, commons/proper/math/trunk/src/mantissa/tests-src/org/spaceroots/mantissa/quadrature/scalar/GaussLegendreIntegratorTest.java)
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/test/org/apache/commons/math/analysis/integration/LegendreGaussIntegratorTest.java?p2=commons/proper/math/trunk/src/test/org/apache/commons/math/analysis/integration/LegendreGaussIntegratorTest.java&p1=commons/proper/math/trunk/src/mantissa/tests-src/org/spaceroots/mantissa/quadrature/scalar/GaussLegendreIntegratorTest.java&r1=735530&r2=735876&rev=735876&view=diff
==============================================================================
--- commons/proper/math/trunk/src/mantissa/tests-src/org/spaceroots/mantissa/quadrature/scalar/GaussLegendreIntegratorTest.java
(original)
+++ commons/proper/math/trunk/src/test/org/apache/commons/math/analysis/integration/LegendreGaussIntegratorTest.java
Mon Jan 19 15:40:14 2009
@@ -1,102 +1,117 @@
-// Licensed to the Apache Software Foundation (ASF) under one
-// or more contributor license agreements.  See the NOTICE file
-// distributed with this work for additional information
-// regarding copyright ownership.  The ASF licenses this file
-// to you under the Apache License, Version 2.0 (the
-// "License"); you may not use this file except in compliance
-// with the License.  You may obtain a copy of the License at
-// 
-//   http://www.apache.org/licenses/LICENSE-2.0
-// 
-// Unless required by applicable law or agreed to in writing,
-// software distributed under the License is distributed on an
-// "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY
-// KIND, either express or implied.  See the License for the
-// specific language governing permissions and limitations
-// under the License.
-
-package org.spaceroots.mantissa.quadrature.scalar;
-
-import org.spaceroots.mantissa.functions.scalar.ComputableFunction;
-import org.spaceroots.mantissa.functions.FunctionException;
+/*
+ * Licensed to the Apache Software Foundation (ASF) under one or more
+ * contributor license agreements.  See the NOTICE file distributed with
+ * this work for additional information regarding copyright ownership.
+ * The ASF licenses this file to You under the Apache License, Version 2.0
+ * (the "License"); you may not use this file except in compliance with
+ * the License.  You may obtain a copy of the License at
+ *
+ *      http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+package org.apache.commons.math.analysis.integration;
 
 import java.util.Random;
 
-import junit.framework.*;
+import org.apache.commons.math.ConvergenceException;
+import org.apache.commons.math.FunctionEvaluationException;
+import org.apache.commons.math.MathException;
+import org.apache.commons.math.analysis.QuinticFunction;
+import org.apache.commons.math.analysis.SinFunction;
+import org.apache.commons.math.analysis.UnivariateRealFunction;
+import org.apache.commons.math.analysis.polynomials.PolynomialFunction;
 
-public class GaussLegendreIntegratorTest
-  extends TestCase {
+import junit.framework.*;
 
-  public GaussLegendreIntegratorTest(String name) {
-    super(name);
-  }
-
-  public void testExactIntegration()
-    throws FunctionException {
-    Random random = new Random(86343623467878363l);
-    int order = 0;
-    while (true) {
-      GaussLegendreIntegrator integrator = new GaussLegendreIntegrator(order,
-                                                                       7.0);
-      int availableOrder = integrator.getEvaluationsPerStep();
-      if (availableOrder < order) {
-        // we have tested all available orders
-        return;
-      }
-
-      // an order n Gauss-Legendre integrator integrates
-      // 2n-1 degree polynoms exactly
-      for (int degree = 0; degree <= 2 * availableOrder - 1; ++degree) {
-        for (int i = 0; i < 10; ++i) {
-          Polynom p = new Polynom(degree, random, 100.0);
-          double s0 = integrator.integrate(p, -5.0, 15.0);
-          double s1 = p.exactIntegration(-5.0, 15.0);
-          assertTrue(Math.abs(s0 - s1) < 1.0e-12 * (1.0 + Math.abs(s0)));
-        }
-      }
+public class LegendreGaussIntegratorTest
+extends TestCase {
 
-      ++order;
+    public LegendreGaussIntegratorTest(String name) {
+        super(name);
+    }
 
+    public void testSinFunction() throws MathException {
+        UnivariateRealFunction f = new SinFunction();
+        UnivariateRealIntegrator integrator = new LegendreGaussIntegrator(5, 64);
+        integrator.setAbsoluteAccuracy(1.0e-10);
+        integrator.setRelativeAccuracy(1.0e-14);
+        integrator.setMinimalIterationCount(2);
+        integrator.setMaximalIterationCount(15);
+        double min, max, expected, result, tolerance;
+
+        min = 0; max = Math.PI; expected = 2;
+        tolerance = Math.max(integrator.getAbsoluteAccuracy(),
+                             Math.abs(expected * integrator.getRelativeAccuracy()));
+        result = integrator.integrate(f, min, max);
+        assertEquals(expected, result, tolerance);
+
+        min = -Math.PI/3; max = 0; expected = -0.5;
+        tolerance = Math.max(integrator.getAbsoluteAccuracy(),
+                Math.abs(expected * integrator.getRelativeAccuracy()));
+        result = integrator.integrate(f, min, max);
+        assertEquals(expected, result, tolerance);
     }
-  }
 
-  public static Test suite() {
-    return new TestSuite(GaussLegendreIntegratorTest.class);
-  }
+    public void testQuinticFunction() throws MathException {
+        UnivariateRealFunction f = new QuinticFunction();
+        UnivariateRealIntegrator integrator = new LegendreGaussIntegrator(3, 64);
+        double min, max, expected, result;
+
+        min = 0; max = 1; expected = -1.0/48;
+        result = integrator.integrate(f, min, max);
+        assertEquals(expected, result, 1.0e-16);
+
+        min = 0; max = 0.5; expected = 11.0/768;
+        result = integrator.integrate(f, min, max);
+        assertEquals(expected, result, 1.0e-16);
+
+        min = -1; max = 4; expected = 2048/3.0 - 78 + 1.0/48;
+        result = integrator.integrate(f, min, max);
+        assertEquals(expected, result, 1.0e-16);
+    }
 
-  private static class Polynom implements ComputableFunction {
+    public void testExactIntegration()
+        throws ConvergenceException, FunctionEvaluationException {
+        Random random = new Random(86343623467878363l);
+        for (int n = 2; n < 6; ++n) {
+            LegendreGaussIntegrator integrator =
+                new LegendreGaussIntegrator(n, 64);
+
+            // an n points Gauss-Legendre integrator integrates 2n-1 degree polynoms exactly
+            for (int degree = 0; degree <= 2 * n - 1; ++degree) {
+                for (int i = 0; i < 10; ++i) {
+                    double[] coeff = new double[degree + 1];
+                    for (int k = 0; k < coeff.length; ++k) {
+                        coeff[k] = 2 * random.nextDouble() - 1;
+                    }
+                    PolynomialFunction p = new PolynomialFunction(coeff);
+                    double result    = integrator.integrate(p, -5.0, 15.0);
+                    double reference = exactIntegration(p, -5.0, 15.0);
+                    assertEquals(n + " " + degree + " " + i, reference, result, 1.0e-12 *
(1.0 + Math.abs(reference)));
+                }
+            }
 
-    public Polynom(int degree, Random random, double max) {
-      coeffs = new double[degree + 1];
-      for (int i = 0; i <= degree; ++i) {
-        coeffs[i] = 2.0 * max * (random.nextDouble() - 0.5);
-      }
+        }
     }
 
-    public double valueAt(double t)
-      throws FunctionException {
-      double y = coeffs[coeffs.length - 1];
-      for (int i = coeffs.length - 2; i >= 0; --i) {
-        y = y * t + coeffs[i];
-      }
-      return y;
+    private double exactIntegration(PolynomialFunction p, double a, double b) {
+        final double[] coeffs = p.getCoefficients();
+        double yb = coeffs[coeffs.length - 1] / coeffs.length;
+        double ya = yb;
+        for (int i = coeffs.length - 2; i >= 0; --i) {
+            yb = yb * b + coeffs[i] / (i + 1);
+            ya = ya * a + coeffs[i] / (i + 1);
+        }
+        return yb * b - ya * a;
     }
 
-    public double exactIntegration(double a, double b)
-      throws FunctionException {
-      double yb = coeffs[coeffs.length - 1] / coeffs.length;
-      double ya = yb;
-      for (int i = coeffs.length - 2; i >= 0; --i) {
-        yb = yb * b + coeffs[i] / (i + 1);
-        ya = ya * a + coeffs[i] / (i + 1);
-      }
-      return yb * b - ya * a;
+    public static Test suite() {
+        return new TestSuite(LegendreGaussIntegratorTest.class);
     }
 
-    private double[] coeffs;
-
-     private static final long serialVersionUID = -7304282612679254557L;
-
-  }
-
 }

Propchange: commons/proper/math/trunk/src/test/org/apache/commons/math/analysis/integration/LegendreGaussIntegratorTest.java
------------------------------------------------------------------------------
    svn:eol-style = native

Propchange: commons/proper/math/trunk/src/test/org/apache/commons/math/analysis/integration/LegendreGaussIntegratorTest.java
------------------------------------------------------------------------------
    svn:keywords = Author Date Id Revision

Propchange: commons/proper/math/trunk/src/test/org/apache/commons/math/analysis/integration/LegendreGaussIntegratorTest.java
------------------------------------------------------------------------------
    svn:mergeinfo = 



Mime
View raw message