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.math4.legacy.analysis.integration.gauss; 18 19 import org.apache.commons.math4.legacy.analysis.UnivariateFunction; 20 import org.apache.commons.math4.legacy.exception.DimensionMismatchException; 21 import org.apache.commons.math4.legacy.core.MathArrays; 22 import org.apache.commons.math4.legacy.core.Pair; 23 24 /** 25 * Class that implements the Gaussian rule for 26 * {@link #integrate(UnivariateFunction) integrating} a weighted 27 * function. 28 * 29 * @since 3.1 30 */ 31 public class GaussIntegrator { 32 /** Nodes. */ 33 private final double[] points; 34 /** Nodes weights. */ 35 private final double[] weights; 36 37 /** 38 * Creates an integrator from the given {@code points} and {@code weights}. 39 * The integration interval is defined by the first and last value of 40 * {@code points} which must be sorted in increasing order. 41 * 42 * @param points Integration points. 43 * @param weights Weights of the corresponding integration nodes. 44 * @throws org.apache.commons.math4.legacy.exception.NonMonotonicSequenceException if the {@code points} are not 45 * sorted in increasing order. 46 * @throws DimensionMismatchException if points and weights don't have the same length 47 */ 48 public GaussIntegrator(double[] points, 49 double[] weights) { 50 if (points.length != weights.length) { 51 throw new DimensionMismatchException(points.length, 52 weights.length); 53 } 54 55 MathArrays.checkOrder(points, MathArrays.OrderDirection.INCREASING, true, true); 56 57 this.points = points.clone(); 58 this.weights = weights.clone(); 59 } 60 61 /** 62 * Creates an integrator from the given pair of points (first element of 63 * the pair) and weights (second element of the pair. 64 * 65 * @param pointsAndWeights Integration points and corresponding weights. 66 * @throws org.apache.commons.math4.legacy.exception.NonMonotonicSequenceException if the {@code points} are not 67 * sorted in increasing order. 68 * 69 * @see #GaussIntegrator(double[], double[]) 70 */ 71 public GaussIntegrator(Pair<double[], double[]> pointsAndWeights) { 72 this(pointsAndWeights.getFirst(), pointsAndWeights.getSecond()); 73 } 74 75 /** 76 * Returns an estimate of the integral of {@code f(x) * w(x)}, 77 * where {@code w} is a weight function that depends on the actual 78 * flavor of the Gauss integration scheme. 79 * The algorithm uses the points and associated weights, as passed 80 * to the {@link #GaussIntegrator(double[],double[]) constructor}. 81 * 82 * @param f Function to integrate. 83 * @return the integral of the weighted function. 84 */ 85 public double integrate(UnivariateFunction f) { 86 double s = 0; 87 double c = 0; 88 for (int i = 0; i < points.length; i++) { 89 final double x = points[i]; 90 final double w = weights[i]; 91 final double y = w * f.value(x) - c; 92 final double t = s + y; 93 c = (t - s) - y; 94 s = t; 95 } 96 return s; 97 } 98 99 /** 100 * @return the order of the integration rule (the number of integration 101 * points). 102 */ 103 public int getNumberOfPoints() { 104 return points.length; 105 } 106 107 /** 108 * Gets the integration point at the given index. 109 * The index must be in the valid range but no check is performed. 110 * @param index index of the integration point 111 * @return the integration point. 112 */ 113 public double getPoint(int index) { 114 return points[index]; 115 } 116 117 /** 118 * Gets the weight of the integration point at the given index. 119 * The index must be in the valid range but no check is performed. 120 * @param index index of the integration point 121 * @return the weight. 122 */ 123 public double getWeight(int index) { 124 return weights[index]; 125 } 126 }