MATLAB中圆柱坐标系下2D各向同性谐振子哈密顿矩阵构建问题

我正在尝试用MATLAB求解圆柱坐标系下二维各向同性谐振子的哈密顿矩阵，但遇到了问题。需要将角动量为$l$的项0.5*((l^2)/(rho^2))添加到哈密顿矩阵的对角线上，寻求帮助。相关公式如下：

现有MATLAB代码（已翻译） #

% 计算权重因子mu
mu= sqrt(1/2)*sqrt(rho1.^2 - rho2.^2);
mu_i_plus_1= mu(2:end);
mu_i_min_1= [mu(1)];
mu_i_min_1= [mu_i_min_1, mu(1:end)];
% 步长 h = rho(2) - rho(1);
% 计算哈密顿量的对角元和非对角元
H_diag = (-rho1-rho2)./(h.*mu.^2) - (2*l ./ rho(1:end-1).^2) - rho(1:end-1).^2;
H_off = rho1(1:end-1)./(mu(1:end-1).*mu_i_plus_1.*h);
H_0 = -0.5*(diag(H_diag)+ diag(H_off(1: length(H_diag)-1),1) + diag(H_off(1: length(H_diag)-1),-1)); 

错误修正方案 #

当前代码的核心错误在于角动量项的实现：

1. 原公式是正的势能项0.5*(l²/ρ²)，但代码里写成了负的-(2*l ./ rho(1:end-1).^2)，系数和符号完全错误。
2. 角量子数应该取平方l^2，而非2*l。

修正后的H_diag计算行：

H_diag = (-rho1-rho2)./(h.*mu.^2) - rho(1:end-1).^2 + 0.5*(l^2 ./ rho(1:end-1).^2);

修改说明：

• 移除错误的-(2*l ./ ...)项，替换为正确的角动量势能项+0.5*(l^2 ./ ...)，确保符号和形式匹配公式要求。
• 角动量项是势能的一部分，需要加到哈密顿量的对角元（包含动能和谐振子势能）中。

内容的提问来源于stack exchange，提问作者helene de witte