积分

符号积分

积分与求导的关系：

\frac{d}{dx} F(x) = f(x)\Rightarrow F(x) = \int f(x) dx

符号运算可以用 sympy 模块完成。

先导入 init_printing 模块方便其显示：

In [1]:

  1. from sympy import init_printing
  2. init_printing()

In [2]:

  1. from sympy import symbols, integrate
  2. import sympy

产生 x 和 y 两个符号变量，并进行运算：

In [3]:

  1. x, y = symbols('x y')
  2. sympy.sqrt(x ** 2 + y ** 2)

Out[3]:

\sqrt{x^{2} + y^{2}}

对于生成的符号变量 z，我们将其中的 x 利用 subs 方法替换为 3：

In [4]:

  1. z = sympy.sqrt(x ** 2 + y ** 2)
  2. z.subs(x, 3)

Out[4]:

\sqrt{y^{2} + 9}

再替换 y：

In [5]:

  1. z.subs(x, 3).subs(y, 4)

Out[5]:

5

还可以从 sympy.abc 中导入现成的符号变量：

In [6]:

  1. from sympy.abc import theta
  2. y = sympy.sin(theta) ** 2
  3. y

Out[6]:

\sin^{2}{\left (\theta \right )}

对 y 进行积分：

In [7]:

  1. Y = integrate(y)
  2. Y

Out[7]:

\frac{\theta}{2} - \frac{1}{2} \sin{\left (\theta \right )} \cos{\left (\theta \right )}

计算 $Y(\pi) - Y(0)$：

In [8]:

  1. import numpy as np
  2. np.set_printoptions(precision=3)
  3.  
  4. Y.subs(theta, np.pi) - Y.subs(theta, 0)

Out[8]:

1.5707963267949

计算 $\int_0^\pi y d\theta$ ：

In [9]:

  1. integrate(y, (theta, 0, sympy.pi))

Out[9]:

\frac{\pi}{2}

显示的是字符表达式，查看具体数值可以使用 evalf() 方法，或者传入 numpy.pi，而不是 sympy.pi ：

In [10]:

  1. integrate(y, (theta, 0, sympy.pi)).evalf()

Out[10]:

1.5707963267949

In [11]:

  1. integrate(y, (theta, 0, np.pi))

Out[11]:

1.5707963267949

根据牛顿莱布尼兹公式，这两个数值应该相等。

产生不定积分对象：

In [12]:

  1. Y_indef = sympy.Integral(y)
  2. Y_indef

Out[12]:

\int \sin^{2}{\left (\theta \right )}\, d\theta

In [13]:

  1. print type(Y_indef)
  1. <class 'sympy.integrals.integrals.Integral'>

定积分：

In [14]:

  1. Y_def = sympy.Integral(y, (theta, 0, sympy.pi))
  2. Y_def

Out[14]:

\int_{0}^{\pi} \sin^{2}{\left (\theta \right )}\, d\theta

产生函数 $Y(x) = \int_0^x sin^2(\theta) d\theta$，并将其向量化：

In [15]:

  1. Y_raw = lambda x: integrate(y, (theta, 0, x))
  2. Y = np.vectorize(Y_raw)

In [16]:

  1. %matplotlib inline
  2. import matplotlib.pyplot as plt
  3.  
  4. x = np.linspace(0, 2 * np.pi)
  5. p = plt.plot(x, Y(x))
  6. t = plt.title(r'$Y(x) = \int_0^x sin^2(\theta) d\theta$')

04.06 积分 - 图1

数值积分

数值积分：

F(x) = \lim{n \rightarrow \infty} \sum{i=0}^{n-1} f(xi)(x{i+1}-xi) \Rightarrow F(x) = \int{x_0}^{x_n} f(x) dx

导入贝塞尔函数：

In [17]:

  1. from scipy.special import jv

In [18]:

  1. def f(x):
  2.     return jv(2.5, x)

In [19]:

  1. x = np.linspace(0, 10)
  2. p = plt.plot(x, f(x), 'k-')

04.06 积分 - 图2

quad 函数

Quadrature 积分的原理参见：

http://en.wikipedia.org/wiki/Numerical_integration#Quadrature_rules_based_on_interpolating_functions

quad 返回一个 (积分值，误差) 组成的元组：

In [20]:

  1. from scipy.integrate import quad
  2. interval = [0, 6.5]
  3. value, max_err = quad(f, *interval)

积分值：

In [21]:

  1. print value
  1. 1.28474297234

最大误差：

In [22]:

  1. print max_err
  1. 2.34181853668e-09

积分区间图示，蓝色为正，红色为负：

In [23]:

  1. print "integral = {:.9f}".format(value)
  2. print "upper bound on error: {:.2e}".format(max_err)
  3. x = np.linspace(0, 10, 100)
  4. p = plt.plot(x, f(x), 'k-')
  5. x = np.linspace(0, 6.5, 45)
  6. p = plt.fill_between(x, f(x), where=f(x)>0, color="blue")
  7. p = plt.fill_between(x, f(x), where=f(x)<0, color="red", interpolate=True)
  1. integral = 1.284742972
  2. upper bound on error: 2.34e-09

04.06 积分 - 图3

积分到无穷

In [24]:

  1. from numpy import inf
  2. interval = [0., inf]
  3.  
  4. def g(x):
  5.     return np.exp(-x ** 1/2)

In [25]:

  1. value, max_err = quad(g, *interval)
  2. x = np.linspace(0, 10, 50)
  3. fig = plt.figure(figsize=(10,3))
  4. p = plt.plot(x, g(x), 'k-')
  5. p = plt.fill_between(x, g(x))
  6. plt.annotate(r"$\int_0^{\infty}e^{-x^1/2}dx = $" + "{}".format(value), (4, 0.6),
  7.          fontsize=16)
  8. print "upper bound on error: {:.1e}".format(max_err)
  1. upper bound on error: 7.2e-11

04.06 积分 - 图4

双重积分

假设我们要进行如下的积分：

I_n = \int \limits_0^{\infty} \int \limits_1^{\infty} \frac{e^{-xt}}{t^n}dt dx = \frac{1}{n}

In [26]:

  1. def h(x, t, n):
  2.     """core function, takes x, t, n"""
  3.     return np.exp(-x * t) / (t ** n)

一种方式是调用两次 quad 函数，不过这里 quad 的返回值不能向量化，所以使用了修饰符 vectorize 将其向量化：

In [27]:

  1. from numpy import vectorize
  2. @vectorize
  3. def int_h_dx(t, n):
  4.     """Time integrand of h(x)."""
  5.     return quad(h, 0, np.inf, args=(t, n))[0]

In [28]:

  1. @vectorize
  2. def I_n(n):
  3.     return quad(int_h_dx, 1, np.inf, args=(n))

In [29]:

  1. I_n([0.5, 1.0, 2.0, 5])

Out[29]:

  1. (array([ 1.97,  1.  ,  0.5 ,  0.2 ]),
  2.  array([  9.804e-13,   1.110e-14,   5.551e-15,   2.220e-15]))

或者直接调用 dblquad 函数，并将积分参数传入，传入方式有多种，后传入的先进行积分：

In [30]:

  1. from scipy.integrate import dblquad
  2. @vectorize
  3. def I(n):
  4.     """Same as I_n, but using the built-in dblquad"""
  5.     x_lower = 0
  6.     x_upper = np.inf
  7.     return dblquad(h,
  8.                    lambda t_lower: 1, lambda t_upper: np.inf,
  9.                    x_lower, x_upper, args=(n,))

In [31]:

  1. I_n([0.5, 1.0, 2.0, 5])

Out[31]:

  1. (array([ 1.97,  1.  ,  0.5 ,  0.2 ]),
  2.  array([  9.804e-13,   1.110e-14,   5.551e-15,   2.220e-15]))

采样点积分

trapz 方法 和 simps 方法

In [32]:

  1. from scipy.integrate import trapz, simps

sin 函数， 100 个采样点和 5 个采样点：

In [33]:

  1. x_s = np.linspace(0, np.pi, 5)
  2. y_s = np.sin(x_s)
  3. x = np.linspace(0, np.pi, 100)
  4. y = np.sin(x)

In [34]:

  1. p = plt.plot(x, y, 'k:')
  2. p = plt.plot(x_s, y_s, 'k+-')
  3. p = plt.fill_between(x_s, y_s, color="gray")

04.06 积分 - 图5

采用 trapezoidal 方法 和 simpson 方法 对这些采样点进行积分（函数积分为 2）：

In [35]:

  1. result_s = trapz(y_s, x_s)
  2. result_s_s = simps(y_s, x_s)
  3. result = trapz(y, x)
  4. print "Trapezoidal Integration over 5 points : {:.3f}".format(result_s)
  5. print "Simpson Integration over 5 points : {:.3f}".format(result_s_s)
  6. print "Trapezoidal Integration over 100 points : {:.3f}".format(result)
  1. Trapezoidal Integration over 5 points : 1.896
  2. Simpson Integration over 5 points : 2.005
  3. Trapezoidal Integration over 100 points : 2.000

使用 ufunc 进行积分

Numpy 中有很多 ufunc 对象：

In [36]:

  1. type(np.add)

Out[36]:

  1. numpy.ufunc

In [37]:

  1. np.info(np.add.accumulate)
  1. accumulate(array, axis=0, dtype=None, out=None)
  2.  
  3. Accumulate the result of applying the operator to all elements.
  4.  
  5. For a one-dimensional array, accumulate produces results equivalent to::
  6.  
  7.   r = np.empty(len(A))
  8.   t = op.identity        # op = the ufunc being applied to A's  elements
  9.   for i in range(len(A)):
  10.       t = op(t, A[i])
  11.       r[i] = t
  12.   return r
  13.  
  14. For example, add.accumulate() is equivalent to np.cumsum().
  15.  
  16. For a multi-dimensional array, accumulate is applied along only one
  17. axis (axis zero by default; see Examples below) so repeated use is
  18. necessary if one wants to accumulate over multiple axes.
  19.  
  20. Parameters
  21. ----------
  22. array : array_like
  23.     The array to act on.
  24. axis : int, optional
  25.     The axis along which to apply the accumulation; default is zero.
  26. dtype : data-type code, optional
  27.     The data-type used to represent the intermediate results. Defaults
  28.     to the data-type of the output array if such is provided, or the
  29.     the data-type of the input array if no output array is provided.
  30. out : ndarray, optional
  31.     A location into which the result is stored. If not provided a
  32.     freshly-allocated array is returned.
  33.  
  34. Returns
  35. -------
  36. r : ndarray
  37.     The accumulated values. If `out` was supplied, `r` is a reference to
  38.     `out`.
  39.  
  40. Examples
  41. --------
  42. 1-D array examples:
  43.  
  44. >>> np.add.accumulate([2, 3, 5])
  45. array([ 2,  5, 10])
  46. >>> np.multiply.accumulate([2, 3, 5])
  47. array([ 2,  6, 30])
  48.  
  49. 2-D array examples:
  50.  
  51. >>> I = np.eye(2)
  52. >>> I
  53. array([[ 1.,  0.],
  54.        [ 0.,  1.]])
  55.  
  56. Accumulate along axis 0 (rows), down columns:
  57.  
  58. >>> np.add.accumulate(I, 0)
  59. array([[ 1.,  0.],
  60.        [ 1.,  1.]])
  61. >>> np.add.accumulate(I) # no axis specified = axis zero
  62. array([[ 1.,  0.],
  63.        [ 1.,  1.]])
  64.  
  65. Accumulate along axis 1 (columns), through rows:
  66.  
  67. >>> np.add.accumulate(I, 1)
  68. array([[ 1.,  1.],
  69.        [ 0.,  1.]])

np.add.accumulate 相当于 cumsum ：

In [38]:

  1. result_np = np.add.accumulate(y) * (x[1] - x[0]) - (x[1] - x[0]) / 2

In [39]:

  1. p = plt.plot(x, - np.cos(x) + np.cos(0), 'rx')
  2. p = plt.plot(x, result_np)

04.06 积分 - 图6

速度比较

计算积分：\int_0^x sin \theta d\theta

In [40]:

  1. import sympy
  2. from sympy.abc import x, theta
  3. sympy_x = x

In [41]:

  1. x = np.linspace(0, 20 * np.pi, 1e+4)
  2. y = np.sin(x)
  3. sympy_y = vectorize(lambda x: sympy.integrate(sympy.sin(theta), (theta, 0, x)))

numpy 方法：

In [42]:

  1. %timeit np.add.accumulate(y) * (x[1] - x[0])
  2. y0 = np.add.accumulate(y) * (x[1] - x[0])
  3. print y0[-1] 
  1. The slowest run took 4.32 times longer than the fastest. This could mean that an intermediate result is being cached 
  2. 10000 loops, best of 3: 56.2 µs per loop
  3. -2.34138044756e-17

quad 方法：

In [43]:

  1. %timeit quad(np.sin, 0, 20 * np.pi)
  2. y2 = quad(np.sin, 0, 20 * np.pi, full_output=True)
  3. print "result = ", y2[0]
  4. print "number of evaluations", y2[-1]['neval']
  1. 10000 loops, best of 3: 40.5 µs per loop
  2. result =  3.43781337153e-15
  3. number of evaluations 21

trapz 方法：

In [44]:

  1. %timeit trapz(y, x)
  2. y1 = trapz(y, x)
  3. print y1
  1. 10000 loops, best of 3: 105 µs per loop
  2. -4.4408920985e-16

simps 方法：

In [45]:

  1. %timeit simps(y, x)
  2. y3 = simps(y, x)
  3. print y3
  1. 1000 loops, best of 3: 801 µs per loop
  2. 3.28428554968e-16

sympy 积分方法：

In [46]:

  1. %timeit sympy_y(20 * np.pi)
  2. y4 = sympy_y(20 * np.pi)
  3. print y4
  1. 100 loops, best of 3: 6.86 ms per loop
  2. 0

原文: https://nbviewer.jupyter.org/github/lijin-THU/notes-python/blob/master/04-scipy/04.06-integration-in-python.ipynb