Kansas方法在非线性薛定谔方程（NLSE）的数值求解
“Kansas方法”在非线性薛定谔方程（NLSE）的数值求解中通常指一种分裂步长法（Split-Step Method），也称为分步傅里叶法（Split-Step Fourier Method, SSFM）。该方法将方程中的线性部分（拉普拉斯项）和非线性部分分开处理，利用傅里叶变换高效计算线性部分的演化，而对非线性部分采用局部近似（线性化）。时间方向上的推进可以采用有限差分（如Crank-Nicolson）或直接指数积分。
一、二维非线性薛定谔方程（NLSE）
其中：
- (\(\psi(x,y,t)\))：波函数（复值）
- (\(\nabla^2 = \partial^2/\partial x^2 + \partial^2/\partial y^2\))
- (\(V(\mathbf{r})\))：外势（如谐振子势）
- (\(g\))：非线性系数（(\(g>0\)) 聚焦，(\(g<0\)) 散焦）
二、Kansas方法（分裂步长法）原理
将方程写为：
其中线性算子 (\(\hat{L} = -\frac{1}{2}\nabla^2 + V\))，非线性算子 (\(\hat{N} = g|\psi|^2\))。
分裂步长的思想是将时间步长 (\(\Delta t\)) 分成若干子步，交替处理线性和非线性部分。最常用的是Strang分裂（二阶精度）：
- 半步非线性演化：(\(\psi^* = e^{-i\Delta t \hat{N}/2}\,\psi^n\))
- 整步线性演化：(\(\psi^{**} = e^{-i\Delta t \hat{L}}\,\psi^*\))
- 半步非线性演化：(\(\psi^{n+1} = e^{-i\Delta t \hat{N}/2}\,\psi^{**}\))
由于 (\(\hat{L}\)) 在傅里叶空间中是对角的，线性部分可以精确计算：
但实际上 (V) 是空间的函数，不能直接与拉普拉斯项同时对角化。因此通常将 (\hat{L}) 进一步分裂为动能项 (\hat{T} = -\frac{1}{2}\nabla^2) 和势能项 (\hat{V})，但这会引入额外的分裂误差。更常见的做法是采用对称分裂：
- 半步势能：(\(\psi^* = e^{-i\Delta t V/2}\,\psi^n\))
- 整步动能：(\(\psi^{**} = \mathcal{F}^{-1}\left[e^{-i\Delta t (k^2/2)}\mathcal{F}[\psi^*]\right]\))
- 半步势能：(\(\psi^{***} = e^{-i\Delta t V/2}\,\psi^{**}\))
- 非线性部分（线性化处理）...
注意：非线性项 (\(g|\psi|^2\psi\)) 依赖于 (\(\psi\)) 本身，在分裂法中通常将其视为一个“势” (\(U_{\text{nl}}=g|\psi|^2\))，并在每个子步中认为它保持不变（即线性化）。这就是您所说的对非线性项进行线性处理。
三、非线性项的线性化处理
在时间步 (\([t_n, t_{n+1}]\)) 内，我们近似 (\(|\psi(t)|^2 \approx |\psi(t_n)|^2\))，从而非线性项变为一个已知的“势”：
\(N(\psi) \approx g\,|\psi^n|^2\,\psi\)
这样在子步中求解线性方程：
\(i\frac{\partial \psi}{\partial t} = g|\psi^n|^2\psi \quad\Rightarrow\quad \psi(t_{n+1}) = e^{-i\Delta t\,g|\psi^n|^2}\psi(t_n)\)
这种线性化是显式的，稳定性受限于非线性强度和时间步长。为了提高稳定性，可以采用半隐式或迭代线性化（如固定点迭代），但会增加计算量。
四、MATLAB 完整实现
MATLAB 代码，求解二维 NLSE，采用 Strang 分裂 + 线性化非线性项。
%% Kansas方法（分裂步长法）求解二维非线性薛定谔方程
clear; clc; close all;
%% 参数设置
L = 10; % 空间半宽度
Nx = 128; % 网格点数（2的幂有利于FFT）
Ny = 128;
dx = 2*L/Nx;
dy = 2*L/Ny;
x = linspace(-L, L-dx, Nx);
y = linspace(-L, L-dy, Ny);
[X, Y] = meshgrid(x, y);
dt = 0.001; % 时间步长
T = 1; % 总时间
Nt = round(T/dt); % 时间步数
% 物理参数
V = 0.5*(X.^2 + Y.^2); % 谐振子势
g = 1.0; % 非线性系数（聚焦）
% 初始波包
sigma = 0.5;
x0 = -2; y0 = -2;
kx0 = 2; ky0 = 2;
psi = exp(-((X-x0).^2 + (Y-y0).^2)/(2*sigma^2)) .* exp(1i*(kx0*X + ky0*Y));
psi = psi / sqrt(sum(abs(psi(:)).^2)*dx*dy); % 归一化
%% 傅里叶波数
kx = 2*pi/L * [0:Nx/2-1, -Nx/2:-1];
ky = 2*pi/L * [0:Ny/2-1, -Ny/2:-1];
[KX, KY] = meshgrid(kx, ky);
K2 = KX.^2 + KY.^2; % 拉普拉斯算子对应波数平方
%% 分裂步长法主循环
psi_save = cell(1, 5); % 保存几个时刻的结果
save_times = round(linspace(1, Nt, 5));
cnt = 1;
for n = 1:Nt
% --- 第一步：半步势能 ---
psi = psi .* exp(-1i*dt/2 * V);
% --- 第二步：整步动能（傅里叶空间） ---
psi_hat = fft2(psi);
psi_hat = psi_hat .* exp(-1i*dt/2 * K2); % 注意：这里用了 dt/2？实际上Strang分裂中动能步应为整步
% 修正：Strang分裂中动能步应该是整步，但为了与势能配合，通常写作：
% psi_hat = psi_hat .* exp(-1i*dt * K2/2); 但这里我们已经写了dt/2，需要调整
% 更标准的写法见下方注释
psi = ifft2(psi_hat);
% --- 第三步：半步势能 ---
psi = psi .* exp(-1i*dt/2 * V);
% --- 第四步：非线性项（线性化） ---
% 注意：非线性项应该在势能之后处理，但顺序可调，这里放在最后
psi = psi .* exp(-1i*dt * g * abs(psi).^2);
% 保存结果
if any(n == save_times)
psi_save{cnt} = psi;
cnt = cnt + 1;
end
end
%% 绘图
figure('Position', [100 100 1200 450]);
titles = {'t=0', 't=0.25T', 't=0.5T', 't=0.75T', 't=T'};
for k = 1:5
subplot(1,5,k);
imagesc(x, y, abs(psi_save{k}).^2);
axis xy equal tight; colorbar;
title(titles{k}); xlabel('x'); ylabel('y');
end
sgtitle('概率密度演化（Kansas分裂步长法）');
重要修正：上述代码中的分裂顺序并非严格的 Strang 分裂，因为非线性项被单独处理。更准确的 Strang 分裂（二阶）应该如下安排：
- 半步非线性：(\(\psi^* = e^{-i\Delta t\,g|\psi|^2/2}\,\psi^n\))
- 整步线性（势能+动能）：(\(\psi^{**} = e^{-i\Delta t(\hat{T}+V)}\psi^*\))
- 半步非线性：(\(\psi^{n+1} = e^{-i\Delta t\,g|\psi^{**}|^2/2}\,\psi^{**}\))
但由于线性部分 (\(\hat{T}+V\)) 不易联合对角化，通常进一步分裂为：
- 半步势能：(\(\psi^{(1)} = e^{-i\Delta t V/2}\,\psi^n\))
- 整步动能：(\(\psi^{(2)} = \mathcal{F}^{-1}[e^{-i\Delta t k^2/2}\mathcal{F}[\psi^{(1)}]]\))
- 半步势能：(\(\psi^{(3)} = e^{-i\Delta t V/2}\,\psi^{(2)}\))
- 整步非线性：(\(\psi^{n+1} = e^{-i\Delta t\,g|\psi^{(3)}|^2}\,\psi^{(3)}\)) （或分两步）
这样的分裂仍然是二阶精度的。下面是修正后的标准代码：
%% 修正版：标准Strang分裂（二阶）
for n = 1:Nt
% 1. 半步势能
psi = psi .* exp(-1i*dt/2 * V);
% 2. 整步动能（傅里叶空间）
psi_hat = fft2(psi);
psi_hat = psi_hat .* exp(-1i*dt * K2/2); % 注意：这里 dt 整步，K2/2 对应 -1/2 ∇²
psi = ifft2(psi_hat);
% 3. 半步势能
psi = psi .* exp(-1i*dt/2 * V);
% 4. 整步非线性（线性化）
psi = psi .* exp(-1i*dt * g * abs(psi).^2);
% 保存...
end
五、结果与讨论
运行上述代码，您将看到高斯波包在谐振子势中演化的概率密度图像。由于非线性项的存在，波包可能会发生聚焦或散焦，并与势场相互作用。
线性化处理的有效性：当时间步长 (\(\Delta t\)) 足够小时，(\(|\psi|^2\)) 的变化不大，线性化近似是合理的。但如果 (\(g\)) 很大或 (\(\Delta t\)) 较大，可能会出现数值不稳定。此时可采用隐式-显式（IMEX） 方法或迭代线性化（如 Picard 迭代）。
参考代码 利用Kansas方法求解二维薛定谔方程，时间方向利用有限差分方法，其中对非线性项进行了线性处理 www.youwenfan.com/contentcnv/81149.html
六、与有限差分的关系
您提到“时间方向利用有限差分方法”，但在分裂步长法中，时间演化是通过指数算子直接进行的，并没有使用有限差分。如果您希望在时间方向上使用有限差分（如 Crank-Nicolson），则需要将整个方程视为一个整体，并对空间导数进行离散，得到一个大型线性系统。此时非线性项的处理通常采用固定点迭代或牛顿法。例如：
其中 (\(\psi^{n+1/2} = (\psi^n + \psi^{n+1})/2\))，这样得到关于 (\(\psi^{n+1}\)) 的线性系统（因为非线性项被冻结在 (\(n\)) 时刻）。这种方法称为半隐式 Crank-Nicolson，也是一种常见的线性化处理。
浙公网安备 33010602011771号