1   /*
2    * Licensed to the Apache Software Foundation (ASF) under one or more
3    * contributor license agreements.  See the NOTICE file distributed with
4    * this work for additional information regarding copyright ownership.
5    * The ASF licenses this file to You under the Apache License, Version 2.0
6    * (the "License"); you may not use this file except in compliance with
7    * the License.  You may obtain a copy of the License at
8    *
9    *      http://www.apache.org/licenses/LICENSE-2.0
10   *
11   * Unless required by applicable law or agreed to in writing, software
12   * distributed under the License is distributed on an "AS IS" BASIS,
13   * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
14   * See the License for the specific language governing permissions and
15   * limitations under the License.
16   */
17  package org.apache.commons.math.analysis.integration;
18  
19  import java.util.Random;
20  
21  import org.apache.commons.math.ConvergenceException;
22  import org.apache.commons.math.FunctionEvaluationException;
23  import org.apache.commons.math.MathException;
24  import org.apache.commons.math.analysis.QuinticFunction;
25  import org.apache.commons.math.analysis.SinFunction;
26  import org.apache.commons.math.analysis.UnivariateRealFunction;
27  import org.apache.commons.math.analysis.polynomials.PolynomialFunction;
28  
29  import junit.framework.*;
30  
31  public class LegendreGaussIntegratorTest
32  extends TestCase {
33  
34      public LegendreGaussIntegratorTest(String name) {
35          super(name);
36      }
37  
38      public void testSinFunction() throws MathException {
39          UnivariateRealFunction f = new SinFunction();
40          UnivariateRealIntegrator integrator = new LegendreGaussIntegrator(5, 64);
41          integrator.setAbsoluteAccuracy(1.0e-10);
42          integrator.setRelativeAccuracy(1.0e-14);
43          integrator.setMinimalIterationCount(2);
44          integrator.setMaximalIterationCount(15);
45          double min, max, expected, result, tolerance;
46  
47          min = 0; max = Math.PI; expected = 2;
48          tolerance = Math.max(integrator.getAbsoluteAccuracy(),
49                               Math.abs(expected * integrator.getRelativeAccuracy()));
50          result = integrator.integrate(f, min, max);
51          assertEquals(expected, result, tolerance);
52  
53          min = -Math.PI/3; max = 0; expected = -0.5;
54          tolerance = Math.max(integrator.getAbsoluteAccuracy(),
55                  Math.abs(expected * integrator.getRelativeAccuracy()));
56          result = integrator.integrate(f, min, max);
57          assertEquals(expected, result, tolerance);
58      }
59  
60      public void testQuinticFunction() throws MathException {
61          UnivariateRealFunction f = new QuinticFunction();
62          UnivariateRealIntegrator integrator = new LegendreGaussIntegrator(3, 64);
63          double min, max, expected, result;
64  
65          min = 0; max = 1; expected = -1.0/48;
66          result = integrator.integrate(f, min, max);
67          assertEquals(expected, result, 1.0e-16);
68  
69          min = 0; max = 0.5; expected = 11.0/768;
70          result = integrator.integrate(f, min, max);
71          assertEquals(expected, result, 1.0e-16);
72  
73          min = -1; max = 4; expected = 2048/3.0 - 78 + 1.0/48;
74          result = integrator.integrate(f, min, max);
75          assertEquals(expected, result, 1.0e-16);
76      }
77  
78      public void testExactIntegration()
79          throws ConvergenceException, FunctionEvaluationException {
80          Random random = new Random(86343623467878363l);
81          for (int n = 2; n < 6; ++n) {
82              LegendreGaussIntegrator integrator =
83                  new LegendreGaussIntegrator(n, 64);
84  
85              // an n points Gauss-Legendre integrator integrates 2n-1 degree polynoms exactly
86              for (int degree = 0; degree <= 2 * n - 1; ++degree) {
87                  for (int i = 0; i < 10; ++i) {
88                      double[] coeff = new double[degree + 1];
89                      for (int k = 0; k < coeff.length; ++k) {
90                          coeff[k] = 2 * random.nextDouble() - 1;
91                      }
92                      PolynomialFunction p = new PolynomialFunction(coeff);
93                      double result    = integrator.integrate(p, -5.0, 15.0);
94                      double reference = exactIntegration(p, -5.0, 15.0);
95                      assertEquals(n + " " + degree + " " + i, reference, result, 1.0e-12 * (1.0 + Math.abs(reference)));
96                  }
97              }
98  
99          }
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 }