/
hermite.go
315 lines (282 loc) · 9.92 KB
/
hermite.go
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
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
// Copyright ©2016 The Gonum Authors. All rights reserved.
// Use of this source code is governed by a BSD-style
// license that can be found in the LICENSE file.
package quad
import (
"math"
"gonum.org/v1/gonum/floats"
"gonum.org/v1/gonum/mathext"
)
// Hermite generates sample locations and weights for performing quadrature with
// a squared-exponential weight
//
// int_-inf^inf e^(-x^2) f(x) dx .
type Hermite struct{}
func (h Hermite) FixedLocations(x, weight []float64, min, max float64) {
// TODO(btracey): Implement the case where x > 20, x < 200 so that we don't
// need to store all of that data.
// Algorithm adapted from Chebfun http://www.chebfun.org/.
//
// References:
// Algorithm:
// G. H. Golub and J. A. Welsch, "Calculation of Gauss quadrature rules",
// Math. Comp. 23:221-230, 1969.
// A. Glaser, X. Liu and V. Rokhlin, "A fast algorithm for the
// calculation of the roots of special functions", SIAM Journal
// on Scientific Computing", 29(4):1420-1438:, 2007.
// A. Townsend, T. Trogdon, and S.Olver, Fast computation of Gauss quadrature
// nodes and weights on the whole real line, IMA J. Numer. Anal., 36: 337–358,
// 2016. http://arxiv.org/abs/1410.5286
if len(x) != len(weight) {
panic("hermite: slice length mismatch")
}
if min >= max {
panic("hermite: min >= max")
}
if !math.IsInf(min, -1) || !math.IsInf(max, 1) {
panic("hermite: non-infinite bound")
}
h.locations(x, weight)
}
func (h Hermite) locations(x, weights []float64) {
n := len(x)
switch {
case 0 < n && n <= 200:
copy(x, xCacheHermite[n-1])
copy(weights, wCacheHermite[n-1])
case n > 200:
h.locationsAsy(x, weights)
}
}
// Algorithm adapted from Chebfun http://www.chebfun.org/. Specific code
// https://github.com/chebfun/chebfun/blob/development/hermpts.m.
// Original Copyright Notice:
/*
Copyright (c) 2015, The Chancellor, Masters and Scholars of the University
of Oxford, and the Chebfun Developers. All rights reserved.
Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions are met:
* Redistributions of source code must retain the above copyright
notice, this list of conditions and the following disclaimer.
* Redistributions in binary form must reproduce the above copyright
notice, this list of conditions and the following disclaimer in the
documentation and/or other materials provided with the distribution.
* Neither the name of the University of Oxford nor the names of its
contributors may be used to endorse or promote products derived from
this software without specific prior written permission.
THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR
ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
(INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
*/
// locationAsy returns the node locations and weights of a Hermite quadrature rule
// with len(x) points.
func (h Hermite) locationsAsy(x, w []float64) {
// A. Townsend, T. Trogdon, and S.Olver, Fast computation of Gauss quadrature
// nodes and weights the whole real line, IMA J. Numer. Anal.,
// 36: 337–358, 2016. http://arxiv.org/abs/1410.5286
// Find the positive locations and weights.
n := len(x)
l := n / 2
xa := x[l:]
wa := w[l:]
for i := range xa {
xa[i], wa[i] = h.locationsAsy0(i, n)
}
// Flip around zero -- copy the negative x locations with the corresponding
// weights.
if n%2 == 0 {
l--
}
for i, v := range xa {
x[l-i] = -v
}
for i, v := range wa {
w[l-i] = v
}
sumW := floats.Sum(w)
c := math.SqrtPi / sumW
floats.Scale(c, w)
}
// locationsAsy0 returns the location and weight for location i in an n-point
// quadrature rule. The rule is symmetric, so i should be <= n/2 + n%2.
func (h Hermite) locationsAsy0(i, n int) (x, w float64) {
const convTol = 1e-16
const convIter = 20
theta0 := h.hermiteInitialGuess(i, n)
t0 := theta0 / math.Sqrt(2*float64(n)+1)
theta0 = math.Acos(t0)
sqrt2np1 := math.Sqrt(2*float64(n) + 1)
var vali, dvali float64
for k := 0; k < convIter; k++ {
vali, dvali = h.hermpolyAsyAiry(i, n, theta0)
dt := -vali / (math.Sqrt2 * sqrt2np1 * dvali * math.Sin(theta0))
theta0 -= dt
if math.Abs(theta0) < convTol {
break
}
}
x = sqrt2np1 * math.Cos(theta0)
ders := x*vali + math.Sqrt2*dvali
w = math.Exp(-x*x) / (ders * ders)
return x, w
}
// hermpolyAsyAiry evaluates the Hermite polynomials using the Airy asymptotic
// formula in theta-space.
func (h Hermite) hermpolyAsyAiry(i, n int, t float64) (valVec, dvalVec float64) {
musq := 2*float64(n) + 1
cosT := math.Cos(t)
sinT := math.Sin(t)
sin2T := 2 * cosT * sinT
eta := 0.5*t - 0.25*sin2T
chi := -math.Pow(3*eta/2, 2.0/3)
phi := math.Pow(-chi/(sinT*sinT), 1.0/4)
cnst := 2 * math.SqrtPi * math.Pow(musq, 1.0/6) * phi
airy0 := real(mathext.AiryAi(complex(math.Pow(musq, 2.0/3)*chi, 0)))
airy1 := real(mathext.AiryAiDeriv(complex(math.Pow(musq, 2.0/3)*chi, 0)))
// Terms in 12.10.43:
const (
a1 = 15.0 / 144
b1 = -7.0 / 5 * a1
a2 = 5.0 * 7 * 9 * 11.0 / 2.0 / 144.0 / 144.0
b2 = -13.0 / 11 * a2
a3 = 7.0 * 9 * 11 * 13 * 15 * 17 / 6.0 / 144.0 / 144.0 / 144.0
b3 = -19.0 / 17 * a3
)
// Pre-compute terms.
cos2T := cosT * cosT
cos3T := cos2T * cosT
cos4T := cos3T * cosT
cos5T := cos4T * cosT
cos7T := cos5T * cos2T
cos9T := cos7T * cos2T
chi2 := chi * chi
chi3 := chi2 * chi
chi4 := chi3 * chi
chi5 := chi4 * chi
phi6 := math.Pow(phi, 6)
phi12 := phi6 * phi6
phi18 := phi12 * phi6
// u polynomials in 12.10.9.
u1 := (cos3T - 6*cosT) / 24.0
u2 := (-9*cos4T + 249*cos2T + 145) / 1152.0
u3 := (-4042*cos9T + 18189*cos7T - 28287*cos5T - 151995*cos3T - 259290*cosT) / 414720.0
val := airy0
B0 := -(phi6*u1 + a1) / chi2
val += B0 * airy1 / math.Pow(musq, 4.0/3)
A1 := (phi12*u2 + b1*phi6*u1 + b2) / chi3
val += A1 * airy0 / (musq * musq)
B1 := -(phi18*u3 + a1*phi12*u2 + a2*phi6*u1 + a3) / chi5
val += B1 * airy1 / math.Pow(musq, 4.0/3+2)
val *= cnst
// Derivative.
eta = 0.5*t - 0.25*sin2T
chi = -math.Pow(3*eta/2, 2.0/3)
phi = math.Pow(-chi/(sinT*sinT), 1.0/4)
cnst = math.Sqrt2 * math.SqrtPi * math.Pow(musq, 1.0/3) / phi
// v polynomials in 12.10.10.
v1 := (cos3T + 6*cosT) / 24
v2 := (15*cos4T - 327*cos2T - 143) / 1152
v3 := (259290*cosT + 238425*cos3T - 36387*cos5T + 18189*cos7T - 4042*cos9T) / 414720
C0 := -(phi6*v1 + b1) / chi
dval := C0 * airy0 / math.Pow(musq, 2.0/3)
dval += airy1
C1 := -(phi18*v3 + b1*phi12*v2 + b2*phi6*v1 + b3) / chi4
dval += C1 * airy0 / math.Pow(musq, 2.0/3+2)
D1 := (phi12*v2 + a1*phi6*v1 + a2) / chi3
dval += D1 * airy1 / (musq * musq)
dval *= cnst
return val, dval
}
// hermiteInitialGuess returns the initial guess for node i in an n-point Hermite
// quadrature rule. The rule is symmetric, so i should be <= n/2 + n%2.
func (h Hermite) hermiteInitialGuess(i, n int) float64 {
// There are two different formulas for the initial guesses of the hermite
// quadrature locations. The first uses the Gatteschi formula and is good
// near x = sqrt(n+0.5)
// [1] L. Gatteschi, Asymptotics and bounds for the zeros of Laguerre
// polynomials: a survey, J. Comput. Appl. Math., 144 (2002), pp. 7-27.
// The second is the Tricomi initial guesses, good near x = 0. This is
// equation 2.1 in [1] and is originally from
// [2] F. G. Tricomi, Sugli zeri delle funzioni di cui si conosce una
// rappresentazione asintotica, Ann. Mat. Pura Appl. 26 (1947), pp. 283-300.
// If the number of points is odd, there is a quadrature point at 1, which
// has an initial guess of 0.
if n%2 == 1 {
if i == 0 {
return 0
}
i--
}
m := n / 2
a := -0.5
if n%2 == 1 {
a = 0.5
}
nu := 4*float64(m) + 2*a + 2
// Find the split between Gatteschi guesses and Tricomi guesses.
p := 0.4985 + math.SmallestNonzeroFloat64
pidx := int(math.Floor(p * float64(n)))
// Use the Tricomi initial guesses in the first half where x is nearer to zero.
// Note: zeros of besselj(+/-.5,x) are integer and half-integer multiples of pi.
if i < pidx {
rhs := math.Pi * (4*float64(m) - 4*(float64(i)+1) + 3) / nu
tnk := math.Pi / 2
for k := 0; k < 7; k++ {
val := tnk - math.Sin(tnk) - rhs
dval := 1 - math.Cos(tnk)
dTnk := val / dval
tnk -= dTnk
if math.Abs(dTnk) < 1e-14 {
break
}
}
vc := math.Cos(tnk / 2)
t := vc * vc
return math.Sqrt(nu*t - (5.0/(4.0*(1-t)*(1-t))-1.0/(1-t)-1+3*a*a)/3/nu)
}
// Use Gatteschi guesses in the second half where x is nearer to sqrt(n+0.5)
i = m - (i + 1)
var ar float64
if i < len(airyRtsExact) {
ar = airyRtsExact[i]
} else {
t := 3.0 / 8 * math.Pi * (4*(float64(i)+1) - 1)
ar = math.Pow(t, 2.0/3) * (1 +
5.0/48*math.Pow(t, -2) -
5.0/36*math.Pow(t, -4) +
77125.0/82944*math.Pow(t, -6) -
108056875.0/6967296*math.Pow(t, -8) +
162375596875.0/334430208*math.Pow(t, -10))
}
r := nu + math.Pow(2, 2.0/3)*ar*math.Pow(nu, 1.0/3) +
0.2*math.Pow(2, 4.0/3)*ar*ar*math.Pow(nu, -1.0/3) +
(11.0/35-a*a-12.0/175*ar*ar*ar)/nu +
(16.0/1575*ar+92.0/7875*math.Pow(ar, 4))*math.Pow(2, 2.0/3)*math.Pow(nu, -5.0/3) -
(15152.0/3031875*math.Pow(ar, 5)+1088.0/121275*ar*ar)*math.Pow(2, 1.0/3)*math.Pow(nu, -7.0/3)
if r < 0 {
ar = 0
} else {
ar = math.Sqrt(r)
}
return ar
}
// airyRtsExact are the first airy roots.
var airyRtsExact = []float64{
-2.338107410459762,
-4.087949444130970,
-5.520559828095555,
-6.786708090071765,
-7.944133587120863,
-9.022650853340979,
-10.040174341558084,
-11.008524303733260,
-11.936015563236262,
-12.828776752865757,
}