1 /*
2 * Copyright (C) 2023 Alberto Irurueta Carro (alberto@irurueta.com)
3 *
4 * Licensed under the Apache License, Version 2.0 (the "License");
5 * you may not use this file except in compliance with the License.
6 * You may obtain a copy of the License at
7 *
8 * http://www.apache.org/licenses/LICENSE-2.0
9 *
10 * Unless required by applicable law or agreed to in writing, software
11 * distributed under the License is distributed on an "AS IS" BASIS,
12 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
13 * See the License for the specific language governing permissions and
14 * limitations under the License.
15 */
16 package com.irurueta.numerical.integration;
17
18 import com.irurueta.numerical.EvaluationException;
19 import com.irurueta.numerical.SingleDimensionFunctionEvaluatorListener;
20
21 /**
22 * Integrates functions given a quadrature implementation up to desired accuracy.
23 * If assumptions can be made about the smoothness of a function other implementations such as
24 * Simpon's or Romberg's are more efficient and require less function evaluations. Otherwise, this
25 * is the simplest integrator that can be used for general purpose integrations when no assumptions
26 * can be made.
27 *
28 * @param <T> a quadrature.
29 */
30 public abstract class QuadratureIntegrator<T extends Quadrature> extends Integrator {
31
32 /**
33 * Default accuracy.
34 */
35 public static final double EPS = 1e-10;
36
37 /**
38 * Minimum required number of steps.
39 */
40 private static final int JMIN = 5;
41
42 /**
43 * Maximum number of allowed steps.
44 */
45 private static final int JMAX = 35;
46
47 /**
48 * Quadrature used for integration.
49 */
50 private final T q;
51
52 /**
53 * Required accuracy.
54 */
55 private final double eps;
56
57 /**
58 * Constructor.
59 *
60 * @param q Quadrature used for integration.
61 * @param eps Required accuracy.
62 */
63 protected QuadratureIntegrator(final T q, final double eps) {
64 this.q = q;
65 this.eps = eps;
66 }
67
68 /**
69 * Integrates function between provided lower and upper limits.
70 *
71 * @return result of integration.
72 * @throws IntegrationException if integration fails for numerical reasons.
73 */
74 @Override
75 public double integrate() throws IntegrationException {
76 try {
77 double s;
78 // Initial value of olds is arbitrary.
79 var olds = 0.0;
80
81 for (var j = 0; j < JMAX; j++) {
82 s = q.next();
83 if (j > JMIN && (Math.abs(s - olds) < eps * Math.abs(olds) || (s == 0.0 && olds == 0.0))) {
84 // Avoid spurious early convergence.
85 return s;
86 }
87 olds = s;
88 }
89 } catch (final EvaluationException e) {
90 throw new IntegrationException(e);
91 }
92
93 // Too many steps
94 throw new IntegrationException();
95 }
96
97 /**
98 * Gets type of integrator.
99 *
100 * @return type of integrator.
101 */
102 @Override
103 public IntegratorType getIntegratorType() {
104 return IntegratorType.QUADRATURE;
105 }
106
107 /**
108 * Creates a quadrature integrator.
109 *
110 * @param a Lower limit of integration.
111 * @param b Upper limit of integration.
112 * @param listener listener to evaluate a single dimension function at required points.
113 * @param eps required accuracy.
114 * @param quadratureType quadrature type.
115 * @return created integrator.
116 * @throws IllegalArgumentException if provided quadrature type is not supported.
117 */
118 public static QuadratureIntegrator<Quadrature> create(
119 final double a, final double b, final SingleDimensionFunctionEvaluatorListener listener, final double eps,
120 final QuadratureType quadratureType) {
121 return switch (quadratureType) {
122 case TRAPEZOIDAL -> cast(new TrapezoidalQuadratureIntegrator(a, b, listener, eps));
123 case MID_POINT -> cast(new MidPointQuadratureIntegrator(a, b, listener, eps));
124 case INFINITY_MID_POINT -> cast(new InfinityMidPointQuadratureIntegrator(a, b, listener, eps));
125 case LOWER_SQUARE_ROOT_MID_POINT ->
126 cast(new LowerSquareRootMidPointQuadratureIntegrator(a, b, listener, eps));
127 case UPPER_SQUARE_ROOT_MID_POINT ->
128 cast(new UpperSquareRootMidPointQuadratureIntegrator(a, b, listener, eps));
129 case DOUBLE_EXPONENTIAL_RULE -> cast(new DoubleExponentialRuleQuadratureIntegrator(a, b, listener, eps));
130 default -> throw new IllegalArgumentException();
131 };
132 }
133
134 /**
135 * Creates a quadrature integrator with default accuracy.
136 *
137 * @param a Lower limit of integration.
138 * @param b Upper limit of integration.
139 * @param listener listener to evaluate a single dimension function at required points.
140 * @param quadratureType quadrature type.
141 * @return created integrator.
142 * @throws IllegalArgumentException if provided quadrature type is not supported.
143 */
144 public static QuadratureIntegrator<Quadrature> create(
145 final double a, final double b, final SingleDimensionFunctionEvaluatorListener listener,
146 final QuadratureType quadratureType) {
147 return switch (quadratureType) {
148 case TRAPEZOIDAL -> cast(new TrapezoidalQuadratureIntegrator(a, b, listener));
149 case MID_POINT -> cast(new MidPointQuadratureIntegrator(a, b, listener));
150 case INFINITY_MID_POINT -> cast(new InfinityMidPointQuadratureIntegrator(a, b, listener));
151 case LOWER_SQUARE_ROOT_MID_POINT -> cast(new LowerSquareRootMidPointQuadratureIntegrator(a, b, listener));
152 case UPPER_SQUARE_ROOT_MID_POINT -> cast(new UpperSquareRootMidPointQuadratureIntegrator(a, b, listener));
153 case DOUBLE_EXPONENTIAL_RULE -> cast(new DoubleExponentialRuleQuadratureIntegrator(a, b, listener));
154 default -> throw new IllegalArgumentException();
155 };
156 }
157
158 /**
159 * Creates a quadrature integrator using default quadrature type.
160 *
161 * @param a Lower limit of integration.
162 * @param b Upper limit of integration.
163 * @param listener listener to evaluate a single dimension function at required points.
164 * @param eps required accuracy.
165 * @return created integrator.
166 */
167 public static QuadratureIntegrator<Quadrature> create(
168 final double a, final double b, final SingleDimensionFunctionEvaluatorListener listener, final double eps) {
169 return create(a, b, listener, eps, DEFAULT_QUADRATURE_TYPE);
170 }
171
172 /**
173 * Creates a quadrature integrator using default accuracy and quadrature type.
174 *
175 * @param a Lower limit of integration.
176 * @param b Upper limit of integration.
177 * @param listener listener to evaluate a single dimension function at required points.
178 * @return created integrator.
179 */
180 public static QuadratureIntegrator<Quadrature> create(
181 final double a, final double b, final SingleDimensionFunctionEvaluatorListener listener) {
182 return create(a, b, listener, DEFAULT_QUADRATURE_TYPE);
183 }
184
185 /**
186 * Casts integrator to a quadrature integrator without wildcard parameter.
187 * .
188 * @param integrator integrator to be cast.
189 * @return cast integrator.
190 */
191 private static QuadratureIntegrator<Quadrature> cast(QuadratureIntegrator<?> integrator) {
192 //noinspection unchecked
193 return (QuadratureIntegrator<Quadrature>) integrator;
194 }
195 }