「SymPy」实战之Maxwell分布律分子最概然、均方根与平均速率

admin • • 人工智能

目录


在这里插入图片描述

0 导言

Python之符号运算库SymPy的讲解系列已经包含了大部分常用的符号运算和推导方法，现在我们拿分子动理论——Maxwell平衡态速度分布律来练练手，巩固一下知识。本讲会用到以下知识，忘记的小伙伴可以先复习一下：

传送链接：

「SymPy」符号运算(1) 简介/符号/变量/函数/表达式/等式/不等式/运算符

「SymPy」符号运算(2) 各种形式输出、表达式的化简合并与展开

「SymPy」符号运算(3) (非)线性方程(组)求解、数列求和、连乘、求极限

「SymPy」符号运算(4) 微积分与有限差分

1 Maxwell速度分布律

Maxwell于1859年推导出了分子速率的分布公式，后来Boltzmann用统计力学的方法也得到了相同的公式1 。对于做分子动理论和介观数值方法的工作者来说，Maxwell平衡态速度分布是相当重要的，掌握其推导、演变有利于深入理解气体动理学本质和拓展模型方法。话不多说，上公式：

f

(

v

)

=

4

π

(

m

2

π

k

T

)

1.5

exp

⁡

(

−

m

v

2

2

k

T

)

v

2

f(v) = 4pileft( dfrac{m}{2pi 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

)

d

v

=

d

N

v

/

N

f(v)mathrm d v = mathrm d N_v / N

f(v)dv=dNv​/N，其中分子速率

v

=

v

x

2

+

v

y

2

+

v

z

2

v = 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{}

LATE​X数学公式，如果无法正常使用的话可以直接使用SymPy的pprint()函数用字符串输出公式，效果会差一点。也可以使用print(sympy.latex(<expr>))打印输出

LaTeX

LaTeX{}

LATE​X公式的源代码，放进支持

LaTeX

LaTeX{}

LATE​X编译的编辑器中查看。

plot函数是SymPy内的函数绘图库，后端默认是以matplotlib库为基础的2，可以直接以符号作图，免去了离散数据的步骤，自动赋予坐标轴名称，但是功能相对简单，这里作演示用还是很方便的。

Rational之前没有讲过，该函数是获得有理数的方法，在计算过程中不会将这种有理数进行小数展开运算，使输出结果更简洁清晰。

2.2 初始化输出

sympy.init_printing()

可以根据控制台的背景优化SymPy输出的

LaTeX

LaTeX{}

LATE​X渲染公式。

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常数、温度、分子质量和速度的符号，并且对正负进行了限制。

代码以氧气分子举例，摩尔质量取

32

g

/

m

o

l

32mathrm{g/mol}

32g/mol。

2.4 设置符号替换

在演算过程中，我们既希望获得通用的符号表达式，又希望能过针对具体分子获得具体的值，这时定义一个符号与值的替换字典，结合SymPy表达式的.subs()函数，可以方便地做到这一点：

## Parameters substitution dict
subsdict = {m: M / 6.02214076E+23 * 1E-3,     # g -> kg
            Av: 6.02214076E+23,               # Avogadro constant
            kB: 1.380649E-23,                 # Boltzmann constant
            T: 298,                           # Temperature
            }

这里定义了温度为

298

K

298mathrm K

298K。

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。

输出结果：

2

m

3

2

v

2

e

−

m

v

2

2

T

k

B

π

T

3

2

k

B

3

2

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}}}

π

​T23​kB23​​2

​m23​v2e−2TkB​mv2​​
结果比较清晰，但是SymPy不会根据指数的基进行自动合并。

值得一提的是，在之前讲指数式合并和化简时，提到了一个函数.powsimp(<expr>)，它有一个参数设置combine = ‘base’，并且进行强制化简force = True时，会将指数相同的基进行合并：

a

b

c

b

=

(

a

c

)

b

a^bc^b = (ac)^b

abcb=(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 最概然速率

所谓最概然速率（最可几速率）就是在给定条件下分子最有可能出现的速率，或者一堆分子中速率相同的分子数最多的那个速率，就是求

v

m

v_m

vm​使

f

(

v

m

)

=

max

⁡

{

f

(

v

)

}

f(v_m) = max{f(v)}

f(vm​)=max{f(v)}，就是求

d

f

(

v

)

/

d

v

=

0

mathrm d f(v) / mathrm d v = 0

df(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.()函数表示控制输出的精度。

输出最概然速率的符号表达式：

2

T

k

B

m

frac{sqrt{2} sqrt{T} sqrt{k_{B}}}{sqrt{m}}

m

​2

​T

​kB​

​​
简单变化，就是我们常见的形式

2

k

B

T

m

sqrt{dfrac{2k_B T}{m}}

m2kB​T​

​
最概然速率的值：

393.52

m

/

s

393.52 mathrm{m/s}

393.52m/s，对应的概率

f

(

v

)

=

0.0021097

f(v) = 0.0021097

f(v)=0.0021097

2.8 数学平均速率

平均速率的求法：

v

a

=

∫

0

∞

v

f

(

v

)

d

v

v_a = int_0^{infty} vf(v) mathrm d v

va​=∫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")

输出：

2

2

T

k

B

π

m

frac{2 sqrt{2} sqrt{T} sqrt{k_{B}}}{sqrt{pi} sqrt{m}}

π

​m

​22

​T

​kB​

​​
简单手动化简：

8

k

B

T

π

m

sqrt{dfrac{8k_B T}{pi m}}

πm8kB​T​

​
数值为：

444.04

m

/

s

444.04 mathrm{m/s}

444.04m/s

2.9 均方根速率

根据压力的微观统计表达（

p

=

1

3

n

m

u

2

p = frac{1}{3}nmu^2

p=31​nmu2）是可以推出均方根速率的，也是Maxwell推导过程中的重要一环。根据Maxwell倒退均方根速率：

u

=

∫

0

∞

v

2

f

(

v

)

d

v

u = 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")

输出：

3

T

k

B

m

frac{sqrt{3} sqrt{T} sqrt{k_{B}}}{sqrt{m}}

m

​3

​T

​kB​

​​
简单化简：

3

k

B

T

m

sqrt{dfrac{3k_B T}{m}}

m3kB​T​

​
数值：

481.96

m

/

s

481.96 mathrm{m/s}

481.96m/s

2.A 做个积分

由于Maxwell速率分布是一种概率分布，因此对

f

(

v

)

f(v)

f(v)从

0

0

0到

+

∞

+infty

+∞积分应该为

1

1

1：

# 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.0

1.0

1.0

3

N

2

mathrm N_2

N2​与

H

2

mathrm H_2

H2​的速率分布

复现一下《物理化学》(第五版, 傅献彩, 高等教育出版社) 中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, oo

sympy.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 -> kg
            Av: 6.02214076E+23,               # Avogadro constant
            kB: 1.380649E-23,                 # Boltzmann constant
            T: 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}")

N

2

mathrm N_2

N2​与

H

2

mathrm H_2

H2​的速率分布绘图的代码块：

'''
以运行上一个Cell为基础
复现《物理化学》(第五版, 傅献彩, 高等教育出版社)中28页的图1.8
分子速率分布曲线与温度及分子质量的关系图
'''

# N2, 100K
subsdict = {m: 28 / 6.02214076E+23 * 1E-3,     # g -> kg
            Av: 6.02214076E+23,                # Avogadro constant
            kB: 1.380649E-23,                  # Boltzmann constant
            T: 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 -> kg
            Av: 6.02214076E+23,                # Avogadro constant
            kB: 1.380649E-23,                  # Boltzmann constant
            T: 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 -> kg
            Av: 6.02214076E+23,               # Avogadro constant
            kB: 1.380649E-23,                 # Boltzmann constant
            T: 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 -> kg
            Av: 6.02214076E+23,               # Avogadro constant
            kB: 1.380649E-23,                 # Boltzmann constant
            T: 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()

  1. 傅献彩, 沈文霞，姚天扬等. 物理化学 (第五版).上册[M]. 人民教育出版社, 2005. ↩︎

  2. 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 ↩︎

本图文内容来源于网友网络收集整理提供，作为学习参考使用，版权属于原作者。
THE END
分享
二维码
< <上一篇
下一篇>>