几年前,我参加了数学建模竞赛,当时流行从matlab和python中选择一个编程。由于导师使用的是matlab,我也用这个参加了竞赛,主要做数值计算和结果展示。

        但是慢慢地,我发现了matlab这个语言的局限性——比如有一次,我希望将经纬坐标标注在真实地图上,matlab很难完成这点需求,但是用python就可以调库实现。包括这个在内,很多很炫酷新奇的东西都需要用python实现。所以,我一直希望能用好python。

        然而真正使用了之后,我又觉得python有点“死板”,我习惯了matlab的工作区,变量可以一直保存下来,随时修改,重复计算,很适合快速验证自己的想法(虽然代码规范很差,但是探索性计算还挺有意思),相比之下python没有这种能力。直到我前几天尝试了Jupyter,怎么有这么好用的东西,既有matlab的快捷方便,又有python的丰富生态。于是我想记录一下用Jupyter使我的计算从matlab转向python的过程。

        对我来说,matlab就是一个超大号计算器,做数值比较方便,所以我接下来也主要写一些数值计算的问题。对于环境搭建,就暂且略过。

目录

1.python常用库

2.Jupyter、Kernel 与“工作区”

3.NumPy 数组基础

4.索引和切片

5.矩阵基础运算

6.一维向量、形状与转置

7. 按维度处理:axis,keepdims:

8.numpy广播

9.引用、视图与复制

10.SymPy 符号计算

11.NumPy 数值求导

12.SciPy 数值算法

13.Matplotlib 二维作图

14.Matplotlib 三维作图

15.随机数与可复现


1.python常用库

matlab把工具都集成完毕,python需要调用各种库。

大家都很眼熟,但要问每个库具体是干什么的,可能不太清楚,下面简单整理了一下

MATLAB 中的功能 Python 中对应库
矩阵和数组运算 NumPy
数值积分、优化、ODE、插值 SciPy
符号计算 SymPy
plotsurfscatter3 Matplotlib
table、CSV、Excel pandas
稀疏线性代数 SciPy sparse

2.Jupyter、Kernel 与“工作区”

Notebook 文件保存代码、Markdown、输出和图像;真正保存变量状态的是背后的 Kernel

只要 Kernel 没重启,前面单元创建的变量就可以在后面继续使用,不必从头运行。这种体验很接近 MATLAB Workspace。

import numpy as np
import matplotlib.pyplot as plt
import sympy as sp

a = 10
x = np.linspace(0, 2 * np.pi, 100)
y = np.sin(x)
A = np.array([[1, 2],
              [3, 4]])

print('变量已创建。')

%whos查看当前“工作区”

删除单个变量:del variable
清空用户变量:%reset -f
重启 Kernel:变量全部消失。

但是似乎直接编辑变量不如matlab那么方便?不过这个本来就是坏习惯,没有这个也是好事。

3.NumPy 数组基础

Python 原生 list 是通用容器;科学计算通常使用 NumPy 的

a = [1, 2, 3]
b = [4, 5, 6]

print(a + b)

 list 的结果是[1, 2, 3, 4, 5, 6],是通用的拼接,而非数学的计算

下面是 NumPy 创建数组的一些基本函数

MATLAB NumPy
1:5 np.arange(1, 6)
linspace(0,1,6) np.linspace(0,1,6)
zeros(2,3) np.zeros((2,3))
ones(2,3) np.ones((2,3))
eye(3) np.eye(3)
[1 2;3 4] np.array([[1,2],[3,4]])

以及查看数组属性的一些基本函数

MATLAB NumPy
数组形状 size(B) B.shape
数组维数 ndims(B) B.ndim
元素总数 numel(B) B.size
数据类型 class(B) B.dtype
转换成双精度浮点 double(B) B.astype(float)
查看转换后的类型 class(double(B)) B.astype(float).dtype

4.索引和切片

Python 与 MATLAB 的主要区别:

  1. Python 下标从 0 开始
  2. 切片右端点不包含
  3. 使用方括号 []
  4. 负数索引从末尾倒数
MATLAB Python
第一个元素 x(1) x[0]
最后一个元素 x(end) x[-1]
第2到第4个 x(2:4) x[1:4]
第2列 A(:,2) A[:,1]
前两行 A(1:2,:) A[:2,:]

5.矩阵基础运算

MATLAB NumPy
矩阵乘法 A*B A@B
元素乘法 A.*B A*B
元素除法 A./B A/B
元素幂 A.^2 A**2
转置 A' A.T

线性代数:求解 Ax=b 时优先使用 solve,不要先求逆。

A = np.array([[3.0, 1.0], [1.0, 2.0]])
b = np.array([9.0, 8.0])
solution = np.linalg.solve(A, b)
det_A = np.linalg.det(A)
eigenvalues, eigenvectors = np.linalg.eig(A)
print('解:', solution)
print('行列式:', det_A)
print('特征值:', eigenvalues)
print('特征向量:\n', eigenvectors)
解: [2. 3.]
行列式: 5.000000000000001
特征值: [3.61803399 1.38196601]
特征向量:
 [[ 0.85065081 -0.52573111]
 [ 0.52573111  0.85065081]]

6.一维向量、形状与转置

np.array([1,2,3]) 的形状是 (3,),既不是严格行向量,也不是严格列向量;对它使用 .T 不会改变形状。

x_vector = np.array([1, 2, 3])
print(x_vector.shape)
print(x_vector.T.shape)

x_column = x_vector[:, None]
x_row = x_vector[None, :]
print('列向量:\n', x_column, x_column.shape)
print('行向量:\n', x_row, x_row.shape)
(3,)
(3,)
列向量:
 [[1]
 [2]
 [3]] (3, 1)
行向量:
 [[1 2 3]] (1, 3)

7. 按维度处理:axis,keepdims:

matlab的axis是用来调整绘图的,而在python里axis=n 可理解为:沿第 n 个维度压缩运算,该维度默认消失。如果使用keepdims=True,则维度不消失。

A = np.array([[1, 2, 3], [4, 5, 6],[7,8,9]])
print('所有元素:', A.sum())
print('axis=0,每列:', A.sum(axis=0,keepdims=False))
print('axis=1,每行:', A.sum(axis=1,keepdims=False))
print(A.sum(axis=0).shape, A.sum(axis=1).shape)
所有元素: 45
axis=0,每列: [12 15 18]
axis=1,每行: [ 6 15 24]
(3,) (3,)

8.numpy广播

广播允许不同形状数组直接运算,不必手动复制。从最后一维开始比较,只要两维相同、其中一个为 1,或某个数组缺少这一维,就可能广播。如果是一个(n,:),也就是matlab的一维数组,会在前面填入若干个为1的维度,与另一数组对齐。

9.引用、视图与复制

B=A 不复制数组,只是让两个变量指向同一个对象。

A = np.array([1, 2, 3])
B = A
B[0] = 100
print('A =', A)
print('A is B:', A is B)
A = [100   2   3]
A is B: True

真正复制:B=A.copy()。切片通常返回共享底层数据的视图。

A = np.array([1, 2, 3, 4])
view = A[:2]
view[0] = 100
print('视图修改后 A =', A)

A = np.array([1, 2, 3, 4])
copy_part = A[:2].copy()
copy_part[0] = 100
print('副本修改后 A =', A)
视图修改后 A = [100   2   3   4]
副本修改后 A = [1 2 3 4]

10.SymPy 符号计算

SymPy 中的变量是数学符号,最接近 MATLAB Symbolic Math Toolbox。

x = sp.symbols('x', real=True)
f = sp.exp(-x**2) * sp.sin(x)
df = sp.diff(f, x)
d2f = sp.diff(f, x, 2)
display(f)
display(df)
display(d2f)


还可以进行更丰富的计算

display(sp.integrate(x**2, x))
print('定积分:', sp.integrate(x**2, (x, 0, 1)))
print('极限:', sp.limit(sp.sin(x) / x, x, 0))
print('方程根:', sp.solve(sp.Eq(x**2 - 5*x + 6, 0), x))
display(sp.expand((x + 1)**3))
display(sp.factor(x**2 - 5*x + 6))
display(sp.simplify(sp.sin(x)**2 + sp.cos(x)**2))

我感觉这真的很疯狂,究竟是怎么实现的呢,如果以后有机会,我可能会另外写一篇博客

lambdify:把符号表达式转换成 NumPy 可计算函数

f_numpy = sp.lambdify(x, f, 'numpy')
df_numpy = sp.lambdify(x, df, 'numpy')
x_grid = np.linspace(-3, 3, 400)

fig, ax = plt.subplots(figsize=(7, 4))
ax.plot(x_grid, f_numpy(x_grid), label='f(x)')
ax.plot(x_grid, df_numpy(x_grid), label="f'(x)")
ax.axhline(0, linewidth=0.8)
ax.set_xlabel('x'); ax.set_ylabel('y')
ax.set_title('Function and symbolic derivative')
ax.grid(True); ax.legend()
plt.show()

11.NumPy 数值求导

只有离散数据、没有解析公式时,可用 np.gradient。只是一个方便的默认数值求导工具,并不是能自由选择所有差分格式。内部点默认使用二阶精度的中心差分,不能通过参数改成前向差分、五点四阶中心差分等其他格式。但它也支持非均匀网格坐标。

x_data = np.linspace(0, 2 * np.pi, 200)
y_data = np.sin(x_data)
dy_dx = np.gradient(y_data, x_data)

fig, ax = plt.subplots(figsize=(7, 4))
ax.plot(x_data, dy_dx, label='Numerical derivative')
ax.plot(x_data, np.cos(x_data), '--', label='Exact cos(x)')
ax.set_xlabel('x'); ax.set_ylabel('dy/dx')
ax.set_title('Numerical differentiation')
ax.grid(True); ax.legend()
plt.show()

12.SciPy 数值算法

展开写太漫长了,手头也暂时找不到什么特别好的实例,就此带过。

如果发出去真的有人读的话,我会回来完善的。

功能 模块
积分和 ODE scipy.integrate
优化和求根 scipy.optimize
插值 scipy.interpolate
信号处理 scipy.signal
统计分布 scipy.stats
稀疏矩阵 scipy.sparse
高级线性代数 scipy.linalg
MATLAB 文件 scipy.io

13.Matplotlib 二维作图

图片真的很重要啊

rng = np.random.default_rng(42)
x_scatter = rng.normal(size=100)
y_scatter = 2*x_scatter + rng.normal(scale=0.8, size=100)

fig, ax = plt.subplots(figsize=(7, 4))
ax.scatter(x_scatter, y_scatter)
ax.set_title('Scatter plot'); ax.grid(True)
plt.show()

fig, axes = plt.subplots(1, 2, figsize=(10, 4))
axes[0].plot(x, np.sin(x)); axes[0].set_title('sin(x)'); axes[0].grid(True)
axes[1].plot(x, np.cos(x)); axes[1].set_title('cos(x)'); axes[1].grid(True)
fig.tight_layout(); plt.show()

保存fig.savefig('save.png', dpi=200, bbox_inches='tight')

14.Matplotlib 三维作图

三维曲面

x = np.linspace(-6, 6, 150)
y = np.linspace(-6, 6, 150)
X, Y = np.meshgrid(x, y)
R = np.sqrt(X**2 + Y**2)
Z = np.sin(R)

fig = plt.figure(figsize=(8, 6))
ax = fig.add_subplot(projection='3d')
surface = ax.plot_surface(X, Y, Z, cmap='viridis', linewidth=0)
ax.set_xlabel('x'); ax.set_ylabel('y'); ax.set_zlabel('z')
ax.set_title(r'$z=\sin(\sqrt{x^2+y^2})$')
fig.colorbar(surface, ax=ax, shrink=0.65, label='z')
plt.show()

线框和三维散点

fig = plt.figure(figsize=(8, 6))
ax = fig.add_subplot(projection='3d')
ax.plot_wireframe(X, Y, Z, rstride=6, cstride=6)
ax.set_title('3D wireframe')
plt.show()

rng = np.random.default_rng(42)
x3 = rng.normal(size=150); y3 = rng.normal(size=150); z3 = x3**2 - y3**2
fig = plt.figure(figsize=(8, 6))
ax = fig.add_subplot(projection='3d')
points = ax.scatter(x3, y3, z3, c=z3, cmap='coolwarm')
ax.set_title(r'$z=x^2-y^2$')
fig.colorbar(points, ax=ax, shrink=0.65, label='z')
plt.show()

15.随机数与可复现

rng = np.random.default_rng(42)
print(rng.uniform(0, 1, size=5))
print(rng.normal(0, 1, size=5))
print(rng.random((2, 3)))

固定种子后,从头运行会得到相同随机结果,便于复现与做展示。

哇写个文章真的累,我居然写完了。写的内容也比较宽泛,感觉到最后反而像是numpy介绍?主要是把之前用matlab做数值方面一些工具在python里整理了下,有很多写得不足的地方,还请多指点。

不过真的会有人看吗?

Logo

这里是“一人公司”的成长家园。我们提供从产品曝光、技术变现到法律财税的全栈内容,并连接云服务、办公空间等稀缺资源,助你专注创造,无忧运营。

更多推荐