4. Scipy Tutorial-积分函数

在Scipy的integrate子包里提供了很多和积分相关的函数，例如quad可以对一般一元函数进行积分。

• 利用quad求一重积分 $$ I = \int_0^1 (3x^2 + 2x + 1)dx = (x^3 + x^2 + x)|_0^1 = 3 + c $$

import numpy as np
from scipy.integrate import quad
f = np.poly1d([3, 2, 1])
print f
ret = quad(f, 0, 1)
print ret

变量f是求积分的多项式$\ f(x) = 3x^2 + 2x + 1$，quad函数则是求定积分函数，计算结果应该是3。

在积分里$\pm\infty$用$\pm np.inf$表示。例如求下面的积分： $$ I = \int_{-\infty}^{-1} \frac{1}{x^2}dx = -\frac{1}{x}|^{-1}_{-\infty} = 1 - \frac{1}{-\infty} = 1 + \frac{1}{\infty} = 1 + 0 = 1 $$

import numpy as np
from scipy.integrate import quad
f = np.poly1d([3, 2, 1])
print f
ret = quad(f, 0, 1)
print ret
ret = quad(lambda x : 1 / x ** 2, -np.inf, -1)
print ret
import matplotlib.pyplot as plt
x = np.arange(-10,0.1,0.1)
print x
def f(x):
    return 1 / x ** 2
y = f(x)
print y
plt.ylim(-0.1, 5)
plt.plot(x, y)
plt.show()

程序执行结果：

   2
3 x + 2 x + 1
(3.0, 3.3306690738754696e-14)
(1.0, 1.1102230246251565e-14)

如果积分的函数$f$带参数，即： $$ I(a, b) = \int_0^1 (ax^2 + b)dx $$ 那么a和b可以通过args传入quad函数。

from scipy.integrate import quad
def f(x, a, b):
    return a * (x ** 2) + b
ret = quad(f, 0, 1, args=(3, 1))
print ret

• 利用dblquad 求二重积分和tplquad求三重积分。

scipy.integrate.dblquad(func, a, b, gfun, hfun, args=(), ...)

这里$x \in [a, b]$而$gfun(x) < y < hfun(x)$,即a是x的下限，gfun(x)是y的下限。 $$ I = \int_a^b[\int_{gfun(x)}^{hfun(x)}f(x,y)dy]dx $$ 下面求公式： $$ I = \int_{y = 0}^{\frac{1}{2}}\int_{x = 0}^{1 - 2 x}xydydx $$ 那么$gfun(x) = lambda\quad x : 0$、$hfun(x) = lambda\quad x :1 - 2*x$、$a = 0$、$b = 0.5$。

from scipy.integrate import dblquad
ret = dblquad(lambda x,y : x * y, 0.0, 0.5, lambda x : 0, lambda x : 1 - 2 * x)
print ret
print 1.0 / 96