【数据分析】基于能量离散化的耦合分数阶Gross-Pitaevskii方程能量守恒格式的MATLAB代码

✅作者简介热爱科研的Matlab仿真开发者擅长毕业设计辅导、数学建模、数据处理、程序设计科研仿真。完整代码获取 定制创新 论文复现点击Matlab科研工作室 关注我领取海量matlab电子书和数学建模资料个人信条做科研博学之、审问之、慎思之、明辨之、笃行之是为博学慎思明辨笃行。 内容介绍一、引言耦合分数阶 Gross - Pitaevskii 方程CFGPE在描述诸如玻色 - 爱因斯坦凝聚BEC、非线性光学等物理现象中具有重要应用。然而数值求解 CFGPE 时保证能量守恒对于准确模拟物理过程至关重要。基于能量离散化构建能量守恒格式能够有效维持系统能量特性为研究相关物理现象提供可靠的数值方法。二、耦合分数阶 Gross - Pitaevskii 方程概述一方程形式二物理意义该方程描述了多分量系统中粒子间的相互作用以及外部势场对波函数的影响。分数阶拉普拉斯算子的引入使得方程能够更准确地描述具有长程相互作用或非局部特性的物理系统如在一些复杂的凝聚态物理和光学系统中这种非局部特性起着关键作用。三、能量离散化方法一能量表达式推导⛳️ 运行结果 部分代码%%Computes the **mass and energy errors** for the Coupled Fractional Gross-Pitaevskii Equations under different numerical settings.function [ErrM, ErrE] GP_ErrM_ErrE(alpha_u,alpha_v,tau,h)% Initialize problem domaina -10; b 10; % Spatial domain [a, b]T0.5; % Final simulation timeNfix(T/tau); % Number of time stepsMfix((b-a)/h); % Number of spatial stepst0:tau:T; % Time vectorxa:h:b; % Spatial grid% System parametersD0.2;beta_111;beta_221;beta_120.4;lambda-0.3;V(1/2*x(1:M).^2).;% Initial conditions (Gaussian functions)phi_1 1/(sqrt(pi))*exp(-x.^2);phi_2 1/(sqrt(pi))*exp(-x.^2);u_tempphi_1(1:M).; % Initial condition for uv_tempphi_2(1:M).; % Initial condition for vuv_temp[u_temp;v_temp]; % Combined state vector% Construct A (Fractional Laplacian)rones(1,M);for j1:Mr(j)1/(2*M)*(abs(-M/2*(2*pi/(b-a))))^alpha_u*exp(1i*(-M/2*(2*pi/(b-a)))*((j-1)*h))...1/(2*M)*(abs(M/2*(2*pi/(b-a))))^alpha_u*exp(1i*(M/2*(2*pi/(b-a)))*((j-1)*h));for k-M/21:M/2-1r(j)r(j)1/M*(abs(k*(2*pi/(b-a))))^alpha_u*exp(1i*(k*(2*pi/(b-a)))*(j-1)*h);endendA_utoeplitz(r);for j1:Mr(j)1/(2*M)*(abs(-M/2*(2*pi/(b-a))))^alpha_v*exp(1i*(-M/2*(2*pi/(b-a)))*((j-1)*h))...1/(2*M)*(abs(M/2*(2*pi/(b-a))))^alpha_v*exp(1i*(M/2*(2*pi/(b-a)))*((j-1)*h));for k-M/21:M/2-1r(j)r(j)1/M*(abs(k*(2*pi/(b-a))))^alpha_v*exp(1i*(k*(2*pi/(b-a)))*(j-1)*h);endendA_vtoeplitz(r);% Construct D1 (Fractional derivative operator)dones(1,M);for j1:Md(j)1/(2*M)*(abs(-M/2*(2*pi/(b-a))))^(alpha_u/2)*exp(1i*(-M/2*(2*pi/(b-a)))*((j-1)*h))...1/(2*M)*(abs(M/2*(2*pi/(b-a))))^(alpha_u/2)*exp(1i*(M/2*(2*pi/(b-a)))*((j-1)*h));for k-M/21:M/2-1d(j)d(j)1/M*(abs(k*(2*pi/(b-a))))^(alpha_u/2)*exp(1i*(k*(2*pi/(b-a)))*(j-1)*h);endendD1_utoeplitz(d);for j1:Md(j)1/(2*M)*(abs(-M/2*(2*pi/(b-a))))^(alpha_v/2)*exp(1i*(-M/2*(2*pi/(b-a)))*((j-1)*h))...1/(2*M)*(abs(M/2*(2*pi/(b-a))))^(alpha_v/2)*exp(1i*(M/2*(2*pi/(b-a)))*((j-1)*h));for k-M/21:M/2-1d(j)d(j)1/M*(abs(k*(2*pi/(b-a))))^(alpha_v/2)*exp(1i*(k*(2*pi/(b-a)))*(j-1)*h);endendD1_vtoeplitz(d);%First Time Step Calculationu_Iu_temp-1i*tau*(1/2*A_u*u_temp)-1i*tau*(Vones(M,1).*Dbeta_11*abs(u_temp).^2beta_12*abs(v_temp).^2).*u_temp...-1i*tau*lambda*v_temp;v_Iv_temp-1i*tau*(1/2*A_v*v_temp)-1i*tau*(Vbeta_12*abs(u_temp).^2beta_22*abs(v_temp).^2).*v_temp...-1i*tau*lambda*u_temp;uv_I[u_I;v_I];% Compute mass at t 0 and t 1M1zeros(1,N1);M1(1)h*(1/2*(abs(u_temp(1))^2abs(v_temp(1))^2)sum(abs(u_temp(2:M-1)).^2)sum(abs(v_temp(2:M-1)).^2)1/2*(abs(u_temp(M))^2abs(v_temp(M))^2));M1(2)h*(1/2*(abs(u_I(1))^2abs(v_I(1))^2)sum(abs(u_I(2:M-1)).^2)sum(abs(v_I(2:M-1)).^2)1/2*(abs(u_I(M))^2abs(v_I(M))^2));% Compute mass conservation errorErrM(1)abs(M1(1)-M1(1));ErrM(2)abs(M1(2)-M1(1));% Compute energy at t 0 and t 1Gzeros(1,M);G1/2*abs(D1_u*u_temp).^21/2*abs(D1_v*v_temp).^2V.*(abs(u_temp).^2abs(v_temp).^2)...D.*abs(u_temp).^21/2*beta_11*abs(u_temp).^41/2*beta_22*abs(v_temp).^4beta_12*(abs(u_temp).^2).*(abs(v_temp).^2)...lambda*(u_temp.*conj(v_temp)v_temp.*conj(u_temp));E(1)h*(1/2*G(1)sum(G(2:M-1))1/2*G(M));% Compute energy conservation errorErrE(1)abs(E(1)-E(1));G1/2*abs(D1_u*u_I).^21/2*abs(D1_v*v_I).^2V.*(abs(u_I).^2abs(v_I).^2)...D.*abs(u_I).^21/2*beta_11*abs(u_I).^41/2*beta_22*abs(v_I).^4beta_12*(abs(u_I).^2).*(abs(v_I).^2)...lambda*(u_I.*conj(v_I)v_I.*conj(u_I));E(2)h*(1/2*G(1)sum(G(2:M-1))1/2*G(M));ErrE(2)abs(E(2)-E(1));% Construct the third time step calculationIones(M,1)*(1i/(2*tau));for n2:Na1I-1/4*V-1/4*ones(M,1)*D-1/4*beta_11*abs(u_I).^2-beta_12/4*(abs(v_I).^2);a21/2*V1/2*ones(M,1)*D1/2*beta_11*abs(u_I).^2beta_12/2*(abs(v_I).^2);a3I1/4*V1/4*ones(M,1)*D1/4*beta_11*abs(u_I).^2beta_12/4*(abs(v_I).^2);A1diag(a1,0)-1/8.*A_u;A2diag(a2,0)1/4.*A_u;A3diag(a3,0)1/8.*A_u;A4diag(ones(M,1)*lambda,0);b1I-1/4*V-1/4*beta_22*abs(v_I).^2-beta_12/4*(abs(u_I).^2);b21/2*V1/2*beta_22*abs(v_I).^2beta_12/2*(abs(u_I).^2);b3I1/4*V1/4*beta_22*abs(v_I).^2beta_12/4*(abs(u_I).^2);B1diag(b1,0)-1/8.*A_v;B2diag(b2,0)1/4.*A_v;B3diag(b3,0)1/8.*A_v;AB1[A1,A4.*(-1/4);A4.*(-1/4),B1];AB2[A2,A4.*(1/2);A4.*(1/2),B2];AB3[A3,A4.*(1/4);A4.*(1/4),B3];uv_IIAB1\( AB2*uv_IAB3*uv_temp);% Mass and energy updateM1(n1)h*(1/2*(abs(uv_II(1))^2abs(uv_II(M1))^2)sum(abs(uv_II(2:M-1)).^2)sum(abs(uv_II(M2:2*M-1)).^2)1/2*(abs(uv_II(M))^2abs(uv_II(2*M))^2));ErrM(n1)abs(M1(n1)-M1(1));G1/2*abs(D1_u*uv_II(1:M)).^21/2*abs(D1_v*uv_II(M1:2*M)).^2 ...V.*(abs(uv_II(1:M)).^2abs(uv_II(M1:2*M)).^2)...D.*abs(uv_II(1:M)).^21/2*beta_11*abs(uv_II(1:M)).^41/2*beta_22*abs(uv_II(M1:2*M)).^4beta_12*(abs(uv_II(1:M)).^2).*(abs(uv_II(M1:2*M)).^2)...lambda*(uv_II(1:M).*conj(uv_II(M1:2*M))uv_II(M1:2*M).*conj(uv_II(1:M)));E(n1)h*(1/2*G(1)sum(G(2:M-1))1/2*G(M));ErrE(n1)abs(E(n1)-E(1));% Update time layersuv_tempuv_I;uv_Iuv_II;u_Iuv_I(1:M); v_Iuv_I(M1:2*M);endErrMErrM(N1);ErrEErrE(N1); 参考文献更多免费数学建模和仿真教程关注领取