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 * Base integrator for implementations based on Simpson's method.
23 * Simpson's method is an optimization of Trapezoidal quadrature integrator.
24 * Implementations of this class will in general be more efficient than
25 * Trapezoidal quadrature integrators (i.e., require fewer function evaluations) when the function
26 * to be integrated has a finite fourth derivative (i.e., a continuous third derivative).
27 *
28 * @param <T> a quadrature.
29 */
30 public abstract class SimpsonIntegrator<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 = 20;
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 SimpsonIntegrator(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 double st;
79 var ost = 0.0;
80 var os = 0.0;
81 for (var j = 0; j < JMAX; j++) {
82 st = q.next();
83 s = (4.0 * st - ost) / 3.0;
84 if (j > JMIN && (Math.abs(s - os) < eps * Math.abs(os) || (s == 0.0 && os == 0.0))) {
85 return s;
86 }
87 os = s;
88 ost = st;
89 }
90 } catch (final EvaluationException e) {
91 throw new IntegrationException(e);
92 }
93
94 // Too many steps
95 throw new IntegrationException();
96 }
97
98 /**
99 * Gets type of integrator.
100 *
101 * @return type of integrator.
102 */
103 @Override
104 public IntegratorType getIntegratorType() {
105 return IntegratorType.SIMPSON;
106 }
107
108 /**
109 * Creates an integrator using Simpson's method.
110 *
111 * @param a Lower limit of integration.
112 * @param b Upper limit of integration.
113 * @param listener listener to evaluate a single dimension function at required points.
114 * @param eps required accuracy.
115 * @param quadratureType quadrature type.
116 * @return created integrator.
117 * @throws IllegalArgumentException if provided quadrature type is not supported.
118 */
119 public static SimpsonIntegrator<Quadrature> create(
120 final double a, final double b, final SingleDimensionFunctionEvaluatorListener listener, final double eps,
121 final QuadratureType quadratureType) {
122 return switch (quadratureType) {
123 case TRAPEZOIDAL -> cast(new SimpsonTrapezoidalQuadratureIntegrator(a, b, listener, eps));
124 case MID_POINT -> cast(new SimpsonMidPointQuadratureIntegrator(a, b, listener, eps));
125 case INFINITY_MID_POINT -> cast(new SimpsonInfinityMidPointQuadratureIntegrator(a, b, listener, eps));
126 case LOWER_SQUARE_ROOT_MID_POINT ->
127 cast(new SimpsonLowerSquareRootMidPointQuadratureIntegrator(a, b, listener, eps));
128 case UPPER_SQUARE_ROOT_MID_POINT ->
129 cast(new SimpsonUpperSquareRootMidPointQuadratureIntegrator(a, b, listener, eps));
130 case DOUBLE_EXPONENTIAL_RULE ->
131 cast(new SimpsonDoubleExponentialRuleQuadratureIntegrator(a, b, listener, eps));
132 default -> throw new IllegalArgumentException();
133 };
134 }
135
136 /**
137 * Creates an integrator using Simpson's method and having default accuracy.
138 *
139 * @param a Lower limit of integration.
140 * @param b Upper limit of integration.
141 * @param listener listener to evaluate a single dimension function at required points.
142 * @param quadratureType quadrature type.
143 * @return created integrator.
144 * @throws IllegalArgumentException if provided quadrature type is not supported.
145 */
146 public static SimpsonIntegrator<Quadrature> create(
147 final double a, final double b, final SingleDimensionFunctionEvaluatorListener listener,
148 final QuadratureType quadratureType) {
149 return switch (quadratureType) {
150 case TRAPEZOIDAL -> cast(new SimpsonTrapezoidalQuadratureIntegrator(a, b, listener));
151 case MID_POINT -> cast(new SimpsonMidPointQuadratureIntegrator(a, b, listener));
152 case INFINITY_MID_POINT -> cast(new SimpsonInfinityMidPointQuadratureIntegrator(a, b, listener));
153 case LOWER_SQUARE_ROOT_MID_POINT ->
154 cast(new SimpsonLowerSquareRootMidPointQuadratureIntegrator(a, b, listener));
155 case UPPER_SQUARE_ROOT_MID_POINT ->
156 cast(new SimpsonUpperSquareRootMidPointQuadratureIntegrator(a, b, listener));
157 case DOUBLE_EXPONENTIAL_RULE -> cast(new SimpsonDoubleExponentialRuleQuadratureIntegrator(a, b, listener));
158 default -> throw new IllegalArgumentException();
159 };
160 }
161
162 /**
163 * Creates an integrator using Simpson's method and default quadrature type.
164 *
165 * @param a Lower limit of integration.
166 * @param b Upper limit of integration.
167 * @param listener listener to evaluate a single dimension function at required points.
168 * @param eps required accuracy.
169 * @return created integrator.
170 */
171 public static SimpsonIntegrator<Quadrature> create(
172 final double a, final double b, final SingleDimensionFunctionEvaluatorListener listener, final double eps) {
173 return create(a, b, listener, eps, DEFAULT_QUADRATURE_TYPE);
174 }
175
176 /**
177 * Creates an integrator using Simpson's method and having default accuracy and quadrature type.
178 *
179 * @param a Lower limit of integration.
180 * @param b Upper limit of integration.
181 * @param listener listener to evaluate a single dimension function at required points.
182 * @return created integrator.
183 */
184 public static SimpsonIntegrator<Quadrature> create(
185 final double a, final double b, final SingleDimensionFunctionEvaluatorListener listener) {
186 return create(a, b, listener, DEFAULT_QUADRATURE_TYPE);
187 }
188
189 /**
190 * Casts integrator to a quadrature integrator without wildcard parameter.
191 *
192 * @param integrator integrator to be cast.
193 * @return cast integrator.
194 */
195 private static SimpsonIntegrator<Quadrature> cast(final SimpsonIntegrator<?> integrator) {
196 //noinspection unchecked
197 return (SimpsonIntegrator<Quadrature>) integrator;
198 }
199 }