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.analysis.integration.gauss;
018
019import java.util.Arrays;
020
021import org.apache.commons.math4.analysis.polynomials.PolynomialFunction;
022import org.apache.commons.math4.analysis.polynomials.PolynomialsUtils;
023import org.apache.commons.math4.linear.EigenDecomposition;
024import org.apache.commons.math4.linear.MatrixUtils;
025import org.apache.commons.math4.linear.RealMatrix;
026import org.apache.commons.math4.util.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 * 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}