「SymPy」实战之Maxwell分布律分子最概然、均方根与平均速率
目录
- 0 导言
- 1 Maxwell速度分布律
- 2 SymPy推导特征速度
- 2.1 导入库和方法
- 2.2 初始化输出
- 2.3 定义符号
- 2.4 设置符号替换
- 2.5 定义分布律
- 2.6 分布律作图
- 2.7 最概然速率
- 2.8 数学平均速率
- 2.9 均方根速率
- 2.A 做个积分
- 3 N2\mathrm N_2N2与H2\mathrm H_2H2的速率分布
- 4 小结
- 5 源代码
0 导言
Python之符号运算库SymPy的讲解系列已经包含了大部分常用的符号运算和推导方法，现在我们拿分子动理论——Maxwell平衡态速度分布律来练练手，巩固一下知识。本讲会用到以下知识，忘记的小伙伴可以先复习一下：
传送链接：
「SymPy」符号运算(1) 简介/符号/变量/函数/表达式/等式/不等式/运算符
「SymPy」符号运算(2) 各种形式输出、表达式的化简合并与展开
「SymPy」符号运算(3) (非)线性方程(组)求解、数列求和、连乘、求极限
「SymPy」符号运算(4) 微积分与有限差分
1 Maxwell速度分布律
Maxwell于1859年推导出了分子速率的分布公式，后来Boltzmann用统计力学的方法也得到了相同的公式1 。对于做分子动理论和介观数值方法的工作者来说，Maxwell平衡态速度分布是相当重要的，掌握其推导、演变有利于深入理解气体动理学本质和拓展模型方法。话不多说，上公式：
f(v)=4π(m2πkT)1.5exp(−mv22kT)v2f(v) = 4\pi\left( \dfrac{m}{2\pi kT} \right)^{1.5} \exp{\left(\dfrac{-mv^2}{2kT} \right)} v^2 f(v)=4π(2πkTm)1.5exp(2kT−mv2)v2
Maxwell速率分布函数描述了气体分子关于速率的分布概率，f(v)f(v)f(v)可以看作数密度，f(v)dv=dNv/Nf(v)\mathrm d v = \mathrm d N_v / Nf(v)dv=dNv/N，其中分子速率v=vx2+vy2+vz2v = \sqrt{v_x^2 + v_y^2 + v_z^2}v=vx2+vy2+vz2，根据该速率分布函数，可以推导出分子的数学平均速率、均方根速率和最概然（最可几）速率，我们接下来就用Python的符号运算库SymPy进行推导。
先让大家看看输出效果（Anaconda→Jupyter编辑器）
2 SymPy推导特征速度
2.1 导入库和方法
import sympy
from IPython.display import display, Math
from sympy.plotting import plot
from sympy import Rational, pi, pprint, oo
其中display函数是为了在控制台打印经过渲染后的LaTeX\LaTeX{}LATEX数学公式，如果无法正常使用的话可以直接使用SymPy的pprint()函数用字符串输出公式，效果会差一点。也可以使用print(sympy.latex(<expr>))打印输出LaTeX\LaTeX{}LATEX公式的源代码，放进支持LaTeX\LaTeX{}LATEX编译的编辑器中查看。
plot函数是SymPy内的函数绘图库，后端默认是以matplotlib库为基础的2，可以直接以符号作图，免去了离散数据的步骤，自动赋予坐标轴名称，但是功能相对简单，这里作演示用还是很方便的。
Rational之前没有讲过，该函数是获得有理数的方法，在计算过程中不会将这种有理数进行小数展开运算，使输出结果更简洁清晰。
2.2 初始化输出
sympy.init_printing()
可以根据控制台的背景优化SymPy输出的LaTeX\LaTeX{}LATEX渲染公式。
2.3 定义符号
## Define constants/variables
Av, kB = sympy.symbols('A_v k_B', positive = True)
T = sympy.symbols('T', positive = True)
m = sympy.symbols('m', positive = True) # molecular mass
v = sympy.symbols('v', positive = True) # \sqrt{v_x^2 + v_y^2 + v_z^2}
M = 32 # Mole mass, take O2 (oxygen) for example
上面的代码定义了Avogadro常数、Boltzmann常数、温度、分子质量和速度的符号，并且对正负进行了限制。
代码以氧气分子举例，摩尔质量取32g/mol32\mathrm{g/mol}32g/mol。
2.4 设置符号替换
在演算过程中，我们既希望获得通用的符号表达式，又希望能过针对具体分子获得具体的值，这时定义一个符号与值的替换字典，结合SymPy表达式的.subs()函数，可以方便地做到这一点：
## Parameters substitution dict
subsdict = {m: M / 6.02214076E+23 * 1E-3, # g -> kgAv: 6.02214076E+23, # Avogadro constantkB: 1.380649E-23, # Boltzmann constantT: 298, # Temperature}
这里定义了温度为298K298\mathrm K298K。
2.5 定义分布律
直接给出Maxwell分布律，暂时不管分布律是怎么来的，先拿来用。
## Maxwell distribution
lam = m / (2 * kB * T)
f = 4*pi * (lam / pi)**Rational(1.5) * sympy.exp(- lam * v**2) * v**2
fs = f.subs(subsdict)
print("** Maxwell分布:")
# pprint(f) # 字符串显示模式
display(f) # 控制台渲染模式
pprint(f)和display(f)是两种输出方式，选用即可。
fs是根据刚才定义的符号与值的替换字典进行替换，只保留了未知的符号v。
输出结果：
2m32v2e−mv22TkBπT32kB32\frac{\sqrt{2} m^{\frac{3}{2}} v^{2} e^{- \frac{m v^{2}}{2 T k_{B}}}}{\sqrt{\pi} T^{\frac{3}{2}} k_{B}^{\frac{3}{2}}} πT23kB232m23v2e−2TkBmv2
结果比较清晰，但是SymPy不会根据指数的基进行自动合并。
值得一提的是，在之前讲指数式合并和化简时，提到了一个函数.powsimp(<expr>)，它有一个参数设置combine = ‘base’，并且进行强制化简force = True时，会将指数相同的基进行合并：abcb=(ac)ba^bc^b = (ac)^babcb=(ac)b，但是现在这里加入了分式，表达式更加复杂，按基进行指数式合并不起作用（还会继续探索，有知道的大佬望不吝赐教）。
2.6 分布律作图
## Draw velocity distribution
p1 = plot(fs, (v, 0, 50), show = False)
p1.show()
绘制结果：
可以看到SymPy在符号绘图时，图片的采样密度比较低，这时可以通过设置plot里的参数adaptive = False, nb_of_points=500使之取消自适应采样，手动设置采样点数，曲线会更加光滑。
2.7 最概然速率
所谓最概然速率（最可几速率）就是在给定条件下分子最有可能出现的速率，或者一堆分子中速率相同的分子数最多的那个速率，就是求vmv_mvm使f(vm)=max{f(v)}f(v_m) = \max\{f(v)\}f(vm)=max{f(v)}，就是求df(v)/dv=0\mathrm d f(v) / \mathrm d v = 0df(v)/dv=0的零点。
# most probable rate
vmst = sympy.solve(sympy.Eq(sympy.diff(f, v), 0), v)
print("\n** 最概然速率(表达式): ")
# pprint(vmst[0])
display(vmst[0])
print(f"\n** 最概然速率(值): {vmst[0].subs(subsdict).n(5)} m/s")
其中n.()函数表示控制输出的精度。
输出最概然速率的符号表达式：
2TkBm\frac{\sqrt{2} \sqrt{T} \sqrt{k_{B}}}{\sqrt{m}} m2TkB
简单变化，就是我们常见的形式
2kBTm\sqrt{\dfrac{2k_B T}{m}} m2kBT
最概然速率的值：393.52m/s393.52 \mathrm{m/s}393.52m/s，对应的概率f(v)=0.0021097f(v) = 0.0021097f(v)=0.0021097
2.8 数学平均速率
平均速率的求法：va=∫0∞vf(v)dvv_a = \int_0^{\infty} vf(v) \mathrm d vva=∫0∞vf(v)dv：
# average velocity
vave = sympy.integrate(v * f, (v, 0, oo))
print("\n** 平均速率(表达式): ")
# pprint(vave)
display(vave)
print(f"\n** 平均速率(值): {vave.subs(subsdict).n(5)} m/s")
输出：
22TkBπm\frac{2 \sqrt{2} \sqrt{T} \sqrt{k_{B}}}{\sqrt{\pi} \sqrt{m}} πm22TkB
简单手动化简：
8kBTπm\sqrt{\dfrac{8k_B T}{\pi m}} πm8kBT
数值为：444.04m/s444.04 \mathrm{m/s}444.04m/s
2.9 均方根速率
根据压力的微观统计表达（p=13nmu2p = \frac{1}{3}nmu^2p=31nmu2）是可以推出均方根速率的，也是Maxwell推导过程中的重要一环。根据Maxwell倒退均方根速率：u=∫0∞v2f(v)dvu = \sqrt{\int_0^{\infty} v^2 f(v) \mathrm dv}u=∫0∞v2f(v)dv
# velocity of mean square root
vrms = sympy.sqrt(sympy.integrate(v**2 * f, (v, 0, oo)))
print("\n** 均方根速率(表达式): ")
# pprint(vrms)
display(vrms)
print(f"\n** 均方根速率(值): {vrms.subs(subsdict).n(5)} m/s")
输出：
3TkBm\frac{\sqrt{3} \sqrt{T} \sqrt{k_{B}}}{\sqrt{m}} m3TkB
简单化简：
3kBTm\sqrt{\dfrac{3k_B T}{m}} m3kBT
数值：481.96m/s481.96 \mathrm{m/s}481.96m/s
2.A 做个积分
由于Maxwell速率分布是一种概率分布，因此对f(v)f(v)f(v)从000到+∞+\infty+∞积分应该为111：
# integration of f(v) over (0, +oo)
Integ = sympy.integrate(f, (v, 0, oo))
print(f"\n** f(v)在(0, +oo)积分为: {Integ.n(2)}")
输出：1.01.01.0
3 N2\mathrm N_2N2与H2\mathrm H_2H2的速率分布
复现一下《物理化学》(第五版, 傅献彩, 高等教育出版社) 中28页的图1.8：分子速率分布曲线与温度及分子质量的关系图。
只需要简单重复地定一下替换字典subsdict里面的参数就可以了，直接上结果：
变化趋势对得上，但是在概率极值和速率分布上面和书里有一点微妙的区别。暂时没有做深究，发现问题的大佬望不吝赐教~
4 小结
代码看上去很复杂，但大部分是为了美化输出，核心求解的代码并不多。用好SymPy还是能对自己推导公式进行辅助和验证，对科研有很大助力的。之后我还会根据科研中的学习和探索，充分发挥SymPy的功能，继续应用SymPy做一些验证和推导，敬请关注~
5 源代码
以下是程序源代码，在Anacona中的Spyder编译器可以正常运行，供大家参考批评。
'''
利用SymPy符号运算，根据Maxwell速率分布律，推导最概然、均方根、数学平均速率
针对特定气体，计算(低压、较高温度)下分子的最概然、均方根、数学平均速率数值
'''import sympy
from IPython.display import display, Math
from sympy.plotting import plot
from sympy import Rational, pi, pprint, oosympy.init_printing()
print(f"{' For Oxygen (M=32 g/mol) ':-^80}")#------------------------------------------------------------------------------
## Define constants/variables
Av, kB = sympy.symbols('A_v k_B', positive = True)
T = sympy.symbols('T', positive = True)
m = sympy.symbols('m', positive = True) # molecular mass
v = sympy.symbols('v', positive = True) # \sqrt{v_x^2 + v_y^2 + v_z^2}
M = 32 # Mole mass, take O2 (oxygen) for example#------------------------------------------------------------------------------
## Parameters substitution dict
subsdict = {m: M / 6.02214076E+23 * 1E-3, # g -> kgAv: 6.02214076E+23, # Avogadro constantkB: 1.380649E-23, # Boltzmann constantT: 298, # Temperature}#------------------------------------------------------------------------------
## Maxwell distribution
lam = m / (2 * kB * T)
f = 4*pi * (lam / pi)**Rational(1.5) * sympy.exp(- lam * v**2) * v**2
fs = f.subs(subsdict)
print("** Maxwell分布:")
# pprint(f) # 字符串显示模式
display(f) # 控制台渲染模式#------------------------------------------------------------------------------
## Quantative relations
# most probable rate
vmst = sympy.solve(sympy.Eq(sympy.diff(f, v), 0), v)
print("\n** 最概然速率(表达式): ")
# pprint(vmst[0])
display(vmst[0])
print(f"\n** 最概然速率(值): {vmst[0].subs(subsdict).n(5)} m/s")fmax = sympy.maximum(f, v, sympy.Interval.open(0, oo))
print("\n** 最概然速率对应概率(表达式):")
# pprint(fmax)
display(fmax)
print(f"\n** 最概然速率对应概率(值): {fmax.subs(subsdict).n(5)}")# average velocity
vave = sympy.integrate(v * f, (v, 0, oo))
print("\n** 平均速率(表达式): ")
# pprint(vave)
display(vave)
print(f"\n** 平均速率(值): {vave.subs(subsdict).n(5)} m/s")# velocity of mean square root
vrms = sympy.sqrt(sympy.integrate(v**2 * f, (v, 0, oo)))
print("\n** 均方根速率(表达式): ")
# pprint(vrms)
display(vrms)
print(f"\n** 均方根速率(值): {vrms.subs(subsdict).n(5)} m/s")# integration of f(v) over (0, +oo)
Integ = sympy.integrate(f, (v, 0, oo))
print(f"\n** f(v)在(0, +oo)积分为: {Integ.n(2)}")#------------------------------------------------------------------------------
## Draw velocity distribution
p1 = plot(fs, (v, 0, 50), show = False)
# p1.show()#------------------------------------------------------------------------------
print(f"\n{' End ':-^80}")
N2\mathrm N_2N2与H2\mathrm H_2H2的速率分布绘图的代码块：
'''
以运行上一个Cell为基础
复现《物理化学》(第五版, 傅献彩, 高等教育出版社)中28页的图1.8
分子速率分布曲线与温度及分子质量的关系图
'''# N2, 100K
subsdict = {m: 28 / 6.02214076E+23 * 1E-3, # g -> kgAv: 6.02214076E+23, # Avogadro constantkB: 1.380649E-23, # Boltzmann constantT: 100, # Temperature}
p1 = plot(f.subs(subsdict), (v, 0, 2500), show = False, label = "N2(100K)", legend = True, size = (6, 4),adaptive = False, nb_of_points=500)# N2, 300K
subsdict = {m: 28 / 6.02214076E+23 * 1E-3, # g -> kgAv: 6.02214076E+23, # Avogadro constantkB: 1.380649E-23, # Boltzmann constantT: 300, # Temperature}
p2 = plot(f.subs(subsdict), (v, 0, 2500), show = False, label = "N2(300K)",adaptive = False, nb_of_points=500)# H2, 100K
subsdict = {m: 2 / 6.02214076E+23 * 1E-3, # g -> kgAv: 6.02214076E+23, # Avogadro constantkB: 1.380649E-23, # Boltzmann constantT: 100, # Temperature}
p3 = plot(f.subs(subsdict), (v, 0, 2500), show = False, label = "H2(100K)",adaptive = False, nb_of_points=500)# H2, 300K
subsdict = {m: 2 / 6.02214076E+23 * 1E-3, # g -> kgAv: 6.02214076E+23, # Avogadro constantkB: 1.380649E-23, # Boltzmann constantT: 300, # Temperature}
p4 = plot(f.subs(subsdict), (v, 0, 2500), show = False, label = "H2(300K)",adaptive = False, nb_of_points=500)p1.append(p2[0])
p1.append(p3[0])
p1.append(p4[0])p1.show()
傅献彩, 沈文霞，姚天扬等. 物理化学 (第五版).上册[M]. 人民教育出版社, 2005. ↩︎
Meurer A, Smith CP, Paprocki M, Čertík O, Kirpichev SB, Rocklin M, Kumar A, Ivanov S, Moore JK, Singh S, Rathnayake T, Vig S, Granger BE, Muller RP, Bonazzi F, Gupta H, Vats S, Johansson F, Pedregosa F, Curry MJ, Terrel AR, Roučka Š, Saboo A, Fernando I, Kulal S, Cimrman R, Scopatz A. (2017) SymPy: symbolic computing in Python. PeerJ Computer Science 3:e103 https://doi.org/10.7717/peerj-cs.103 ↩︎
「SymPy」实战之Maxwell分布律分子最概然、均方根与平均速率相关推荐
- 「SymPy」符号运算(5) Vector向量及运算
目录 导言 坐标系和向量 笛卡尔坐标系与向量 坐标原点 向量运算 四则运算 点乘.叉乘 并矢/外积(Dyadic) ∇ \nabla ∇算子(Hamiltonian算子) 梯度 旋度 散度 向量操作 ...
- 「SymPy」符号运算(4) 微积分与有限差分
目录 导言 积分 不定积分 定积分 多重积分 求导 一阶导数 高阶导数 偏导数 有限差分 常微分差分 差分系数 高阶差分 偏微分差分 导言 在前几篇中,我们学习了SymPy的基本语法.方程求解等基础知 ...
- 「SymPy」符号运算(2) 各种形式输出、表达式的化简合并与展开
目录 导言 输出 替换.演化 化简.合并与展开 化简 展开 合并 `cancel`函数 `apart`函数 `rewrite`函数 `expand_func`函数 导言 在前一篇文章中,我们简单学习了 ...
- 「SymPy」符号运算(6) 矩阵Matrix及基础运算
目录 导言 创建矩阵 列表初始化 行向量 列向量 维度和数集 二元函数 `lambda`函数 特殊矩阵 基本操作 索引 增删 基础运算 向量运算 导言 在前几篇文章中,我们学习了SymPy基础/高级用 ...
- 「SymPy」符号运算(3) (非)线性方程(组)求解、数列求和、连乘、求极限
目录 导言 解方程(组) solve函数 solveset函数 求和 ∑ \sum ∑ 连乘 ∏ \prod ∏ 求函数极限 求数列极限 导言 在前两篇文章中,我们学习了SymPy的输入输出.基本符号 ...
- 「SymPy」符号运算(7) Matrix进阶运算与线性代数
目录 0 导言 1 矩阵一般属性 秩 逆 迹 转置 共轭 伴随 行列式 特征值 特征向量 特征多项式 2 矩阵运算与线性代数 范数 标准化 条件数 矩阵分解 黑塞矩阵 雅克比矩阵 Jordan标准型 ...
- 「SymPy」符号运算(1) 简介/符号/变量/函数/表达式/等式/不等式/运算符
目录 1 简介 2 导入库 3 定义符号 4 定义函数 5 表达式 6 等式/不等式 7 `SymPy`假设与限制 8 运算符/函数 常用运算符/函数 数学常数 三角函数 复杂函数 指数运算 1 简介 ...
- Java应用诊断工具】「BTrace」基本概念和原理的介绍（1）
Java应用诊断工具]「BTrace」基本概念和原理的介绍(1) 系列文章 「BTrace」基本概念和初步介绍(1) 「BTrace」安装介绍和使用原理(2)- 未完成 「BTrace」实战代码进行调 ...
- 利用Spring扩展点模拟MyBatis的注解编程「知识点多多」「扩展点实战系列」- 第448篇
历史文章(文章累计440+) <国内最全的Spring Boot系列之一> <国内最全的Spring Boot系列之二> <国内最全的Spring Boot系列之三> ...
最新文章
- Mysql 事务中Update 会锁表吗？
- bind日志配置详解
- 调整屏幕亮度，调整字体大小
- tomact+apache实现web网页动静结合
- 模拟退火算法解决TSP（python实现 110+行代码）【gif生成】
- 1135 Is It A Red-Black Tree (30 分)【难度: 难 / 知识点: 红黑树 未完成】
- 图解10大机器学习算法
- [网文摘录]云计算平台管理
- Win7系统删除微软拼音
- mac brew 测速 软件_敏捷过程中的软件持续建模
- js 通过图片URL地址将图片转为可操作的File文件对象
- LM393 电压比较器及其典型电路介绍
- 【工具篇】使用OpenCV播放视频并截取图片
- 景联文科技提供步态数据采集服务、提供21000个id步态视频训练数据集
- 小白入门之海康威视摄像机的二次开发
- 像西方知识分子那样登场
- 英文歌曲：my heart will go on(我心永恒)
- Redis数据莫名其妙全部丢失
- 使用了未经检查或不安全的操作
- English--基础知识点--9--used