From f211265a9a44b26ee94d80c124d47b9a68eb4098 Mon Sep 17 00:00:00 2001 From: Per A Brodtkorb Date: Wed, 22 Feb 2017 02:54:17 +0100 Subject: [PATCH] pep8 --- wafo/integrate.py | 620 +++++++++++++++++++++++----------------------- 1 file changed, 305 insertions(+), 315 deletions(-) diff --git a/wafo/integrate.py b/wafo/integrate.py index c69382a..a92e9bd 100644 --- a/wafo/integrate.py +++ b/wafo/integrate.py @@ -7,14 +7,14 @@ import scipy.special.orthogonal as ort from scipy import special as sp from scipy.integrate import simps, trapz -from .plotbackend import plotbackend as plt -from .demos import humps -from .misc import dea3 -from .dctpack import dct +from wafo.plotbackend import plotbackend as plt +from wafo.demos import humps +from wafo.misc import dea3 +from wafo.dctpack import dct # from pychebfun import Chebfun _EPS = np.finfo(float).eps -_POINTS_AND_WEIGHTS = {} +_NODES_AND_WEIGHTS = {} __all__ = ['dea3', 'clencurt', 'romberg', 'h_roots', 'j_roots', 'la_roots', 'p_roots', 'qrule', @@ -31,7 +31,7 @@ def _assert_warn(cond, msg): warnings.warn(msg) -def clencurt(fun, a, b, n0=5, trace=False): +def clencurt(fun, a, b, n=5, trace=False): """ Numerical evaluation of an integral, Clenshaw-Curtis method. @@ -78,31 +78,27 @@ def clencurt(fun, a, b, n0=5, trace=False): [2] Clenshaw, C.W. and Curtis, A.R. (1960), Numerische Matematik, Vol. 2, pp. 197--205 """ - - # make sure n is even - n = 2 * int(n0) + # make sure n_2 is even + n_2 = 2 * int(n) a, b = np.atleast_1d(a, b) a_shape = a.shape - af = a.ravel() - bf = b.ravel() + a = a.ravel() + b = b.ravel() - na = np.prod(a_shape) + a_size = np.prod(a_shape) - s = np.r_[0:n + 1] - s2 = np.r_[0:n + 1:2] - s2.shape = (-1, 1) - x1 = np.cos(np.pi * s / n) - x1.shape = (-1, 1) - x = x1 * (bf - af) / 2. + (bf + af) / 2 + s = np.c_[0:n_2 + 1:1] + s_2 = np.c_[0:n_2 + 1:2] + x = np.cos(np.pi * s / n_2) * (b - a) / 2. + (b + a) / 2 if hasattr(fun, '__call__'): f = fun(x) else: - x0 = np.flipud(fun[:, 0]) - n = len(x0) - 1 - if abs(x - x0) > 1e-8: + x_0 = np.flipud(fun[:, 0]) + n_2 = len(x_0) - 1 + if abs(x - x_0) > 1e-8: raise ValueError( - 'Input vector x must equal cos(pi*s/n)*(b-a)/2+(b+a)/2') + 'Input vector x must equal cos(pi*s/n_2)*(b-a)/2+(b+a)/2') f = np.flipud(fun[:, 1::]) @@ -112,27 +108,27 @@ def clencurt(fun, a, b, n0=5, trace=False): # using a Gauss-Lobatto variant, i.e., first and last # term f(a) and f(b) is multiplied with 0.5 f[0, :] = f[0, :] / 2 - f[n, :] = f[n, :] / 2 + f[n_2, :] = f[n_2, :] / 2 - # x = cos(pi*0:n/n) + # x = cos(pi*0:n_2/n_2) # f = f(x) # # N+1 # c(k) = (2/N) sum f''(n)*cos(pi*(2*k-2)*(n-1)/N), 1 <= k <= N/2+1. # n=1 - n0 = n // 2 + n = n_2 // 2 fft = np.fft.fft - tmp = np.real(fft(f[:n, :], axis=0)) - c = 2 / n * (tmp[0:n0 + 1, :] + np.cos(np.pi * s2) * f[n, :]) + tmp = np.real(fft(f[:n_2, :], axis=0)) + c = 2 / n_2 * (tmp[0:n + 1, :] + np.cos(np.pi * s_2) * f[n_2, :]) c[0, :] = c[0, :] / 2 - c[n0, :] = c[n0, :] / 2 + c[n, :] = c[n, :] / 2 - c = c[0:n0 + 1, :] / ((s2 - 1) * (s2 + 1)) - q_val = (af - bf) * np.sum(c, axis=0) + c = c[0:n + 1, :] / ((s_2 - 1) * (s_2 + 1)) + q_val = (a - b) * np.sum(c, axis=0) - abserr = (bf - af) * np.abs(c[n0, :]) + abserr = (b - a) * np.abs(c[n, :]) - if na > 1: + if a_size > 1: abserr = np.reshape(abserr, a_shape) q_val = np.reshape(q_val, a_shape) return q_val, abserr @@ -181,51 +177,54 @@ def romberg(fun, a, b, releps=1e-3, abseps=1e-3): rom[0, 0] = h * (fun(a) + fun(b)) / 2 ipower = 1 - fp = ones(table_limit) * 4 + f_p = ones(table_limit) * 4 - # ih1 = 0 - ih2 = 0. - ih4 = rom[0, 0] - abserr = ih4 + # q_val1 = 0 + q_val2 = 0. + q_val4 = rom[0, 0] + abserr = q_val4 # epstab = zeros(1,decdigs+7) # newflg = 1 - # [res,abserr,epstab,newflg] = dea(newflg,ih4,abserr,epstab) + # [res,abserr,epstab,newflg] = dea(newflg,q_val4,abserr,epstab) two = 1 one = 0 + converged = False for i in range(1, table_limit): h *= 0.5 - un5 = np.sum(fun(a + np.arange(1, 2 * ipower, 2) * h)) * h + u_n5 = np.sum(fun(a + np.arange(1, 2 * ipower, 2) * h)) * h # trapezoidal approximations - # T2n = 0.5 * (Tn + Un) = 0.5*Tn + un5 - rom[two, 0] = 0.5 * rom[one, 0] + un5 + # T2n = 0.5 * (Tn + Un) = 0.5*Tn + u_n5 + rom[two, 0] = 0.5 * rom[one, 0] + u_n5 - fp[i] = 4 * fp[i - 1] + f_p[i] = 4 * f_p[i - 1] # Richardson extrapolation for k in range(i): rom[two, k + 1] = (rom[two, k] + - (rom[two, k] - rom[one, k]) / (fp[k] - 1)) + (rom[two, k] - rom[one, k]) / (f_p[k] - 1)) - ih1 = ih2 - ih2 = ih4 - ih4 = rom[two, i] + q_val1 = q_val2 + q_val2 = q_val4 + q_val4 = rom[two, i] if 2 <= i: - res, abserr = dea3(ih1, ih2, ih4) - # ih4 = res - if abserr <= max(abseps, releps * abs(res)): + res, abserr = dea3(q_val1, q_val2, q_val4) + # q_val4 = res + converged = abserr <= max(abseps, releps * abs(res)) + if converged: break # rom(1,1:i) = rom(2,1:i) two = one one = (one + 1) % 2 ipower *= 2 + _assert(converged, "Integral did not converge to the required accuracy!") return res, abserr def _h_roots_newton(n, releps=3e-14, max_iter=10): - # PIM4=0.7511255444649425 - PIM4 = np.pi ** (-1. / 4) + # pim4=0.7511255444649425 + pim4 = np.pi ** (-1. / 4) # The roots are symmetric about the origin, so we have to # find only half of them. @@ -237,49 +236,42 @@ def _h_roots_newton(n, releps=3e-14, max_iter=10): theta = _get_theta(rhs) z = sqrt(anu) * np.cos(theta) - L = zeros((3, len(z))) - k0 = 0 - kp1 = 1 + p = zeros((3, len(z))) + k_0 = 0 + k_p1 = 1 - for i in range(max_iter): # @UnusedVariable + for _i in range(max_iter): # Newtons method carried out simultaneously on the roots. - L[k0, :] = 0 - L[kp1, :] = PIM4 + p[k_0, :] = 0 + p[k_p1, :] = pim4 for j in range(1, n + 1): # Loop up the recurrence relation to get the Hermite # polynomials evaluated at z. - km1 = k0 - k0 = kp1 - kp1 = np.mod(kp1 + 1, 3) + k_m1 = k_0 + k_0 = k_p1 + k_p1 = np.mod(k_p1 + 1, 3) - L[kp1, :] = (z * sqrt(2 / j) * L[k0, :] - - np.sqrt((j - 1) / j) * L[km1, :]) + p[k_p1, :] = (z * sqrt(2 / j) * p[k_0, :] - + sqrt((j - 1) / j) * p[k_m1, :]) - # L now contains the desired Hermite polynomials. - # We next compute pp, the derivatives, + # p now contains the desired Hermite polynomials. + # We next compute p_deriv, the derivatives, # by the relation (4.5.21) using p2, the polynomials # of one lower order. - pp = sqrt(2 * n) * L[k0, :] - dz = L[kp1, :] / pp + p_deriv = sqrt(2 * n) * p[k_0, :] + d_z = p[k_p1, :] / p_deriv - z = z - dz # Newtons formula. + z = z - d_z # Newtons formula. - converged = not np.any(abs(dz) > releps) + converged = not np.any(abs(d_z) > releps) if converged: break _assert_warn(converged, 'Newton iteration did not converge!') - w = 2. / pp ** 2 - return _expand_roots(z, w, n, m) -# x = np.empty(n) -# w = np.empty(n) -# x[0:m] = z # Store the root -# x[n - 1:n - m - 1:-1] = -z # and its symmetric counterpart. -# w[0:m] = 2. / pp ** 2 # Compute the weight -# w[n - 1:n - m - 1:-1] = w[0:m] # and its symmetric counterpart. -# return x, w + weights = 2. / p_deriv ** 2 + return _expand_roots(z, weights, n, m) def h_roots(n, method='newton'): @@ -337,48 +329,50 @@ def _j_roots_newton(n, alpha, beta, releps=3e-14, max_iter=10): z = np.cos(np.pi * (np.arange(1, n + 1) - 0.25 + 0.5 * alpha) / (n + 0.5 * (alfbet + 1))) - L = zeros((3, len(z))) - k0 = 0 - kp1 = 1 - for i in range(max_iter): + p = zeros((3, len(z))) + k_0 = 0 + k_p1 = 1 + for _i in range(max_iter): # Newton's method carried out simultaneously on the roots. tmp = 2 + alfbet - L[k0, :] = 1 - L[kp1, :] = (alpha - beta + tmp * z) / 2 + p[k_0, :] = 1 + p[k_p1, :] = (alpha - beta + tmp * z) / 2 for j in range(2, n + 1): # Loop up the recurrence relation to get the Jacobi # polynomials evaluated at z. - km1 = k0 - k0 = kp1 - kp1 = np.mod(kp1 + 1, 3) + k_m1 = k_0 + k_0 = k_p1 + k_p1 = np.mod(k_p1 + 1, 3) a = 2. * j * (j + alfbet) * tmp tmp = tmp + 2 c = 2 * (j - 1 + alpha) * (j - 1 + beta) * tmp b = (tmp - 1) * (alpha ** 2 - beta ** 2 + tmp * (tmp - 2) * z) - L[kp1, :] = (b * L[k0, :] - c * L[km1, :]) / a + p[k_p1, :] = (b * p[k_0, :] - c * p[k_m1, :]) / a - # L now contains the desired Jacobi polynomials. - # We next compute pp, the derivatives with a standard + # p now contains the desired Jacobi polynomials. + # We next compute p_deriv, the derivatives with a standard # relation involving the polynomials of one lower order. - pp = ((n * (alpha - beta - tmp * z) * L[kp1, :] + - 2 * (n + alpha) * (n + beta) * L[k0, :]) / (tmp * (1 - z ** 2))) - dz = L[kp1, :] / pp - z = z - dz # Newton's formula. + p_deriv = ((n * (alpha - beta - tmp * z) * p[k_p1, :] + + 2 * (n + alpha) * (n + beta) * p[k_0, :]) / + (tmp * (1 - z ** 2))) + d_z = p[k_p1, :] / p_deriv + z = z - d_z # Newton's formula. - if not any(abs(dz) > releps * abs(z)): + converged = not any(abs(d_z) > releps * abs(z)) + if converged: break - else: - warnings.warn('too many iterations in jrule') + + _assert_warn(converged, 'too many iterations in jrule') x = z # Store the root and the weight. f = (sp.gammaln(alpha + n) + sp.gammaln(beta + n) - sp.gammaln(n + 1) - sp.gammaln(alpha + beta + n + 1)) - w = (np.exp(f) * tmp * 2 ** alfbet / (pp * L[k0, :])) - return x, w + weights = (np.exp(f) * tmp * 2 ** alfbet / (p_deriv * p[k_0, :])) + return x, weights def j_roots(n, alpha, beta, method='newton'): @@ -483,10 +477,10 @@ def la_roots(n, alpha=0, method='newton'): def _get_theta(rhs): - r3 = rhs ** (1. / 3) - r2 = r3 ** 2 - C = [9.084064e-01, 5.214976e-02, 2.579930e-03, 3.986126e-03] - theta = r3 * (C[0] + r2 * (C[1] + r2 * (C[2] + r2 * C[3]))) + r_3 = rhs ** (1. / 3) + r_2 = r_3 ** 2 + c = [9.084064e-01, 5.214976e-02, 2.579930e-03, 3.986126e-03] + theta = r_3 * (c[0] + r_2 * (c[1] + r_2 * (c[2] + r_2 * c[3]))) return theta @@ -498,57 +492,57 @@ def _la_roots_newton(n, alpha, releps=3e-14, max_iter=10): theta = _get_theta(rhs) z = anu * np.cos(theta) ** 2 - dz = zeros(len(z)) - L = zeros((3, len(z))) - Lp = zeros((1, len(z))) - pp = zeros((1, len(z))) - k0 = 0 - kp1 = 1 + d_z = zeros(len(z)) + p = zeros((3, len(z))) + p_previous = zeros((1, len(z))) + p_deriv = zeros((1, len(z))) + k_0 = 0 + k_p1 = 1 k = slice(len(z)) - for _its in range(max_iter): + for _i in range(max_iter): # Newton's method carried out simultaneously on the roots. - L[k0, k] = 0. - L[kp1, k] = 1. + p[k_0, k] = 0. + p[k_p1, k] = 1. - for jj in range(1, n + 1): + for j in range(1, n + 1): # Loop up the recurrence relation to get the Laguerre # polynomials evaluated at z. - km1 = k0 - k0 = kp1 - kp1 = np.mod(kp1 + 1, 3) + km1 = k_0 + k_0 = k_p1 + k_p1 = np.mod(k_p1 + 1, 3) - L[kp1, k] = ((2 * jj - 1 + alpha - z[k]) * L[k0, k] - - (jj - 1 + alpha) * L[km1, k]) / jj + p[k_p1, k] = ((2 * j - 1 + alpha - z[k]) * p[k_0, k] - + (j - 1 + alpha) * p[km1, k]) / j # end - # L now contains the desired Laguerre polynomials. - # We next compute pp, the derivatives with a standard - # relation involving the polynomials of one lower order. + # p now contains the desired Laguerre polynomials. + # We next compute p_deriv, the derivatives with a standard + # relation involving the polynomials of one lower order. - Lp[k] = L[k0, k] - pp[k] = (n * L[kp1, k] - (n + alpha) * Lp[k]) / z[k] + p_previous[k] = p[k_0, k] + p_deriv[k] = (n * p[k_p1, k] - (n + alpha) * p_previous[k]) / z[k] - dz[k] = L[kp1, k] / pp[k] - z[k] = z[k] - dz[k] # % Newton?s formula. - # k = find((abs(dz) > releps.*z)) + d_z[k] = p[k_p1, k] / p_deriv[k] + z[k] = z[k] - d_z[k] # Newton?s formula. + # k = find((abs(d_z) > releps.*z)) - converged = not np.any(abs(dz) > releps) + converged = not np.any(abs(d_z) > releps) if converged: break - _assert_warn('too many iterations!') + _assert_warn(converged, 'too many iterations!') - x = z - w = -np.exp(sp.gammaln(alpha + n) - sp.gammaln(n)) / (pp * n * Lp) - return x, w + nodes = z + weights = -np.exp(sp.gammaln(alpha + n) - + sp.gammaln(n)) / (p_deriv * n * p_previous) + return nodes, weights def _p_roots_newton_start(n): m = int(np.fix((n + 1) / 2)) - mm = 4 * m - 1 - t = (np.pi / (4 * n + 2)) * np.arange(3, mm + 1, 4) - nn = 1 - (1 - 1 / n) / (8 * n * n) - xo = nn * np.cos(t) - return m, xo + t = (np.pi / (4 * n + 2)) * np.arange(3, 4 * m, 4) + a = 1 - (1 - 1 / n) / (8 * n * n) + x = a * np.cos(t) + return m, x def _p_roots_newton(n, ): @@ -556,52 +550,52 @@ def _p_roots_newton(n, ): Algorithm given by Davis and Rabinowitz in 'Methods of Numerical Integration', page 365, Academic Press, 1975. """ - m, xo = _p_roots_newton_start(n) + m, x = _p_roots_newton_start(n) - e1 = n * (n + 1) + e_1 = n * (n + 1) for _j in range(2): - pkm1 = 1 - pk = xo + p_km1 = 1 + p_k = x for k in range(2, n + 1): - t1 = xo * pk - pkp1 = t1 - pkm1 - (t1 - pkm1) / k + t1 - pkm1 = pk - pk = pkp1 - - den = 1. - xo * xo - d1 = n * (pkm1 - xo * pk) - dpn = d1 / den - d2pn = (2. * xo * dpn - e1 * pk) / den - d3pn = (4. * xo * d2pn + (2 - e1) * dpn) / den - d4pn = (6. * xo * d3pn + (6 - e1) * d2pn) / den - u = pk / dpn - v = d2pn / dpn - h = -u * (1 + (.5 * u) * (v + u * (v * v - u * d3pn / (3 * dpn)))) - p = pk + h * (dpn + (.5 * h) * (d2pn + - (h / 3) * (d3pn + .25 * h * d4pn))) - dp = dpn + h * (d2pn + (.5 * h) * (d3pn + h * d4pn / 3)) - h = h - p / dp - xo = xo + h - - x = -xo - h - fx = d1 - h * e1 * (pk + (h / 2) * (dpn + (h / 3) * - (d2pn + (h / 4) * - (d3pn + (.2 * h) * d4pn)))) - w = 2 * (1 - x ** 2) / (fx ** 2) - return _expand_roots(x, w, n, m) + t_1 = x * p_k + p_kp1 = t_1 - p_km1 - (t_1 - p_km1) / k + t_1 + p_km1 = p_k + p_k = p_kp1 + + den = 1. - x * x + d_1 = n * (p_km1 - x * p_k) + d_pn = d_1 / den + d_2pn = (2. * x * d_pn - e_1 * p_k) / den + d_3pn = (4. * x * d_2pn + (2 - e_1) * d_pn) / den + d_4pn = (6. * x * d_3pn + (6 - e_1) * d_2pn) / den + u = p_k / d_pn + v = d_2pn / d_pn + h = -u * (1 + (.5 * u) * (v + u * (v * v - u * d_3pn / (3 * d_pn)))) + p = p_k + h * (d_pn + (.5 * h) * (d_2pn + (h / 3) * + (d_3pn + .25 * h * d_4pn))) + d_p = d_pn + h * (d_2pn + (.5 * h) * (d_3pn + h * d_4pn / 3)) + h = h - p / d_p + x = x + h + + nodes = -x - h + f_x = d_1 - h * e_1 * (p_k + (h / 2) * (d_pn + (h / 3) * + (d_2pn + (h / 4) * + (d_3pn + (.2 * h) * d_4pn)))) + weights = 2 * (1 - nodes ** 2) / (f_x ** 2) + return _expand_roots(nodes, weights, n, m) def _p_roots_newton1(n, releps=1e-15, max_iter=100): - m, xo = _p_roots_newton_start(n) + m, x = _p_roots_newton_start(n) # Compute the zeros of the N+1 Legendre Polynomial # using the recursion relation and the Newton-Raphson method # Legendre-Gauss Polynomials - L = zeros((3, m)) + p = zeros((3, m)) # Derivative of LGP - Lp = zeros((m,)) - dx = zeros((m,)) + p_deriv = zeros((m,)) + d_x = zeros((m,)) # Compute the zeros of the N+1 Legendre Polynomial # using the recursion relation and the Newton-Raphson method @@ -609,32 +603,33 @@ def _p_roots_newton1(n, releps=1e-15, max_iter=100): # Iterate until new points are uniformly within epsilon of old # points k = slice(m) - k0 = 0 - kp1 = 1 + k_0 = 0 + k_p1 = 1 for _ix in range(max_iter): - L[k0, k] = 1 - L[kp1, k] = xo[k] - - for jj in range(2, n + 1): - km1 = k0 - k0 = kp1 - kp1 = np.mod(k0 + 1, 3) - L[kp1, k] = ((2 * jj - 1) * xo[k] * L[k0, k] - - (jj - 1) * L[km1, k]) / jj - - Lp[k] = n * (L[k0, k] - xo[k] * L[kp1, k]) / (1 - xo[k] ** 2) - - dx[k] = L[kp1, k] / Lp[k] - xo[k] = xo[k] - dx[k] - k, = np.nonzero((abs(dx) > releps * np.abs(xo))) - if len(k) == 0: + p[k_0, k] = 1 + p[k_p1, k] = x[k] + + for j in range(2, n + 1): + k_m1 = k_0 + k_0 = k_p1 + k_p1 = np.mod(k_0 + 1, 3) + p[k_p1, k] = ((2 * j - 1) * x[k] * p[k_0, k] - + (j - 1) * p[k_m1, k]) / j + + p_deriv[k] = n * (p[k_0, k] - x[k] * p[k_p1, k]) / (1 - x[k] ** 2) + + d_x[k] = p[k_p1, k] / p_deriv[k] + x[k] = x[k] - d_x[k] + k, = np.nonzero((abs(d_x) > releps * np.abs(x))) + converged = len(k) == 0 + if converged: break - else: - warnings.warn('Too many iterations!') - x = -xo - w = 2. / ((1 - x ** 2) * (Lp ** 2)) - return _expand_roots(x, w, n, m) + _assert(converged, 'Too many iterations!') + + nodes = -x + weights = 2. / ((1 - nodes ** 2) * (p_deriv ** 2)) + return _expand_roots(nodes, weights, n, m) def _expand_roots(x, w, n, m): @@ -650,7 +645,7 @@ def _expand_roots(x, w, n, m): def p_roots(n, method='newton', a=-1, b=1): """ Returns the roots (x) of the nth order Legendre polynomial, P_n(x), - and weights (w) to use in Gaussian Quadrature over [-1,1] with weighting + and weights to use in Gaussian Quadrature over [-1,1] with weighting function 1. Parameters @@ -664,9 +659,9 @@ def p_roots(n, method='newton', a=-1, b=1): Returns ------- - x : ndarray + nodes : ndarray roots - w : ndarray + weights : ndarray weights @@ -674,8 +669,8 @@ def p_roots(n, method='newton', a=-1, b=1): ------- Integral of exp(x) from a = 0 to b = 3 is: exp(3)-exp(0)= >>> import numpy as np - >>> x, w = p_roots(11, a=0, b=3) - >>> np.allclose(np.sum(np.exp(x)*w), 19.085536923187668) + >>> nodes, weights = p_roots(11, a=0, b=3) + >>> np.allclose(np.sum(np.exp(nodes) * weights), 19.085536923187668) True See also @@ -697,74 +692,74 @@ def p_roots(n, method='newton', a=-1, b=1): """ if not method.startswith('n'): - x, w = ort.p_roots(n) + nodes, weights = ort.p_roots(n) else: if method.endswith('1'): - x, w = _p_roots_newton1(n) + nodes, weights = _p_roots_newton1(n) else: - x, w = _p_roots_newton(n) + nodes, weights = _p_roots_newton(n) if (a != -1) | (b != 1): # Linear map from[-1,1] to [a,b] - dh = (b - a) / 2 - x = dh * (x + 1) + a - w = w * dh + d_h = (b - a) / 2 + nodes = d_h * (nodes + 1) + a + weights = weights * d_h - return x, w + return nodes, weights def q5_roots(n): """ 5 : p(x) = 1/sqrt((x-a)*(b-x)), a =-1, b = 1 Chebyshev 1'st kind """ - jj = np.arange(1, n + 1) - wf = ones(n) * np.pi / n - bp = np.cos((2 * jj - 1) * np.pi / (2 * n)) - return bp, wf + j = np.arange(1, n + 1) + weights = ones(n) * np.pi / n + nodes = np.cos((2 * j - 1) * np.pi / (2 * n)) + return nodes, weights def q6_roots(n): """ 6 : p(x) = sqrt((x-a)*(b-x)), a =-1, b = 1 Chebyshev 2'nd kind """ - jj = np.arange(1, n + 1) - xj = jj * np.pi / (n + 1) - wf = np.pi / (n + 1) * np.sin(xj) ** 2 - bp = np.cos(xj) - return bp, wf + j = np.arange(1, n + 1) + x_j = j * np.pi / (n + 1) + weights = np.pi / (n + 1) * np.sin(x_j) ** 2 + nodes = np.cos(x_j) + return nodes, weights def q7_roots(n): """ 7 : p(x) = sqrt((x-a)/(b-x)), a = 0, b = 1 """ - jj = np.arange(1, n + 1) - xj = (jj - 0.5) * pi / (2 * n + 1) - bp = np.cos(xj) ** 2 - wf = 2 * np.pi * bp / (2 * n + 1) - return bp, wf + j = np.arange(1, n + 1) + x_j = (j - 0.5) * pi / (2 * n + 1) + nodes = np.cos(x_j) ** 2 + weights = 2 * np.pi * nodes / (2 * n + 1) + return nodes, weights def q8_roots(n): """ 8 : p(x) = 1/sqrt(b-x), a = 0, b = 1 """ - bp1, wf1 = p_roots(2 * n) - k, = np.where(0 <= bp1) - wf = 2 * wf1[k] - bp = 1 - bp1[k] ** 2 - return bp, wf + nodes_1, weights_1 = p_roots(2 * n) + k, = np.where(0 <= nodes_1) + weights = 2 * weights_1[k] + nodes = 1 - nodes_1[k] ** 2 + return nodes, weights def q9_roots(n): """ 9 : p(x) = sqrt(b-x), a = 0, b = 1 """ - bp1, wf1 = p_roots(2 * n + 1) - k, = np.where(0 < bp1) - wf = 2 * bp1[k] ** 2 * wf1[k] - bp = 1 - bp1[k] ** 2 - return bp, wf + nodes_1, weights_1 = p_roots(2 * n + 1) + k, = np.where(0 < nodes_1) + weights = 2 * nodes_1[k] ** 2 * weights_1[k] + nodes = 1 - nodes_1[k] ** 2 + return nodes, weights def qrule(n, wfun=1, alpha=0, beta=0): @@ -805,18 +800,18 @@ def qrule(n, wfun=1, alpha=0, beta=0): >>> import numpy as np # integral of x^2 from a = -1 to b = 1 - >>> [bp,wf] = qrule(10) + >>> bp, wf = qrule(10) >>> np.allclose(sum(bp**2*wf), 0.66666666666666641) True - # integral of exp(-x.^2)*x.^2 from a = -inf to b = inf - >>> [bp,wf] = qrule(10,2) - >>> np.allclose(sum(bp**2*wf), 0.88622692545275772) + # integral of exp(-x**2)*x**2 from a = -inf to b = inf + >>> bp, wf = qrule(10,2) + >>> np.allclose(sum(bp ** 2 * wf), 0.88622692545275772) True - # integral of (x+1)*(1-x)^2 from a = -1 to b = 1 - >>> [bp,wf] = qrule(10,4,1,2) - >>> np.allclose((bp*wf).sum(), 0.26666666666666755) + # integral of (x+1)*(1-x)**2 from a = -1 to b = 1 + >>> bp, wf = qrule(10,4,1,2) + >>> np.allclose((bp * wf).sum(), 0.26666666666666755) True See also @@ -928,22 +923,22 @@ class _Gaussq(object): def fun1(x): return x if wfun == 4: - dx = jacob ** (alpha + beta + 1) + d_x = jacob ** (alpha + beta + 1) else: - dx = [None, fun1, fun1, fun1, None, lambda x: ones(np.shape(x)), - lambda x: x ** 2, fun1, lambda x: sqrt(x), - lambda x: sqrt(x) ** 3][wfun](jacob) - return dx.ravel() + d_x = [None, fun1, fun1, fun1, None, lambda x: ones(np.shape(x)), + lambda x: x ** 2, fun1, sqrt, + lambda x: sqrt(x) ** 3][wfun](jacob) + return d_x.ravel() @staticmethod - def _points_and_weights(gn, wfun, alpha, beta): - global _POINTS_AND_WEIGHTS - name = 'wfun{:d}_{:d}_{:g}_{:g}'.format(wfun, gn, alpha, beta) - x_and_w = _POINTS_AND_WEIGHTS.setdefault(name, []) - if len(x_and_w) == 0: - x_and_w.extend(qrule(gn, wfun, alpha, beta)) - xn, w = x_and_w - return xn, w + def _nodes_and_weights(num_nodes, wfun, alpha, beta): + global _NODES_AND_WEIGHTS + name = 'wfun{:d}_{:d}_{:g}_{:g}'.format(wfun, num_nodes, alpha, beta) + nodes_and_weights = _NODES_AND_WEIGHTS.setdefault(name, []) + if len(nodes_and_weights) == 0: + nodes_and_weights.extend(qrule(num_nodes, wfun, alpha, beta)) + nodes, weights = nodes_and_weights + return nodes, weights def _initialize_trace(self, max_iter): if self.trace: @@ -963,24 +958,23 @@ class _Gaussq(object): plt.plot(np.hstack(self.x_trace), np.hstack(self.y_trace), '+') @staticmethod - def _get_jacob(wfun, A, B): + def _get_jacob(wfun, a, b): if wfun in [2, 3]: - nk = np.size(A) - jacob = ones((nk, 1)) + jacob = ones((np.size(a), 1)) else: - jacob = (B - A) * 0.5 + jacob = (b - a) * 0.5 if wfun in [7, 8, 9]: - jacob = jacob * 2 + jacob *= 2 return jacob @staticmethod def _warn(k, a_shape): - nk = len(k) - if nk > 1: - if nk == np.prod(a_shape): + n = len(k) + if n > 1: + if n == np.prod(a_shape): tmptxt = 'All integrals did not converge' else: - tmptxt = '%d integrals did not converge' % (nk, ) + tmptxt = '%d integrals did not converge' % (n, ) tmptxt = tmptxt + '--singularities likely!' else: tmptxt = 'Integral did not converge--singularity likely!' @@ -1000,42 +994,43 @@ class _Gaussq(object): def __call__(self, fun, a, b, releps=1e-3, abseps=1e-3, alpha=0, beta=0, wfun=1, trace=False, args=(), max_iter=11): self.trace = trace - gn = 2 + num_nodes = 2 - aa, bb, args, a_shape = self._initialize(wfun, a, b, args) + a_0, b_0, args, a_shape = self._initialize(wfun, a, b, args) - jacob = self._get_jacob(wfun, aa, bb) + jacob = self._get_jacob(wfun, a_0, b_0) shift = int(wfun in [1, 4, 5, 6]) - dx = self._get_dx(wfun, jacob, alpha, beta) + d_x = self._get_dx(wfun, jacob, alpha, beta) self._initialize_trace(max_iter) # Break out of the iteration loop for three reasons: # 1) the last update is very small (compared to int and to releps) # 2) There are more than 11 iterations. This should NEVER happen. - dtype = np.result_type(fun((aa+bb)*0.5, *args)) - nk = np.prod(a_shape) # # of integrals we have to compute - k = np.arange(nk) - opts = (nk, dtype) + dtype = np.result_type(fun((a_0+b_0)*0.5, *args)) + n_k = np.prod(a_shape) # # of integrals we have to compute + k = np.arange(n_k) + opts = (n_k, dtype) val, val_old, abserr = zeros(*opts), ones(*opts), zeros(*opts) + nodes_and_weights = self._nodes_and_weights for i in range(max_iter): - xn, w = self._points_and_weights(gn, wfun, alpha, beta) - x = (xn + shift) * jacob[k, :] + aa[k, :] + x_n, weights = nodes_and_weights(num_nodes, wfun, alpha, beta) + x = (x_n + shift) * jacob[k, :] + a_0[k, :] - pari = [xi[k, :] for xi in args] - y = fun(x, *pari) + params = [xi[k, :] for xi in args] + y = fun(x, *params) self._plot_trace(x, y) - val[k] = np.sum(w * y, axis=1) * dx[k] # do the integration + val[k] = np.sum(weights * y, axis=1) * d_x[k] # do the integration if any(np.isnan(val)): val[np.isnan(val)] = val_old[np.isnan(val)] if 1 < i: abserr[k] = abs(val_old[k] - val[k]) # absolute tolerance k, = np.where(abserr > np.maximum(abs(releps * val), abseps)) - nk = len(k) # of integrals we have to compute again - if nk == 0: + n_k = len(k) # of integrals we have to compute again + if n_k == 0: break val_old[k] = val[k] - gn *= 2 # double the # of basepoints and weights + num_nodes *= 2 # double the # of basepoints and weights else: self._warn(k, a_shape) @@ -1048,14 +1043,13 @@ class _Gaussq(object): gaussq = _Gaussq() -def richardson(Q, k): +def richardson(q_val, k): # license BSD # Richardson extrapolation with parameter estimation - c = np.real((Q[k - 1] - Q[k - 2]) / (Q[k] - Q[k - 1])) - 1. + c = np.real((q_val[k - 1] - q_val[k - 2]) / (q_val[k] - q_val[k - 1])) - 1. # The lower bound 0.07 admits the singularity x.^-0.9 c = max(c, 0.07) - R = Q[k] + (Q[k] - Q[k - 1]) / c - return R + return q_val[k] + (q_val[k] - q_val[k - 1]) / c class _Quadgr(object): @@ -1083,15 +1077,15 @@ class _Quadgr(object): @staticmethod def _nodes_and_weights(): # Gauss-Legendre quadrature (12-point) - xq = np.asarray( + x = np.asarray( [0.12523340851146894, 0.36783149899818018, 0.58731795428661748, 0.76990267419430469, 0.9041172563704748, 0.98156063424671924]) - wq = np.asarray( + w = np.asarray( [0.24914704581340288, 0.23349253653835478, 0.20316742672306584, 0.16007832854334636, 0.10693932599531818, 0.047175336386511842]) - xq = np.hstack((xq, -xq)) - wq = np.hstack((wq, wq)) - return xq, wq + nodes = np.hstack((x, -x)) + weights = np.hstack((w, w)) + return nodes, weights @staticmethod def _get_best_estimate(vals0, vals1, vals2, k): @@ -1121,24 +1115,23 @@ class _Quadgr(object): val1 = zeros(max_iter, dtype=dtype) # First Richardson extrapolation val2 = zeros(max_iter, dtype=dtype) # Second Richardson extrapolation - x0, w0 = self._nodes_and_weights() - nx0 = len(x0) + x_n, weights = self._nodes_and_weights() + n = len(x_n) # One interval - hh = (b - a) / 2 # Half interval length - x = (a + b) / 2 + hh * x0 # Nodes + d_x = (b - a) / 2 # Half interval length + x = (a + b) / 2 + d_x * x_n # Nodes # Quadrature - val0[0] = hh * np.sum(w0 * fun(x), axis=0) + val0[0] = d_x * np.sum(weights * fun(x), axis=0) # Successive bisection of intervals for k in range(1, max_iter): # Interval bisection - hh = hh / 2 + d_x = d_x / 2 x = np.hstack([x + a, x + b]) / 2 # Quadrature - val0[k] = hh * np.sum(w0 * np.sum(np.reshape(fun(x), (-1, nx0)), - axis=0), - axis=0) + val0[k] = np.sum(np.sum(np.reshape(fun(x), (-1, n)), axis=0) * + weights, axis=0) * d_x # Richardson extrapolation if k >= 5: @@ -1147,22 +1140,20 @@ class _Quadgr(object): elif k >= 3: val1[k] = richardson(val0, k) - Q, err = self._get_best_estimate(val0, val1, val2, k) - - # Convergence + q_val, err = self._get_best_estimate(val0, val1, val2, k) - converged = (err < abseps) | ~np.isfinite(Q) + converged = (err < abseps) | ~np.isfinite(q_val) if converged: break _assert_warn(converged, 'Max number of iterations reached without ' 'convergence.') - _assert_warn(np.isfinite(Q), + _assert_warn(np.isfinite(q_val), 'Integral approximation is Infinite or NaN.') # The error estimate should not be zero - err = err + 2 * np.finfo(Q).eps - return Q, err + err = err + 2 * np.finfo(q_val).eps + return q_val, err @staticmethod def _order_limits(a, b): @@ -1174,9 +1165,9 @@ class _Quadgr(object): """ Gauss-Legendre quadrature with Richardson extrapolation. - [Q,ERR] = QUADGR(FUN,A,B,TOL) approximates the integral of a function + [q_val,ERR] = QUADGR(FUN,A,B,TOL) approximates the integral of a function FUN from A to B with an absolute error tolerance TOL. FUN is a function - handle and must accept vector arguments. TOL is 1e-6 by default. Q is + handle and must accept vector arguments. TOL is 1e-6 by default. q_val is the integral approximation and ERR is an estimate of the absolute error. @@ -1188,7 +1179,7 @@ class _Quadgr(object): Examples -------- >>> import numpy as np - >>> Q, err = quadgr(np.log,0,1) + >>> q_val, err = quadgr(np.log,0,1) >>> q, err = quadgr(np.exp,0,9999*1j*np.pi) >>> np.allclose(q, -2.0000000000122662), err < 1.0e-08 (True, True) @@ -1223,25 +1214,24 @@ class _Quadgr(object): a = np.asarray(a) b = np.asarray(b) if a == b: - Q = b - a + q_val = b - a err = b - a - return Q, err + return q_val, err a, b, reverse = self._order_limits(a, b) - # Infinite limits improper_integral = np.isinf(a) | np.isinf(b) - if improper_integral: - Q, err = self._change_variable_and_integrate(fun, a, b, abseps, - max_iter) + if improper_integral: # Infinite limits + q_val, err = self._change_variable_and_integrate(fun, a, b, abseps, + max_iter) else: - Q, err = self._integrate(fun, a, b, abseps, max_iter) + q_val, err = self._integrate(fun, a, b, abseps, max_iter) # Reverse direction if reverse: - Q = -Q + q_val = -q_val - return Q, err + return q_val, err quadgr = _Quadgr() @@ -1380,13 +1370,13 @@ def qdemo(f, a, b, kmax=9, plot_error=False): err_dic.setdefault(name, []).append(abs(q[0] - true_val)) name = 'Chebychev' - ck = np.polynomial.chebyshev.chebfit(x, y, deg=min(n-1, 36)) - cki = np.polynomial.chebyshev.chebint(ck) - q = np.polynomial.chebyshev.chebval(x[-1], cki) + c_k = np.polynomial.chebyshev.chebfit(x, y, deg=min(n-1, 36)) + c_ki = np.polynomial.chebyshev.chebint(c_k) + q = np.polynomial.chebyshev.chebval(x[-1], c_ki) vals_dic.setdefault(name, []).append(q) err_dic.setdefault(name, []).append(abs(q - true_val)) - # ck = chebfit(f,n,a,b) - # q = chebval(b,chebint(ck,a,b),a,b) + # c_k = chebfit(f,n,a,b) + # q = chebval(b,chebint(c_k,a,b),a,b) # qc2[k] = q; ec2[k] = abs(q - true) name = 'Gauss-Legendre' # quadrature