【实验楼-Python 科学计算】SciPy - 科学计算库（上）

SciPy 库建立在 Numpy 库之上，提供了大量科学算法，主要包括这些主题：

· 特殊函数 (scipy.special)

· 积分(scipy.integrate)

· 最优化 (scipy.optimize)

· 插值(scipy.interpolate)

· 傅立叶变换 (scipy.fftpack)

· 信号处理 (scipy.signal)

· 线性代数 (scipy.linalg)

· 稀疏特征值 (scipy.sparse)

· 统计(scipy.stats)

· 多维图像处理 (scipy.ndimage)

· 文件IO (scipy.io)

• 特定函数
• 在计算科学问题时，常常会用到很多特定的函数，SciPy 提供了一个非常广泛的特定函数集合。函数列表可参考： http://docs.scipy.org/doc/scipy/reference/special.html#module-scipy.special
• 为了演示特定函数的一般用法我们拿贝塞尔函数举例：
• #
• # The scipy.special module includes alarge number of Bessel-functions
• # Here we will use the functions jnand yn, which are the Bessel functions
• # of the first and second kind andreal-valued order. We also include the
• # function jn_zeros and yn_zeros thatgives the zeroes of the functions jn
• # and yn.
• #
• %matplotlib qt
• from scipy.special import jn, yn,jn_zeros, yn_zeros
• importmatplotlib.pyplot as plt
• n = 0 # order
• x = 0.0
• # Bessel function of first kind
• print"J_%d(%f)= %f" % (n, x, jn(n, x))
• x = 1.0
• # Bessel function of second kind
• print"Y_%d(%f)= %f" % (n, x, yn(n, x))
• => J_0(0.000000) = 1.000000
• Y_0(1.000000) = 0.088257
• x = linspace(0, 10, 100)
• fig, ax = plt.subplots()
• for n in range(4):
• ax.plot(x, jn(n, x), label=r"$J_%d(x)$" % n)
• ax.legend();
• fig

• # zeros of Bessel functions
• n = 0# order
• m = 4# number ofroots to compute
• jn_zeros(n, m)
• => array([ 2.40482556, 5.52007811, 8.65372791, 11.79153444])
• 积分
• 数值积分: 求积

• 被称作 数值求积，Scipy提供了一些列不同类型的求积函数，像是 quad, dblquad 还有 tplquad 分别对应单积分，双重积分，三重积分。
• fromscipy.integrate import quad, dblquad, tplquad
• quad 函数有许多参数选项来调整该函数的行为（详情见help(quad)）。
• 一般用法如下:
• # define a simple function for theintegrand
• def f(x):
• return x
• x_lower = 0# the lowerlimit of x
• x_upper = 1# the upperlimit of x
• val, abserr = quad(f, x_lower,x_upper)
• print"integralvalue =", val, ", absolute error =", abserr
• => integral value = 0.5 , absoluteerror = 5.55111512313e-15
• 如果我们需要传递额外的参数，可以使用 args 关键字：
• def integrand(x, n):
• """
• Bessel function of first kind and order n.
• """
• return jn(n, x)
• x_lower = 0 # the lower limit of x
• x_upper = 10# the upperlimit of x
• val, abserr = quad(integrand, x_lower,x_upper, args=(3,))
• print val, abserr
• => 0.7366751370819.38925687719e-13
• 对于简单的函数我们可以直接使用匿名函数：
• val, abserr = quad(lambda x: exp(-x ** 2), -Inf, Inf)
• print"numerical =", val, abserr
• analytical = sqrt(pi)
• print"analytical=", analytical
• => numerical = 1.772453850911.42026367809e-08
• analytical = 1.77245385091
• 如例子所示，'Inf' 与 '-Inf' 可以表示数值极限。
• 高阶积分用法类似:
• def integrand(x, y):
• return exp(-x**2-y**2)
• x_lower = 0
• x_upper = 10
• y_lower = 0
• y_upper = 10
• val, abserr = dblquad(integrand,x_lower, x_upper, lambda x : y_lower, lambda x: y_upper)
• print val, abserr
• => 0.7853981633971.63822994214e-13
• 注意到我们为y积分的边界传参的方式，这样写是因为y可能是关于x的函数。
• 常微分方程 (ODEs)
• SciPy 提供了两种方式来求解常微分方程：基于函数 odeint 的API与基于 ode 类的面相对象的API。通常 odeint 更好上手一些，而 ode 类更灵活一些。
• 这里我们将使用 odeint 函数，首先让我们载入它：
• fromscipy.integrate import odeint, ode
• 常微分方程组的标准形式如下:

• 当

• 为了求解常微分方程我们需要知道方程

与初始条件

。 注意到高阶常微分方程常常写成引入新的变量作为中间导数的形式。 一旦我们定义了函数 f 与数组 y_0 我们可以使用 odeint 函数：

• y_t = odeint(f, y_0,t)
• 我们将会在下面的例子中看到 Python 代码是如何实现 f 与 y_0 。
• 示例:双摆
• 让我们思考一个物理学上的例子：双摆
• 关于双摆，参考：http://en.wikipedia.org/wiki/Double_pendulum
• Image(url='http://upload.wikimedia.org/wikipedia/commons/c/c9/Double-compound-pendulum-dimensioned.svg')
• 维基上已给出双摆的运动方程:

• 为了使 Python 代码更容易实现，让我们介绍新的变量名与向量表示法：

• g = 9.82
• L = 0.5
• m = 0.1
• def dx(x, t):
• """
• The right-hand side of the pendulum ODE
• """
• x1, x2, x3, x4 = x[0], x[1], x[2], x[3]
• dx1 = 6.0/(m*L**2) * (2 * x3 - 3 * cos(x1-x2) *x4)/(16 - 9 * cos(x1-x2)**2)
• dx2 = 6.0/(m*L**2) * (8 * x4 - 3 * cos(x1-x2) *x3)/(16 - 9 * cos(x1-x2)**2)
• dx3 = -0.5 * m * L**2 * ( dx1 * dx2* sin(x1-x2) + 3 * (g/L) * sin(x1))
• dx4 = -0.5 * m * L**2 * (-dx1 * dx2* sin(x1-x2) + (g/L) * sin(x2))
• return [dx1, dx2,dx3, dx4]
• # choose an initial state
• x0 = [pi/4, pi/2, 0, 0]
• # time coodinate to solve the ODE for:from 0 to 10 seconds
• t = linspace(0, 10, 250)
• # solve the ODE problem
• x = odeint(dx, x0, t)
• # plot the angles as a function oftime
• fig, axes = plt.subplots(1,2, figsize=(12,4))
• axes[0].plot(t, x[:, 0], 'r', label="theta1")
• axes[0].plot(t, x[:, 1], 'b', label="theta2")
• x1 = + L * sin(x[:, 0])
• y1 = - L * cos(x[:, 0])
• x2 = x1 + L * sin(x[:, 1])
• y2 = y1 - L * cos(x[:, 1])
• axes[1].plot(x1, y1, 'r', label="pendulum1")
• axes[1].plot(x2, y2, 'b', label="pendulum2")
• axes[1].set_ylim([-1, 0])
• axes[1].set_xlim([1, -1]);
• fig

• 我们将在第四节课看到如何做出更好的演示动画。
• from IPython.display importclear_output
• import time
• fig, ax = plt.subplots(figsize=(4,4))
• for t_idx, tt inenumerate(t[:200]):
• x1 = + L * sin(x[t_idx, 0])
• y1 = - L * cos(x[t_idx, 0])
• x2 = x1 + L * sin(x[t_idx, 1])
• y2 = y1 - L * cos(x[t_idx, 1])
• ax.cla()
• ax.plot([0, x1], [0, y1], 'r.-')
• ax.plot([x1, x2], [y1, y2], 'b.-')
• ax.set_ylim([-1.5, 0.5])
• ax.set_xlim([1, -1])
• display(fig)
• clear_output()
• time.sleep(0.1)
• fig

• 示例：阻尼谐震子
• 常微分方程问题在计算物理学中非常重要，所以我们接下来要看另一个例子：阻尼谐震子。wiki地址：http://en.wikipedia.org/wiki/Damping
• 阻尼震子的运动公式：

• 在这个例子的实现中，我们会加上额外的参数到 RHS 方程中：
• def dy(y, t, zeta,w0):
• """
• The right-hand side of the dampedoscillator ODE
• """
• x, p = y[0], y[1]
• dx = p
• dp = -2 * zeta * w0 *p - w0**2 * x
• return [dx, dp]
• # initial state:
• y0 = [1.0, 0.0]
• # time coodinate to solve the ODE for
• t = linspace(0, 10, 1000)
• w0 = 2*pi*1.0
• # solve the ODE problem for threedifferent values of the damping ratio
• y1 = odeint(dy, y0, t, args=(0.0, w0)) # undamped
• y2 = odeint(dy, y0, t, args=(0.2, w0)) # under damped
• y3 = odeint(dy, y0, t, args=(1.0, w0)) # critialdamping
• y4 = odeint(dy, y0, t, args=(5.0, w0)) # over damped
• fig, ax = plt.subplots()
• ax.plot(t, y1[:,0], 'k', label="undamped", linewidth=0.25)
• ax.plot(t, y2[:,0], 'r', label="underdamped")
• ax.plot(t, y3[:,0], 'b', label=r"criticaldamping")
• ax.plot(t, y4[:,0], 'g', label="overdamped")
• ax.legend();
• fig

• 傅立叶变换
• 傅立叶变换是计算物理学所用到的通用工具之一。Scipy 提供了使用 NetLib FFTPACK 库的接口，它是用FORTRAN写的。Scipy 还另外提供了很多便捷的函数。不过大致上接口都与 NetLib 的接口差不多。
• 让我们加载它：
• from scipy.fftpack import *
• 下面演示快速傅立叶变换，例子使用上节阻尼谐震子的例子：
• N = len(t)
• dt = t[1]-t[0]
• # calculate the fast fourier transform
• # y2 is the solution to theunder-damped oscillator from the previous section
• F = fft(y2[:,0])
• # calculate the frequencies for thecomponents in F
• w = fftfreq(N, dt)
• fig, ax = plt.subplots(figsize=(9,3))
• ax.plot(w, abs(F));
• fig

• 既然信号是实数，同时频谱是对称的。那么我们只需要画出正频率所对应部分的图：
• indices = where(w >0) # select only indices for elements that corresponds to positive frequencies
• w_pos = w[indices]
• F_pos = F[indices]
• fig, ax = subplots(figsize=(9,3))
• ax.plot(w_pos, abs(F_pos))
• ax.set_xlim(0, 5);
• fig

• 正如预期的那样，我们可以看到频谱的峰值在1处。1就是我们在上节例子中所选的频率。
• 求土豪给知识的搬运工一点奖励：