001/* 002 * Licensed to the Apache Software Foundation (ASF) under one or more 003 * contributor license agreements. See the NOTICE file distributed with 004 * this work for additional information regarding copyright ownership. 005 * The ASF licenses this file to You under the Apache License, Version 2.0 006 * (the "License"); you may not use this file except in compliance with 007 * the License. You may obtain a copy of the License at 008 * 009 * http://www.apache.org/licenses/LICENSE-2.0 010 * 011 * Unless required by applicable law or agreed to in writing, software 012 * distributed under the License is distributed on an "AS IS" BASIS, 013 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. 014 * See the License for the specific language governing permissions and 015 * limitations under the License. 016 */ 017package org.apache.commons.math4.legacy.analysis.integration.gauss; 018 019import org.apache.commons.math4.legacy.analysis.UnivariateFunction; 020import org.apache.commons.math4.legacy.exception.DimensionMismatchException; 021import org.apache.commons.math4.legacy.core.MathArrays; 022import org.apache.commons.math4.legacy.core.Pair; 023 024/** 025 * Class that implements the Gaussian rule for 026 * {@link #integrate(UnivariateFunction) integrating} a weighted 027 * function. 028 * 029 * @since 3.1 030 */ 031public class GaussIntegrator { 032 /** Nodes. */ 033 private final double[] points; 034 /** Nodes weights. */ 035 private final double[] weights; 036 037 /** 038 * Creates an integrator from the given {@code points} and {@code weights}. 039 * The integration interval is defined by the first and last value of 040 * {@code points} which must be sorted in increasing order. 041 * 042 * @param points Integration points. 043 * @param weights Weights of the corresponding integration nodes. 044 * @throws org.apache.commons.math4.legacy.exception.NonMonotonicSequenceException if the {@code points} are not 045 * sorted in increasing order. 046 * @throws DimensionMismatchException if points and weights don't have the same length 047 */ 048 public GaussIntegrator(double[] points, 049 double[] weights) { 050 if (points.length != weights.length) { 051 throw new DimensionMismatchException(points.length, 052 weights.length); 053 } 054 055 MathArrays.checkOrder(points, MathArrays.OrderDirection.INCREASING, true, true); 056 057 this.points = points.clone(); 058 this.weights = weights.clone(); 059 } 060 061 /** 062 * Creates an integrator from the given pair of points (first element of 063 * the pair) and weights (second element of the pair. 064 * 065 * @param pointsAndWeights Integration points and corresponding weights. 066 * @throws org.apache.commons.math4.legacy.exception.NonMonotonicSequenceException if the {@code points} are not 067 * sorted in increasing order. 068 * 069 * @see #GaussIntegrator(double[], double[]) 070 */ 071 public GaussIntegrator(Pair<double[], double[]> pointsAndWeights) { 072 this(pointsAndWeights.getFirst(), pointsAndWeights.getSecond()); 073 } 074 075 /** 076 * Returns an estimate of the integral of {@code f(x) * w(x)}, 077 * where {@code w} is a weight function that depends on the actual 078 * flavor of the Gauss integration scheme. 079 * The algorithm uses the points and associated weights, as passed 080 * to the {@link #GaussIntegrator(double[],double[]) constructor}. 081 * 082 * @param f Function to integrate. 083 * @return the integral of the weighted function. 084 */ 085 public double integrate(UnivariateFunction f) { 086 double s = 0; 087 double c = 0; 088 for (int i = 0; i < points.length; i++) { 089 final double x = points[i]; 090 final double w = weights[i]; 091 final double y = w * f.value(x) - c; 092 final double t = s + y; 093 c = (t - s) - y; 094 s = t; 095 } 096 return s; 097 } 098 099 /** 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}