专注旅游网站网站开发,wordpress saml,河南手机网站建设公司排名,国风网页设计欣赏要点 Python物理学差分数值和符号计算 热动力学计算#xff1a;统计力学#xff0c;分子动力学模型 Python寻找弹性物体的运动#xff0c;LAMMPS 分子动力学模拟器模拟2D气体分子#xff0c;Python原子模拟绘图#xff0c;Python数值计算原子平衡性#xff0c;Python绘制…要点 Python物理学差分数值和符号计算 热动力学计算统计力学分子动力学模型 Python寻找弹性物体的运动LAMMPS 分子动力学模拟器模拟2D气体分子Python原子模拟绘图Python数值计算原子平衡性Python绘制平衡时原子波动。 MATLAB随机速度原子晶格编辑读写lammps轨迹文件函数。使用LAMMPSMATLAB和Python 二维模拟 Lennard-Jones 系统。 Python模拟理想气体热动力学结果体积压力和温度。 LAMMPS模拟固体原子数据Python基于模拟测量模型左侧粒子的动能。Python计算绘制爱因斯坦晶体宏态的多重性。MATLAB模拟爱因斯坦晶体热运动状态。Python计算两个耦合的晶体熵。MATLAB计算绘制二维晶体热分子运动。 Python数值计算绘图爱因斯坦晶体在微规范系统中谐振子运动概率。Python蒙特卡洛方法模拟亥姆霍兹自由能。Python使用 命中和错过 算法符号计算蒙特卡洛积分估计。MATLAB对比Python 计算热浴的蒙特卡洛伊辛模型。 LAMMPS 模拟低温下液相和气相之间的聚结和相共存Python基于模拟数值计算。 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(x) k(x)它是薛定谔方程的重新设计版本用于求解变量 k k k以及我们的 Ψ \Psi Ψ 方程它将由 psione ( x ) (x) (x) 和 psitwo ( x ) (x) (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(rk1/2,xh/2)k3 h*deriv(rk2/2,xh/2)k4 h*deriv(rk3,xh)r (k12*k22*k3k4)/6x h 在这里循环几乎涵盖了龙格-库塔的整个过程。通过使用由 k k k 值定义的斜率近似值每个 k k k 值有助于近似下一个斜率使我们更接近求解 f ( x ) f(x) f(x)。此外在获得每个斜率之后我们获得加权平均值并使用这些新值更新我们的数组为下一次迭代做好准备。这个过程将在 x \mathrm{x} 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(rk1/2,xh/2)k3 h*deriv(rk2/2,xh/2)k4 h*deriv(rk3,xh)r (k12*k22*k3k4)/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; v02; 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()参阅一计算思维
参阅二亚图跨际
