模型可视化>> quadpy>> 返回
项目作者： nschloe

项目描述 ：
Numerical integration (quadrature, cubature) in Python
高级语言： Python
项目地址: git://github.com/nschloe/quadpy.git
创建时间： 2016-10-03T11:35:39Z
项目社区：https://github.com/nschloe/quadpy
开源协议：GNU General Public License v3.0
下载



quadpy

Your one-stop shop for numerical integration in Python.


PyPi Version
PyPI pyversions
GitHub stars
PyPi downloads

Discord
awesome

More than 1500 numerical integration schemes for
line segments,
circles,
disks,
triangles,
quadrilaterals,
spheres,
balls,
tetrahedra,
hexahedra,
wedges,
pyramids,
n-spheres,
n-balls,
n-cubes,
n-simplices,
the 1D half-space with weight functions exp(-r),
the 2D space with weight functions exp(-r),
the 3D space with weight functions exp(-r),
the nD space with weight functions exp(-r),
the 1D space with weight functions exp(-r2),
the 2D space with weight functions exp(-r2),
the 3D space with weight functions exp(-r2),
and
the nD space with weight functions exp(-r2),
for fast integration of real-, complex-, and vector-valued functions.

Installation

Install quadpy from PyPI with

  1. pip install quadpy

See here on how to get a license.

Using quadpy

Quadpy provides integration schemes for many different 1D, 2D, even nD domains.

To start off easy: If you’d numerically integrate any function over any given
1D interval, do

  1. import numpy as np
  2. import quadpy
  5. def f(x):
  6.     return np.sin(x) - x
  9. val, err = quadpy.quad(f, 0.0, 6.0)

This is like
scipy
with the addition that quadpy handles complex-, vector-, matrix-valued integrands,
and “intervals” in spaces of arbitrary dimension.

To integrate over a triangle, do

  1. import numpy as np
  2. import quadpy
  5. def f(x):
  6.     return np.sin(x[0]) * np.sin(x[1])
  9. triangle = np.array([[0.0, 0.0], [1.0, 0.0], [0.7, 0.5]])
  11. # get a "good" scheme of degree 10
  12. scheme = quadpy.t2.get_good_scheme(10)
  13. val = scheme.integrate(f, triangle)

Most domains have get_good_scheme(degree). If you would like to use a
particular scheme, you can pick one from the dictionary quadpy.t2.schemes.

All schemes have

  1. scheme.points
  2. scheme.weights
  3. scheme.degree
  4. scheme.source
  5. scheme.test_tolerance
  7. scheme.show()
  8. scheme.integrate(
  9.     # ...
  10. )

and many have

  1. scheme.points_symbolic
  2. scheme.weights_symbolic

You can explore schemes on the command line with, e.g.,

  1. quadpy info s2 rabinowitz_richter_3
  1. <quadrature scheme for S2>
  2.   name:                 Rabinowitz-Richter 2
  3.   source:               Perfectly Symmetric Two-Dimensional Integration Formulas with Minimal Numbers of Points
  4.                         Philip Rabinowitz, Nira Richter
  5.                         Mathematics of Computation, vol. 23, no. 108, pp. 765-779, 1969
  6.                         https://doi.org/10.1090/S0025-5718-1969-0258281-4
  7.   degree:               9
  8.   num points/weights:   21
  9.   max/min weight ratio: 7.632e+01
  10.   test tolerance:       9.417e-15
  11.   point position:       outside
  12.   all weights positive: True

Also try quadpy show!

quadpy is fully vectorized, so if you like to compute the integral of a function on many
domains at once, you can provide them all in one integrate() call, e.g.,

  1. # shape (3, 5, 2), i.e., (corners, num_triangles, xy_coords)
  2. triangles = np.stack(
  3.     [
  4.         [[0.0, 0.0], [1.0, 0.0], [0.0, 1.0]],
  5.         [[1.2, 0.6], [1.3, 0.7], [1.4, 0.8]],
  6.         [[26.0, 31.0], [24.0, 27.0], [33.0, 28]],
  7.         [[0.1, 0.3], [0.4, 0.4], [0.7, 0.1]],
  8.         [[8.6, 6.0], [9.4, 5.6], [7.5, 7.4]],
  9.     ],
  10.     axis=-2,
  11. )

The same goes for functions with vectorized output, e.g.,

  1. def f(x):
  2.     return [np.sin(x[0]), np.sin(x[1])]

More examples under test/examples_test.py.

Read more about the dimensionality of the input/output arrays in the
wiki.

Advanced topics:

Schemes

Line segment (C1)

  • Chebyshev-Gauss (type 1 and 2, arbitrary degree)
  • Clenshaw-Curtis (arbitrary degree)
  • Fejér (type 1 and 2, arbitrary degree)
  • Gauss-Jacobi (arbitrary degree)
  • Gauss-Legendre (arbitrary degree)
  • Gauss-Lobatto (arbitrary degree)
  • Gauss-Kronrod (arbitrary degree)
  • Gauss-Patterson (9 nested schemes up to degree 767)
  • Gauss-Radau (arbitrary degree)
  • Newton-Cotes (open and closed, arbitrary degree)

See
here
for how to generate Gauss formulas for your own weight functions.

Example:

  1. import numpy as np
  2. import quadpy
  4. scheme = quadpy.c1.gauss_patterson(5)
  5. scheme.show()
  6. val = scheme.integrate(lambda x: np.exp(x), [0.0, 1.0])

1D half-space with weight function exp(-r) (E1r)

Example:

  1. import quadpy
  3. scheme = quadpy.e1r.gauss_laguerre(5, alpha=0)
  4. scheme.show()
  5. val = scheme.integrate(lambda x: x**2)

1D space with weight function exp(-r2) (E1r2)

  • Gauss-Hermite (arbitrary degree)
  • Genz-Keister (1996, 8 nested schemes up to degree 67)

Example:

  1. import quadpy
  3. scheme = quadpy.e1r2.gauss_hermite(5)
  4. scheme.show()
  5. val = scheme.integrate(lambda x: x**2)

Circle (U2)

  • Krylov (1959, arbitrary degree)

Example:

  1. import numpy as np
  2. import quadpy
  4. scheme = quadpy.u2.get_good_scheme(7)
  5. scheme.show()
  6. val = scheme.integrate(lambda x: np.exp(x[0]), [0.0, 0.0], 1.0)

Triangle (T2)

Apart from the classical centroid, vertex, and seven-point schemes we have

Example:

  1. import numpy as np
  2. import quadpy
  4. scheme = quadpy.t2.get_good_scheme(12)
  5. scheme.show()
  6. val = scheme.integrate(lambda x: np.exp(x[0]), [[0.0, 0.0], [1.0, 0.0], [0.5, 0.7]])

Disk (S2)

  • Radon (1948, degree 5)
  • Peirce (1956, 3 schemes up to degree 11)
  • Peirce (1957, arbitrary degree)
  • Albrecht-Collatz (1958, degree 3)
  • Hammer-Stroud (1958, 8 schemes up to degree 15)
  • Albrecht (1960, 8 schemes up to degree 17)
  • Mysovskih (1964, 3 schemes up to degree 15)
  • Rabinowitz-Richter (1969, 6 schemes up to degree 15)
  • Lether (1971, arbitrary degree)
  • Piessens-Haegemans (1975, 1 scheme of degree 9)
  • Haegemans-Piessens (1977, degree 9)
  • Cools-Haegemans (1985, 4 schemes up to degree 13)
  • Wissmann-Becker (1986, 3 schemes up to degree 8)
  • Kim-Song (1997, 15 schemes up to degree 17)
  • Cools-Kim (2000, 3 schemes up to degree 21)
  • Luo-Meng (2007, 6 schemes up to degree 17)
  • Takaki-Forbes-Rolland (2022, 19 schemes up to degree 77)
  • all schemes from the n-ball

Example:

  1. import numpy as np
  2. import quadpy
  4. scheme = quadpy.s2.get_good_scheme(6)
  5. scheme.show()
  6. val = scheme.integrate(lambda x: np.exp(x[0]), [0.0, 0.0], 1.0)

Quadrilateral (C2)

Example:

  1. import numpy as np
  2. import quadpy
  4. scheme = quadpy.c2.get_good_scheme(7)
  5. val = scheme.integrate(
  6.     lambda x: np.exp(x[0]),
  7.     [[[0.0, 0.0], [1.0, 0.0]], [[0.0, 1.0], [1.0, 1.0]]],
  8. )

The points are specified in an array of shape (2, 2, …) such that arr[0][0]
is the lower left corner, arr[1][1] the upper right. If your c2
has its sides aligned with the coordinate axes, you can use the convenience
function

  1. quadpy.c2.rectangle_points([x0, x1], [y0, y1])

to generate the array.

2D space with weight function exp(-r) (E2r)

Example:

  1. import quadpy
  3. scheme = quadpy.e2r.get_good_scheme(5)
  4. scheme.show()
  5. val = scheme.integrate(lambda x: x[0] ** 2)

2D space with weight function exp(-r2) (E2r2)

Example:

  1. import quadpy
  3. scheme = quadpy.e2r2.get_good_scheme(3)
  4. scheme.show()
  5. val = scheme.integrate(lambda x: x[0] ** 2)

Sphere (U3)

Example:

  1. import numpy as np
  2. import quadpy
  4. scheme = quadpy.u3.get_good_scheme(19)
  5. # scheme.show()
  6. val = scheme.integrate(lambda x: np.exp(x[0]), [0.0, 0.0, 0.0], 1.0)

Integration on the sphere can also be done for functions defined in spherical
coordinates:

  1. import numpy as np
  2. import quadpy
  5. def f(theta_phi):
  6.     theta, phi = theta_phi
  7.     return np.sin(phi) ** 2 * np.sin(theta)
  10. scheme = quadpy.u3.get_good_scheme(19)
  11. val = scheme.integrate_spherical(f)

Ball (S3)

Example:

  1. import numpy as np
  2. import quadpy
  4. scheme = quadpy.s3.get_good_scheme(4)
  5. # scheme.show()
  6. val = scheme.integrate(lambda x: np.exp(x[0]), [0.0, 0.0, 0.0], 1.0)

Tetrahedron (T3)

Example:

  1. import numpy as np
  2. import quadpy
  4. scheme = quadpy.t3.get_good_scheme(5)
  5. # scheme.show()
  6. val = scheme.integrate(
  7.     lambda x: np.exp(x[0]),
  8.     [[0.0, 0.0, 0.0], [1.0, 0.0, 0.0], [0.5, 0.7, 0.0], [0.3, 0.9, 1.0]],
  9. )

Hexahedron (C3)

Example:

  1. import numpy as np
  2. import quadpy
  4. scheme = quadpy.c3.product(quadpy.c1.newton_cotes_closed(3))
  5. # scheme.show()
  6. val = scheme.integrate(
  7.     lambda x: np.exp(x[0]),
  8.     quadpy.c3.cube_points([0.0, 1.0], [-0.3, 0.4], [1.0, 2.1]),
  9. )

Pyramid (P3)

  • Felippa (2004, 9 schemes up to degree 5)

Example:

  1. import numpy as np
  2. import quadpy
  4. scheme = quadpy.p3.felippa_5()
  6. val = scheme.integrate(
  7.     lambda x: np.exp(x[0]),
  8.     [
  9.         [0.0, 0.0, 0.0],
  10.         [1.0, 0.0, 0.0],
  11.         [0.5, 0.7, 0.0],
  12.         [0.3, 0.9, 0.0],
  13.         [0.0, 0.1, 1.0],
  14.     ],
  15. )

Wedge (W3)

Example:

  1. import numpy as np
  2. import quadpy
  4. scheme = quadpy.w3.felippa_3()
  5. val = scheme.integrate(
  6.     lambda x: np.exp(x[0]),
  7.     [
  8.         [[0.0, 0.0, 0.0], [1.0, 0.0, 0.0], [0.5, 0.7, 0.0]],
  9.         [[0.0, 0.0, 1.0], [1.0, 0.0, 1.0], [0.5, 0.7, 1.0]],
  10.     ],
  11. )

3D space with weight function exp(-r) (E3r)

Example:

  1. import quadpy
  3. scheme = quadpy.e3r.get_good_scheme(5)
  4. # scheme.show()
  5. val = scheme.integrate(lambda x: x[0] ** 2)

3D space with weight function exp(-r2) (E3r2)

Example:

  1. import quadpy
  3. scheme = quadpy.e3r2.get_good_scheme(6)
  4. # scheme.show()
  5. val = scheme.integrate(lambda x: x[0] ** 2)

n-Simplex (Tn)

Example:

  1. import numpy as np
  2. import quadpy
  4. dim = 4
  5. scheme = quadpy.tn.grundmann_moeller(dim, 3)
  6. val = scheme.integrate(
  7.     lambda x: np.exp(x[0]),
  8.     np.array(
  9.         [
  10.             [0.0, 0.0, 0.0, 0.0],
  11.             [1.0, 2.0, 0.0, 0.0],
  12.             [0.0, 1.0, 0.0, 0.0],
  13.             [0.0, 3.0, 1.0, 0.0],
  14.             [0.0, 0.0, 4.0, 1.0],
  15.         ]
  16.     ),
  17. )

n-Sphere (Un)

Example:

  1. import numpy as np
  2. import quadpy
  4. dim = 4
  5. scheme = quadpy.un.dobrodeev_1978(dim)
  6. val = scheme.integrate(lambda x: np.exp(x[0]), np.zeros(dim), 1.0)

n-Ball (Sn)

Example:

  1. import numpy as np
  2. import quadpy
  4. dim = 4
  5. scheme = quadpy.sn.dobrodeev_1970(dim)
  6. val = scheme.integrate(lambda x: np.exp(x[0]), np.zeros(dim), 1.0)

n-Cube (Cn)

Example:

  1. import numpy as np
  2. import quadpy
  4. dim = 4
  5. scheme = quadpy.cn.stroud_cn_3_3(dim)
  6. val = scheme.integrate(
  7.     lambda x: np.exp(x[0]),
  8.     quadpy.cn.ncube_points([0.0, 1.0], [0.1, 0.9], [-1.0, 1.0], [-1.0, -0.5]),
  9. )

nD space with weight function exp(-r) (Enr)

Example:

  1. import quadpy
  3. dim = 4
  4. scheme = quadpy.enr.stroud_enr_5_4(dim)
  5. val = scheme.integrate(lambda x: x[0] ** 2)

nD space with weight function exp(-r2) (Enr2)

Example:

  1. import quadpy
  3. dim = 4
  4. scheme = quadpy.enr2.stroud_enr2_5_2(dim)
  5. val = scheme.integrate(lambda x: x[0] ** 2)