You are viewing a plain text version of this content. The canonical link for it is here.
Posted to commits@commons.apache.org by er...@apache.org on 2012/07/12 16:43:19 UTC
svn commit: r1360706 - in /commons/proper/math/trunk/src:
main/java/org/apache/commons/math3/analysis/integration/gauss/
test/java/org/apache/commons/math3/analysis/integration/gauss/
Author: erans
Date: Thu Jul 12 14:43:18 2012
New Revision: 1360706
URL: http://svn.apache.org/viewvc?rev=1360706&view=rev
Log:
MATH-797
Initial version for Gauss-Legendre quadrature rules: the integration is
performed on the whole interval using a single rule. [Whereas the approach
used in class "analysis.integration.LegendreGaussIntegrator" is to divide
iteratively into sub-intervals (over which the integration rule is used)
until some covergence criterion is met.]
Adapted from an original code donated by Sébastien Brisard.
In the current implementation, the Gauss-Legendre rules are computed in
double precision in "LegendreRuleFactory" and high precision (using
"java.math.BigDecimal") in "LegendreHighPrecisionRuleFactory". However,
the "GaussIntegrator" class performs the integration using "double"s
whatever the precision of the rule.
The framework of "BaseRuleFactory" enables the addition of other quadrature
schemes (by overriding the "computeRule" method). [Sébastien's code already
provides Gauss-Chebyshev and Gauss-Hermite schemes (in double precision).]
Added:
commons/proper/math/trunk/src/main/java/org/apache/commons/math3/analysis/integration/gauss/
commons/proper/math/trunk/src/main/java/org/apache/commons/math3/analysis/integration/gauss/BaseRuleFactory.java (with props)
commons/proper/math/trunk/src/main/java/org/apache/commons/math3/analysis/integration/gauss/GaussIntegrator.java (with props)
commons/proper/math/trunk/src/main/java/org/apache/commons/math3/analysis/integration/gauss/GaussIntegratorFactory.java (with props)
commons/proper/math/trunk/src/main/java/org/apache/commons/math3/analysis/integration/gauss/LegendreHighPrecisionRuleFactory.java (with props)
commons/proper/math/trunk/src/main/java/org/apache/commons/math3/analysis/integration/gauss/LegendreRuleFactory.java (with props)
commons/proper/math/trunk/src/test/java/org/apache/commons/math3/analysis/integration/gauss/
commons/proper/math/trunk/src/test/java/org/apache/commons/math3/analysis/integration/gauss/GaussianQuadratureAbstractTest.java (with props)
commons/proper/math/trunk/src/test/java/org/apache/commons/math3/analysis/integration/gauss/LegendreHighPrecisionParametricTest.java (with props)
commons/proper/math/trunk/src/test/java/org/apache/commons/math3/analysis/integration/gauss/LegendreHighPrecisionTest.java (with props)
commons/proper/math/trunk/src/test/java/org/apache/commons/math3/analysis/integration/gauss/LegendreParametricTest.java (with props)
commons/proper/math/trunk/src/test/java/org/apache/commons/math3/analysis/integration/gauss/LegendreTest.java (with props)
Added: commons/proper/math/trunk/src/main/java/org/apache/commons/math3/analysis/integration/gauss/BaseRuleFactory.java
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/main/java/org/apache/commons/math3/analysis/integration/gauss/BaseRuleFactory.java?rev=1360706&view=auto
==============================================================================
--- commons/proper/math/trunk/src/main/java/org/apache/commons/math3/analysis/integration/gauss/BaseRuleFactory.java (added)
+++ commons/proper/math/trunk/src/main/java/org/apache/commons/math3/analysis/integration/gauss/BaseRuleFactory.java Thu Jul 12 14:43:18 2012
@@ -0,0 +1,113 @@
+/*
+ * 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.math3.analysis.integration.gauss;
+
+import java.util.Map;
+import java.util.TreeMap;
+import org.apache.commons.math3.util.Pair;
+import org.apache.commons.math3.exception.DimensionMismatchException;
+
+/**
+ * Base class for rules that determines the integration nodes and their
+ * weights.
+ * Subclasses must implement the {@link #computeRule(int) computeRule} method.
+ *
+ * <T> Type of the number used to represent the points and weights of the
+ * quadrature rules.
+ *
+ * @version $Id$
+ * @since 3.1
+ */
+public abstract class BaseRuleFactory<T extends Number> {
+ /** List of points and weights, indexed by the order of the rule. */
+ private final Map<Integer, Pair<T[], T[]>> pointsAndWeights
+ = new TreeMap<Integer, Pair<T[], T[]>>();
+
+ /**
+ * Gets a copy of the quadrature rule with given number of integration points.
+ *
+ * @param numberOfPoints Number of integration points.
+ * @return a copy of the integration rule.
+ */
+ public Pair<double[], double[]> getRule(int numberOfPoints) {
+ return convertToDouble(getRuleInternal(numberOfPoints));
+ }
+
+ /**
+ * Gets a rule.
+ * Rules are computed once, and cached.
+ * The returned rule is a reference into the cache.
+ *
+ * @param numberOfPoints Order of the rule to be retrieved.
+ * @return the points and weights corresponding to the given order.
+ */
+ protected Pair<T[], T[]> getRuleInternal(int numberOfPoints) {
+ final Pair<T[], T[]> rule = pointsAndWeights.get(numberOfPoints);
+ if (rule == null) {
+ addRule(computeRule(numberOfPoints));
+ // The rule should be available now.
+ return getRuleInternal(numberOfPoints);
+ }
+ return rule;
+ }
+
+ /**
+ * Stores a rule.
+ *
+ * @param rule Rule to be stored.
+ * @throws DimensionMismatchException if the elements of the pair do not
+ * have the same length.
+ */
+ protected void addRule(Pair<T[], T[]> rule) {
+ if (rule.getFirst().length != rule.getSecond().length) {
+ throw new DimensionMismatchException(rule.getFirst().length,
+ rule.getSecond().length);
+ }
+
+ pointsAndWeights.put(rule.getFirst().length, rule);
+ }
+
+ /**
+ * Computes the rule for the given order.
+ *
+ * @param numberOfPoints Order of the rule to be computed.
+ * @return the computed rule.
+ */
+ protected abstract Pair<T[], T[]> computeRule(int numberOfPoints);
+
+ /**
+ * Converts the from the actual {@code Number} type to {@code double}
+ *
+ * @param rule Points and weights.
+ * @return points and weights as {@code double}s.
+ */
+ private static <T extends Number> Pair<double[], double[]> convertToDouble(Pair<T[], T[]> rule) {
+ final T[] pT = rule.getFirst();
+ final T[] wT = rule.getSecond();
+
+ final int len = pT.length;
+ final double[] pD = new double[len];
+ final double[] wD = new double[len];
+
+ for (int i = 0; i < len; i++) {
+ pD[i] = pT[i].doubleValue();
+ wD[i] = wT[i].doubleValue();
+ }
+
+ return new Pair<double[], double[]>(pD, wD);
+ }
+}
Propchange: commons/proper/math/trunk/src/main/java/org/apache/commons/math3/analysis/integration/gauss/BaseRuleFactory.java
------------------------------------------------------------------------------
svn:eol-style = native
Added: commons/proper/math/trunk/src/main/java/org/apache/commons/math3/analysis/integration/gauss/GaussIntegrator.java
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/main/java/org/apache/commons/math3/analysis/integration/gauss/GaussIntegrator.java?rev=1360706&view=auto
==============================================================================
--- commons/proper/math/trunk/src/main/java/org/apache/commons/math3/analysis/integration/gauss/GaussIntegrator.java (added)
+++ commons/proper/math/trunk/src/main/java/org/apache/commons/math3/analysis/integration/gauss/GaussIntegrator.java Thu Jul 12 14:43:18 2012
@@ -0,0 +1,106 @@
+/*
+ * 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.math3.analysis.integration.gauss;
+
+import org.apache.commons.math3.analysis.UnivariateFunction;
+import org.apache.commons.math3.exception.DimensionMismatchException;
+import org.apache.commons.math3.util.MathArrays;
+import org.apache.commons.math3.util.Pair;
+
+/**
+ * Class that implements the Gaussian rule for
+ * {@link #integrate(UnivariateFunction) integrating} a weighted
+ * function.
+ *
+ * @version $Id$
+ * @since 3.1
+ */
+public class GaussIntegrator {
+ /** Nodes. */
+ private final double[] points;
+ /** Nodes weights. */
+ private final double[] weights;
+
+ /**
+ * Creates an integrator from the given {@code points} and {@code weights}.
+ * The integration interval is defined by the first and last value of
+ * {@code points} which must be sorted in increasing order.
+ *
+ * @param points Integration points.
+ * @param weights Weights of the corresponding integration nodes.
+ * @throws org.apache.commons.math3.exception.NonMonotonicSequenceException
+ * if the {@code points} are not sorted in increasing order.
+ */
+ public GaussIntegrator(double[] points,
+ double[] weights) {
+ if (points.length != weights.length) {
+ throw new DimensionMismatchException(points.length,
+ weights.length);
+ }
+
+ MathArrays.checkOrder(points, MathArrays.OrderDirection.INCREASING, true, true);
+
+ this.points = points.clone();
+ this.weights = weights.clone();
+ }
+
+ /**
+ * Creates an integrator from the given pair of points (first element of
+ * the pair) and weights (second element of the pair.
+ *
+ * @param pointsAndWeights Integration points and corresponding weights.
+ * @throws org.apache.commons.math3.exception.NonMonotonicSequenceException
+ * if the {@code points} are not sorted in increasing order.
+ *
+ * @see #GaussIntegrator(double[], double[])
+ */
+ public GaussIntegrator(Pair<double[], double[]> pointsAndWeights) {
+ this(pointsAndWeights.getFirst(), pointsAndWeights.getSecond());
+ }
+
+ /**
+ * Returns an estimate of the integral of {@code f(x) * w(x)},
+ * where {@code w} is a weight function that depends on the actual
+ * flavor of the Gauss integration scheme.
+ * The algorithm uses the points and associated weights, as passed
+ * to the {@link #GaussIntegrator(double[],double[]) constructor}.
+ *
+ * @param f Function to integrate.
+ * @return the integral of the weighted function.
+ */
+ public double integrate(UnivariateFunction f) {
+ double s = 0;
+ double c = 0;
+ for (int i = 0; i < points.length; i++) {
+ final double x = points[i];
+ final double w = weights[i];
+ final double y = w * f.value(x) - c;
+ final double t = s + y;
+ c = (t - s) - y;
+ s = t;
+ }
+ return s;
+ }
+
+ /**
+ * @return the order of the integration rule (the number of integration
+ * points).
+ */
+ public int getNumberOfPoints() {
+ return points.length;
+ }
+}
Propchange: commons/proper/math/trunk/src/main/java/org/apache/commons/math3/analysis/integration/gauss/GaussIntegrator.java
------------------------------------------------------------------------------
svn:eol-style = native
Added: commons/proper/math/trunk/src/main/java/org/apache/commons/math3/analysis/integration/gauss/GaussIntegratorFactory.java
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/main/java/org/apache/commons/math3/analysis/integration/gauss/GaussIntegratorFactory.java?rev=1360706&view=auto
==============================================================================
--- commons/proper/math/trunk/src/main/java/org/apache/commons/math3/analysis/integration/gauss/GaussIntegratorFactory.java (added)
+++ commons/proper/math/trunk/src/main/java/org/apache/commons/math3/analysis/integration/gauss/GaussIntegratorFactory.java Thu Jul 12 14:43:18 2012
@@ -0,0 +1,134 @@
+/*
+ * 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.math3.analysis.integration.gauss;
+
+import java.util.Map;
+import java.util.TreeMap;
+import org.apache.commons.math3.util.Pair;
+import org.apache.commons.math3.exception.DimensionMismatchException;
+
+/**
+ * Class that provides different ways to compute the nodes and weights to be
+ * used by the {@link GaussIntegrator Gaussian integration rule}.
+ *
+ * @version $Id$
+ * @since 3.1
+ */
+public class GaussIntegratorFactory {
+ /** Generator of Gauss-Legendre integrators. */
+ private final BaseRuleFactory legendre = new LegendreRuleFactory();
+ /** Generator of Gauss-Legendre integrators. */
+ private final BaseRuleFactory legendreHighPrecision = new LegendreHighPrecisionRuleFactory();
+
+ /**
+ * Creates an integrator of the given order, and whose call to the
+ * {@link GaussIntegrator#integrate(org.apache.commons.math3.analysis.UnivariateFunction)
+ * integrate} method will perform an integration on the natural interval
+ * {@code [-1 , 1]}.
+ *
+ * @param numberOfPoints Order of the integration rule.
+ * @return a Gauss-Legendre integrator.
+ */
+ public GaussIntegrator legendre(int numberOfPoints) {
+ return new GaussIntegrator(getRule(legendre, numberOfPoints));
+ }
+
+ /**
+ * Creates an integrator of the given order, and whose call to the
+ * {@link GaussIntegrator#integrate(org.apache.commons.math3.analysis.UnivariateFunction)
+ * integrate} method will perform an integration on the given interval.
+ *
+ * @param numberOfPoints Order of the integration rule.
+ * @param lowerBound Lower bound of the integration interval.
+ * @param upperBound Upper bound of the integration interval.
+ * @return a Gauss-Legendre integrator.
+ */
+ public GaussIntegrator legendre(int numberOfPoints,
+ double lowerBound,
+ double upperBound) {
+ return new GaussIntegrator(transform(getRule(legendre, numberOfPoints),
+ lowerBound, upperBound));
+ }
+
+ /**
+ * Creates an integrator of the given order, and whose call to the
+ * {@link GaussIntegrator#integrate(org.apache.commons.math3.analysis.UnivariateFunction)
+ * integrate} method will perform an integration on the natural interval
+ * {@code [-1 , 1]}.
+ *
+ * @param numberOfPoints Order of the integration rule.
+ * @return a Gauss-Legendre integrator.
+ */
+ public GaussIntegrator legendreHighPrecision(int numberOfPoints) {
+ return new GaussIntegrator(getRule(legendreHighPrecision, numberOfPoints));
+ }
+
+ /**
+ * Creates an integrator of the given order, and whose call to the
+ * {@link GaussIntegrator#integrate(org.apache.commons.math3.analysis.UnivariateFunction)
+ * integrate} method will perform an integration on the given interval.
+ *
+ * @param numberOfPoints Order of the integration rule.
+ * @param lowerBound Lower bound of the integration interval.
+ * @param upperBound Upper bound of the integration interval.
+ * @return a Gauss-Legendre integrator.
+ */
+ public GaussIntegrator legendreHighPrecision(int numberOfPoints,
+ double lowerBound,
+ double upperBound) {
+ return new GaussIntegrator(transform(getRule(legendreHighPrecision, numberOfPoints),
+ lowerBound, upperBound));
+ }
+
+ /**
+ * @param factory Integration rule factory.
+ * @param numberOfPoints Order of the integration rule.
+ * @return the integration nodes and weights.
+ */
+ private static Pair<double[], double[]> getRule(BaseRuleFactory factory,
+ int numberOfPoints) {
+ return factory.getRule(numberOfPoints);
+ }
+
+ /**
+ * Performs a change of variable so that the integration can be performed
+ * on an arbitrary interval {@code [a, b]}.
+ * It is assumed that the natural interval is {@code [-1, 1]}.
+ *
+ * @param rule Original points and weights.
+ * @param a Lower bound of the integration interval.
+ * @param b Lower bound of the integration interval.
+ * @return the points and weights adapted to the new interval.
+ */
+ private static Pair<double[], double[]> transform(Pair<double[], double[]> rule,
+ double a,
+ double b) {
+ final double[] points = rule.getFirst();
+ final double[] weights = rule.getSecond();
+
+ // Scaling
+ final double scale = (b - a) / 2;
+ final double shift = a + scale;
+
+ for (int i = 0; i < points.length; i++) {
+ points[i] = points[i] * scale + shift;
+ weights[i] *= scale;
+ }
+
+ return new Pair<double[], double[]>(points, weights);
+ }
+}
Propchange: commons/proper/math/trunk/src/main/java/org/apache/commons/math3/analysis/integration/gauss/GaussIntegratorFactory.java
------------------------------------------------------------------------------
svn:eol-style = native
Added: commons/proper/math/trunk/src/main/java/org/apache/commons/math3/analysis/integration/gauss/LegendreHighPrecisionRuleFactory.java
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/main/java/org/apache/commons/math3/analysis/integration/gauss/LegendreHighPrecisionRuleFactory.java?rev=1360706&view=auto
==============================================================================
--- commons/proper/math/trunk/src/main/java/org/apache/commons/math3/analysis/integration/gauss/LegendreHighPrecisionRuleFactory.java (added)
+++ commons/proper/math/trunk/src/main/java/org/apache/commons/math3/analysis/integration/gauss/LegendreHighPrecisionRuleFactory.java Thu Jul 12 14:43:18 2012
@@ -0,0 +1,210 @@
+/*
+ * 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.math3.analysis.integration.gauss;
+
+import java.math.MathContext;
+import java.math.BigDecimal;
+import org.apache.commons.math3.util.Pair;
+
+/**
+ * Factory that creates Gauss-type quadrature rule using Legendre polynomials.
+ * In this implementation, the lower and upper bounds of the natural interval
+ * of integration are -1 and 1, respectively.
+ * The Legendre polynomials are evaluated using the recurrence relation
+ * presented in <a href="http://en.wikipedia.org/wiki/Abramowitz_and_Stegun"
+ * Abramowitz and Stegun, 1964</a>.
+ *
+ * @version $Id$
+ * @since 3.1
+ */
+public class LegendreHighPrecisionRuleFactory extends BaseRuleFactory<BigDecimal> {
+ private final MathContext mContext;
+ /** The number {@code 2}. */
+ private final BigDecimal two;
+ /** The number {@code -1}. */
+ private final BigDecimal minusOne;
+ /** The number {@code 0.5}. */
+ private final BigDecimal oneHalf;
+
+ /**
+ * Default precision is {@link MathContext#DECIMAL128 DECIMAL128}.
+ */
+ public LegendreHighPrecisionRuleFactory() {
+ this(MathContext.DECIMAL128);
+ }
+
+ /**
+ * @param mContext Precision setting for computing the quadrature rules.
+ */
+ public LegendreHighPrecisionRuleFactory(MathContext mContext) {
+ this.mContext = mContext;
+ two = new BigDecimal("2", mContext);
+ minusOne = new BigDecimal("-1", mContext);
+ oneHalf = new BigDecimal("0.5", mContext);
+ }
+
+ /**
+ * {@inheritDoc}
+ */
+ @Override
+ protected Pair<BigDecimal[], BigDecimal[]> computeRule(int numberOfPoints) {
+ if (numberOfPoints == 1) {
+ // Break recursion.
+ return new Pair<BigDecimal[], BigDecimal[]>(new BigDecimal[] { BigDecimal.ZERO },
+ new BigDecimal[] { two });
+ }
+
+ // Get previous rule.
+ // If it has not been computed yet it will trigger a recursive call
+ // to this method.
+ final BigDecimal[] previousPoints = getRuleInternal(numberOfPoints - 1).getFirst();
+
+ // Compute next rule.
+ final BigDecimal[] points = new BigDecimal[numberOfPoints];
+ final BigDecimal[] weights = new BigDecimal[numberOfPoints];
+
+ // Find i-th root of P[n+1] by bracketing.
+ final int iMax = numberOfPoints / 2;
+ for (int i = 0; i < iMax; i++) {
+ // Lower-bound of the interval.
+ BigDecimal a = (i == 0) ? minusOne : previousPoints[i - 1];
+ // Upper-bound of the interval.
+ BigDecimal b = (iMax == 1) ? BigDecimal.ONE : previousPoints[i];
+ // P[j-1](a)
+ BigDecimal pma = BigDecimal.ONE;
+ // P[j](a)
+ BigDecimal pa = a;
+ // P[j-1](b)
+ BigDecimal pmb = BigDecimal.ONE;
+ // P[j](b)
+ BigDecimal pb = b;
+ for (int j = 1; j < numberOfPoints; j++) {
+ final BigDecimal b_two_j_p_1 = new BigDecimal(2 * j + 1, mContext);
+ final BigDecimal b_j = new BigDecimal(j, mContext);
+ final BigDecimal b_j_p_1 = new BigDecimal(j + 1, mContext);
+
+ // Compute P[j+1](a)
+ // ppa = ((2 * j + 1) * a * pa - j * pma) / (j + 1);
+
+ BigDecimal tmp1 = a.multiply(b_two_j_p_1, mContext);
+ tmp1 = pa.multiply(tmp1, mContext);
+ BigDecimal tmp2 = pma.multiply(b_j, mContext);
+ // P[j+1](a)
+ BigDecimal ppa = tmp1.subtract(tmp2, mContext);
+ ppa = ppa.divide(b_j_p_1, mContext);
+
+ // Compute P[j+1](b)
+ // ppb = ((2 * j + 1) * b * pb - j * pmb) / (j + 1);
+
+ tmp1 = b.multiply(b_two_j_p_1, mContext);
+ tmp1 = pb.multiply(tmp1, mContext);
+ tmp2 = pmb.multiply(b_j, mContext);
+ // P[j+1](b)
+ BigDecimal ppb = tmp1.subtract(tmp2, mContext);
+ ppb = ppb.divide(b_j_p_1, mContext);
+
+ pma = pa;
+ pa = ppa;
+ pmb = pb;
+ pb = ppb;
+ }
+ // Now pa = P[n+1](a), and pma = P[n](a). Same holds for b.
+ // Middle of the interval.
+ BigDecimal c = a.add(b, mContext).multiply(oneHalf, mContext);
+ // P[j-1](c)
+ BigDecimal pmc = BigDecimal.ONE;
+ // P[j](c)
+ BigDecimal pc = c;
+ boolean done = false;
+ while (!done) {
+ BigDecimal tmp1 = b.subtract(a, mContext);
+ BigDecimal tmp2 = c.ulp().multiply(BigDecimal.TEN, mContext);
+ done = tmp1.compareTo(tmp2) <= 0;
+ pmc = BigDecimal.ONE;
+ pc = c;
+ for (int j = 1; j < numberOfPoints; j++) {
+ final BigDecimal b_two_j_p_1 = new BigDecimal(2 * j + 1, mContext);
+ final BigDecimal b_j = new BigDecimal(j, mContext);
+ final BigDecimal b_j_p_1 = new BigDecimal(j + 1, mContext);
+
+ // Compute P[j+1](c)
+ tmp1 = c.multiply(b_two_j_p_1, mContext);
+ tmp1 = pc.multiply(tmp1, mContext);
+ tmp2 = pmc.multiply(b_j, mContext);
+ // P[j+1](c)
+ BigDecimal ppc = tmp1.subtract(tmp2, mContext);
+ ppc = ppc.divide(b_j_p_1, mContext);
+
+ pmc = pc;
+ pc = ppc;
+ }
+ // Now pc = P[n+1](c) and pmc = P[n](c).
+ if (!done) {
+ if (pa.signum() * pc.signum() <= 0) {
+ b = c;
+ pmb = pmc;
+ pb = pc;
+ } else {
+ a = c;
+ pma = pmc;
+ pa = pc;
+ }
+ c = a.add(b, mContext).multiply(oneHalf, mContext);
+ }
+ }
+ final BigDecimal nP = new BigDecimal(numberOfPoints, mContext);
+ BigDecimal tmp1 = pmc.subtract(c.multiply(pc, mContext), mContext);
+ tmp1 = tmp1.multiply(nP);
+ tmp1 = tmp1.pow(2, mContext);
+ BigDecimal tmp2 = c.pow(2, mContext);
+ tmp2 = BigDecimal.ONE.subtract(tmp2, mContext);
+ tmp2 = tmp2.multiply(two, mContext);
+ tmp2 = tmp2.divide(tmp1, mContext);
+
+ points[i] = c;
+ weights[i] = tmp2;
+
+ final int idx = numberOfPoints - i - 1;
+ points[idx] = c.negate(mContext);
+ weights[idx] = tmp2;
+ }
+ // If "numberOfPoints" is odd, 0 is a root.
+ if (numberOfPoints % 2 == 1) {
+ BigDecimal pmc = BigDecimal.ONE;
+ for (int j = 1; j < numberOfPoints; j += 2) {
+ final BigDecimal b_j = new BigDecimal(j, mContext);
+ final BigDecimal b_j_p_1 = new BigDecimal(j + 1, mContext);
+
+ // pmc = -j * pmc / (j + 1);
+ pmc = pmc.multiply(b_j, mContext);
+ pmc = pmc.divide(b_j_p_1, mContext);
+ pmc = pmc.negate(mContext);
+ }
+
+ // 2 / pow(numberOfPoints * pmc, 2);
+ final BigDecimal nP = new BigDecimal(numberOfPoints, mContext);
+ BigDecimal tmp1 = pmc.multiply(nP, mContext);
+ tmp1 = tmp1.pow(2, mContext);
+ BigDecimal tmp2 = two.divide(tmp1, mContext);
+
+ points[iMax] = BigDecimal.ZERO;
+ weights[iMax] = tmp2;
+ }
+
+ return new Pair<BigDecimal[], BigDecimal[]>(points, weights);
+ }
+}
Propchange: commons/proper/math/trunk/src/main/java/org/apache/commons/math3/analysis/integration/gauss/LegendreHighPrecisionRuleFactory.java
------------------------------------------------------------------------------
svn:eol-style = native
Added: commons/proper/math/trunk/src/main/java/org/apache/commons/math3/analysis/integration/gauss/LegendreRuleFactory.java
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/main/java/org/apache/commons/math3/analysis/integration/gauss/LegendreRuleFactory.java?rev=1360706&view=auto
==============================================================================
--- commons/proper/math/trunk/src/main/java/org/apache/commons/math3/analysis/integration/gauss/LegendreRuleFactory.java (added)
+++ commons/proper/math/trunk/src/main/java/org/apache/commons/math3/analysis/integration/gauss/LegendreRuleFactory.java Thu Jul 12 14:43:18 2012
@@ -0,0 +1,139 @@
+/*
+ * 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.math3.analysis.integration.gauss;
+
+import java.math.MathContext;
+import java.math.BigDecimal;
+import org.apache.commons.math3.util.Pair;
+
+/**
+ * Factory that creates Gauss-type quadrature rule using Legendre polynomials.
+ * In this implementation, the lower and upper bounds of the natural interval
+ * of integration are -1 and 1, respectively.
+ * The Legendre polynomials are evaluated using the recurrence relation
+ * presented in <a href="http://en.wikipedia.org/wiki/Abramowitz_and_Stegun"
+ * Abramowitz and Stegun, 1964</a>.
+ *
+ * @version $Id$
+ * @since 3.1
+ */
+public class LegendreRuleFactory extends BaseRuleFactory<Double> {
+ /**
+ * {@inheritDoc}
+ */
+ @Override
+ protected Pair<Double[], Double[]> computeRule(int numberOfPoints) {
+ if (numberOfPoints == 1) {
+ // Break recursion.
+ return new Pair<Double[], Double[]>(new Double[] { 0d },
+ new Double[] { 2d });
+ }
+
+ // Get previous rule.
+ // If it has not been computed yet it will trigger a recursive call
+ // to this method.
+ final Double[] previousPoints = getRuleInternal(numberOfPoints - 1).getFirst();
+
+ // Compute next rule.
+ final Double[] points = new Double[numberOfPoints];
+ final Double[] weights = new Double[numberOfPoints];
+
+ // Find i-th root of P[n+1] by bracketing.
+ final int iMax = numberOfPoints / 2;
+ for (int i = 0; i < iMax; i++) {
+ // Lower-bound of the interval.
+ double a = (i == 0) ? -1 : previousPoints[i - 1].doubleValue();
+ // Upper-bound of the interval.
+ double b = (iMax == 1) ? 1 : previousPoints[i].doubleValue();
+ // P[j-1](a)
+ double pma = 1;
+ // P[j](a)
+ double pa = a;
+ // P[j-1](b)
+ double pmb = 1;
+ // P[j](b)
+ double pb = b;
+ for (int j = 1; j < numberOfPoints; j++) {
+ final int two_j_p_1 = 2 * j + 1;
+ final int j_p_1 = j + 1;
+ // P[j+1](a)
+ final double ppa = (two_j_p_1 * a * pa - j * pma) / j_p_1;
+ // P[j+1](b)
+ final double ppb = (two_j_p_1 * b * pb - j * pmb) / j_p_1;
+ pma = pa;
+ pa = ppa;
+ pmb = pb;
+ pb = ppb;
+ }
+ // Now pa = P[n+1](a), and pma = P[n](a) (same holds for b).
+ // Middle of the interval.
+ double c = 0.5 * (a + b);
+ // P[j-1](c)
+ double pmc = 1;
+ // P[j](c)
+ double pc = c;
+ boolean done = false;
+ while (!done) {
+ done = b - a <= Math.ulp(c);
+ pmc = 1;
+ pc = c;
+ for (int j = 1; j < numberOfPoints; j++) {
+ // P[j+1](c)
+ final double ppc = ((2 * j + 1) * c * pc - j * pmc) / (j + 1);
+ pmc = pc;
+ pc = ppc;
+ }
+ // Now pc = P[n+1](c) and pmc = P[n](c).
+ if (!done) {
+ if (pa * pc <= 0) {
+ b = c;
+ pmb = pmc;
+ pb = pc;
+ } else {
+ a = c;
+ pma = pmc;
+ pa = pc;
+ }
+ c = 0.5 * (a + b);
+ }
+ }
+ final double d = numberOfPoints * (pmc - c * pc);
+ final double w = 2 * (1 - c * c) / (d * d);
+
+ points[i] = c;
+ weights[i] = w;
+
+ final int idx = numberOfPoints - i - 1;
+ points[idx] = -c;
+ weights[idx] = w;
+ }
+ // If "numberOfPoints" is odd, 0 is a root.
+ if (numberOfPoints % 2 == 1) {
+ double pmc = 1;
+ for (int j = 1; j < numberOfPoints; j += 2) {
+ pmc = -j * pmc / (j + 1);
+ }
+ final double d = numberOfPoints * pmc;
+ final double w = 2 / (d * d);
+
+ points[iMax] = 0d;
+ weights[iMax] = w;
+ }
+
+ return new Pair<Double[], Double[]>(points, weights);
+ }
+}
Propchange: commons/proper/math/trunk/src/main/java/org/apache/commons/math3/analysis/integration/gauss/LegendreRuleFactory.java
------------------------------------------------------------------------------
svn:eol-style = native
Added: commons/proper/math/trunk/src/test/java/org/apache/commons/math3/analysis/integration/gauss/GaussianQuadratureAbstractTest.java
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/test/java/org/apache/commons/math3/analysis/integration/gauss/GaussianQuadratureAbstractTest.java?rev=1360706&view=auto
==============================================================================
--- commons/proper/math/trunk/src/test/java/org/apache/commons/math3/analysis/integration/gauss/GaussianQuadratureAbstractTest.java (added)
+++ commons/proper/math/trunk/src/test/java/org/apache/commons/math3/analysis/integration/gauss/GaussianQuadratureAbstractTest.java Thu Jul 12 14:43:18 2012
@@ -0,0 +1,100 @@
+package org.apache.commons.math3.analysis.integration.gauss;
+
+import org.apache.commons.math3.analysis.function.Power;
+import org.junit.Test;
+import org.junit.Assert;
+
+/**
+ * Base class for standard testing of Gaussian quadrature rules,
+ * which are exact for polynomials up to a certain degree. In this test, each
+ * monomial in turn is tested against the specified quadrature rule.
+ *
+ * @version $Id$
+ */
+public abstract class GaussianQuadratureAbstractTest {
+ /**
+ * The maximum absolute error (for zero testing).
+ */
+ private final double eps;
+ /**
+ * The maximum relative error (in ulps).
+ */
+ private final double numUlps;
+ /**
+ * The quadrature rule under test.
+ */
+ private final GaussIntegrator integrator;
+ /**
+ * Maximum degree of monomials to be tested.
+ */
+ private final int maxDegree;
+
+ /**
+ * Creates a new instance of this abstract test with the specified
+ * quadrature rule.
+ * If the expected value is non-zero, equality of actual and expected values
+ * is checked in the relative sense <center>
+ * |x<sub>act</sub> - x<sub>exp</sub>| ≤ n
+ * <code>Math.ulp(</code>x<sub>exp</sub><code>)</code>, </center> where n is
+ * the maximum relative error (in ulps). If the expected value is zero, the
+ * test checks that <center> |x<sub>act</sub>| ≤ ε,
+ * </center> where ε is the maximum absolute error.
+ *
+ * @param integrator Quadrature rule under test.
+ * @param maxDegree Maximum degree of monomials to be tested.
+ * @param eps ε.
+ * @param numUlps Value of the maximum relative error (in ulps).
+ */
+ public GaussianQuadratureAbstractTest(GaussIntegrator integrator,
+ int maxDegree,
+ double eps,
+ double numUlps) {
+ this.integrator = integrator;
+ this.maxDegree = maxDegree;
+ this.eps = eps;
+ this.numUlps = numUlps;
+ }
+
+ /**
+ * Returns the expected value of the integral of the specified monomial.
+ * The integration is carried out on the natural interval of the quadrature
+ * rule under test.
+ *
+ * @param n Degree of the monomial.
+ * @return the expected value of the integral of x<sup>n</sup>.
+ */
+ public abstract double getExpectedValue(final int n);
+
+ /**
+ * Checks that the value of the integral of each monomial
+ * <code>x<sup>0</sup>, ... , x<sup>p</sup></code>
+ * returned by the quadrature rule under test conforms with the expected
+ * value.
+ * Here {@code p} denotes the degree of the highest polynomial for which
+ * exactness is to be expected.
+ */
+ @Test
+ public void testAllMonomials() {
+ for (int n = 0; n <= maxDegree; n++) {
+ final double expected = getExpectedValue(n);
+
+ final Power monomial = new Power(n);
+ final double actual = integrator.integrate(monomial);
+
+ // System.out.println(n + "/" + maxDegree + " " + integrator.getNumberOfPoints()
+ // + " " + expected + " " + actual + " " + Math.ulp(expected));
+ if (expected == 0) {
+ Assert.assertEquals("while integrating monomial x**" + n +
+ " with a " +
+ integrator.getNumberOfPoints() + "-point quadrature rule",
+ expected, actual, eps);
+ } else {
+ double err = Math.abs(actual - expected) / Math.ulp(expected);
+ Assert.assertEquals("while integrating monomial x**" + n + " with a " +
+ + integrator.getNumberOfPoints() + "-point quadrature rule, " +
+ " error was " + err + " ulps",
+ expected, actual, Math.ulp(expected) * numUlps);
+ }
+ }
+ }
+}
Propchange: commons/proper/math/trunk/src/test/java/org/apache/commons/math3/analysis/integration/gauss/GaussianQuadratureAbstractTest.java
------------------------------------------------------------------------------
svn:eol-style = native
Added: commons/proper/math/trunk/src/test/java/org/apache/commons/math3/analysis/integration/gauss/LegendreHighPrecisionParametricTest.java
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/test/java/org/apache/commons/math3/analysis/integration/gauss/LegendreHighPrecisionParametricTest.java?rev=1360706&view=auto
==============================================================================
--- commons/proper/math/trunk/src/test/java/org/apache/commons/math3/analysis/integration/gauss/LegendreHighPrecisionParametricTest.java (added)
+++ commons/proper/math/trunk/src/test/java/org/apache/commons/math3/analysis/integration/gauss/LegendreHighPrecisionParametricTest.java Thu Jul 12 14:43:18 2012
@@ -0,0 +1,69 @@
+package org.apache.commons.math3.analysis.integration.gauss;
+
+import java.util.ArrayList;
+import java.util.Collection;
+
+import org.junit.runner.RunWith;
+import org.junit.runners.Parameterized;
+import org.junit.runners.Parameterized.Parameters;
+
+/**
+ * Test of the {@link LegendreHighPrecisionRuleFactory}.
+ * This parameterized test extends the standard test for Gaussian quadrature
+ * rule, where each monomial is tested in turn.
+ * Parametrization allows to test automatically 0, 1, ... , {@link #MAX_NUM_POINTS}
+ * quadrature rules.
+ *
+ * @version $Id$
+ */
+@RunWith(value=Parameterized.class)
+public class LegendreHighPrecisionParametricTest extends GaussianQuadratureAbstractTest {
+ private static GaussIntegratorFactory factory = new GaussIntegratorFactory();
+
+ /**
+ * The highest order quadrature rule to be tested.
+ */
+ public static final int MAX_NUM_POINTS = 30;
+
+ /**
+ * Creates a new instance of this test, with the specified number of nodes
+ * for the Gauss-Legendre quadrature rule.
+ *
+ * @param numberOfPoints Order of integration rule.
+ * @param maxDegree Maximum degree of monomials to be tested.
+ * @param eps Value of ε.
+ * @param numUlps Value of the maximum relative error (in ulps).
+ */
+ public LegendreHighPrecisionParametricTest(int numberOfPoints,
+ int maxDegree,
+ double eps,
+ double numUlps) {
+ super(factory.legendreHighPrecision(numberOfPoints),
+ maxDegree, eps, numUlps);
+ }
+
+ /**
+ * Returns the collection of parameters to be passed to the constructor of
+ * this class.
+ * Gauss-Legendre quadrature rules of order 1, ..., {@link #MAX_NUM_POINTS}
+ * will be constructed.
+ *
+ * @return the collection of parameters for this parameterized test.
+ */
+ @Parameters
+ public static Collection<Object[]> getParameters() {
+ final ArrayList<Object[]> parameters = new ArrayList<Object[]>();
+ for (int k = 1; k <= MAX_NUM_POINTS; k++) {
+ parameters.add(new Object[] { k, 2 * k - 1, Math.ulp(1d), 13d });
+ }
+ return parameters;
+ }
+
+ @Override
+ public double getExpectedValue(final int n) {
+ if (n % 2 == 1) {
+ return 0;
+ }
+ return 2d / (n + 1);
+ }
+}
Propchange: commons/proper/math/trunk/src/test/java/org/apache/commons/math3/analysis/integration/gauss/LegendreHighPrecisionParametricTest.java
------------------------------------------------------------------------------
svn:eol-style = native
Added: commons/proper/math/trunk/src/test/java/org/apache/commons/math3/analysis/integration/gauss/LegendreHighPrecisionTest.java
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/test/java/org/apache/commons/math3/analysis/integration/gauss/LegendreHighPrecisionTest.java?rev=1360706&view=auto
==============================================================================
--- commons/proper/math/trunk/src/test/java/org/apache/commons/math3/analysis/integration/gauss/LegendreHighPrecisionTest.java (added)
+++ commons/proper/math/trunk/src/test/java/org/apache/commons/math3/analysis/integration/gauss/LegendreHighPrecisionTest.java Thu Jul 12 14:43:18 2012
@@ -0,0 +1,43 @@
+package org.apache.commons.math3.analysis.integration.gauss;
+
+import org.apache.commons.math3.analysis.function.Cos;
+import org.apache.commons.math3.analysis.function.Inverse;
+import org.apache.commons.math3.analysis.function.Log;
+import org.apache.commons.math3.analysis.UnivariateFunction;
+import org.junit.Test;
+import org.junit.Assert;
+
+/**
+ * Test of the {@link LegendreHighPrecisionRuleFactory}.
+ *
+ * @version $Id$
+ */
+public class LegendreHighPrecisionTest {
+ private static GaussIntegratorFactory factory = new GaussIntegratorFactory();
+
+ @Test
+ public void testCos() {
+ final UnivariateFunction cos = new Cos();
+
+ final GaussIntegrator integrator = factory.legendreHighPrecision(7, 0, Math.PI / 2);
+ final double s = integrator.integrate(cos);
+ // System.out.println("s=" + s + " e=" + 1);
+ Assert.assertEquals(1, s, Math.ulp(1d));
+ }
+
+
+ @Test
+ public void testInverse() {
+ final UnivariateFunction inv = new Inverse();
+ final UnivariateFunction log = new Log();
+
+ final double lo = 12.34;
+ final double hi = 456.78;
+
+ final GaussIntegrator integrator = factory.legendreHighPrecision(60, lo, hi);
+ final double s = integrator.integrate(inv);
+ final double expected = log.value(hi) - log.value(lo);
+ // System.out.println("s=" + s + " e=" + expected);
+ Assert.assertEquals(expected, s, 1e-15);
+ }
+}
Propchange: commons/proper/math/trunk/src/test/java/org/apache/commons/math3/analysis/integration/gauss/LegendreHighPrecisionTest.java
------------------------------------------------------------------------------
svn:eol-style = native
Added: commons/proper/math/trunk/src/test/java/org/apache/commons/math3/analysis/integration/gauss/LegendreParametricTest.java
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/test/java/org/apache/commons/math3/analysis/integration/gauss/LegendreParametricTest.java?rev=1360706&view=auto
==============================================================================
--- commons/proper/math/trunk/src/test/java/org/apache/commons/math3/analysis/integration/gauss/LegendreParametricTest.java (added)
+++ commons/proper/math/trunk/src/test/java/org/apache/commons/math3/analysis/integration/gauss/LegendreParametricTest.java Thu Jul 12 14:43:18 2012
@@ -0,0 +1,69 @@
+package org.apache.commons.math3.analysis.integration.gauss;
+
+import java.util.ArrayList;
+import java.util.Collection;
+
+import org.junit.runner.RunWith;
+import org.junit.runners.Parameterized;
+import org.junit.runners.Parameterized.Parameters;
+
+/**
+ * Test of the {@link LegendreRuleFactory}.
+ * This parameterized test extends the standard test for Gaussian quadrature
+ * rule, where each monomial is tested in turn.
+ * Parametrization allows to test automatically 0, 1, ... , {@link #MAX_NUM_POINTS}
+ * quadrature rules.
+ *
+ * @version $Id$
+ */
+@RunWith(value=Parameterized.class)
+public class LegendreParametricTest extends GaussianQuadratureAbstractTest {
+ private static GaussIntegratorFactory factory = new GaussIntegratorFactory();
+
+ /**
+ * The highest order quadrature rule to be tested.
+ */
+ public static final int MAX_NUM_POINTS = 30;
+
+ /**
+ * Creates a new instance of this test, with the specified number of nodes
+ * for the Gauss-Legendre quadrature rule.
+ *
+ * @param numberOfPoints Order of integration rule.
+ * @param maxDegree Maximum degree of monomials to be tested.
+ * @param eps Value of ε.
+ * @param numUlps Value of the maximum relative error (in ulps).
+ */
+ public LegendreParametricTest(int numberOfPoints,
+ int maxDegree,
+ double eps,
+ double numUlps) {
+ super(factory.legendre(numberOfPoints),
+ maxDegree, eps, numUlps);
+ }
+
+ /**
+ * Returns the collection of parameters to be passed to the constructor of
+ * this class.
+ * Gauss-Legendre quadrature rules of order 1, ..., {@link #MAX_NUM_POINTS}
+ * will be constructed.
+ *
+ * @return the collection of parameters for this parameterized test.
+ */
+ @Parameters
+ public static Collection<Object[]> getParameters() {
+ final ArrayList<Object[]> parameters = new ArrayList<Object[]>();
+ for (int k = 1; k <= MAX_NUM_POINTS; k++) {
+ parameters.add(new Object[] { k, 2 * k - 1, Math.ulp(1d), 91d });
+ }
+ return parameters;
+ }
+
+ @Override
+ public double getExpectedValue(final int n) {
+ if (n % 2 == 1) {
+ return 0;
+ }
+ return 2d / (n + 1);
+ }
+}
Propchange: commons/proper/math/trunk/src/test/java/org/apache/commons/math3/analysis/integration/gauss/LegendreParametricTest.java
------------------------------------------------------------------------------
svn:eol-style = native
Added: commons/proper/math/trunk/src/test/java/org/apache/commons/math3/analysis/integration/gauss/LegendreTest.java
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/test/java/org/apache/commons/math3/analysis/integration/gauss/LegendreTest.java?rev=1360706&view=auto
==============================================================================
--- commons/proper/math/trunk/src/test/java/org/apache/commons/math3/analysis/integration/gauss/LegendreTest.java (added)
+++ commons/proper/math/trunk/src/test/java/org/apache/commons/math3/analysis/integration/gauss/LegendreTest.java Thu Jul 12 14:43:18 2012
@@ -0,0 +1,43 @@
+package org.apache.commons.math3.analysis.integration.gauss;
+
+import org.apache.commons.math3.analysis.function.Cos;
+import org.apache.commons.math3.analysis.function.Inverse;
+import org.apache.commons.math3.analysis.function.Log;
+import org.apache.commons.math3.analysis.UnivariateFunction;
+import org.junit.Test;
+import org.junit.Assert;
+
+/**
+ * Test of the {@link LegendreRuleFactory}.
+ *
+ * @version $Id$
+ */
+public class LegendreTest {
+ private static GaussIntegratorFactory factory = new GaussIntegratorFactory();
+
+ @Test
+ public void testCos() {
+ final UnivariateFunction cos = new Cos();
+
+ final GaussIntegrator integrator = factory.legendre(7, 0, Math.PI / 2);
+ final double s = integrator.integrate(cos);
+ // System.out.println("s=" + s + " e=" + 1);
+ Assert.assertEquals(1, s, Math.ulp(1d));
+ }
+
+
+ @Test
+ public void testInverse() {
+ final UnivariateFunction inv = new Inverse();
+ final UnivariateFunction log = new Log();
+
+ final double lo = 12.34;
+ final double hi = 456.78;
+
+ final GaussIntegrator integrator = factory.legendre(60, lo, hi);
+ final double s = integrator.integrate(inv);
+ final double expected = log.value(hi) - log.value(lo);
+ // System.out.println("s=" + s + " e=" + expected);
+ Assert.assertEquals(expected, s, 1e-14);
+ }
+}
Propchange: commons/proper/math/trunk/src/test/java/org/apache/commons/math3/analysis/integration/gauss/LegendreTest.java
------------------------------------------------------------------------------
svn:eol-style = native