c_ngeshi.m

用C_N格式求解耦合非线性薛定谔方程,matlab程序

clear
%定义时间步长，空间步长
t=0.05/2/2/2/2;h=0.4/2/2/2/2;r=t/(h^2); J=90/h;dt=t;%
epc=1e-6;%迭代精度
u0=zeros(1,J+1);
arf=0.2;v=0.5;baita=1;%各参数取值
u0=sqrt(arf)*sech(sqrt(arf)*x).*exp(i*v*x);
plot3(x,ones(1,J+1)*0,abs(u0));
hold on;
v0=u0;%初值相同
Q0=sum(abs(u0).*abs(u0))*h;%初始电荷
s2=0;
   for m=1:J
       au2=abs(u0(m))^2;av2=abs(v0(m))^2;
       s2=s2+1/(h*h)*abs(u0(m+1)-u0(m))*abs(u0(m+1)-u0(m))+1/(h*h)*abs(v0(m+1)-v0(m))*abs(v0(m+1)-v0(m))-(au2^2+av2^2)/2-au2*av2;
   end
   s2=s2+abs(u0(J+1))*abs(u0(J+1))/(h*h)+abs(v0(J+1))*abs(v0(J+1))/(h*h)-...
          (abs(u0(J+1))^4+abs(v0(J+1))^4)/2-abs(u0(J+1))^2*abs(v0(J+1))^2 ;
s2=s2*h/2;
E0=s2;%初始能量
us=u0;vs=v0;%v的初值跟一样
UU=[];T=0.5;diedaicishu=zeros(1,T/t);
%以下for循环为时间层上的剖分
time=cputime;
for n=1:T/t%N不断变化，可以不断的更新N   
    e=1;
    while(e>=epc)%对每一时间层上，用迭代法逐渐逼近解
        %以下为追赶法求解差分方程组
        uvnorm2u=norm(u0(1))^2+norm(v0(1))^2+norm(us(1))^2+norm(vs(1))^2;
        uvnorm2v=norm(u0(1))^2+norm(v0(1))^2+norm(us(1))^2+norm(vs(1))^2;
        du=-r/2*(-2*u0(1)+u0(2))+i*u0(1)-t*uvnorm2u*(us(1)+u0(1))/4;
        dv=-r/2*(-2*v0(1)+v0(2))+i*v0(1)-t*uvnorm2v*(vs(1)+v0(1))/4;
        p=b;q(1)=-c/b;
        u1(1)=du/p;
        v1(1)=dv/p;
        %%%%%%%%%%%%%%%%%%%%%%%以下求解u%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        for j=2:J  %追赶法过程
            p=a*q(j-1)+b;
            q(j)=-c/p;
            uvnorm2u=norm(u0(j))^2+norm(v0(j))^2+norm(us(j))^2+norm(vs(j))^2;
            uvnorm2v=norm(u0(j))^2+norm(v0(j))^2+norm(us(j))^2+norm(vs(j))^2;
            du=-r/2*(u0(j-1)-2*u0(j)+u0(j+1))+i*u0(j)-t/4*uvnorm2u*(us(j)+u0(j));
            dv=-r/2*(v0(j-1)-2*v0(j)+v0(j+1))+i*v0(j)-t/4*uvnorm2v*(vs(j)+v0(j));
            u1(j)=(du-a*u1(j-1))/p;
            v1(j)=(dv-a*v1(j-1))/p;
        end
        p=a*q(J)+b;
        q(J+1)=-c/p;
        uvnorm2u=norm(u0(J+1))^2+norm(v0(J+1))^2+norm(us(J+1))^2+norm(vs(J+1))^2;
        uvnorm2v=norm(u0(J+1))^2+norm(v0(J+1))^2+norm(us(J+1))^2+norm(vs(J+1))^2;
        du=-r/2*(u0(J)-2*u0(J+1))+i*u0(J+1)-t/4*uvnorm2u*(us(J+1)+u0(J+1));
        dv=-r/2*(v0(J)-2*v0(J+1))+i*v0(J+1)-t/4*uvnorm2v*(vs(J+1)+v0(J+1));
        u1(J+1)=(du-a*u1(J))/p;
        v1(J+1)=(dv-a*v1(J))/p;     
        for j=J:-1:1
            u1(j)=q(j)*u1(j+1)+u1(j);
            v1(j)=q(j)*v1(j+1)+v1(j);
        end
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        e=norm(us-u1)^2+norm(vs-v1)^2;
        us=u1; vs=v1;diedaicishu(n)=diedaicishu(n)+1;
    end
    u0=u1;v0=v1; 
    exactu=sqrt(arf)*sech(sqrt(arf)*(x-2*v*n*t)).*exp(i*(v*x-(v^2-arf)*n*t));
    uerror(n)=max(abs((exactu-u1)));
    UU(:,n)=u0';
    VV(:,n)=v0';
end
err=uerror(n);%最后时刻的误差
finaltime=cputime-time;%计算总时间