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 ∫ Li<sup>2</sup> where Li (x) =
+ * ∏ (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 =
