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.math3.analysis.integration.gauss; 018 019import org.apache.commons.math3.analysis.UnivariateFunction; 020import org.apache.commons.math3.exception.DimensionMismatchException; 021import org.apache.commons.math3.exception.NonMonotonicSequenceException; 022import org.apache.commons.math3.util.MathArrays; 023import org.apache.commons.math3.util.Pair; 024 025/** 026 * Class that implements the Gaussian rule for 027 * {@link #integrate(UnivariateFunction) integrating} a weighted 028 * function. 029 * 030 * @since 3.1 031 */ 032public class GaussIntegrator { 033 /** Nodes. */ 034 private final double[] points; 035 /** Nodes weights. */ 036 private final double[] weights; 037 038 /** 039 * Creates an integrator from the given {@code points} and {@code weights}. 040 * The integration interval is defined by the first and last value of 041 * {@code points} which must be sorted in increasing order. 042 * 043 * @param points Integration points. 044 * @param weights Weights of the corresponding integration nodes. 045 * @throws NonMonotonicSequenceException if the {@code points} are not 046 * sorted in increasing order. 047 * @throws DimensionMismatchException if points and weights don't have the same length 048 */ 049 public GaussIntegrator(double[] points, 050 double[] weights) 051 throws NonMonotonicSequenceException, DimensionMismatchException { 052 if (points.length != weights.length) { 053 throw new DimensionMismatchException(points.length, 054 weights.length); 055 } 056 057 MathArrays.checkOrder(points, MathArrays.OrderDirection.INCREASING, true, true); 058 059 this.points = points.clone(); 060 this.weights = weights.clone(); 061 } 062 063 /** 064 * Creates an integrator from the given pair of points (first element of 065 * the pair) and weights (second element of the pair. 066 * 067 * @param pointsAndWeights Integration points and corresponding weights. 068 * @throws NonMonotonicSequenceException if the {@code points} are not 069 * sorted in increasing order. 070 * 071 * @see #GaussIntegrator(double[], double[]) 072 */ 073 public GaussIntegrator(Pair<double[], double[]> pointsAndWeights) 074 throws NonMonotonicSequenceException { 075 this(pointsAndWeights.getFirst(), pointsAndWeights.getSecond()); 076 } 077 078 /** 079 * Returns an estimate of the integral of {@code f(x) * w(x)}, 080 * where {@code w} is a weight function that depends on the actual 081 * flavor of the Gauss integration scheme. 082 * The algorithm uses the points and associated weights, as passed 083 * to the {@link #GaussIntegrator(double[],double[]) constructor}. 084 * 085 * @param f Function to integrate. 086 * @return the integral of the weighted function. 087 */ 088 public double integrate(UnivariateFunction f) { 089 double s = 0; 090 double c = 0; 091 for (int i = 0; i < points.length; i++) { 092 final double x = points[i]; 093 final double w = weights[i]; 094 final double y = w * f.value(x) - c; 095 final double t = s + y; 096 c = (t - s) - y; 097 s = t; 098 } 099 return s; 100 } 101 102 /** 103 * @return the order of the integration rule (the number of integration 104 * points). 105 */ 106 public int getNumberOfPoints() { 107 return points.length; 108 } 109 110 /** 111 * Gets the integration point at the given index. 112 * The index must be in the valid range but no check is performed. 113 * @param index index of the integration point 114 * @return the integration point. 115 */ 116 public double getPoint(int index) { 117 return points[index]; 118 } 119 120 /** 121 * Gets the weight of the integration point at the given index. 122 * The index must be in the valid range but no check is performed. 123 * @param index index of the integration point 124 * @return the weight. 125 */ 126 public double getWeight(int index) { 127 return weights[index]; 128 } 129}