From 5b00715f84012809904f9ecb5c7fb33b8e9cf0f7 Mon Sep 17 00:00:00 2001 From: guanjihuan Date: Wed, 27 Aug 2025 17:24:32 +0800 Subject: [PATCH] Create one-dimensional_infinite_square_well_using_finite_difference_method.py --- ...are_well_using_finite_difference_method.py | 117 ++++++++++++++++++ 1 file changed, 117 insertions(+) create mode 100644 2025.08.27_one-dimensional_infinite_square_well_using_finite_difference_method/one-dimensional_infinite_square_well_using_finite_difference_method.py diff --git a/2025.08.27_one-dimensional_infinite_square_well_using_finite_difference_method/one-dimensional_infinite_square_well_using_finite_difference_method.py b/2025.08.27_one-dimensional_infinite_square_well_using_finite_difference_method/one-dimensional_infinite_square_well_using_finite_difference_method.py new file mode 100644 index 0000000..2a4f4c7 --- /dev/null +++ b/2025.08.27_one-dimensional_infinite_square_well_using_finite_difference_method/one-dimensional_infinite_square_well_using_finite_difference_method.py @@ -0,0 +1,117 @@ +""" +This code is supported by the website: https://www.guanjihuan.com +The newest version of this code is on the web page: https://www.guanjihuan.com/archives/47514 +""" + +import numpy as np +import matplotlib.pyplot as plt +import scipy + +plt.rcParams['font.sans-serif'] = ['SimHei'] # 在画图中显示中文 +plt.rcParams['axes.unicode_minus'] = False # 显示中文后，加上这个才可以正常显示负号 + +# 物理常数定义 +hbar = 1.0545718e-34 # 约化普朗克常数 (J·s) +m = 9.109e-31 # 电子质量 (kg) +a = 1e-9 # 势阱宽度 1nm (m) + +# 解析波函数 +def infinite_well_wavefunction(n, x, a): + """返回第n能级的波函数""" + return np.sqrt(2/a) * np.sin(n * np.pi * x / a) + +# 解析能级 +def infinite_well_energy(n, a, m, hbar): + """返回第n能级的能量""" + return (n**2 * np.pi**2 * hbar**2) / (2 * m * a**2) + +def solve_infinite_well_numerically(a=1e-9, N=3000, num_states=5): + """ + 有限差分数值求解一维无限深势阱的薛定谔方程 + 参数: + a: 势阱宽度 + N: 网格点数 + num_states: 要计算的本征态数量 + """ + # 离散化参数 + x_array = np.linspace(0, a, N) + dx = x_array[1] - x_array[0] + + # 构建哈密顿矩阵：二阶导数的离散化（三点中心差分） + coefficient = -hbar**2 / (2 * m) + main_diag = -2 * np.ones(N) * coefficient / dx**2 + off_diag = np.ones(N-1) * coefficient / dx**2 + H = scipy.sparse.diags([main_diag, off_diag, off_diag], [0, -1, 1], format='csc') + + # 对于N=5个网格点的例子： + # H = coefficient/dx² × + # ⎡ -2 1 0 0 0 ⎤ + # ⎢ 1 -2 1 0 0 ⎥ + # ⎢ 0 1 -2 1 0 ⎥ + # ⎢ 0 0 1 -2 1 ⎥ + # ⎣ 0 0 0 1 -2 ⎦ + + # 求解本征值和本征向量 + eigenvalues, eigenvectors = scipy.sparse.linalg.eigs(H, k=num_states, which='SM') + eigenvalues = np.real(eigenvalues) + + # 排序本征值和本征向量 + idx = eigenvalues.argsort() + eigenvalues = eigenvalues[idx] + eigenvectors = eigenvectors[:, idx] + + return x_array, eigenvalues, eigenvectors + +# 数值求解 +x_array, E_num, psi_num = solve_infinite_well_numerically() +E_num_ev = E_num / 1.602e-19 # J 单位转换为 eV 单位 + +print("数值求解结果:") +for i, E in enumerate(E_num_ev): + print(f"状态 {i+1}: E = {E:.3f} eV") + +print("解析求解结果:") +for i in range(5): + n = i + 1 + E_analytical = infinite_well_energy(n, a, m, hbar) / 1.602e-19 + error = abs(E_num_ev[i] - E_analytical) / E_analytical * 100 + print(f"状态 {i+1}: 解析解 E = {E_analytical:.4f} eV, 数值求解误差 Error = {error:.2f}%") + +# 绘制数值求解的能级图 +plt.figure(figsize=(8, 6)) +for n, E in enumerate(E_num_ev[:6], 1): + plt.axhline(y=E, color='r', linestyle='-', linewidth=2) # 主要的画图命令 + plt.text(a*1e9/2, E + 0.1, f'n={n}', ha='center', va='bottom', fontsize=12) +plt.xlabel('势阱位置') +plt.ylabel('能量 (eV)') +plt.title('一维无限深势阱的能级（数值解）') +plt.grid(True, alpha=0.3) +plt.ylim(0, max(E_num_ev[:6]) + 1) +plt.show() + +# 绘制数值解与解析解的比较 +plt.figure(figsize=(12, 8)) +for i in range(4): + # 数值解 + psi_numerical = np.real(psi_num[:, i]) + norm = np.sqrt(np.trapz(psi_numerical**2, x_array)) + psi_numerical /= norm # 归一化数值解 + + # 确保符号一致（本征向量符号可能任意） + if np.dot(psi_numerical, infinite_well_wavefunction(i+1, x_array, a)) < 0: + psi_numerical *= -1 + + # 解析解 + psi_analytical = infinite_well_wavefunction(i+1, x_array, a) + + plt.subplot(2, 2, i+1) + plt.plot(x_array*1e9, psi_numerical, 'b-', label='数值解', linewidth=2) + plt.plot(x_array*1e9, psi_analytical, 'r--', label='解析解', linewidth=2) + plt.xlabel('位置 x (nm)') + plt.ylabel('波函数') + plt.title(f'n = {i+1}, E_num = {E_num_ev[i]:.2f} eV') + plt.legend() + plt.grid(True, alpha=0.3) + +plt.tight_layout() +plt.show() \ No newline at end of file