SymmetricGaussIntegrator.java

  1. /*
  2.  * Licensed to the Apache Software Foundation (ASF) under one or more
  3.  * contributor license agreements.  See the NOTICE file distributed with
  4.  * this work for additional information regarding copyright ownership.
  5.  * The ASF licenses this file to You under the Apache License, Version 2.0
  6.  * (the "License"); you may not use this file except in compliance with
  7.  * the License.  You may obtain a copy of the License at
  8.  *
  9.  *      http://www.apache.org/licenses/LICENSE-2.0
  10.  *
  11.  * Unless required by applicable law or agreed to in writing, software
  12.  * distributed under the License is distributed on an "AS IS" BASIS,
  13.  * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
  14.  * See the License for the specific language governing permissions and
  15.  * limitations under the License.
  16.  */
  17. package org.apache.commons.math3.analysis.integration.gauss;

  18. import org.apache.commons.math3.analysis.UnivariateFunction;
  19. import org.apache.commons.math3.exception.DimensionMismatchException;
  20. import org.apache.commons.math3.exception.NonMonotonicSequenceException;
  21. import org.apache.commons.math3.util.Pair;

  22. /**
  23.  * This class's implements {@link #integrate(UnivariateFunction) integrate}
  24.  * method assuming that the integral is symmetric about 0.
  25.  * This allows to reduce numerical errors.
  26.  *
  27.  * @since 3.3
  28.  */
  29. public class SymmetricGaussIntegrator extends GaussIntegrator {
  30.     /**
  31.      * Creates an integrator from the given {@code points} and {@code weights}.
  32.      * The integration interval is defined by the first and last value of
  33.      * {@code points} which must be sorted in increasing order.
  34.      *
  35.      * @param points Integration points.
  36.      * @param weights Weights of the corresponding integration nodes.
  37.      * @throws NonMonotonicSequenceException if the {@code points} are not
  38.      * sorted in increasing order.
  39.      * @throws DimensionMismatchException if points and weights don't have the same length
  40.      */
  41.     public SymmetricGaussIntegrator(double[] points,
  42.                                     double[] weights)
  43.         throws NonMonotonicSequenceException, DimensionMismatchException {
  44.         super(points, weights);
  45.     }

  46.     /**
  47.      * Creates an integrator from the given pair of points (first element of
  48.      * the pair) and weights (second element of the pair.
  49.      *
  50.      * @param pointsAndWeights Integration points and corresponding weights.
  51.      * @throws NonMonotonicSequenceException if the {@code points} are not
  52.      * sorted in increasing order.
  53.      *
  54.      * @see #SymmetricGaussIntegrator(double[], double[])
  55.      */
  56.     public SymmetricGaussIntegrator(Pair<double[], double[]> pointsAndWeights)
  57.         throws NonMonotonicSequenceException {
  58.         this(pointsAndWeights.getFirst(), pointsAndWeights.getSecond());
  59.     }

  60.     /**
  61.      * {@inheritDoc}
  62.      */
  63.     @Override
  64.     public double integrate(UnivariateFunction f) {
  65.         final int ruleLength = getNumberOfPoints();

  66.         if (ruleLength == 1) {
  67.             return getWeight(0) * f.value(0d);
  68.         }

  69.         final int iMax = ruleLength / 2;
  70.         double s = 0;
  71.         double c = 0;
  72.         for (int i = 0; i < iMax; i++) {
  73.             final double p = getPoint(i);
  74.             final double w = getWeight(i);

  75.             final double f1 = f.value(p);
  76.             final double f2 = f.value(-p);

  77.             final double y = w * (f1 + f2) - c;
  78.             final double t = s + y;

  79.             c = (t - s) - y;
  80.             s = t;
  81.         }

  82.         if (ruleLength % 2 != 0) {
  83.             final double w = getWeight(iMax);

  84.             final double y = w * f.value(0d) - c;
  85.             final double t = s + y;

  86.             s = t;
  87.         }

  88.         return s;
  89.     }
  90. }