利用 numpy 实现物理计算–物理向量符号与numpy数组的对应 ( jupyterlab 例子)

文章目录

*
– 前言
– 运行环境和依赖包
– 实现一个向量
– 简化公式表达
– 与符号运算一起使用
– 构建交互模拟试验
– 小结

前言

通过一个简单的例子，聊一聊数理代码设计入门。
使用 py 函数、numpy数组和 Pandas 进行位移计算，逐步演示如何将数理符号的数学形式化公式，转换成易读和易用的py代码，最终给出了一个在 jupyterlab 或 jupyter notebook 中的交互式模型的实现。

运行环境和依赖包

运行环境：jupyter notebook 或 jupyterlab
依赖包：numpy pandas matplotlib sympy ipywidgets

实现一个向量

我们经常使用 python 来进行物理计算，这是由于 py 采用伪代码风格，利于清晰的表达物理模型的业务逻辑，同时，代码又具备一定的工程实现能力。今天想分享一个小例子，演示一下将形式化的数理公式，转换成易读、易⽤的py代码。

使用 numpy 的关键是具有数列(或向量)的编程思维，其实就是如何使用表的思维，设一个向量r ⃗ \vec{r}r 具有三个分量，即：r ⃗ = ( x , y , z ) \vec{r}= (x,y,z)r =(x ,y ,z ) ，那么该向量在不同时刻的变化，可以记录到如下表格：

xyz000111222333

如果我们定义这个表格每一行代表一个时刻，那么每一列就是r ⃗ \vec{r}r其中一个向量在每一个时刻的值。当我们用numpy来实现这个向量时，则可以采用如下代码：

import numpy as np

r = np.array([[0,0,0],[1,1,1],[2,2,2],[3,3,3]],dtype='float16')
r.shape

假设我们定义第一行为 t 0 t_0 t 0 ​ 时刻的向量取值，那么可以用 numpy 的切片功能轻松获取该值：

r[0]

同理，其他时刻的值，只需要选择不同的行即可。那么有了这个观点，实现物理计算就容易多了。我个人觉得 numpy 实现数理计算的优势在于代码的表达，首先代码本来也是形式化的，如果代码接近于原始物理模型，那么易读性自然也就具备了。

简化公式表达

来一个简单的二维运动的例子，熟悉一下 numpy 的实现过程，看看代码能不能对小白也很好懂。

例子: 有一只兔子跑过一个停车场，非常奇怪的是停车场上正巧有个坐标系，兔子的位置坐标对于时间的函数有下式给出：
x = − 0.31 t 2 + 7.2 t + 28 y = 0.22 t 2 − 9.1 t + 30 \begin{aligned} x &= -0.31t^2+7.2t+28 \ y &= 0.22t^2-9.1t+30 \end{aligned}x y ​=−0 .3 1 t 2 +7 .2 t +2 8 =0 .2 2 t 2 −9 .1 t +3 0 ​
a) 计算并绘制出兔子在 t = [ 0 s , 15 s ] t =[0s, 15s]t =[0 s ,1 5 s ] 的所有位矢，并求出其大小和角度。
b) 计算并绘制出兔子在 t = [ 0 s , 15 s ] t =[0s, 15s]t =[0 s ,1 5 s ] 的所有速度矢量，并求出其大小和角度。

首先，需要设计一个数据结构，我们直接使用np的数组，也就是前面说过的线性表来表示。我们假设对时间进行离散化 t ( n ) = t 0 + n Δ t t(n) =t_0 +n \Delta t t (n )=t 0 ​+n Δt，那么可以得到：

r ⃗ [ n ] ≡ [ [ x 1 , y 1 ] , [ x 2 , y 2 ] , . . . , [ x n , y n ] ] \vec {r}[n] \equiv [[x_1,y_1],[x_2,y_2],…,[x_n,y_n]]r [n ]≡[[x 1 ​,y 1 ​],[x 2 ​,y 2 ​],…,[x n ​,y n ​]]

上式中，我们使用 ≡ \equiv ≡ 而不是用= ==。其含义是上式左值的物理量r ⃗ \vec {r}r 逻辑等价与右值 np 数组。因此，我们可以建立一个线性表，来记录每个时刻的位矢变化情况：

txy0
x 1 x_1 x 1 ​y 1 y_1 y 1 ​x 2 x_2 x 2 ​y 2 y_2 y 2 ​x 3 x_3 x 3 ​y 3 y_3 y 3 ​

………n
x n x_n x n ​y n y_n y n ​

那么我们只需要编写一个生成上表的py函数，即可完成该任务。

def r(t) -> tuple[float, float]:
    """计算t时刻位置矢量。

    前条件：t 为某时刻值。
    后条件: x,y 为坐标值。

"""
    x = -0.31 * t**2 + 7.2*t + 28
    y = 0.22 * t**2 - 9.1*t + 30
    return x, y

r(t) 函数实际就是兔子位置函数的公式，按照py语法进行的翻译。为了表示位移，还需要用一个调用 numpy 对象的函数来模拟兔子的位移：

def timestep_pos(n, t0, dt) -> np.ndarray:
    """计算兔子在n个位置的位移向量

    前条件：n 为位移的总步长，t0 为初始时刻，dt 为时间增加步长
    后条件: 利用 np.ndarray 记录每个时刻的位矢变化

"""
    return np.array(
        [r(t0+i*dt) for i in range(n+1)], dtype='float16'
    )

那么兔子每个时刻的位矢变化情况表就可以按如下代码计算：

positions = timestep_pos(15, 0, 1)
positions.shape

这样名为 positions 的 numpy 数组就计算机内存中，保存在该例子所需的表格中。可以进行下一步计算，也可以直接输出。这里我们打算利用 numpy 的另一个好搭档 pandas 来对数据进行显性操作，以计算其每一步的位矢大小和角度。

pandas 库为 numpy 提供了一个表的便捷且高级的操作方式，使得我们可以利用py字典方式的方便的定义数据，同时进行各种图形化的输出，首先我们要导入该工具包，并初始化DF对象：

import pandas as pd
pos = pd.DataFrame(positions, columns=['x', 'y'])
pos.iloc[0]

我们知道，位矢大小可以根据坐标公式进行计算：

∣ r ∣ = x 2 + y 2 |r| = \sqrt{x^2+y^2}∣r ∣=x 2 +y 2 ​

对应的代码可以为：

pos['distance/m'] = np.sqrt(pos['x']**2 + pos['y']**2)

位矢的角度计算可以用:

θ = a r c t a n y x \theta = arctan \frac{y}{x}θ=a r c t a n x y ​

对应的代码可以为：

pos['theta/deg'] = np.rad2deg(np.arctan(pos['y']/pos['x']))

四舍五入以后，输出：

pos = pos.round()
pos

结果为：

与符号运算一起使用

如果你还听说过sympy这个库，那么你甚至可以快速进行对这个r ( t ) r(t)r (t ) 物理模型进行微分变换。下面，如果我们需要在某个时刻的速度向量v ⃗ \vec {v}v, 那么可以利用sympy 进行直接计算:

v ⃗ = [ d x d t , d y d t ] \vec {v} = [\frac{dx}{dt},\frac{dy}{dt}]v =[d t d x ​,d t d y ​]

那么其速度方程可以为：

from sympy import *

t = symbols("t")

x = Function('x')
x = -0.31 * t**2 + 7.2*t + 28

y = Function('t')
y = 0.22 * t**2 - 9.1*t + 30

v_x = diff(x)
ccode(v_x)

v_y = diff(y)
ccode(v_y)

将v_x 和v_y 的输出拷贝粘贴，再简单修正下，根据前面的r(t)函数的套路，我们可以直接写出 v(t) 函数及其b) 的实现代码：

def v(t):
    v_x = 7.2 - 0.62*t
    v_y = 0.44*t - 9.1
    return v_x, v_y

def timestep_v(n, t0, dt) -> np.ndarray:
    return np.array(
        [v(t0+i*dt) for i in range(n+1)], dtype='float16'
    )

velocities = timestep_v(15, 0, 1)
velocities = pd.DataFrame(velocities, columns=['x', 'y'])

同理，我们计算出幅值和角度的：

velocities['magnitude'] = np.sqrt(velocities['x']**2 + velocities['y']**2)
velocities['theta'] = np.rad2deg(np.arctan(velocities['y']/pos['x']))
velocities

b) 的结果为：

构建交互模拟试验

观察计算代码可以发现，其实我们得到了计算兔子在任意时刻的运动位移和速度计算代码，只要再做一点升级，就可以形成一个具有交互效果的模拟试验。

在这里我们使用 ipywidgets，假设我们要通过交互的形式来展示从0时刻开始，不同时间情况下，兔子位置轨迹，那么可以用以下方法获得。

%matplotlib widget

import matplotlib.pyplot as plt
import ipywidgets as widgets

fig, ax = plt.subplots(figsize=(6, 4), dpi=96)
ax.grid(True)
ax.set_xlim([-500, 150])
ax.set_ylim([-150, 150])
ax.set_xlabel('x/m')
ax.set_ylabel('y/m')
ax.set_title('the displacement diagram of a rabbit')

def simulate_rabbit(t1):
    """模拟实验
"""
    t0 = 0
    dt = 1
    n = int((t1-t0)/dt)

    return timestep_pos(n, t0, dt)

@widgets.interact(t1=(0, 100))
def update(t1=30):
    [l.remove() for l in ax.lines]
    pos = simulate_mov(t1)
    x = pos[:, 0]
    y = pos[:, 1]
    ax.plot(x, y, color='C0')

结果为：

如果你在jupyterlab 的环境实现了以上代码，只要滑动 t1 即可轻松演示这只奇葩兔子的位移轨迹。

小结

本例给出了使用 py 函数、numpy数组和Pandas 进行位移计算的基本例子，其目标是为了演示如何将数理符号的数学形式化公式方便的转换成易读和易用的py代码，最终给出了一个在 jupyterlab 或 jupyter notebook 中的交互式模型的实现。

真实的数理模拟虽然复杂，但是，其业务逻辑其实也和这个简单例子相近，数理模型 -> 代码实现 -> 可视化展示，一篇博文肯定是无法详解的。

但我们可以看到，使用 py 的显著优势，极大的减低工作量。在可读性方面，代码甚至也可以形如散文。同时，提升输出和演示效果。

如果你仔细阅读你会发现以上计算，对于三维向量的情况也是适用的。只需要在r ⃗ ( t ) \vec {r}(t)r (t ) 中增加一个z分量即可。

后期如有时间，将逐步给出 CP 计算中经常使用的数值计算、可视化优化、误差分析和算法分析，以及 OOP 风格的 CP 设计代码。

扩展阅读： 聊一聊 python 程序的理想化设计流程

如果以上代码对您有一定的益处，请多多点赞；若有不同意见，欢迎来怼：）

Original: https://blog.csdn.net/viola_aoitech/article/details/127895453
Author: viola_aoitech
Title: 利用 numpy 实现物理计算–物理向量符号与numpy数组的对应 ( jupyterlab 例子)