2026/9/16 14:25:32
FEATURED · 精选文章 
简介基于Matlab的一维量子和经典谐振子仿真资源包主要面向物理、电子信息工程、数学等专业的高校学生及相关科研人员适用于课程设计、期末大作业与毕业设计也适合想借助数值计算理解量子力学模型的初学者。资源共包含8个文件其中有5个脚本文件、2个说明文档和1张结果示意图压缩包仅约87KB轻量易用已有306人学习下载。程序采用参数化编程关键参数可灵活修改注释明细、结构清晰并附可直接运行的案例数据利于快速复现和二次开发。脚本涵盖厄米多项式系数计算、经典概率密度函数、量子谐振子波函数与概率分布绘制等模块可直观对比一维量子和经典谐振子的行为差异帮助读者将Matlab数值计算与物理理论结合起来。调整参数即可观察不同状态下的波形变化是教学演示和实践演练的实用工具。1. 一维谐振子用量子和经典两条路线同时读懂波函数一个粒子被束缚在势场里经典力学会告诉你它按正弦规律来回摆动量子力学则只肯给出波函数和概率分布告诉你它在某个位置出现的概率。同样是谐振子经典给出的是确定轨迹量子给出的是离散能级和波函数干涉条纹。把两者算出来、画在同一张图上立刻就能看到著名的玻尔对应原理——高量子数下量子概率分布逐步收敛为经典概率分布。本文就以“一维量子和经典谐振子的波函数和概率分布”为对象直接用 Matlab 做解析计算与绘制不改薛定谔方程求解器不做蒙特卡洛纯靠埃尔米特多项式和经典相空间密度就能在几分钟内看到双峰、零点和概率密度包络。适合正在学量子力学或数值方法、需要快速可复现波形图的本科生、研究生和科研入门人员。2. 从物理模型到 Matlab 计算方案为什么解析式比数值打靶更快2.1 量子端固有频率、自然单位与定态解一维量子谐振子的哈密顿量写作H -ℏ²/(2m) · d²/dx² (1/2)mω²x²标准的定态解是厄米多项式与高斯函数的乘积ψₙ(x) (1/√(2ⁿn!√π)) · Hₙ(√(mω/ℏ)x) · exp(-mωx²/(2ℏ))能量本征值为 Eₙ ℏω(n1/2)。很多教材把 ℏ m ω 1 定为自然单位这时波函数形式最简洁Matlab 代码也只跟无量纲的 x 打交道。实际绘图时我推荐直接采用这种约定坐标轴标成“自然单位”显得专业且能绕开繁杂的量纲转换。x linspace(-8, 8, 2000); % 位置网格自然单位 n 0; % 量子数从基态开始 % 埃尔米特多项式 H_n(x) Hn hermiteH(n, x); % 归一化常数 norm_const 1 / sqrt(2^n * factorial(n) * sqrt(pi)); psi_n norm_const * Hn .* exp(-x.^2 / 2); prob abs(psi_n).^2;hermiteH是符号数学工具箱里的函数也可以换成polyval手动求埃尔米特多项式的系数再计算核心好处是不用自己推导递推公式。norm_const的写法对应自然单位下 mωℏ1其它单位制下需要在 x 前乘上 √(mω/ℏ)。2.2 经典端轨迹的概率密度是什么经典谐振子的位置轨迹是x(t) A·cos(ωt φ)若给它一个随机初始相位 φ粒子长时间停留的位置分布就不再是轨迹本身而是一个概率密度。物理上著名的结论是粒子在两端点停留时间最长位置概率密度反比于速度ρ_cl(x) 1 / (π·√(A² - x²))这个密度在 x±A 处发散经典的“粒子最常待在边界附近”和量子“基态最常在中心”形成非常强的认知冲击。A 的选取不能随便给通常将量子能级 Eₙ 对应到经典振幅Eₙ (1/2)mω²A²自然单位下 A √(2n1)。这是把量子数 n 和经典轨迹联系起来的桥梁玻尔对应原理的图景就靠这个 A 落地。A sqrt(2 * n 1); % 对应能级 n 的经典振幅 x_cl linspace(-A, A, 1000); rho_cl 1 ./ (pi * sqrt(A^2 - x_cl.^2));2.3 为什么不去数值解薛定谔方程一维定态问题当然可以用有限差分法配合特征值分解去解但这样做有三个麻烦一是边界截断位置需要反复试二是对高激发态数值特征向量误差会变大三是无法直接解释波函数里振荡节数的物理来源。解析解直接给出节点数和归一化常数实现成本几乎为零。对于需要展示物理图像、做课程作业或教学 demo 的场景我一般直接解析计算。只有研究非谐势场、含时演化时才切换成 split-operator 或 Crank-Nicolson 数值方法。3. 用 Matlab 计算并绘制一维量子谐振子的波函数与概率分布3.1 最小可运行代码从基态到第 4 激发态实际绘图时玻尔对应原理和波函数正交性都需要一组量子数共同展示而不是孤立画一条曲线。下面这段代码画出 n0,1,2,3,4 五个态的波函数与概率分布并自动分两行子图呈现。x linspace(-8, 8, 2001); n_list 0:4; figure(Color, w); for idx 1:length(n_list) n n_list(idx); Hn hermiteH(n, x); psi_n Hn .* exp(-x.^2/2) / sqrt(2^n * factorial(n) * sqrt(pi)); subplot(2, 3, idx); plot(x, psi_n, LineWidth, 1.5); hold on; plot(x, abs(psi_n).^2, LineWidth, 1.2, Color, [0.85 0.33 0.1]); yline(0, --, Color, [0.5 0.5 0.5]); grid on; title(sprintf(n %d, E %.2f \\hbar\\omega, n, n 0.5)); legend(ψ_n(x), |ψ_n(x)|², Location, northeast); xlim([-6 6]); end sgtitle(一维量子谐振子波函数与概率分布自然单位);3.2 参数说明与调整逻辑上面代码里x取 [-8,8] 是为了给 n4 甚至 n10 的态留出足够空间。若只画基态可以缩窄到 [-5,5]网格点数也能降到 1000。hermiteH直接接受向量输入所以Hn是一个与x等长的数组。注意Hn .* exp(-x.^2/2)用了逐元素乘法如果漏掉点乘Matlab 会尝试矩阵乘法并报维度错误这是最频繁的报错点。legend(ψ_n(x), |ψ_n(x)|², ...)里的上标和希腊字母直接写入字符串Matlab 默认支持 TeX 语法\hbar和\omega都会被正确渲染。用yline(0,--)画零轴参考线比grid on更直观能立刻看出波函数正负振荡的节数。3.3 经典谐振子概率分布的同图绘制把经典密度和量子概率放在同一坐标系时经典密度在端点处发散直接用plot会画出极高尖峰掩盖量子包络形态。因此我通常对经典密度做采样并绘制为直方图或者直接画密度曲线但用截断方式处理端点n 10; A sqrt(2*n1); x linspace(-8, 8, 2001); psi_q hermiteH(n, x) .* exp(-x.^2/2) / sqrt(2^n * factorial(n) * sqrt(pi)); prob_q abs(psi_q).^2; x_cl linspace(-A0.02, A-0.02, 500); % 避开奇点 rho_cl 1 ./ (pi * sqrt(A^2 - x_cl.^2)); figure(Color, w); plot(x, prob_q, b, LineWidth, 1.4); hold on; plot(x_cl, rho_cl, r--, LineWidth, 1.4); legend(量子 |ψ(x)|², 经典 \rho(x), Location, north); xlabel(x (自然单位)); ylabel(概率密度); title(sprintf(一维谐振子量子与经典概率分布对比 (n%d), n)); xlim([-8 8]);经典振幅 A 采用能级匹配方式计算即把量子能量 Eₙ n1/2 换算成振幅 √(2n1)。端点处密度发散所以从 ±A±0.02 开始采样避免除以零。这个截断值不是物理量只是为了绘图美观实际数值计算中若要积分必须解析处理该奇点。4. 量子概率分布向经典分布过渡玻尔对应原理的可视化验证4.1 高量子数下量子概率分布的振荡态量子概率分布不是一座光滑的山包而是带有 n1 个峰谷的振荡函数。经典密度则是一个以中心对称的 U 形曲线两端无限高。直觉上两者差得很远但随着 n 变大量子峰的间距变小包络逐渐逼向经典曲线。要验证这一点只需要循环画不同 n 的对比图。n_list [5, 15, 30, 60]; figure(Color, w); for idx 1:4 n n_list(idx); A sqrt(2*n1); x linspace(-A-2, A2, 5001); psi_q hermiteH(n, x) .* exp(-x.^2/2) / sqrt(2^n*factorial(n)*sqrt(pi)); subplot(2, 2, idx); plot(x, abs(psi_q).^2, b); hold on; x_cl linspace(-A0.05, A-0.05, 500); rho_cl 1 ./ (pi * sqrt(A^2 - x_cl.^2)); plot(x_cl, rho_cl, r--, LineWidth, 1.2); ylim([0 0.35]); legend(量子, 经典, Location, north); title(sprintf(n%d, A%.1f, n, A)); end sgtitle(玻尔对应原理高量子数下量子概率向经典分布逼近);n60 时波函数节点数增加到 61相邻节点距离变小量子概率密度出现高频振荡。为了看清包络网格点数已经提升到 5001。若继续用 2000 点节点附近的尖峰会丢失看起来像粗糙毛刺影响对比效果。4.2 为什么高能态的量子曲线会和经典包络重合量子概率密度的包络趋近经典结果在数学上对应斯坦福大学教材里常提的“局部平均”量子概率密度在小区间上的平均值等于经典概率密度在同样区间上的积分平均值。实际绘图时可以直接对量子概率做滑窗平均再与经典曲线对比效果比裸画振荡曲线更好看。prob_sm movmean(abs(psi_q).^2, 201); % 滑动平均窗口宽度 201 点 plot(x, prob_sm, g-, LineWidth, 1.8);movmean窗口宽度取决于网格密度2001 个点分布在 [-8,8]窗口宽度 201 对应约 1.6 个自然长度单位。窗口太小平滑不掉振荡窗口太大包络被压扁。如果只是做数值实验也可以换成filter实现任意形状的滑窗卷积。4.3 经典粒子的数值采样与直方图对比换一种验证思路直接用随机相位模拟经典粒子的位置分布并把它与量子概率密度画成对比直方图。这里的直方图不是用解析概率密度画而是通过蒙特卡洛采样生成样本更贴近实验探测的物理图景。rng(42); num_samples 1e5; phases 2*pi*rand(num_samples, 1); x_samples A * cos(phases); histogram(x_samples, 200, Normalization, pdf, FaceAlpha, 0.4, EdgeColor, none); hold on; plot(x, abs(psi_q).^2, b, LineWidth, 1.5);采样数量取 1e5直方图分箱 200 个概率密度归一化方式要和概率密度曲线保持同一量纲。相位取均匀分布x 分布自然形成端点多、中间少的形状。这里的物理假设是粒子的初始相位随机且遍历足够长时间与解析密度 ρ_cl(x) 互为等价。5. 数值细节与易踩的坑归一化、网格密度与积分自检5.1 归一化自检trapz 或 integral 都必须接近 1量子波函数画的再漂亮如果概率密度积分不为 1物理上就没有意义。解析归一化常数已写在系数里但离散网格采样后数值积分会因为截断而损失部分概率特别是 n 较大时波函数尾部延伸变长。每次绘图前做这个检查norm_check trapz(x, abs(psi_q).^2); fprintf(归一化积分值 %.6f\n, norm_check);如果结果不在 0.999 到 1.001 之间先检查 x 网格范围是否覆盖到波函数尾部再检查网格点数是否足够。n30 时波函数尾部能达到 x≈±10如果网格只取到 ±6数值积分丢掉的概率会超过百分之一图上线段看起来没事积分却会警告你。5.2 高量子数下的振荡分辨原则量子数 n 越大波函数节点越多每个振荡半周期越来越窄。奈奎斯特采样定律在空间域同样生效一个振荡周期至少撞到两个网格点否则波峰会失真。经验法则如下表。量子数 n推荐网格范围推荐最小网格点数备注0–2[-5, 5]801基态波函数平滑粗网格即可3–10[-6, 6]1601中等振荡需保证两端收敛11–30[-8, 8]2001–4001高激发态尾部延伸到 ±8 以上30–80按 A3 动态扩展5000建议用 x linspace(-A-3, A3, N)linspace第三个参数不能一次给太大n1 用 5000 点纯属浪费计算时间。动态生成网格 A3 是更合理的工程习惯避免每改一个 n 就去手动修改 x 区间。5.3 经典密度在端点发散的处理原则经典概率密度 ρ_cl(x) 在 x±A 处有可积奇点数值上直接计算会产生 1/0 警告和画图尖峰。处理方法有两种一种是把 x_cl 起始点改成 ±A±ε另一种是绘制直方图代替连续曲线。若需要精确积分经典概率密度则用变量替换 xA·cosθ把积分变为对 θ 的均匀积分这是最稳妥的路径。不要试图用很大的数值替换无限值那只会让 Y 轴刻度失衡。5.4 常见报错与误用对照hermiteH要求输入符号或数值输入若 n 是浮点数 1.5会报“期望整数”类错误。另一个高频问题是忘记在乘法处加点号Hn * exp(-x.^2/2)会触发矩阵乘法错误报错信息是“Inner matrix dimensions must agree”解决办法是改成.*。还有一类场景是用户把prob psi_n.^2当作概率密度忘了取模平方abs(psi_n).^2导致负区域也画成正值看起来波函数和概率分布完全重合实际上方向错了。6. 进阶技巧数值积分验证概率守恒并观察波包振荡六个章节到这里静态波函数已经能画得很熟练。最后一个值得动手的技巧是用含时演化把基态波函数偏移后释放观察概率密度的周期性振荡并从时间演化中提取振荡频率和经典 ω 做对比。### 6.1 快速含时演化实现Crank-Nicolson 一维版本 dx 0.05; x -20:dx:20; N length(x); dt 0.01; V 0.5 * x.^2; % 谐振子势 % 偏移高斯波包 psi0 exp(-(x - 5).^2 / 2) / pi^(1/4); psi0 psi0 / sqrt(trapz(x, abs(psi0).^2)); % 构造三对角矩阵 T -diag(ones(N-1,1), -1) 2*eye(N) - diag(ones(N-1,1), 1); T T / (2*dx^2); % 泡利矩阵形式的半步演化算子 A eye(N) 0.5i*dt*T - 0.5i*dt*diag(V); B eye(N) - 0.5i*dt*T 0.5i*dt*diag(V); psi psi0; num_steps 4000; for step 1:num_steps psi A \ (B * psi); if mod(step, 500) 0 x_mean(step/500) trapz(x, x .* abs(psi).^2); time(step/500) step * dt; end end plot(time, x_mean, o-);这段代码的物理逻辑是先把谐振子基态波函数沿 x 正方向偏移 5 个单位再释放它不再是能量本征态因此会围绕原点做量子相干振荡。Crank-Nicolson 格式天生稳定时间步长 dt 取 0.01 时不会发散但 A 矩阵求逆每次迭代都要做一次4000 步耗时约十几秒属于可接受范围。x_mean的期望位置随时间近似按 cos 振荡其角频率应接近经典值 ω1用fft(x_mean)找主频就能验证数值精度。这套思路进一步可以推广到双阱势场的隧穿问题只需修改V矩阵即可。本文还有配套的精品资源点击获取
一维谐振子波函数与概率分布：量子经典对比的Matlab可视化
简介基于Matlab的一维量子和经典谐振子仿真资源包主要面向物理、电子信息工程、数学等专业的高校学生及相关科研人员适用于课程设计、期末大作业与毕业设计也适合想借助数值计算理解量子力学模型的初学者。资源共包含8个文件其中有5个脚本文件、2个说明文档和1张结果示意图压缩包仅约87KB轻量易用已有306人学习下载。程序采用参数化编程关键参数可灵活修改注释明细、结构清晰并附可直接运行的案例数据利于快速复现和二次开发。脚本涵盖厄米多项式系数计算、经典概率密度函数、量子谐振子波函数与概率分布绘制等模块可直观对比一维量子和经典谐振子的行为差异帮助读者将Matlab数值计算与物理理论结合起来。调整参数即可观察不同状态下的波形变化是教学演示和实践演练的实用工具。1. 一维谐振子用量子和经典两条路线同时读懂波函数一个粒子被束缚在势场里经典力学会告诉你它按正弦规律来回摆动量子力学则只肯给出波函数和概率分布告诉你它在某个位置出现的概率。同样是谐振子经典给出的是确定轨迹量子给出的是离散能级和波函数干涉条纹。把两者算出来、画在同一张图上立刻就能看到著名的玻尔对应原理——高量子数下量子概率分布逐步收敛为经典概率分布。本文就以“一维量子和经典谐振子的波函数和概率分布”为对象直接用 Matlab 做解析计算与绘制不改薛定谔方程求解器不做蒙特卡洛纯靠埃尔米特多项式和经典相空间密度就能在几分钟内看到双峰、零点和概率密度包络。适合正在学量子力学或数值方法、需要快速可复现波形图的本科生、研究生和科研入门人员。2. 从物理模型到 Matlab 计算方案为什么解析式比数值打靶更快2.1 量子端固有频率、自然单位与定态解一维量子谐振子的哈密顿量写作H -ℏ²/(2m) · d²/dx² (1/2)mω²x²标准的定态解是厄米多项式与高斯函数的乘积ψₙ(x) (1/√(2ⁿn!√π)) · Hₙ(√(mω/ℏ)x) · exp(-mωx²/(2ℏ))能量本征值为 Eₙ ℏω(n1/2)。很多教材把 ℏ m ω 1 定为自然单位这时波函数形式最简洁Matlab 代码也只跟无量纲的 x 打交道。实际绘图时我推荐直接采用这种约定坐标轴标成“自然单位”显得专业且能绕开繁杂的量纲转换。x linspace(-8, 8, 2000); % 位置网格自然单位 n 0; % 量子数从基态开始 % 埃尔米特多项式 H_n(x) Hn hermiteH(n, x); % 归一化常数 norm_const 1 / sqrt(2^n * factorial(n) * sqrt(pi)); psi_n norm_const * Hn .* exp(-x.^2 / 2); prob abs(psi_n).^2;hermiteH是符号数学工具箱里的函数也可以换成polyval手动求埃尔米特多项式的系数再计算核心好处是不用自己推导递推公式。norm_const的写法对应自然单位下 mωℏ1其它单位制下需要在 x 前乘上 √(mω/ℏ)。2.2 经典端轨迹的概率密度是什么经典谐振子的位置轨迹是x(t) A·cos(ωt φ)若给它一个随机初始相位 φ粒子长时间停留的位置分布就不再是轨迹本身而是一个概率密度。物理上著名的结论是粒子在两端点停留时间最长位置概率密度反比于速度ρ_cl(x) 1 / (π·√(A² - x²))这个密度在 x±A 处发散经典的“粒子最常待在边界附近”和量子“基态最常在中心”形成非常强的认知冲击。A 的选取不能随便给通常将量子能级 Eₙ 对应到经典振幅Eₙ (1/2)mω²A²自然单位下 A √(2n1)。这是把量子数 n 和经典轨迹联系起来的桥梁玻尔对应原理的图景就靠这个 A 落地。A sqrt(2 * n 1); % 对应能级 n 的经典振幅 x_cl linspace(-A, A, 1000); rho_cl 1 ./ (pi * sqrt(A^2 - x_cl.^2));2.3 为什么不去数值解薛定谔方程一维定态问题当然可以用有限差分法配合特征值分解去解但这样做有三个麻烦一是边界截断位置需要反复试二是对高激发态数值特征向量误差会变大三是无法直接解释波函数里振荡节数的物理来源。解析解直接给出节点数和归一化常数实现成本几乎为零。对于需要展示物理图像、做课程作业或教学 demo 的场景我一般直接解析计算。只有研究非谐势场、含时演化时才切换成 split-operator 或 Crank-Nicolson 数值方法。3. 用 Matlab 计算并绘制一维量子谐振子的波函数与概率分布3.1 最小可运行代码从基态到第 4 激发态实际绘图时玻尔对应原理和波函数正交性都需要一组量子数共同展示而不是孤立画一条曲线。下面这段代码画出 n0,1,2,3,4 五个态的波函数与概率分布并自动分两行子图呈现。x linspace(-8, 8, 2001); n_list 0:4; figure(Color, w); for idx 1:length(n_list) n n_list(idx); Hn hermiteH(n, x); psi_n Hn .* exp(-x.^2/2) / sqrt(2^n * factorial(n) * sqrt(pi)); subplot(2, 3, idx); plot(x, psi_n, LineWidth, 1.5); hold on; plot(x, abs(psi_n).^2, LineWidth, 1.2, Color, [0.85 0.33 0.1]); yline(0, --, Color, [0.5 0.5 0.5]); grid on; title(sprintf(n %d, E %.2f \\hbar\\omega, n, n 0.5)); legend(ψ_n(x), |ψ_n(x)|², Location, northeast); xlim([-6 6]); end sgtitle(一维量子谐振子波函数与概率分布自然单位);3.2 参数说明与调整逻辑上面代码里x取 [-8,8] 是为了给 n4 甚至 n10 的态留出足够空间。若只画基态可以缩窄到 [-5,5]网格点数也能降到 1000。hermiteH直接接受向量输入所以Hn是一个与x等长的数组。注意Hn .* exp(-x.^2/2)用了逐元素乘法如果漏掉点乘Matlab 会尝试矩阵乘法并报维度错误这是最频繁的报错点。legend(ψ_n(x), |ψ_n(x)|², ...)里的上标和希腊字母直接写入字符串Matlab 默认支持 TeX 语法\hbar和\omega都会被正确渲染。用yline(0,--)画零轴参考线比grid on更直观能立刻看出波函数正负振荡的节数。3.3 经典谐振子概率分布的同图绘制把经典密度和量子概率放在同一坐标系时经典密度在端点处发散直接用plot会画出极高尖峰掩盖量子包络形态。因此我通常对经典密度做采样并绘制为直方图或者直接画密度曲线但用截断方式处理端点n 10; A sqrt(2*n1); x linspace(-8, 8, 2001); psi_q hermiteH(n, x) .* exp(-x.^2/2) / sqrt(2^n * factorial(n) * sqrt(pi)); prob_q abs(psi_q).^2; x_cl linspace(-A0.02, A-0.02, 500); % 避开奇点 rho_cl 1 ./ (pi * sqrt(A^2 - x_cl.^2)); figure(Color, w); plot(x, prob_q, b, LineWidth, 1.4); hold on; plot(x_cl, rho_cl, r--, LineWidth, 1.4); legend(量子 |ψ(x)|², 经典 \rho(x), Location, north); xlabel(x (自然单位)); ylabel(概率密度); title(sprintf(一维谐振子量子与经典概率分布对比 (n%d), n)); xlim([-8 8]);经典振幅 A 采用能级匹配方式计算即把量子能量 Eₙ n1/2 换算成振幅 √(2n1)。端点处密度发散所以从 ±A±0.02 开始采样避免除以零。这个截断值不是物理量只是为了绘图美观实际数值计算中若要积分必须解析处理该奇点。4. 量子概率分布向经典分布过渡玻尔对应原理的可视化验证4.1 高量子数下量子概率分布的振荡态量子概率分布不是一座光滑的山包而是带有 n1 个峰谷的振荡函数。经典密度则是一个以中心对称的 U 形曲线两端无限高。直觉上两者差得很远但随着 n 变大量子峰的间距变小包络逐渐逼向经典曲线。要验证这一点只需要循环画不同 n 的对比图。n_list [5, 15, 30, 60]; figure(Color, w); for idx 1:4 n n_list(idx); A sqrt(2*n1); x linspace(-A-2, A2, 5001); psi_q hermiteH(n, x) .* exp(-x.^2/2) / sqrt(2^n*factorial(n)*sqrt(pi)); subplot(2, 2, idx); plot(x, abs(psi_q).^2, b); hold on; x_cl linspace(-A0.05, A-0.05, 500); rho_cl 1 ./ (pi * sqrt(A^2 - x_cl.^2)); plot(x_cl, rho_cl, r--, LineWidth, 1.2); ylim([0 0.35]); legend(量子, 经典, Location, north); title(sprintf(n%d, A%.1f, n, A)); end sgtitle(玻尔对应原理高量子数下量子概率向经典分布逼近);n60 时波函数节点数增加到 61相邻节点距离变小量子概率密度出现高频振荡。为了看清包络网格点数已经提升到 5001。若继续用 2000 点节点附近的尖峰会丢失看起来像粗糙毛刺影响对比效果。4.2 为什么高能态的量子曲线会和经典包络重合量子概率密度的包络趋近经典结果在数学上对应斯坦福大学教材里常提的“局部平均”量子概率密度在小区间上的平均值等于经典概率密度在同样区间上的积分平均值。实际绘图时可以直接对量子概率做滑窗平均再与经典曲线对比效果比裸画振荡曲线更好看。prob_sm movmean(abs(psi_q).^2, 201); % 滑动平均窗口宽度 201 点 plot(x, prob_sm, g-, LineWidth, 1.8);movmean窗口宽度取决于网格密度2001 个点分布在 [-8,8]窗口宽度 201 对应约 1.6 个自然长度单位。窗口太小平滑不掉振荡窗口太大包络被压扁。如果只是做数值实验也可以换成filter实现任意形状的滑窗卷积。4.3 经典粒子的数值采样与直方图对比换一种验证思路直接用随机相位模拟经典粒子的位置分布并把它与量子概率密度画成对比直方图。这里的直方图不是用解析概率密度画而是通过蒙特卡洛采样生成样本更贴近实验探测的物理图景。rng(42); num_samples 1e5; phases 2*pi*rand(num_samples, 1); x_samples A * cos(phases); histogram(x_samples, 200, Normalization, pdf, FaceAlpha, 0.4, EdgeColor, none); hold on; plot(x, abs(psi_q).^2, b, LineWidth, 1.5);采样数量取 1e5直方图分箱 200 个概率密度归一化方式要和概率密度曲线保持同一量纲。相位取均匀分布x 分布自然形成端点多、中间少的形状。这里的物理假设是粒子的初始相位随机且遍历足够长时间与解析密度 ρ_cl(x) 互为等价。5. 数值细节与易踩的坑归一化、网格密度与积分自检5.1 归一化自检trapz 或 integral 都必须接近 1量子波函数画的再漂亮如果概率密度积分不为 1物理上就没有意义。解析归一化常数已写在系数里但离散网格采样后数值积分会因为截断而损失部分概率特别是 n 较大时波函数尾部延伸变长。每次绘图前做这个检查norm_check trapz(x, abs(psi_q).^2); fprintf(归一化积分值 %.6f\n, norm_check);如果结果不在 0.999 到 1.001 之间先检查 x 网格范围是否覆盖到波函数尾部再检查网格点数是否足够。n30 时波函数尾部能达到 x≈±10如果网格只取到 ±6数值积分丢掉的概率会超过百分之一图上线段看起来没事积分却会警告你。5.2 高量子数下的振荡分辨原则量子数 n 越大波函数节点越多每个振荡半周期越来越窄。奈奎斯特采样定律在空间域同样生效一个振荡周期至少撞到两个网格点否则波峰会失真。经验法则如下表。量子数 n推荐网格范围推荐最小网格点数备注0–2[-5, 5]801基态波函数平滑粗网格即可3–10[-6, 6]1601中等振荡需保证两端收敛11–30[-8, 8]2001–4001高激发态尾部延伸到 ±8 以上30–80按 A3 动态扩展5000建议用 x linspace(-A-3, A3, N)linspace第三个参数不能一次给太大n1 用 5000 点纯属浪费计算时间。动态生成网格 A3 是更合理的工程习惯避免每改一个 n 就去手动修改 x 区间。5.3 经典密度在端点发散的处理原则经典概率密度 ρ_cl(x) 在 x±A 处有可积奇点数值上直接计算会产生 1/0 警告和画图尖峰。处理方法有两种一种是把 x_cl 起始点改成 ±A±ε另一种是绘制直方图代替连续曲线。若需要精确积分经典概率密度则用变量替换 xA·cosθ把积分变为对 θ 的均匀积分这是最稳妥的路径。不要试图用很大的数值替换无限值那只会让 Y 轴刻度失衡。5.4 常见报错与误用对照hermiteH要求输入符号或数值输入若 n 是浮点数 1.5会报“期望整数”类错误。另一个高频问题是忘记在乘法处加点号Hn * exp(-x.^2/2)会触发矩阵乘法错误报错信息是“Inner matrix dimensions must agree”解决办法是改成.*。还有一类场景是用户把prob psi_n.^2当作概率密度忘了取模平方abs(psi_n).^2导致负区域也画成正值看起来波函数和概率分布完全重合实际上方向错了。6. 进阶技巧数值积分验证概率守恒并观察波包振荡六个章节到这里静态波函数已经能画得很熟练。最后一个值得动手的技巧是用含时演化把基态波函数偏移后释放观察概率密度的周期性振荡并从时间演化中提取振荡频率和经典 ω 做对比。### 6.1 快速含时演化实现Crank-Nicolson 一维版本 dx 0.05; x -20:dx:20; N length(x); dt 0.01; V 0.5 * x.^2; % 谐振子势 % 偏移高斯波包 psi0 exp(-(x - 5).^2 / 2) / pi^(1/4); psi0 psi0 / sqrt(trapz(x, abs(psi0).^2)); % 构造三对角矩阵 T -diag(ones(N-1,1), -1) 2*eye(N) - diag(ones(N-1,1), 1); T T / (2*dx^2); % 泡利矩阵形式的半步演化算子 A eye(N) 0.5i*dt*T - 0.5i*dt*diag(V); B eye(N) - 0.5i*dt*T 0.5i*dt*diag(V); psi psi0; num_steps 4000; for step 1:num_steps psi A \ (B * psi); if mod(step, 500) 0 x_mean(step/500) trapz(x, x .* abs(psi).^2); time(step/500) step * dt; end end plot(time, x_mean, o-);这段代码的物理逻辑是先把谐振子基态波函数沿 x 正方向偏移 5 个单位再释放它不再是能量本征态因此会围绕原点做量子相干振荡。Crank-Nicolson 格式天生稳定时间步长 dt 取 0.01 时不会发散但 A 矩阵求逆每次迭代都要做一次4000 步耗时约十几秒属于可接受范围。x_mean的期望位置随时间近似按 cos 振荡其角频率应接近经典值 ω1用fft(x_mean)找主频就能验证数值精度。这套思路进一步可以推广到双阱势场的隧穿问题只需修改V矩阵即可。本文还有配套的精品资源点击获取 
