几年前我参加了数学建模竞赛当时流行从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符号计算SymPyplot、surf、scatter3Matplotlibtable、CSV、Excelpandas稀疏线性代数SciPy sparse2.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 创建数组的一些基本函数MATLABNumPy1:5np.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]])以及查看数组属性的一些基本函数MATLABNumPy数组形状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).dtype4.索引和切片Python 与 MATLAB 的主要区别Python 下标从0开始切片右端点不包含使用方括号[]负数索引从末尾倒数MATLABPython第一个元素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.矩阵基础运算MATLABNumPy矩阵乘法A*BAB元素乘法A.*BA*B元素除法A./BA/B元素幂A.^2A**2转置AA.T线性代数求解Axb时优先使用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里axisn可理解为沿第n个维度压缩运算该维度默认消失。如果使用keepdimsTrue则维度不消失。A np.array([[1, 2, 3], [4, 5, 6],[7,8,9]]) print(所有元素, A.sum()) print(axis0每列, A.sum(axis0,keepdimsFalse)) print(axis1每行, A.sum(axis1,keepdimsFalse)) print(A.sum(axis0).shape, A.sum(axis1).shape)所有元素 45 axis0每列 [12 15 18] axis1每行 [ 6 15 24] (3,) (3,)8.numpy广播广播允许不同形状数组直接运算不必手动复制。从最后一维开始比较只要两维相同、其中一个为1或某个数组缺少这一维就可能广播。如果是一个(n,:)也就是matlab的一维数组会在前面填入若干个为1的维度与另一数组对齐。9.引用、视图与复制BA不复制数组只是让两个变量指向同一个对象。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真正复制BA.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, realTrue) 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), labelf(x)) ax.plot(x_grid, df_numpy(x_grid), labelf(x)) ax.axhline(0, linewidth0.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, labelNumerical derivative) ax.plot(x_data, np.cos(x_data), --, labelExact 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 数值算法展开写太漫长了手头也暂时找不到什么特别好的实例就此带过。如果发出去真的有人读的话我会回来完善的。功能模块积分和 ODEscipy.integrate优化和求根scipy.optimize插值scipy.interpolate信号处理scipy.signal统计分布scipy.stats稀疏矩阵scipy.sparse高级线性代数scipy.linalgMATLAB 文件scipy.io13.Matplotlib 二维作图图片真的很重要啊rng np.random.default_rng(42) x_scatter rng.normal(size100) y_scatter 2*x_scatter rng.normal(scale0.8, size100) 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, dpi200, bbox_inchestight)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(projection3d) surface ax.plot_surface(X, Y, Z, cmapviridis, linewidth0) ax.set_xlabel(x); ax.set_ylabel(y); ax.set_zlabel(z) ax.set_title(r$z\sin(\sqrt{x^2y^2})$) fig.colorbar(surface, axax, shrink0.65, labelz) plt.show()线框和三维散点fig plt.figure(figsize(8, 6)) ax fig.add_subplot(projection3d) ax.plot_wireframe(X, Y, Z, rstride6, cstride6) ax.set_title(3D wireframe) plt.show() rng np.random.default_rng(42) x3 rng.normal(size150); y3 rng.normal(size150); z3 x3**2 - y3**2 fig plt.figure(figsize(8, 6)) ax fig.add_subplot(projection3d) points ax.scatter(x3, y3, z3, cz3, cmapcoolwarm) ax.set_title(r$zx^2-y^2$) fig.colorbar(points, axax, shrink0.65, labelz) plt.show()15.随机数与可复现rng np.random.default_rng(42) print(rng.uniform(0, 1, size5)) print(rng.normal(0, 1, size5)) print(rng.random((2, 3)))固定种子后从头运行会得到相同随机结果便于复现与做展示。哇写个文章真的累我居然写完了。写的内容也比较宽泛感觉到最后反而像是numpy介绍主要是把之前用matlab做数值方面一些工具在python里整理了下有很多写得不足的地方还请多指点。不过真的会有人看吗
)
