🥥MATLAB和Python数值和符号计算可视化物理学气体动能和粒子速度

关键词

Python | MATLAB | 数值计算 | 符号计算 | 可视化 | 物理学 | 热力学 | 气体 | 动能 | 粒子 | 速度 | 压力 | 体积 | 平衡性 | 粒子波动 | 原子晶格 | 分子动力学模拟 | 二维 | 晶体熵 | 蒙特卡洛算法 | 亥姆霍兹自由能 | 命中和错过算法 | 伊辛模型 | 玻色-爱因斯坦分布 | 费米-狄拉克分布 | 龙格-库塔方法 | 固体 | 分子

MATLAB和Python数值和符号计算可视化物理学气体动能和粒子速度 | Notion亞圖跨際 on Notion

🏈page指点迷津 | Brief

要点

1. Python物理学差分数值和符号计算

2. 热动力学计算：统计力学，分子动力学模型

   1. Python寻找弹性物体的运动，LAMMPS 分子动力学模拟器模拟2D气体分子，Python原子模拟绘图，Python数值计算原子平衡性，Python绘制平衡时原子波动。

   2. MATLAB随机速度原子晶格，编辑读写lammps轨迹文件函数。使用LAMMPS，MATLAB和Python 二维模拟 Lennard-Jones 系统。

   3. Python模拟理想气体热动力学结果：体积，压力和温度。

   4. LAMMPS模拟固体原子数据，Python基于模拟，测量模型左侧粒子的动能。Python计算绘制爱因斯坦晶体宏态的多重性。MATLAB模拟爱因斯坦晶体热运动状态。Python计算两个耦合的晶体熵。MATLAB计算绘制二维晶体热分子运动。

   5. Python数值计算绘图爱因斯坦晶体在微规范系统中谐振子运动概率。Python蒙特卡洛方法模拟亥姆霍兹自由能。Python使用 命中和错过 算法符号计算蒙特卡洛积分估计。MATLAB对比Python 计算热浴的蒙特卡洛伊辛模型。

   6. LAMMPS 模拟低温下液相和气相之间的聚结和相共存，Python基于模拟数值计算。

   7. MATLAB绘制了玻色-爱因斯坦分布和费米-狄拉克分布。

Python物理学数值计算示例 龙格-库塔方法

鉴于计算效率，以下代码只是演示作用

龙格-库塔方法是常微分方程的数值近似，由 Carl Runge 和 Wilhelm Kutta 开发。 通过使用一个区间内的四个斜率值（不一定落在实际解上）并对斜率进行平均，可以获得一个非常好的近似解。在这个例子中，我们将重点关注四阶龙格-库塔方法来帮助我们解决一维散射问题。

为了开始我们的代码，我们将导入一些包来帮助我们进行数学和可视化。

 import cmath 
 import numpy as np 
 import matplotlib.pyplot as plt 

从这里开始，我们开始定义方程的初始参数。

 mass = 1.0 
 hbar = 1.0 
 v0 = 2.0 
 alpha = 0.5 
 E = 3.0 
 i = 1.0j 
 x = 10.0 
 xf = -10.0 
 h = -.001 
 xaxis = np.array([], float) 
 psi = np.array([], complex) 
 psiprime = np.array([], complex) 

一旦我们有了初始值，我们就开始处理定义我们将要使用的方程的函数。我们的主要方程是 k(x)，它是薛定谔方程的重新设计版本，用于求解变量 k，以及我们的 \Psi 方程，它将由 psione (x) 和 psitwo (x)定义。

 def v(x): 
     return v0/2.0 * (1.0 + np.tanh(x/alpha))
 def k(x): 
     return cmath.sqrt((2*mass/(hbar**2))*(E - v(x)))
 def psione(x): 
     return np.exp(i*k(x)*x)
 def psitwo(x): 
     return i*k(x)*np.exp(i*k(x)*x)

在此，首先，我们需要定义一个包含初始条件波函数的数组。

 r = np.array([psione(x), psitwo(x)]) 

在数组中设置这些方程，我们可以通过下面将定义的龙格-库塔方法迭代这两个方程，并让它们为我们为 psione(x) 和 psitwo(x) 定义的方程提供近似解 。 但在我们到达方程的主要部分之前，我们需要定义一个更重要的函数。

 def deriv(r,x): 
     return np.array([r[1],-(2.0*mass/(hbar**2) * (E - v(x))*r[0])], complex)

deriv 函数是龙格-库塔的输出经过的地方，该函数从数组 r 中获取我们的值，然后将其推入这些条件。 对于返回的第一个值，非常简单，我们的 x 值将被输入到数组的第二个方程中。 然而，第二个值将经历不同的处理。 这次，x 值将经历薛定谔方程的另一次迭代，其中考虑了波函数 psione(x)。

 while (x >= xf ):
         xaxis = np.append(xaxis, x)
         psi = np.append(psi, r[0])
         psiprime = np.append(psiprime, r[1])
         k1 = h*deriv(r,x)
         k2 = h*deriv(r+k1/2,x+h/2)
         k3 = h*deriv(r+k2/2,x+h/2)
         k4 = h*deriv(r+k3,x+h)
         r += (k1+2*k2+2*k3+k4)/6
         x += h 

在这里，循环几乎涵盖了龙格-库塔的整个过程。通过使用由 k 值定义的斜率近似值，每个 k 值有助于近似下一个斜率，使我们更接近求解 f(x)。此外，在获得每个斜率之后，我们获得加权平均值并使用这些新值更新我们的数组，为下一次迭代做好准备。这个过程将在 \mathrm{x}​ 轴上我们定义的范围内继续，这最终将为我们提供绘制即将求解的常微分方程所需的值。

龙格-库塔方法可以很容易地适应许多其他方程，大多数时候我们只需要调整导数函数和我们的初始条件方程。 其他示例包括摆常微分方程和行星运动常微分方程。 在下面，我们现在可以找到完整的代码以及额外的步骤，例如求解反射和透射值的函数，以及如何绘制我们的值。

 import cmath 
 import numpy as np 
 import matplotlib.pyplot as plt 
 mass = 1.0 
 hbar = 1.0 
 v0 = 2.0 
 alpha = 0.5 
 E = 3.0 
 i = 1.0j 
 x = 10.0 
 xf = -10.0 
 h = -.001 
 xaxis = np.array([], float) 
 psi = np.array([], complex) 
 psiprime = np.array([], complex) 
 def v(x): 
     return v0/2.0 * (1.0 + np.tanh(x/alpha))
 def k(x): 
     return cmath.sqrt((2*mass/(hbar**2))*(E - v(x)))
 r = np.array([psione(x), psitwo(x)]) 
 def deriv(r,x): 
     return np.array([r[1],-(2.0*mass/(hbar**2) * (E - v(x))*r[0])], complex)
 while (x >= xf ):
         xaxis = np.append(xaxis, x)
         psi = np.append(psi, r[0])
         psiprime = np.append(psiprime, r[1])
         k1 = h*deriv(r,x)
         k2 = h*deriv(r+k1/2,x+h/2)
         k3 = h*deriv(r+k2/2,x+h/2)
         k4 = h*deriv(r+k3,x+h)
         r += (k1+2*k2+2*k3+k4)/6
         x += h 
 psi1 = psi[20000]; psi2 = psiprime[20000]; x = 10; xf = -10
 def reflection(x, y):
     aa = (psi1 + psi2/(i*k(y)))/(2*np.exp(i*k(y)*y))
     bb = (psi1 - psi2/(i*k(y)))/(2*np.exp(-i*k(y)*y))
     return (np.abs(bb)/np.abs(aa))**2
 def transmission(x,y):
     aa = (psi1 + psi2/(i*k(y)))/(2.0*np.exp(i*k(y)*y))
     return k(x)/k(y) * 1.0/(np.abs(aa))**2
 print('reflection = ',reflection(x,xf))
 print('transmission = ', transmission(x,xf))
 print('r + t = ', reflection(x,xf) + transmission(x,xf))
 fig, ax = plt.subplots(1,2, figsize = (15,5))
 ax[0].plot(xaxis, psi.real, xaxis, psi.imag, xaxis, v(xaxis))
 ax[1].plot(xaxis, psiprime.real, xaxis, psiprime.imag, xaxis, v(xaxis))
 plt.show()

使用Scipy 改写为：

 import cmath
 import numpy as np
 import matplotlib.pyplot as plt
 from scipy.integrate import odeint, solve_ivp
 E = 3; m = 1; h = 1; alpha = .5; v0=2; i = 1.0j; xi = 10; xf = -10
 def v(x): return v0/2.0 * (1.0 + np.tanh(x/alpha))
 def k(x): return cmath.sqrt((2*m/(h**2))*(E - v(x)))
 def psione(x): return np.exp(i*k(x)*x)
 def psitwo(x): return i*k(x)*np.exp(i*k(x)*x)
 def deriv(x, y): return [y[1], -(2.0*m/(h**2.0) * (E - v(x))*y[0])]
 ​
 values = solve_ivp(deriv, [10, -10], [psione(xi), psitwo(xi)], first_step = .001, max_step = .001)
 psi1 = values.y[0,20000]; psi2 = values.y[1,20000]; x = 10; xf = -10
 def reflection(x, y):
     aa = (psi1 + psi2/(i*k(y)))/(2*np.exp(i*k(y)*y))
     bb = (psi1 - psi2/(i*k(y)))/(2*np.exp(-i*k(y)*y))
     return (np.abs(bb)/np.abs(aa))**2
 def transmission(x,y):
     aa = (psi1 + psi2/(i*k(y)))/(2.0*np.exp(i*k(y)*y))
     return k(x)/k(y) * 1.0/(np.abs(aa))**2
 print('reflection = ',reflection(x,xf))
 print('transmission = ', transmission(x,xf))
 print('r + t = ', reflection(x,xf) + transmission(x,xf))
 fig, ax = plt.subplots(1,2, figsize = (15,5))
 ax[0].plot(values.t, values.y[0].real, values.t, values.y[0].imag, values.t, v(values.t))
 ax[1].plot(values.t, values.y[1].real, values.t, values.y[1].imag, values.t, v(values.t))
 plt.show()