mirror of
https://github.com/gonum/gonum.git
synced 2025-11-01 02:52:49 +08:00
315 lines
9.9 KiB
Go
315 lines
9.9 KiB
Go
// 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
|
||
// 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 = i + 1 - m
|
||
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,
|
||
}
|