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 java.util.Arrays; 020 021import org.apache.commons.math4.legacy.analysis.polynomials.PolynomialFunction; 022import org.apache.commons.math4.legacy.analysis.polynomials.PolynomialsUtils; 023import org.apache.commons.math4.legacy.linear.EigenDecomposition; 024import org.apache.commons.math4.legacy.linear.MatrixUtils; 025import org.apache.commons.math4.legacy.linear.RealMatrix; 026import org.apache.commons.math4.legacy.core.Pair; 027 028/** 029 * Factory that creates Gauss-type quadrature rule using Laguerre polynomials. 030 * 031 * @see <a href="http://en.wikipedia.org/wiki/Gauss%E2%80%93Laguerre_quadrature">Gauss-Laguerre quadrature (Wikipedia)</a> 032 * @since 4.0 033 */ 034public class LaguerreRuleFactory extends BaseRuleFactory<Double> { 035 /** {@inheritDoc} */ 036 @Override 037 protected Pair<Double[], Double[]> computeRule(int numberOfPoints) { 038 final RealMatrix companionMatrix = companionMatrix(numberOfPoints); 039 final EigenDecomposition eigen = new EigenDecomposition(companionMatrix); 040 final double[] roots = eigen.getRealEigenvalues(); 041 Arrays.sort(roots); 042 043 final Double[] points = new Double[numberOfPoints]; 044 final Double[] weights = new Double[numberOfPoints]; 045 046 final int n1 = numberOfPoints + 1; 047 final long n1Squared = n1 * (long) n1; 048 final PolynomialFunction laguerreN1 = PolynomialsUtils.createLaguerrePolynomial(n1); 049 for (int i = 0; i < numberOfPoints; i++) { 050 final double xi = roots[i]; 051 points[i] = xi; 052 053 final double val = laguerreN1.value(xi); 054 weights[i] = xi / n1Squared / (val * val); 055 } 056 057 return new Pair<>(points, weights); 058 } 059 060 /** 061 * @param degree Matrix dimension. 062 * @return a square matrix. 063 */ 064 private RealMatrix companionMatrix(final int degree) { 065 final RealMatrix c = MatrixUtils.createRealMatrix(degree, degree); 066 067 for (int i = 0; i < degree; i++) { 068 c.setEntry(i, i, 2.0 * i + 1); 069 if (i + 1 < degree) { 070 // subdiagonal 071 c.setEntry(i+1, i, -(i + 1)); 072 } 073 if (i - 1 >= 0) { 074 // superdiagonal 075 c.setEntry(i-1, i, -i); 076 } 077 } 078 079 return c; 080 } 081}