1.资料合集
scipy相关资料合集：
- https://docs.scipy.org/doc/scipy/reference/index.html 【官网】
- https://www.runoob.com/scipy/scipy-tutorial.html 【菜鸟教程】
- https://www.yggzs.com/technology/article/27319 【妖怪工作室】
- http://www.manongjc.com/detail/54-fucaifqwckblhal.html 【码农教程】
- https://www.osgeo.cn/sagemath/thematic_tutorials/numerical_sage/scipy.html【Sage 中文文档】
2.模块汇总
| 模块 | 功能 | 介绍 |
|---|---|---|
| scipy.special | 特殊函数模块 | 里面有各种特殊的数学函数，可以直接调用，如贝塞尔函数 |
| scipy.integrate | 积分模块 | 可以求多重积分，高斯积分，解常微分方程 |
| scipy.optimize | 优化模块 | 里面有各种优化算法，包括用来求有/ 无约束的多元标量函数最小值算法，最小二乘法，求有/无约束的单变量函数最小值算法，还有解各种复杂方程的算法 |
| scipy.interpolation | 插值模块 | 提供各种一维、二维、N维插值算法，包括B样条插值、径向基函数插值等 |
| scipy.fftpack | FFT（快速傅里叶变换）模块 | 可以进行FFT/ DCT/ DST |
| scipy.signal | 信号处理模块 | 包括样条插值，卷积，差分等滤波方法，还有FIR/IIR滤波，中值、排序、维纳、希尔伯特等滤波器，各种谱分析算法 |
| scipy.linalg | 线代模块 | 提供各种线性代数中的常规操作 |
| scipy.sparse | 稀疏矩阵模块 | 提供了大型稀疏矩阵计算中的各种算法 |
| scipy.spatial | 空间结构模块 | 提供了一些空间相关的数据结构和算法，如Delaunay三角剖分，共面点，凸包，维诺图，Kd树等 |
| scipy.stats | 统计模块 | 提供一些统计学上常用的函数 |
| scipy.ndimage | 多维图像处理模块 | 提供一些多维图像处理上的常用算法 |
| scipy.io | IO模块 | 提供与其他文件的接口，如matlab文件、IDL文件、Wav（音频）文件、ARFF文件 |
上表中的资料来源于妖怪工作室，本文单纯做了资料整理工作。
3.求解非线性方程组
使用fsolve或root求解非线性方程组时，需要先把非线性方程组写成$F(x)=0$，x表示向量；F(x)表示向量函数。
3.1求解方程
求解方程在定初值1.5附近的一个实根：
$$x^{980}-5.01^{979}+7.398^{978}-3.388x^{977}-x^3+5.01x^2-7.398x+3.388=0$$
- 程序实现
from scipy.optimize import fsolve, root
fx = lambda x: x**980-5.01*x**979+7.398*x**978\
-3.388*x**977-x**3+5.01*x**2-7.398*x+3.388
x1 = fsolve(fx, 1.5, maxfev=4000) # maxfev为函数调用次数
x2 = root(fx, 1.5)
print(x1,'\n','-------------'); print(x2)
- 输出结果
[1.21] # 计算结果
------------- # 以下是计算计运算过程
fjac: array([[-1.]])
fun: array([-1.23341756e+69])
message: 'The solution converged.'
nfev: 319
qtf: array([2.00183544e+72])
r: array([2.54210657e+80])
status: 1
success: True
x: array([1.21])
这里设置函数调用次数的原因是对于求解高次方程，没有通式如韦达定理来帮我们求解，即计算机会用本办法不断试错得出结果。关于高次方程的求解科普可以观摩：https://zhuanlan.zhihu.com/p/64354689。
如果设置的迭代次数过少，会提示warning。
3.2求解方程组
求下列方程组的一组数值解。
$\begin{cases}x_1^2 + x_2^2 = 1\x_1 = x_2\end{cases}$
- 程序实现
from scipy.optimize import fsolve, root
fx = lambda x: [x[0]**2+x[1]**2-1, x[0]-x[1]]
s1 = fsolve(fx, [1, 1])
s2 = root(fx, [1, 1])
print(s1,'\n','--------------'); print(s2)
- 输出结果
[0.70710678 0.70710678] # 计算结果
-------------- # 计算过程
fjac: array([[-0.81649679, -0.57734998],
[ 0.57734998, -0.81649679]])
fun: array([4.4408921e-16, 0.0000000e+00])
message: 'The solution converged.'
nfev: 9
qtf: array([-3.64625238e-10, 2.57828785e-10])
r: array([-1.73205167, -0.57735171, 1.63299357])
status: 1
success: True
x: array([0.70710678, 0.70710678])
4.积分
积分主要分为以下两类：
- 对给定函数的数值积分；
- 对给定离散点的数值积分，函数有trapz。
以下为scipy.integrate中的数值积分函数表：
| 函数 | 说明 |
|---|---|
| quad(func, a, b, args) | 计算一重数值积分 |
| dblquad(func, a, b, gfun, hfun, args) | 计算二重数值积分 |
| tplquad(func, a, b, gfun, hfun, qfun, rfun) | 计算三重数值积分 |
| nquad(func, ranges, args) | 计算多变量积分 |
4.1求微积分值
分别计算$a=2,b=1;a=2,b=10$时，$I(a,b)=\int_0^1(ax^2+bx)dx$的值。
- 程序实现
from scipy.integrate import quad
def fun42(x, a, b):
return a*x**2+b*x
I1 = quad(fun42, 0, 1, args=(2, 1))
I2 = quad(fun42, 0, 1, args=(2, 10))
print(I1); print(I2)
- 输出结果
print(I1); print(I2)
(1.1666666666666665, 1.2952601953960159e-14)
(5.666666666666667, 6.291263806209221e-14)
5.最小二乘解
对于非线性方程
$$\begin{cases}f_1(x)=0,\f_2(x)=0,\\vdots\f_x(x)=0\end{cases}$$
其中x为m维向量，一般n>m，求方程的最小二乘解，且方程组是矛盾方程组，所以实际情况在往往求的是以下多元函数的最小值：
$\delta(x) = \sum_{i=1}^nf_i^2(x)$
5.1实例
- 程序实现
from scipy.optimize import least_squares
import numpy as np
a=np.array([[245. , 164. , 192. , 232. ],
[442. , 480. , 281. , 300. ],
[126.2204, 120.7509, 90.1854, 101.4021]])
x0=a[0]; y0=a[1]; d=a[2]
fx=lambda x: np.sqrt((x0-x[0])**2+(y0-x[1])**2)-d # 最小二乘表达式
s=least_squares(fx, np.random.rand(2))
print(s, '\n', '------------', '\n', s.x)
- 输出结果
>>> print(s, '\n', '------------', '\n', s.x)
>>> active_mask: array([0., 0.])
cost: 0.3671505584636804
fun: array([-0.34333429, 0.13598942, -0.49656113, 0.59275341])
grad: array([-2.03507763e-08, 1.28918378e-08])
jac: array([[-0.75860567, -0.65155001],
[-0.11987285, -0.99278925],
[-0.4737608 , 0.88065357],
[-0.80877665, 0.58811592]])
message: '`ftol` termination condition is satisfied.'
nfev: 15
njev: 15
optimality: 2.0350776341881982e-08
status: 2
success: True
x: array([149.50894333, 359.98479551])
------------
[149.50894333 359.98479551]
6.最大模特征值及特征向量
求下列矩阵的最大模特征值及对应的特征向量：
$A = \begin{bmatrix}1&2&3\2&1&3\3&3&6\end{bmatrix}$
- 程序实现
from scipy.sparse.linalg import eigs
import numpy as np
a = np.array([[1, 2, 3], [2, 1, 3], [3, 3, 6]], dtype=float) #必须加float,否则出错
b, c = np.linalg.eig(a)
d, e = eigs(a, 1)
print('最大模特征值为：', d)
print('对应的特征向量为：\n', e)
- 输出结果
最大模特征值为： [9.+0.j]
对应的特征向量为：
[[0.40824829+0.j]
[0.40824829+0.j]
[0.81649658+0.j]]
注：假定矩阵A有n个线性无关的特征向量.n个特征值按模由大到小排列:：│λ1│> =│λ2│> =…> =│λn│，λ1是按模最大特征值。
参考：《python数学建模算法与应用》司守奎 孙玺菁
