com.opengamma.strata.math.impl.integration.GaussianQuadratureIntegrator1D Maven / Gradle / Ivy
Show all versions of strata-math Show documentation
/*
* Copyright (C) 2009 - present by OpenGamma Inc. and the OpenGamma group of companies
*
* Please see distribution for license.
*/
package com.opengamma.strata.math.impl.integration;
import java.util.Objects;
import java.util.function.Function;
import com.opengamma.strata.collect.ArgChecker;
import com.opengamma.strata.math.impl.function.special.OrthogonalPolynomialFunctionGenerator;
/**
* Class that performs integration using Gaussian quadrature.
*
* If a function $f(x)$ can be written as $f(x) = W(x)g(x)$, where $g(x)$ is
* approximately polynomial, then for suitably chosen weights $w_i$ and points
* $x_i$, the integral can be approximated as:
* $$
* \begin{align*}
* \int_{-1}^1 f(x)dx
* &=\int_{-1}^1 W(x)g(x)dx\\
* &\approx \sum_{\i=1}^{n} w_i f(x_i)
* \end{align*}
* $$
* The evaluation points, weights and valid limits of integration depend on the type of orthogonal
* polynomials that are used
* (see {@link OrthogonalPolynomialFunctionGenerator} and {@link GaussLaguerreWeightAndAbscissaFunction}).
*
*/
public abstract class GaussianQuadratureIntegrator1D extends Integrator1D {
private final int size;
private final QuadratureWeightAndAbscissaFunction generator;
private final GaussianQuadratureData quadrature;
/**
* Creates an instance.
*
* @param n The number of sample points to be used in the integration, not negative or zero
* @param generator The generator of weights and abscissas
*/
public GaussianQuadratureIntegrator1D(int n, QuadratureWeightAndAbscissaFunction generator) {
ArgChecker.isTrue(n > 0, "number of intervals must be > 0");
ArgChecker.notNull(generator, "generating function");
this.size = n;
this.generator = generator;
this.quadrature = generator.generate(size);
}
/**
* {@inheritDoc}
*/
@Override
public Double integrate(Function function, Double lower, Double upper) {
ArgChecker.notNull(function, "function");
ArgChecker.notNull(lower, "lower");
ArgChecker.notNull(upper, "upper");
Function integral = getIntegralFunction(function, lower, upper);
return integrateFromPolyFunc(integral);
}
/**
* If a function $g(x)$ can be written as $W(x)f(x)$, where the weight function $W(x)$ corresponds
* to one of the Gaussian quadrature forms, then we may approximate the integral of $g(x)$ over
* a specific range as $\int^b_a g(x) dx =\int^b_a W(x)f(x) dx \approx \sum_{i=0}^{N-1} w_i f(x_i)$,
* were the abscissas $x_i$ and the weights $w_i$ have been precomputed. This is accurate
* if $f(x)$ can be approximated by a polynomial.
*
* @param polyFunction The function $f(x)$ rather than the full function $g(x) = W(x)f(x)$
* This should be well approximated by a polynomial.
* @return The integral
*/
public double integrateFromPolyFunc(Function polyFunction) {
ArgChecker.notNull(polyFunction, "polyFunction");
double[] abscissas = quadrature.getAbscissas();
int n = abscissas.length;
double[] weights = quadrature.getWeights();
double sum = 0;
for (int i = 0; i < n; i++) {
sum += polyFunction.apply(abscissas[i]) * weights[i];
}
return sum;
}
/**
* Gets the limits.
*
* @return The lower and upper limits for which the quadrature is valid
*/
public abstract Double[] getLimits();
/**
* Returns a function that is valid for both the type of quadrature and the limits of integration.
* @param function The function to be integrated, not null
* @param lower The lower integration limit, not null
* @param upper The upper integration limit, not null
* @return A function in the appropriate form for integration
*/
public abstract Function getIntegralFunction(
Function function,
Double lower,
Double upper);
@Override
public int hashCode() {
int prime = 31;
int result = 1;
result = prime * result + generator.hashCode();
result = prime * result + size;
return result;
}
@Override
public boolean equals(Object obj) {
if (this == obj) {
return true;
}
if (obj == null) {
return false;
}
if (getClass() != obj.getClass()) {
return false;
}
GaussianQuadratureIntegrator1D other = (GaussianQuadratureIntegrator1D) obj;
if (this.size != other.size) {
return false;
}
return Objects.equals(this.generator, other.generator);
}
}