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 * Implementation of quadrature using double exponential, which allows integration with a variable
23 * transformation.
24 * This implementation is suitable for improper integrands containing singularities.
25 */
26 public class DoubleExponentialRuleQuadrature extends Quadrature {
27
28 /**
29 * Default transformation of range of integration.
30 */
31 public static final double DEFAULT_HMAX = 3.7;
32
33 /**
34 * Lower limit of integration.
35 */
36 private final double a;
37
38 /**
39 * Upper limit of integration.
40 */
41 private final double b;
42
43 /**
44 * Maximum step size. Determines transformation of range of integration.
45 */
46 private final double hmax;
47
48 /**
49 * Value of the next stage of refinement.
50 */
51 private double s;
52
53 /**
54 * Listener to evaluate single dimension functions at required points.
55 */
56 private final DoubleExponentialSingleDimensionFunctionEvaluatorListener listener;
57
58 /**
59 * Constructor.
60 *
61 * @param listener listener to evaluate function and handle singularities.
62 * @param a Lower limit of integration.
63 * @param b Upper limit of integration.
64 * @param hmax Maximum step size. This quadrature transforms the range of integration to
65 * [-hmax, hmax].
66 */
67 public DoubleExponentialRuleQuadrature(
68 final DoubleExponentialSingleDimensionFunctionEvaluatorListener listener, final double a, final double b,
69 final double hmax) {
70 this.listener = listener;
71 this.a = a;
72 this.b = b;
73 this.hmax = hmax;
74 n = 0;
75 }
76
77 /**
78 * Constructor with default maximum step size.
79 *
80 * @param listener listener to evaluate function and handle singularities.
81 * @param a Lower limit of integration.
82 * @param b Upper limit of integration.
83 */
84 public DoubleExponentialRuleQuadrature(
85 final DoubleExponentialSingleDimensionFunctionEvaluatorListener listener, final double a, final double b) {
86 this(listener, a, b, DEFAULT_HMAX);
87 }
88
89 /**
90 * Constructor.
91 *
92 * @param listener listener to evaluate function if function has mild singularities.
93 * @param a Lower limit of integration.
94 * @param b Upper limit of integration.
95 * @param hmax Maximum step size. This quadrature transforms the range of integration to
96 * [-hmax, hmax].
97 */
98 public DoubleExponentialRuleQuadrature(
99 final SingleDimensionFunctionEvaluatorListener listener, final double a, final double b,
100 final double hmax) {
101 this((x, delta) -> listener.evaluate(x), a, b, hmax);
102 }
103
104 /**
105 * Constructor with default maximum step size.
106 *
107 * @param listener listener to evaluate function if function has mild singularities.
108 * @param a Lower limit of integration.
109 * @param b Upper limit of integration.
110 */
111 public DoubleExponentialRuleQuadrature(
112 final SingleDimensionFunctionEvaluatorListener listener, final double a, final double b) {
113 this(listener, a, b, DEFAULT_HMAX);
114 }
115
116 /**
117 * Returns the value of the integral at the nth stage of refinement.
118 *
119 * @return the value of the integral at the nth stage of refinement.
120 * @throws EvaluationException Raised if something failed during the evaluation.
121 */
122 @SuppressWarnings("Duplicates")
123 @Override
124 public double next() throws EvaluationException {
125 // On the first call to the function next (n = 1), the routine returns the crudest estimate
126 // of integral between a and b of f(x). Subsequent calls to next (n = 2, 3, ...) will
127 // improve the accuracy by adding 2^(n-1) additional interior points.
128 double del;
129 double fact;
130 double q;
131 double sum;
132 double t;
133 final double twoh;
134 int it;
135 int j;
136 n++;
137 if (n == 1) {
138 fact = 0.25;
139 s = hmax * 2.0 * (b - a) * fact * listener.evaluate(0.5 * (b + a), 0.5 * (b - a));
140 } else {
141 for (it = 1, j = 1; j < n - 1; j++) {
142 it <<= 1;
143 }
144
145 // Twice the spacing of the points to be added
146 twoh = hmax / it;
147 t = 0.5 * twoh;
148 for (sum = 0.0, j = 0; j < it; j++) {
149 q = Math.exp(-2.0 * Math.sinh(t));
150 del = (b - a) * q / (1.0 + q);
151 final var value = 1.0 + q;
152 fact = q / (value * value) * Math.cosh(t);
153 sum += fact * (listener.evaluate(a + del, del) + listener.evaluate(b - del, del));
154 t += twoh;
155 }
156
157 // Replace s by its refined value and return.
158 s = 0.5 * s + (b - a) * twoh * sum;
159 }
160 return s;
161 }
162
163 /**
164 * Gets type of quadrature.
165 *
166 * @return type of quadrature.
167 */
168 @Override
169 public QuadratureType getType() {
170 return QuadratureType.DOUBLE_EXPONENTIAL_RULE;
171 }
172 }