1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
|
# Copyright 2017, 2018 Christoph Groth (CEA).
#
# This file is part of Vquad. It is subject to the license terms in the file
# LICENSE.rst found in the top-level directory of this distribution.
import numpy as np
from scipy.linalg import norm
from . import tables as tbls
eps = np.spacing(1)
# If the relative difference between two consecutive approximations is
# lower than this value, the error estimate is considered reliable.
# See section 6.2 of Pedro Gonnet's thesis.
hint = 0.1
# Smallest acceptable relative difference of points in a rule. This was chosen
# such that no artifacts are apparent in plots of (i, log(a_i)), where a_i is
# the sequence of estimates of the integral value of an interval and all its
# ancestors..
min_sep = 16 * eps
min_level = 1
max_level = 4
ndiv_max = 20
_sqrt_one_half = np.sqrt(0.5)
def _calc_coeffs(vals, level):
nans = np.flatnonzero(~np.isfinite(vals))
if nans.size:
# Replace vals by a copy and zero-out non-finite elements.
vals = vals.copy()
vals[nans] = 0
# Prepare things for the loop further down.
b = tbls.newton_coeffs[level].copy()
m = len(b) - 2 # = len(tbls.nodes[level]) - 1
coeffs = tbls.inv_Vs[level] @ vals
# This is a variant of Algorithm 7 from the thesis of Pedro Gonnet where no
# linear system has to be solved explicitly. Instead, Algorithm 5 is used.
for i in nans:
b[m + 1] /= tbls.alpha[m]
x = tbls.nodes[level][i]
b[m] = (b[m] + x * b[m + 1]) / tbls.alpha[m - 1]
for j in range(m - 1, 0, -1):
b[j] = ((b[j] + x * b[j + 1] - tbls.gamma[j + 1] * b[j + 2])
/ tbls.alpha[j - 1])
b = b[1:]
coeffs[:m] -= coeffs[m] / b[m] * b[:m]
coeffs[m] = 0
m -= 1
return coeffs
class DivergentIntegralError(ValueError):
def __init__(self, msg, igral, err, nr_points):
self.igral = igral
self.err = err
self.nr_points = nr_points
super().__init__(msg)
class _Interval:
__slots__ = ['a', 'b', 'coeffs', 'vals', 'igral', 'err', 'level', 'depth',
'ndiv', 'c00', 'unreliable_err']
def __init__(self, a, b, level, depth):
self.a = a
self.b = b
self.level = level
self.depth = depth
def points(self):
a = self.a
b = self.b
return (a + b) / 2 + (b - a) * tbls.nodes[self.level] / 2
def interpolate(self, vals, coeffs_old=None):
self.vals = vals
self.coeffs = coeffs = _calc_coeffs(self.vals, self.level)
if self.level == min_level:
self.c00 = coeffs[0]
if coeffs_old is None:
coeffs_diff = norm(coeffs)
else:
coeffs_diff = np.zeros(max(len(coeffs_old), len(coeffs)))
coeffs_diff[:len(coeffs_old)] = coeffs_old
coeffs_diff[:len(coeffs)] -= coeffs
coeffs_diff = norm(coeffs_diff)
w = self.b - self.a
self.igral = w * coeffs[0] * _sqrt_one_half
self.err = w * coeffs_diff
self.unreliable_err = coeffs_diff > hint * norm(coeffs)
class Vquad:
"""Evaluate an integral using adaptive quadrature.
The algorithm uses Clenshaw-Curtis quadrature rules of increasing
degree in each interval. The error estimate is
sqrt(integrate((f0(x) - f1(x))**2)), where f0 and f1 are two
successive interpolations of the integrand. To fall below the
desired total error, intervals are worked on ranked by their own
absolute error: either the degree of the rule is increased or the
interval is split if either the function does not appear to be
smooth or a rule of maximum degree has been reached.
Reference: "Increasing the Reliability of Adaptive Quadrature
Using Explicit Interpolants", P. Gonnet, ACM Transactions on
Mathematical Software, 37 (3), art. no. 26, 2008.
"""
def __init__(self, f, a, b, level=max_level - 1):
ival = _Interval(a, b, level, 1)
vals = f(ival.points())
ival.interpolate(vals)
ival.c00 = 0.0 # Will go away.
ival.ndiv = 0
self.ivals = [ival] # Active intervals
self.attic = [] # Inactive intervals
self.f = f
self.nr_points = len(vals)
self.igral_excess = 0
self.err_excess = 0
self.i_max = 0
def split(self, ival):
m = (ival.a + ival.b) / 2
f_center = ival.vals[(len(ival.vals) - 1) // 2]
depth = ival.depth + 1
children = [_Interval(ival.a, m, min_level, depth),
_Interval(m, ival.b, min_level, depth)]
points = np.concatenate([child.points()[1:-1] for child in children])
self.nr_points += len(points)
valss = np.empty((2, tbls.sizes[min_level]))
valss[:, 0] = ival.vals[0], f_center
valss[:, -1] = f_center, ival.vals[-1]
valss[:, 1:-1] = self.f(points).reshape((2, -1))
for child, vals, T in zip(children, valss, tbls.Ts):
child.interpolate(vals, T[:, :ival.coeffs.shape[0]] @ ival.coeffs)
child.ndiv = (ival.ndiv
+ (ival.c00 and child.c00 / ival.c00 > 2))
if child.ndiv > ndiv_max and 2*child.ndiv > child.depth:
msg = ('Possibly divergent integral in the interval '
'[{}, {}]! (h={})')
raise DivergentIntegralError(
msg.format(child.a, child.b, child.b - child.a),
child.igral * np.inf, None, self.nr_points)
return children
def refine(self, ival):
"""Increase degree of interval."""
ival.level += 1
points = ival.points()
vals = np.empty(points.shape)
vals[0::2] = ival.vals
vals[1::2] = self.f(points[1::2])
self.nr_points += (len(vals) - 1) // 2
ival.interpolate(vals, ival.coeffs)
return points
def improve(self):
i_max = self.i_max
if self.ivals[i_max].level == max_level:
split = True
else:
points = self.refine(self.ivals[i_max])
split = self.ivals[i_max].unreliable_err
if (points[1] - points[0] < points[0] * min_sep
or points[-1] - points[-2] < points[-2] * min_sep
or (self.ivals[i_max].err
< (abs(self.ivals[i_max].igral) * eps
* tbls.V_cond_nums[self.ivals[i_max].level]))):
# Remove the interval (while remembering the excess integral
# and error), since it is either too narrow, or the estimated
# relative error is already at the limit of numerical accuracy
# and cannot be reduced further.
self.err_excess += self.ivals[i_max].err
self.igral_excess += self.ivals[i_max].igral
self.ivals[i_max] = self.ivals[-1]
self.attic.append(self.ivals.pop())
return
if split:
self.ivals.extend(self.split(self.ivals[i_max]))
self.ivals[i_max] = self.ivals.pop()
def totals(self):
# Compute the total error and new max.
i_max = 0
err = self.err_excess
igral = self.igral_excess
for i in range(len(self.ivals)):
if self.ivals[i].err > self.ivals[i_max].err:
i_max = i
err += self.ivals[i].err
igral += self.ivals[i].igral
self.i_max = i_max
return igral, err
def improve_until(self, tol):
while True:
self.improve()
igral, err = self.totals()
if (err == 0
or err < abs(igral) * tol
or err - self.err_excess < abs(igral) * tol < self.err_excess
or not self.ivals):
return igral, err, self.nr_points
def vquad(f, a, b, tol):
igrator = Vquad(f, a, b)
return igrator.improve_until(tol)
|
