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     */
017    package org.apache.commons.math.analysis.integration;
018    
019    import java.util.Random;
020    
021    import org.apache.commons.math.ConvergenceException;
022    import org.apache.commons.math.FunctionEvaluationException;
023    import org.apache.commons.math.MathException;
024    import org.apache.commons.math.analysis.QuinticFunction;
025    import org.apache.commons.math.analysis.SinFunction;
026    import org.apache.commons.math.analysis.UnivariateRealFunction;
027    import org.apache.commons.math.analysis.polynomials.PolynomialFunction;
028    
029    import junit.framework.*;
030    
031    public class LegendreGaussIntegratorTest
032    extends TestCase {
033    
034        public LegendreGaussIntegratorTest(String name) {
035            super(name);
036        }
037    
038        public void testSinFunction() throws MathException {
039            UnivariateRealFunction f = new SinFunction();
040            UnivariateRealIntegrator integrator = new LegendreGaussIntegrator(5, 64);
041            integrator.setAbsoluteAccuracy(1.0e-10);
042            integrator.setRelativeAccuracy(1.0e-14);
043            integrator.setMinimalIterationCount(2);
044            integrator.setMaximalIterationCount(15);
045            double min, max, expected, result, tolerance;
046    
047            min = 0; max = Math.PI; expected = 2;
048            tolerance = Math.max(integrator.getAbsoluteAccuracy(),
049                                 Math.abs(expected * integrator.getRelativeAccuracy()));
050            result = integrator.integrate(f, min, max);
051            assertEquals(expected, result, tolerance);
052    
053            min = -Math.PI/3; max = 0; expected = -0.5;
054            tolerance = Math.max(integrator.getAbsoluteAccuracy(),
055                    Math.abs(expected * integrator.getRelativeAccuracy()));
056            result = integrator.integrate(f, min, max);
057            assertEquals(expected, result, tolerance);
058        }
059    
060        public void testQuinticFunction() throws MathException {
061            UnivariateRealFunction f = new QuinticFunction();
062            UnivariateRealIntegrator integrator = new LegendreGaussIntegrator(3, 64);
063            double min, max, expected, result;
064    
065            min = 0; max = 1; expected = -1.0/48;
066            result = integrator.integrate(f, min, max);
067            assertEquals(expected, result, 1.0e-16);
068    
069            min = 0; max = 0.5; expected = 11.0/768;
070            result = integrator.integrate(f, min, max);
071            assertEquals(expected, result, 1.0e-16);
072    
073            min = -1; max = 4; expected = 2048/3.0 - 78 + 1.0/48;
074            result = integrator.integrate(f, min, max);
075            assertEquals(expected, result, 1.0e-16);
076        }
077    
078        public void testExactIntegration()
079            throws ConvergenceException, FunctionEvaluationException {
080            Random random = new Random(86343623467878363l);
081            for (int n = 2; n < 6; ++n) {
082                LegendreGaussIntegrator integrator =
083                    new LegendreGaussIntegrator(n, 64);
084    
085                // an n points Gauss-Legendre integrator integrates 2n-1 degree polynoms exactly
086                for (int degree = 0; degree <= 2 * n - 1; ++degree) {
087                    for (int i = 0; i < 10; ++i) {
088                        double[] coeff = new double[degree + 1];
089                        for (int k = 0; k < coeff.length; ++k) {
090                            coeff[k] = 2 * random.nextDouble() - 1;
091                        }
092                        PolynomialFunction p = new PolynomialFunction(coeff);
093                        double result    = integrator.integrate(p, -5.0, 15.0);
094                        double reference = exactIntegration(p, -5.0, 15.0);
095                        assertEquals(n + " " + degree + " " + i, reference, result, 1.0e-12 * (1.0 + Math.abs(reference)));
096                    }
097                }
098    
099            }
100        }
101    
102        private double exactIntegration(PolynomialFunction p, double a, double b) {
103            final double[] coeffs = p.getCoefficients();
104            double yb = coeffs[coeffs.length - 1] / coeffs.length;
105            double ya = yb;
106            for (int i = coeffs.length - 2; i >= 0; --i) {
107                yb = yb * b + coeffs[i] / (i + 1);
108                ya = ya * a + coeffs[i] / (i + 1);
109            }
110            return yb * b - ya * a;
111        }
112    
113        public static Test suite() {
114            return new TestSuite(LegendreGaussIntegratorTest.class);
115        }
116    
117    }