企业官网建设流程全解析

在这种情况下，势只是一个乘法算符，即可以直接指数化的对角矩阵。演化的步骤是：知道t时刻的波函数，将波函数傅立叶变换到动量空间，乘以动量算符，然后逆变换到坐标空间，乘以势能项，接着变换到动量空间，乘以动量算符，最后变换回来，即可得到t+时刻的波函数 。

clear all,clc; % 空间网格 % 空间总长度2457.6a.u.,步长0.15a.u. xmin = -1228.8; xmax = 1228.8; Nx = 16384; dx = (xmax - xmin) / Nx; x = linspace(xmin, xmax, Nx)'; % 时间网格 % 总时间1600a.u.,步长0.1a.u. T = 800; % 总时间 (a.u.) dt = 0.1; % 时间步长 (a.u.) Nt = round(T / dt); % 时间步数 t = (0:Nt-1)' * dt; % 动量空间网格 % 为了分裂算符法中动量算符直接与波函数作用，对动量空间作fftshift dp = 2*pi / (xmax - xmin); p = (-Nx/2:Nx/2-1)' * dp; p = fftshift(p); % 调整为自然序 % 原子参数（软库仑势） a = 1; % 原子序数 b = 1; % 软化参数 Vatom = -a./ sqrt(x.^2 + b); % 激光脉冲 % E = E0·A(t)·cos(wt+ø),取峰值强度1e14W/cm2,中心波长800nm，包络为高斯型 a0 = 5.29177210903e-11; % 玻尔半径 (m) epsilon = 8.854187817e-12; % 真空介电常数(F/m) Eh = 4.3597447222071e-18; % 哈特里能量 (J) c = 299792458; % 光速(m/s) hbar = 1.054571817e-34; % 约化普朗克常数 (J·s) au_time = hbar/Eh; % 原子时间单位 (s) % 激光参数 peak_intensity = 1e14; % 峰值强度 (W/cm²) FWHM = 10e-15; % FWHM wavelength = 800e-9; FWHM = FWHM / au_time; lambda = wavelength / a0; c0 = c/2.1876912633/1e6; w0 = 2*pi*c0 / lambda; % 中心频率(a.u.) E0 = sqrt(2*peak_intensity*1e4/epsilon/c)/5.1422e11; %电场振幅(a.u.) tau = FWHM / (2*sqrt(2*log(2))); E = E0 * exp(-(t-T/2).^2 / (2 * tau^2)*4*log(2)) .* cos(w0 * (t-T/2)); %电场形式 % 吸收边界参数 % 左右各50a.u.范围内快速衰减 x0 = 50; absorb = ones(Nx, 1); idx_left = find(x < x(1) + x0, 1, 'last'); idx_right = find(x > x(end) - x0, 1, 'first'); absorb(1:idx_left) = exp(-(x(1) + x0 - x(1:idx_left)).^2 / (2*x0^2)); absorb(idx_right:end) = exp(-(x(idx_right:end) - (x(end) - x0)).^2 / (2*x0^2)); %% 虚时间演化求基态波函数 max_steps = 400; % 演化时间 dt = 0.1; sigma = 10; % 初始化波函数 psi = exp(-x.^2/(2*sigma^2)); norm_x = sqrt(sum(abs(psi).^2*dx)); psi = (psi/norm_x); % 归一化 % 虚时间演化 for i = 1:max_steps % 动能演化 psi_p = fft(psi); psi_p = psi_p .* exp(-dt*p.^2/4); psi = ifft(psi_p); % 势能演化 psi = psi .* exp(-dt*Vatom); % 动能演化 psi_p = fft(psi); psi_p = psi_p .* exp(-dt*p.^2/4); psi = ifft(psi_p); % 归一化波函数 norm_psi = sqrt(sum(abs(psi).^2)*dx); psi = psi / norm_psi; end %% 计算基态能量 % 取实部并归一化 psi = real(psi); norm_psi = sqrt(sum(psi.^2)*dx); psi = psi / norm_psi; % 势能期望值 V_expect = sum(psi.^2 .* Vatom) * dx; % 动能期望值 psi_p = fft(psi); d2_psi = ifft(-p.^2 .* psi_p); T_expect = -0.5 * real(sum(conj(psi).*d2_psi)) * dx; % 总能量 E_ground = T_expect + V_expect; fprintf('基态能量 = %.4f\n', E_ground); %% % 分裂算符时间演化 % 预计算偶极矩 dipole = zeros(Nt, 1); for n = 1:Nt V_total = Vatom - x * E(n); % 总势能 % 画图，势场变化和波函数实部变化 % subplot(2,1,1) % plot(x,V_total); % xlim([-100,100]);ylim([-3,1]) % pause(0.01); % subplot(2,1,2) % plot(x,real(psi)); % xlim([-10,10]);ylim([-2,1]); pause(0.01); % 时间演化 psi_p = fft(psi); psi_p = psi_p .* exp(-1i * dt/2 * p.^2 / 2); psi = ifft(psi_p); psi = psi .* exp(-1i * dt/2 * V_total); psi_p = fft(psi); psi_p = psi_p .* exp(-1i * dt/2 * p.^2 / 2); psi = ifft(psi_p); % 吸收边界 psi = psi .* absorb; % 偶极矩 <x> = ∫ψ* x ψ dx dipole(n) = trapz(x, conj(psi) .* x .* psi); end % 加窗函数减少频谱泄漏 window = hann(Nt); dipole = dipole .* window; % 傅里叶变换 N_fft = 2^nextpow2(2*Nt); dw = fft(dipole, N_fft); w = 2*pi*(0:N_fft-1)' / (N_fft * dt); % 中心角频率 (a.u.) % 转换为谐波级次 harmonic_order = w / w0; % 结果可视化 % 激光电场 figure('Name', 'Laser Field'); plot(t, E, 'b', 'LineWidth', 1.5); xlabel('Time (a.u.)'); ylabel('Electric Field (a.u.)'); title('Laser Field'); grid on; % 偶极矩 figure('Name', 'Dipole'); plot(t, real(dipole), 'r', 'LineWidth', 1.5); xlabel('Time (a.u.)'); ylabel('Dipole (a.u.)'); title('Dipole'); grid on; % 高次谐波谱 figure('Name', 'HHG Spectrum'); max_harmonic = 100; % 显示的最大谐波级次 final_order = find(harmonic_order <= max_harmonic, 1, 'last'); semilogy(harmonic_order(1:final_order), abs(dw(1:final_order)).^2, 'k', 'LineWidth', 1.5); xlabel('Harmonic Order (w/w0)'); ylabel('Intensity (a.u.)'); title('High-Harmonic Generation Spectrum(I=1e15W/cm2)'); xlim([0 max_harmonic]); grid on;