Scipy——数值计算库

灰太狼 2022-05-18 00:09 308阅读 0赞

**SciPy函数库在NumPy库的基础上增加了众多的数学、科学以及工程计算中常用的库函数。例如线性代数、常微分方程数值求解、信号处理、图像处理、稀疏矩阵等等。**

#### 1.最小二乘拟合 ####

假设有一组实验数据(x\[i\], y\[i\])，我们知道它们之间的函数关系:y = f(x)，通过这些已知信息，需要确定函数中的一些参数项。例如，如果f是一个线型函数f(x) = k\*x+b，那么参数k和b就是我们需要确定的值。如果将这些参数用 p 表示的话，那么我们就是要找到一组 p 值使得如下公式中的S函数最小：  
![这里写图片描述][70]  
这种算法被称之为最小二乘拟合(Least-square fitting)。

#### 2.函数最小值 ####

optimize库提供了几个求函数最小值的算法：fmin, fmin\_powell, fmin\_cg, fmin\_bfg。  
可以使用fmin计算反卷积。

#### 3.非线性方程组 ####

optimize库中的fsolve函数可以用来对非线性方程组进行求解。它的基本调用形式如下：

fsolve(func, x0)

func(x)是计算方程组误差的函数，它的参数x是一个矢量，表示方程组的各个未知数的一组可能解，func返回将x代入方程组之后得到的误差；x0为未知数矢量的初始值。如果要对如下方程组进行求解的话：  
• f1(u1,u2,u3) = 0  
• f2(u1,u2,u3) = 0  
• f3(u1,u2,u3) = 0  
那么func可以如下定义：

def func(x):
    u1,u2,u3 = x
    return [f1(u1,u2,u3), f2(u1,u2,u3), f3(u1,u2,u3)]

在对方程组进行求解时，fsolve会自动计算方程组的雅可比矩阵，如果方程组中的未知数很多，而与每个方程有关的未知数较少时，即雅可比矩阵比较稀疏时，传递一个计算雅可比矩阵的函数将能大幅度提升运算速度。  
（**雅可比矩阵**：是一阶偏导数以一定方式排列的矩阵，它给出了可微分方程与给定点的最优线性逼近，因此类似于多元函数的导数。例如函数f1,f2,f3和未知数u1,u2,u3的雅可比矩阵如下：）  
![这里写图片描述][70 1]

#### 4.B-Spline样条曲线 ####

interpolate库提供了许多对数据进行插值运算的函数。下面是使用直线和B-Spline对正弦波上的点进行插值的例子。

# -*- coding: utf-8 -*-
    import numpy as np
    import pylab as pl
    from scipy import interpolate
    x = np.linspace(0, 2*np.pi+np.pi/4, 10)
    y = np.sin(x)
    x_new = np.linspace(0, 2*np.pi+np.pi/4, 100)
    f_linear = interpolate.interp1d(x, y)
    tck = interpolate.splrep(x, y)
    y_bspline = interpolate.splev(x_new, tck)
    pl.plot(x, y, "o", label=u"原始数据")
    pl.plot(x_new, f_linear(x_new), label=u"线性插值")
    pl.plot(x_new, y_bspline, label=u"B-spline插值")
    pl.legend()pl.show

![这里写图片描述][70 2]  
在这段程序中，通过interp1d函数直接得到一个新的线性插值函数。而B-Spline插值运算需要先使用splrep函数计算出B-Spline曲线的参数，然后将参数传递给splev函数计算出各个取样点的插值结果。

#### 5.数值积分 ####

数值积分是对定积分的数值求解，例如可以利用数值积分计算某个形状的面积。  
多重定积分的求值可以通过多次调用quad函数实现，为了调用方便，integrate库提供了dblquad函数进行二重定积分，tplquad函数进行三重积分。

#### 6.解常微分方程组 ####

scipy.integrate库提供了数值积分和常微分方程组求解算法odeint。用odeint  
计算洛仑兹吸引子的轨迹。洛仑兹吸引子由下面的三个微分方程定义：  
![这里写图片描述][70 3]  
这三个方程定义了三维空间中各个坐标点上的速度矢量。从某个坐标开始沿着速度矢量进行积分，就可以计算出无质量点在此空间中的运动轨迹。其中 σ, ρ, β 为三个常数，不同的参数可以计算出不同的运动轨迹： x(t), y(t), z(t)。 当参数为某些值时，轨迹出现馄饨现象：即微小的初值差别也会显著地影响运动轨迹。下面是洛仑兹吸引子的轨迹计算和绘制程序：

# -*- coding: utf-8 -*-
    from scipy.integrate import odeint
    import numpy as np
    def lorenz(w, t, p, r, b):# 给出位置矢量w，和三个参数p, r, b
    #计算出 dx/dt, dy/dt, dz/dt的值
    x, y, z = w
    # 直接与lorenz的计算公式对应
    return np.array([p*(y-x), x*(r-z)-y, x*y-b*z])
    
    t = np.arange(0, 30, 0.01) # 创建时间点
    # 调用ode对lorenz进行求解, 用两个不同的初始值
    track1 = odeint(lorenz, (0.0, 1.00, 0.0), t, args=(10.0, 28.0, 3.0))
    track2 = odeint(lorenz, (0.0, 1.01, 0.0), t, args=(10.0, 28.0, 3.0)
    # 绘图
    from mpl_toolkits.mplot3d import Axes3D
    import matplotlib.pyplot as plt
    fig = plt.figure()
    ax = Axes3D(fig)
    ax.plot(track1[:,0], track1[:,1], track1[:,2])
    ax.plot(track2[:,0], track2[:,1], track2[:,2])
    plt.show（）

![这里写图片描述][70 4]  
我们看到即使初始值只相差0.01，两条运动轨迹也是完全不同的。  
在程序中先定义一个lorenz函数，它的任务是计算出某个位置的各个方向的微分值，这个计算直接根据洛仑兹吸引子的公式得出。然后调用odeint，对微分方程求解，odeint有许多参数，这里用到的四个参数分别为：  
① lorenz， 它是计算某个位移上的各个方向的速度(位移的微分)  
② (0.0, 1.0, 0.0)，位移初始值。计算常微分方程所需的各个变量的初始值  
③ t， 表示时间的数组，odeint对于此数组中的每个时间点进行求解，得出所有时间点的位置  
④args， 这些参数直接传递给lorenz函数，因此它们都是常量。

#### 7.滤波器设计 ####

scipy.signal库提供了许多信号处理方面的函数。我们可以利用signal库设计滤波器，查看滤波器的频率响应，以及使用滤波器对信号进行滤波。  
\#\#\#\#8.用Weave嵌入C语言  
Python作为动态语言其功能虽然强大，但是在数值计算方面有一个最大的缺点：速度不够快。在Python级别的循环和计算的速度只有C语言程序的百分之一。因此才有了NumPy, SciPy这样的函数库，将高度优化的C、Fortran的函数库进行包装，以供Python程序调用。如果这些高度优化的函数库无法实现我们的算法，必须从头开始写循环、计算的话，那么用Python来做显然是不合适的。因此SciPy提供了快速调用C++语言程序的方法-- Weave。  
用Weave编译的C语言程序比numpy自带的sum函数还要快。而Python的内部函数sum使用数组的迭代器接口进行运算，因此速度是Python语言级别的，只有Weave版本的1/300。

[70]: /images/20220518/fa4ef5a5b6a449709e9e0f1f7a238d74.png
[70 1]: /images/20220518/1c17b607955b4666959322339cae4965.png
[70 2]: /images/20220518/21388aac8c6c4e90b3683adfbfdfeb91.png
[70 3]: /images/20220518/6f3ac459b1ad45f19b89b6994dce5738.png
[70 4]: /images/20220518/c333fb97a8ec449a8fbe098f825d2804.png