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

要点

Python物理学差分数值和符号计算
热动力学计算：统计力学，分子动力学模型
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)/6x += 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)/6x += 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()