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 mid-point algorithm.
23 * Mid-point algorithm is suitable for improper integrations, which consists of the following
24 * situations:
25 * - integrand goes to a finite limiting value at finite upper and lower limits, but cannot be
26 * evaluated right on one of those limits (e.g. sin(x)/x at x = 0)
27 * - the upper limit is positive infinity, or the lower limit is negative infinity.
28 * - integrand has an integrable singularity at either limit (e.g. x^-0.5 at x = 0)
29 * - integrand has an integrable singularity at a known place between its upper and lower limits.
30 * - integrand has an integrable singularity at an unknown place between its upper and lower limits.
31 * <p>
32 * Impossible integrals, such as those being infinite (e.g. integral from 1 to infinity of x^-1) or
33 * those that have no limiting sense (e.g. integral from -infinity to infinity of cos(x)) cannot
34 * be handled by this implementation.
35 */
36 public class MidPointQuadrature extends Quadrature {
37
38 /**
39 * Lower limit of integration.
40 */
41 private final double a;
42
43 /**
44 * Upper limit of integration.
45 */
46 private final double b;
47
48 /**
49 * Current value of integral.
50 */
51 private double s;
52
53 /**
54 * Listener to evaluate single dimension functions at required points.
55 */
56 protected final SingleDimensionFunctionEvaluatorListener listener;
57
58 /**
59 * Constructor.
60 *
61 * @param a Lower limit of integration.
62 * @param b Upper limit of integration.
63 * @param listener listener to evaluate a single dimension function at required points.
64 */
65 public MidPointQuadrature(
66 final double a, final double b, final SingleDimensionFunctionEvaluatorListener listener) {
67 this.n = 0;
68 this.a = a;
69 this.b = b;
70 this.s = 0;
71 this.listener = listener;
72 }
73
74 /**
75 * Gets lower limit of integration.
76 *
77 * @return lower limit of integration.
78 */
79 public double getA() {
80 return a;
81 }
82
83 /**
84 * Gets upper limit of integration.
85 *
86 * @return upper limit of integration.
87 */
88 public double getB() {
89 return b;
90 }
91
92 /**
93 * Gets current value of integral.
94 *
95 * @return current value of integral.
96 */
97 public double getS() {
98 return s;
99 }
100
101 /**
102 * Returns the value of the integral at the nth stage of refinement.
103 *
104 * @return the value of the integral at the nth stage of refinement.
105 * @throws EvaluationException Raised if something failed during the evaluation.
106 */
107 @SuppressWarnings("Duplicates")
108 @Override
109 public double next() throws EvaluationException {
110 int it;
111 int j;
112 double x;
113 double tnm;
114 double sum;
115 double del;
116 double ddel;
117 n++;
118 if (n == 1) {
119 s = (b - a) * func(0.5 * (a + b));
120 return s;
121 } else {
122 for (it = 1, j = 1; j < n - 1; j++) {
123 it *= 3;
124 }
125 tnm = it;
126 del = (b - a) / (3.0 * tnm);
127 // The added points alternate in spacing between del and ddel
128 ddel = del + del;
129 x = a + 0.5 * del;
130 sum = 0.0;
131 for (j = 0; j < it; j++) {
132 sum += func(x);
133 x += ddel;
134 sum += func(x);
135 x += del;
136 }
137 // The new sum is combined with the old integral to give a refined integral
138 s = (s + (b - a) * sum / tnm) / 3.0;
139 return s;
140 }
141 }
142
143 /**
144 * Gets type of quadrature.
145 *
146 * @return type of quadrature.
147 */
148 @Override
149 public QuadratureType getType() {
150 return QuadratureType.MID_POINT;
151 }
152
153 /**
154 * Evaluates function at x.
155 *
156 * @param x point where function is evaluated.
157 * @return result of evaluation.
158 * @throws EvaluationException if evaluation fails.
159 */
160 protected double func(final double x) throws EvaluationException {
161 return listener.evaluate(x);
162 }
163 }