带三次项的非线性四阶薛定谔方程的一个局部能量守恒格式Matlab代码

带三次项的非线性四阶薛定谔方程可以写成如下形式：

$$
i\frac{\partial \psi}{\partial t} = -\frac{1}{2}\nabla^2\psi + V(x,y)\psi + g|\psi|^2\psi + f|\psi|^4\psi
$$

其中 $g$ 和 $f$ 是非线性系数。下面给出一个局部能量守恒格式的 Matlab 代码：

% 定义参数
N = 64;              % 网格数
L = 20;              % 区域长度
h = L/N;             % 网格大小
dt = 0.01;           % 时间步长
T = 10;              % 模拟时间
g = 1;               % 非线性系数 g
f = 1;               % 非线性系数 f

% 初始化波函数和势能
x = linspace(-L/2, L/2, N+1); x(end) = [];
y = x;
[X, Y] = meshgrid(x, y);
V = 0.5*(X.^2 + Y.^2);       % 势能
psi = exp(-X.^2-Y.^2);       % 波函数

% 设定差分算子和局部能量守恒格式中的参数
hbar = 1;
m = 1;
Lap = -1/12*(circshift(psi, [0 -2]) + circshift(psi, [0 2]) ...
            + circshift(psi, [-2 0]) + circshift(psi, [2 0]) ...
            - 8*psi)/h^2;
C = hbar^2/(2*m);
K = exp(-1i*dt*C*Lap);
V = exp(-1i*dt/hbar*V);
psi2 = abs(psi).^2;
E = sum(sum(C*Lap.*conj(psi).*psi + V.*psi2.*(g*abs(psi).^2 + f*psi2)));

% 使用局部能量守恒格式进行时间演化
Nt = round(T/dt);
for n = 1:Nt
    psi = K.*psi;
    psi2 = abs(psi).^2;
    K = exp(-1i*dt*(C*Lap + V.*(g*psi2 + f*psi2.^2 - E)));
end

该局部能量守恒格式中，$Lap$ 代表二阶差分算子，$C$ 为普朗克常数除以 $2m$，$K$ 是时间演化算子，$E$ 是局部能量。在每个时间步长内，首先使用时间演化算子 $K$ 将波函数演化一半时间，然后计算波函数的模平方 $|\psi|^2$，局部能量 $E$ 和新的时间演化算子 $K$，最后使用新的时间演化算子 $K$ 将波函数演化一半时间。