Scipy是一个用于数学、科学、工程领域的常用软件包，可以处理插值、积分、优化、图像处理、常微分方程数值解的求解、信号处理等问题。它用于有效计算Numpy矩阵，使Numpy和Scipy协同工作，高效解决问题。
Scipy 由不同科学计算领域的子模块组成：
| 子模块 | 描述 |
|---|---|
cluster | 聚类算法 |
constants | 物理数学常数 |
fftpack | 快速傅里叶变换 |
integrate | 积分和常微分方程求解 |
interpolate | 插值 |
io | 输入输出 |
linalg | 线性代数 |
odr | 正交距离回归 |
optimize | 优化和求根 |
signal | 信号处理 |
sparse | 稀疏矩阵 |
spatial | 空间数据结构和算法 |
special | 特殊方程 |
stats | 统计分布和函数 |
weave | C/C++ 积分 |
在使用 Scipy 之前，为了方便，假定这些基础的模块已经被导入：
-
import numpy
as np
-
import scipy
as sp
-
import matplotlib
as mpl
-
import matplotlib.pyplot
as plt
由于 Scipy 以 Numpy 为基础，因此很多基础的 Numpy 函数可以在scipy 命名空间中直接调用
1.插值
样条插值法是一种以可变样条来作出一条经过一系列点的光滑曲线的数学方法
-
from scipy.interpolate
import interp1d
-
np.set_printoptions(precision=
2, suppress=
True)
#设置 Numpy 浮点数显示格式
-
#从文本中读入数据
-
data = np.genfromtxt(
"JANAF_CH4.txt",
-
delimiter=
"\t",
# TAB 分隔
-
skiprows=
1,
# 忽略首行
-
names=
True,
# 读入属性
-
missing_values=
"INFINITE",
# 缺失值
-
filling_values=np.inf)
# 填充缺失值
-
ch4_cp = interp1d(data[
'TK'], data[
'Cp'])
#默认情况下，输入值要在插值允许的范围内，否则插值会报错
我们可以通过 kind 参数来调节使用的插值方法，来得到不同的结果:
nearest最近邻插值zero0阶插值linear线性插值quadratic二次插值cubic三次插值4,5,6,7更高阶插值
对于径向基函数，其插值的公式为：
∑njφ(|x−xj|)∑njφ(|x−xj|)
查看具体数值:
-
integrate(y, (theta,
0, sympy.pi))
.evalf()
-
integrate(y, (theta,
0, np.pi))
(2)数值积分
导入贝塞尔函数
-
from scipy.special
import jv
-
def
f(
x):
-
return jv(
2.5, x)
积分到无穷
-
from numpy
import inf
-
interval = [
0., inf]
-
-
def
g(
x):
-
return np.exp(-x **
1/
2)
-
value, max_err = quad(g, *interval)
-
x = np.linspace(
0,
10,
50)
-
fig = plt.figure(figsize=(
10,
3))
-
p = plt.plot(x, g(x),
'k-')
-
p = plt.fill_between(x, g(x))
-
plt.annotate(
r"$\int_0^{\infty}e^{-x^1/2}dx = $" +
"{}".
format(value), (
4,
0.6),
-
fontsize=
16)
-
print
"upper bound on error: {:.1e}".
format(max_err)
6.解微分方程
（1）简单一阶：
-
from scipy.integrate
import odeint
-
def
dy_dt(
y, t):
-
return np.sin(t)
-
t = np.linspace(
0,
2*pi,
100)
-
-
result = odeint(dy_dt,
0, t)
-
-
fig = figure(figsize=(
12,
4))
-
p = plot(t, result,
"rx", label=
r"$\int_{0}^{x}sin(t) dt $")
-
p = plot(t, -cos(t) + cos(
0), label=
r"$cos(0) - cos(t)$")
-
p = plot(t, dy_dt(
0, t),
"g-", label=
r"$\frac{dy}{dt}(t)$")
-
l = legend(loc=
"upper right")
-
xl = xlabel(
"t")
（2）高阶微分方程：
-
def
dy_dt(
y, t):
-
"""Governing equations for projectile motion with drag.
-
y[0] = position
-
y[1] = velocity
-
g = gravity (m/s2)
-
D = drag (1/s) = force/velocity
-
m = mass (kg)
-
"""
-
g = -
9.8
-
D =
0.1
-
m =
0.15
-
dy1 = g - (D/m) * y[
1]
-
dy0 = y[
1]
if y[
0] >=
0
else
0.
-
return [dy0, dy1]
-
-
position_0 =
0.
-
velocity_0 =
100
-
t = linspace(
0,
12,
100)
-
y = odeint(dy_dt, [position_0, velocity_0], t)
-
p = plot(t, y[:,
0])
-
yl = ylabel(
"Height (m)")
-
xl = xlabel(
"Time (s)")
7.稀疏矩阵
稀疏矩阵主要使用 位置 + 值 的方法来存储矩阵的非零元素，根据存储和使用方式的不同，有如下几种类型的稀疏矩阵：
| 类型 | 描述 |
|---|---|
bsr_matrix(arg1[, shape, dtype, copy, blocksize]) | Block Sparse Row matrix |
coo_matrix(arg1[, shape, dtype, copy]) | A sparse matrix in COOrdinate format. |
csc_matrix(arg1[, shape, dtype, copy]) | Compressed Sparse Column matrix |
csr_matrix(arg1[, shape, dtype, copy]) | Compressed Sparse Row matrix |
dia_matrix(arg1[, shape, dtype, copy]) | Sparse matrix with DIAgonal storage |
dok_matrix(arg1[, shape, dtype, copy]) | Dictionary Of Keys based sparse matrix. |
lil_matrix(arg1[, shape, dtype, copy]) | Row-based linked list sparse matrix |
在这些存储格式中：
- COO 格式在构建矩阵时比较高效
- CSC 和 CSR 格式在乘法计算时比较高效
-
from scipy.sparse
import *
-
import numpy
as np
-
-
coo_matrix((
2,
3))
-
#也可以使用一个已有的矩阵或数组或列表中创建新矩阵：
-
A = coo_matrix([[
1,
2,
0],[
0,
0,
3],[
4,
0,
5]])
(0, 0) 1 (0, 1) 2 (1, 2) 3 (2, 0) 4 (2, 2) 5
可以转化为普通矩阵：
C = A.todense()
matrix([[1, 2, 0], [0, 0, 3], [4, 0, 5]])
还可以传入一个 (data, (row, col)) 的元组来构建稀疏矩阵：
-
I
= np.array([
0,3,1,0])
-
J
= np.array([
0,3,1,2])
-
V
= np.array([
4,5,7,9])
-
A
= coo_matrix((V,(I,J)),shape
=(
4,4))
-
print A
(0, 0) 4 (3, 3) 5 (1, 1) 7 (0, 2) 9
8.线性代数
-
import numpy
as np
-
import numpy.linalg
-
import scipy
as sp
-
import scipy.linalg
-
import matplotlib.pyplot
as plt
-
from scipy
import linalg
scipy.linalg 包含 numpy.linalg 中的所有函数，同时还包含了很多 numpy.linalg 中没有的函数，在使用时，我们一般使用 scipy.linalg 而不是 numpy.linalg
（1）求解方程
-
A = np.array([[
1,
3,
5],
-
[
2,
5,
1],
-
[
2,
3,
8]])
-
b = np.array([
10,
8,
3])
-
x = linalg.solve(A, b)
（2）计算行列式
-
A
= np.array([[
1,
3,
5],
-
[
2,
5,
1],
-
[
2,
3,
8]])
-
-
linalg.det(A)
（3）矩阵分解
对于给定的 N×NN×N 矩阵 AA，特征值和特征向量问题相当与寻找标量 λλ 和对应的向量 vv 使得：
Av=λvAv=λv
矩阵的 NN 个特征值（可能相同）可以通过计算特征方程的根得到：
|A−λI|=0|A−λI|=0
然后利用这些特征值求（归一化的）特征向量。
问题求解
linalg.eig(A)- 返回矩阵的特征值与特征向量
linalg.eigvals(A)- 返回矩阵的特征值
linalg.eig(A, B)- 求解 Av=λBvAv=λBv 的问题
（4）对稀疏矩阵
scipy.sparse.linalg.inv- 稀疏矩阵求逆
scipy.sparse.linalg.expm- 求稀疏矩阵的指数函数
scipy.sparse.linalg.norm- 稀疏矩阵求范数
对于特别大的矩阵，原来的方法可能需要太大的内存，考虑使用这两个方法替代：
scipy.sparse.linalg.eigs- 返回前 k 大的特征值和特征向量
scipy.sparse.linalg.svds- 返回前 k 大的奇异值和奇异向量
1万+
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
