查找错误
function [psinew, phinew] = new2D1(h, dt)
% 二维非线性薛定谔方程求解器
L = 10; a = 1; % 空间域半长,波速
T = 1; % 终止时间
M = round(2*L/h); % 每维空间网格数
x = (-L:h:L)';
x = x(1:end-1); % 周期边界处理
[x, y] = meshgrid(x, x); % 二维网格
% 非线性参数
a1 = 0.5; b1 = -2; c1 = -1;
a2 = 0.5; b2 = -3; c2 = -1;
% 2-stage Gauss系数
c1g = 0.5 - sqrt(3)/6;
c2g = 0.5 + sqrt(3)/6;
a11 = 1/4; a12 = 1/4 - sqrt(3)/6;
a21 = 1/4 + sqrt(3)/6; a22 = 1/4;
b1g = 0.5; b2g = 0.5;
%% 计算二维拉普拉斯算子的特征值
kk = 0:M-1;
theta = 2*pi*kk/M;
lambda_A = (150 + 32*cos(theta) - 2*cos(2*theta)) / 180;
coef_S = 1/(60*h^2);
lambda_S = (-126 + 128*cos(theta) - 2*cos(2*theta)) * coef_S;
lambda_L = lambda_S ./ lambda_A;
% 二维特征值矩阵
[KX, KY] = meshgrid(lambda_L, lambda_L);
lambda_B = KX + KY; % 二维拉普拉斯算子特征值
%% 预计算FFT求解器
invD = zeros(4, 4, M, M); % 每个频率点的4x4矩阵
for i = 1:M
for j = 1:M
lam = lambda_B(i, j);
D_k = [1, 0, -a^2*dt*a11*lam, -a^2*dt*a12*lam;
0, 1, -a^2*dt*a21*lam, -a^2*dt*a22*lam;
-dt*a11, -dt*a12, 1, 0;
-dt*a21, -dt*a22, 0, 1];
invD(:, :, i, j) = inv(D_k); % 预计算逆矩阵
end
end
%% 初始条件 - 二维高斯波包
r = sqrt(x.^2 + y.^2);
exact_psi = @(x, y, t) 4*atan(exp(3 - 5*r));
exact_phi = @(x, y, t) atan(exp(3 - 5*r));
psi = exact_psi(x, y, 0);
psi0 = psi;
phi = exact_phi(x, y, 0);
phi0 = phi;
u = zeros(size(psi));
v = zeros(size(phi));
%% 初始质量 M0
M0 = h^2 * sum(sum(abs(psi0).^2 + abs(phi0).^2));
% 初始能量 E0
term1 = h^2 * sum(sum(u.^2));
term2 = -a^2 * h^2 * sum(sum(real(ifft2(lambda_B .* fft2(psi0))) .* psi0));
term3 = a1 * h^2 * sum(sum(psi0.^2));
term4 = h^2 * sum(sum(v.^2));
term5 = -a^2 * h^2 * sum(sum(real(ifft2(lambda_B .* fft2(phi0))) .* phi0));
term6 = a2 * h^2 * sum(sum(phi0.^2));
term7 = (c2*b1/(4*c1)) * h^2 * sum(sum(psi0.^4));
term8 = (b2/4) * h^2 * sum(sum(phi0.^4));
term9 = (c2/2) * h^2 * sum(sum((psi0.*phi0).^2));
E0 = (1/2) * ((c2/c1)*(term1 + term2 + term3) + term4 + term5 + term6) + term7 + term8 + term9;
Nt = round(T/dt);
% 初始化守恒量记录
M_arr = zeros(Nt+1, 1);
E_arr = zeros(Nt+1, 1);
abs_err_M = zeros(Nt+1, 1);
abs_err_E = zeros(Nt+1, 1);
M_arr(1) = M0;
E_arr(1) = E0;
abs_err_M(1) = 0;
abs_err_E(1) = 0;
%% 初始化
prev_u = u;
prev_v = v;
prev_psi = psi;
prev_phi = phi;
prev_psi1 = psi; prev_psi2 = psi;
prev_phi1 = phi; prev_phi2 = phi;
prev_u1 = u; prev_u2 = u;
prev_v1 = v; prev_v2 = v;
prev_t = 0;
%% 时间层循环
for n = 1:Nt
%% 使用预测龙格库塔求高阶预估值
current_t = n*dt;
u = prev_u; u1 = prev_u1; u2 = prev_u2;
v = prev_v; v1 = prev_v1; v2 = prev_v2;
psi = prev_psi; psi1 = prev_psi1; psi2 = prev_psi2;
phi = prev_phi; phi1 = prev_phi1; phi2 = prev_phi2;
% 使用前一时间步的值进行三点插值
t_nodes = [prev_t, prev_t + c1g*dt, prev_t + c2g*dt];
U_nodes = [u, u1, u2];
V_nodes = [v, v1, v2];
psi_nodes = [psi, psi1, psi2];
phi_nodes = [phi, phi1, phi2];
% 当前阶段点时间
t_stage1 = current_t + c1g*dt;
t_stage2 = current_t + c2g*dt;
% 内联拉格朗日插值计算
u1star0 = lagrange3_2D(t_stage1, t_nodes, u, u1, u2);
u2star0 = lagrange3_2D(t_stage2, t_nodes, u, u1, u2);
v1star0 = lagrange3_2D(t_stage1, t_nodes, v, v1, v2);
v2star0 = lagrange3_2D(t_stage2, t_nodes, v, v1, v2);
psi1star0 = lagrange3_2D(t_stage1, t_nodes, psi, psi1, psi2);
psi2star0 = lagrange3_2D(t_stage2, t_nodes, psi, psi1, psi2);
phi1star0 = lagrange3_2D(t_stage1, t_nodes, phi, phi1, phi2);
phi2star0 = lagrange3_2D(t_stage2, t_nodes, phi, phi1, phi2);
% 迭代初值
u1n_s = u1star0; psi1n_s = psi1star0;
u2n_s = u2star0; psi2n_s = psi2star0;
v1n_s = v1star0; phi1n_s = phi1star0;
v2n_s = v2star0; phi2n_s = phi2star0;
iterw = 0; err = inf;
% 迭代循环
while iterw < 50 && err > 1e-14
% 非线性项
p1n_s = a1*psi1n_s + b1*psi1n_s.^3 + c1*psi1n_s.*phi1n_s.^2;
p2n_s = a1*psi2n_s + b1*psi2n_s.^3 + c1*psi2n_s.*phi2n_s.^2;
q1n_s = a2*phi1n_s + b2*phi1n_s.^3 + c2*phi1n_s.*psi1n_s.^2;
q2n_s = a2*phi2n_s + b2*phi2n_s.^3 + c2*phi2n_s.*psi2n_s.^2;
% 列方程组求解
bw = cat(3,u - dt*a11*p1n_s - dt*a12*p2n_s,...
u - dt*a21*p1n_s - dt*a22*p2n_s,...
psi,...
psi);
bo = cat(3, v - dt*a11*q1n_s - dt*a12*q2n_s, ...
v - dt*a21*q1n_s - dt*a22*q2n_s, ...
phi, ...
phi);
% 使用FFT求解器
xw = solve_fft2(bw, invD, M);
xo = solve_fft2(bo, invD, M);
% 更新迭代值
u1n_new = xw(:,:,1);
u2n_new = xw(:,:,2);
psi1n_new = xw(:,:,3);
psi2n_new = xw(:,:,4);
v1n_new = xo(:,:,1);
v2n_new = xo(:,:,2);
phi1n_new = xo(:,:,3);
phi2n_new = xo(:,:,4);
% 计算误差
err_u1n = max(abs(u1n_new(:) - u1n_s(:)));
err_u2n = max(abs(u2n_new(:) - u2n_s(:)));
err_psi1n = max(abs(psi1n_new(:) - psi1n_s(:)));
err_psi2n = max(abs(psi2n_new(:) - psi2n_s(:)));
err_v1n = max(abs(v1n_new(:) - v1n_s(:)));
err_v2n = max(abs(v2n_new(:) - v2n_s(:)));
err_phi1n = max(abs(phi1n_new(:) - phi1n_s(:)));
err_phi2n = max(abs(phi2n_new(:) - phi2n_s(:)));
err = max([err_u1n, err_u2n, err_psi1n, err_psi2n,...
err_v1n, err_v2n, err_phi1n, err_phi2n]);
% 更新迭代值
u1n_s = u1n_new;
u2n_s = u2n_new;
psi1n_s = psi1n_new;
psi2n_s = psi2n_new;
v1n_s = v1n_new;
v2n_s = v2n_new;
phi1n_s = phi1n_new;
phi2n_s = phi2n_new;
iterw = iterw + 1;
end
% 输出最终迭代值
psi1_star = psi1n_s;
psi2_star = psi2n_s;
phi1_star = phi1n_s;
phi2_star = phi2n_s;
%% 计算非线性项
p1n_star = a1*psi1_star + b1*psi1_star.^3 + c1*psi1_star.*phi1_star.^2;
p2n_star = a1*psi2_star + b1*psi2_star.^3 + c1*psi2_star.*phi2_star.^2;
q1n_star = a2*phi1_star + b2*phi1_star.^3 + c2*phi1_star.*psi1_star.^2;
q2n_star = a2*phi2_star + b2*phi2_star.^3 + c2*phi2_star.*psi2_star.^2;
% 右端项
bb0 = cat(3, a^2*real(ifft2(lambda_B.*fft2(psi))), ...
a^2*real(ifft2(lambda_B.*fft2(psi))), ...
u, u);
bb1 = cat(3, -p1n_star, zeros(M), zeros(M), zeros(M));
bb2 = cat(3, zeros(M), -p2n_star, zeros(M), zeros(M));
bb3 = cat(3, -psi1_star, zeros(M), zeros(M), zeros(M));
bb4 = cat(3, zeros(M), -psi2_star, zeros(M), zeros(M));
% 使用FFT求解器
x0 = solve_fft2(bb0, invD, M);
x1 = solve_fft2(bb1, invD, M);
x2 = solve_fft2(bb2, invD, M);
x3 = solve_fft2(bb3, invD, M);
x4 = solve_fft2(bb4, invD, M);
% 右端项
yb0 = cat(3, a^2*real(ifft2(lambda_B.*fft2(phi))), ...
a^2*real(ifft2(lambda_B.*fft2(phi))), ...
v, v);
yb1 = cat(3, -q1n_star, zeros(M), zeros(M), zeros(M));
yb2 = cat(3, zeros(M), -q2n_star, zeros(M), zeros(M));
yb3 = cat(3, -phi1_star, zeros(M), zeros(M), zeros(M));
yb4 = cat(3, zeros(M), -phi2_star, zeros(M), zeros(M));
% 使用FFT求解器
y0 = solve_fft2(yb0, invD, M);
y1 = solve_fft2(yb1, invD, M);
y2 = solve_fft2(yb2, invD, M);
y3 = solve_fft2(yb3, invD, M);
y4 = solve_fft2(yb4, invD, M);
% 提取各部分
k1_1_0 = x0(:,:,1); k1_1_lamda1 = x1(:,:,1); k1_1_lamda2 = x2(:,:,1); k1_1_eita1 = x3(:,:,1); k1_1_eita2 = x4(:,:,1);
k1_2_0 = x0(:,:,2); k1_2_lamda1 = x1(:,:,2); k1_2_lamda2 = x2(:,:,2); k1_2_eita1 = x3(:,:,2); k1_2_eita2 = x4(:,:,2);
k2_1_0 = x0(:,:,3); k2_1_lamda1 = x1(:,:,3); k2_1_lamda2 = x2(:,:,3); k2_1_eita1 = x3(:,:,3); k2_1_eita2 = x4(:,:,3);
k2_2_0 = x0(:,:,4); k2_2_lamda1 = x1(:,:,4); k2_2_lamda2 = x2(:,:,4); k2_2_eita1 = x3(:,:,4); k2_2_eita2 = x4(:,:,4);
k3_1_0 = y0(:,:,1); k3_1_lamda1 = y1(:,:,1); k3_1_lamda2 = y2(:,:,1); k3_1_eita1 = y3(:,:,1); k3_1_eita2 = y4(:,:,1);
k3_2_0 = y0(:,:,2); k3_2_lamda1 = y1(:,:,2); k3_2_lamda2 = y2(:,:,2); k3_2_eita1 = y3(:,:,2); k3_2_eita2 = y4(:,:,2);
k4_1_0 = y0(:,:,3); k4_1_lamda1 = y1(:,:,3); k4_1_lamda2 = y2(:,:,3); k4_1_eita1 = y3(:,:,3); k4_1_eita2 = y4(:,:,3);
k4_2_0 = y0(:,:,4); k4_2_lamda1 = y1(:,:,4); k4_2_lamda2 = y2(:,:,4); k4_2_eita1 = y3(:,:,4); k4_2_eita2 = y4(:,:,4);
%% 求解优化问题(简化版)
% 在实际应用中,这里需要求解一个非线性系统来获得最优参数
% 为简化,我们使用固定值
lamda1 = 1; lamda2 = 1;
eita1 = 0; eita2 = 0;
kesi1 = 1; kesi2 = 1;
%% 组装得到K1(1,2),K2(1,2),K3(1,2),K4(1,2)
k1_1 = k1_1_0 + lamda1.*k1_1_lamda1 + lamda2.*k1_1_lamda2 + eita1.*k1_1_eita1 + eita2.*k1_1_eita2;
k1_2 = k1_2_0 + lamda1.*k1_2_lamda1 + lamda2.*k1_2_lamda2 + eita1.*k1_2_eita1 + eita2.*k1_2_eita2;
k2_1 = k2_1_0 + lamda1.*k2_1_lamda1 + lamda2.*k2_1_lamda2 + eita1.*k2_1_eita1 + eita2.*k2_1_eita2;
k2_2 = k2_2_0 + lamda1.*k2_2_lamda1 + lamda2.*k2_2_lamda2 + eita1.*k2_2_eita1 + eita2.*k2_2_eita2;
k3_1 = k3_1_0 + lamda1.*k3_1_lamda1 + lamda2.*k3_1_lamda2 + eita1.*k3_1_eita1 + eita2.*k3_1_eita2;
k3_2 = k3_2_0 + lamda1.*k3_2_lamda1 + lamda2.*k3_2_lamda2 + eita1.*k3_2_eita1 + eita2.*k3_2_eita2;
k4_1 = k4_1_0 + lamda1.*k4_1_lamda1 + lamda2.*k4_1_lamda2 + eita1.*k4_1_eita1 + eita2.*k4_1_eita2;
k4_2 = k4_2_0 + lamda1.*k4_2_lamda1 + lamda2.*k4_2_lamda2 + eita1.*k4_2_eita1 + eita2.*k4_2_eita2;
%% n+1层的解
unew = u + dt*b1g*k1_1 + dt*b2g*k1_2;
psinew = psi + dt*b1g*k2_1 + dt*b2g*k2_2;
vnew = v + dt*b1g*k3_1 + dt*b2g*k3_2;
phinew = phi + dt*b1g*k4_1 + dt*b2g*k4_2;
% 计算当前时间步的质量和能量
M_current = h^2 * sum(sum(abs(psinew).^2 + abs(phinew).^2));
term1 = h^2 * sum(sum(vnew.^2));
term2 = -a^2 * h^2 * sum(sum(real(ifft2(lambda_B .* fft2(psinew))) .* psinew));
term3 = a1 * h^2 * sum(sum(psinew.^2));
term4 = h^2 * sum(sum(unew.^2));
term5 = -a^2 * h^2 * sum(sum(real(ifft2(lambda_B .* fft2(phinew))) .* phinew));
term6 = a2 * h^2 * sum(sum(phinew.^2));
term7 = (c2*b1/(4*c1)) * h^2 * sum(sum(psinew.^4));
term8 = (b2/4) * h^2 * sum(sum(phinew.^4));
term9 = (c2/2) * h^2 * sum(sum((psinew.*phinew).^2));
E_current = (1/2)*((c2/c1)*(term1 + term2 + term3) + term4 + term5 + term6) + term7 + term8 + term9;
% 记录守恒量及误差
M_arr(n+1) = M_current;
E_arr(n+1) = E_current;
abs_err_M(n+1) = abs(M_current - M0);
abs_err_E(n+1) = abs(E_current - E0);
% 更新历史变量
prev_psi1 = psi + dt*(a11*k2_1 + a12*k2_2);
prev_psi2 = psi + dt*(a21*k2_1 + a22*k2_2);
prev_phi1 = phi + dt*(a11*k4_1 + a12*k4_2);
prev_phi2 = phi + dt*(a21*k4_1 + a22*k4_2);
prev_u1 = u + dt*(a11*k1_1 + a12*k1_2);
prev_u2 = u + dt*(a21*k1_1 + a22*k1_2);
prev_v1 = v + dt*(a11*k3_1 + a12*k3_2);
prev_v2 = v + dt*(a21*k3_1 + a22*k3_2);
% 更新n+1层的值
prev_u = unew;
prev_psi = psinew;
prev_v = vnew;
prev_phi = phinew;
prev_t = n*dt;
% 每10步显示进度
if mod(n, 10) == 0
fprintf('完成时间步:%d/%d, 质量误差:%.4e, 能量误差:%.4e\n', ...
n, Nt, abs_err_M(n+1), abs_err_E(n+1));
end
end
%% 结果可视化
% 绘制守恒量误差图
figure;
semilogy((0:Nt)*dt, abs_err_M, 'b', 'LineWidth', 1.5);
xlabel('t');
ylabel('绝对误差 (log scale)');
title('质量守恒误差');
grid on;
figure;
semilogy((0:Nt)*dt, abs_err_E, 'r', 'LineWidth', 1.5);
xlabel('t');
ylabel('绝对误差 (log scale)');
title('能量守恒误差');
grid on;
% 绘制初始和最终波形
figure;
subplot(2,2,1);
surf(x, y, psi0, 'EdgeColor', 'none');
title('初始 \psi');
xlabel('x'); ylabel('y');
subplot(2,2,2);
surf(x, y, psinew, 'EdgeColor', 'none');
title('最终 \psi');
xlabel('x'); ylabel('y');
subplot(2,2,3);
surf(x, y, phi0, 'EdgeColor', 'none');
title('初始 \phi');
xlabel('x'); ylabel('y');
subplot(2,2,4);
surf(x, y, phinew, 'EdgeColor', 'none');
title('最终 \phi');
xlabel('x'); ylabel('y');
end
%% 二维FFT求解器函数
function x = solve_fft2(b, invD, N)
% 输入b是N×N×4的三维数组
% 输出x是N×N×4的三维数组
% FFT变换
B1 = fft2(b(:,:,1));
B2 = fft2(b(:,:,2));
B3 = fft2(b(:,:,3));
B4 = fft2(b(:,:,4));
% 预分配解的空间
X1 = zeros(N, N);
X2 = zeros(N, N);
X3 = zeros(N, N);
X4 = zeros(N, N);
% 频域求解
for i = 1:N
for j = 1:N
rhs = [B1(i,j); B2(i,j); B3(i,j); B4(i,j)];
sol = invD(:,:,i,j) * rhs;
X1(i,j) = sol(1);
X2(i,j) = sol(2);
X3(i,j) = sol(3);
X4(i,j) = sol(4);
end
end
% IFFT变换回物理空间
x1 = real(ifft2(X1));
x2 = real(ifft2(X2));
x3 = real(ifft2(X3));
x4 = real(ifft2(X4));
% 组合解
x = cat(3, x1, x2, x3, x4);
end
%% 二维三点拉格朗日插值
function y = lagrange3_2D(t, t_nodes, y0, y1, y2)
t0 = t_nodes(1); t1 = t_nodes(2); t2 = t_nodes(3);
L0 = ((t - t1).*(t - t2)) ./ ((t0 - t1).*(t0 - t2));
L1 = ((t - t0).*(t - t2)) ./ ((t1 - t0).*(t1 - t2));
L2 = ((t - t0).*(t - t1)) ./ ((t2 - t0).*(t2 - t1));
y = y0.*L0 + y1.*L1 + y2.*L2;
end
最新发布