doubleuw.m

这是我和同事门一起努力编的的孤立子程序。基本上所以孤立子的情况都考虑到了

clear all;
clc;

UWFlag=1;   %等于1时计算u，其它值时计算w；

deltat=1.0; %时间分辨率
ST=-10.0; 	%开始时间
ET=10.0;	%终止时间

SX=-50.0;	%x开始
EX=50.0;	%x终止
x=SX:1:EX;	%x序列

y=1:3;	%y
[X,Y]=meshgrid(x,y);

lambda1=3;  %参数
r1=-11;
beta1=sqrt(1+lambda1*lambda1);

lambda2=6;  %参数
r2=-15;
beta2=sqrt(1+lambda2*lambda2);

lambda3=10; %参数
r3=-17;
beta3=sqrt(1+lambda3*lambda3);

xi1=beta1*(X+lambda1*ST);   %初值
deltapt1=beta1*lambda1*deltat;

xi2=beta2*(X+lambda2*ST);   %初值
deltapt2=beta2*lambda2*deltat;

xi3=beta3*(X+lambda3*ST);   %初值
deltapt3=beta3*lambda3*deltat;

Delta=zeros(size(X));
DeltaX=Delta;
DeltaA=Delta;
DeltaAX=Delta;

if UWFlag==1
    DeltaT=Delta;
    DeltaTX=Delta;
    DeltaAXX=Delta;
    DeltaXX=Delta;
end 

for t=ST:deltat:ET
    sigma1=beta1*(tanh(xi1)-r1)./(1-r1*tanh(xi1))+lambda1;
    sigma2=beta2*(tanh(xi2)-r2)./(1-r2*tanh(xi2))+lambda2;
    sigma3=beta3*(tanh(xi3)-r3)./(1-r3*tanh(xi3))+lambda3;
    
    for j=1:size(X,2)
        temp1=[1 lambda1 sigma1(1,j);1 lambda2 sigma2(1,j);1 lambda3 sigma3(1,j)];
        Delta(:,j)=det(temp1);
        
        temp1=[1 lambda1 sigma1(1,j)*(2*lambda1-sigma1(1,j));1 lambda2 sigma2(1,j)*(2*lambda2-sigma2(1,j));1 lambda3 sigma3(1,j)*(2*lambda3-sigma3(1,j))];
        DeltaX(:,j)=det(temp1);
        
        temp1=[1 lambda1*(sigma1(1,j)-2*lambda1) sigma1(1,j)*(2*lambda1-sigma1(1,j));1 lambda2*(sigma2(1,j)-2*lambda2) sigma2(1,j)*(2*lambda2-sigma2(1,j));1 lambda3*(sigma3(1,j)-2*lambda3) sigma2(1,j)*(2*lambda2-sigma2(1,j))];
        temp2=[1 lambda1*(1+sigma1(1,j)*(2*lambda1-sigma1(1,j))) sigma1(1,j);1 lambda2*(1+sigma2(1,j)*(2*lambda2-sigma2(1,j)))  sigma2(1,j);1 lambda3*(1+sigma3(1,j)*(2*lambda3-sigma3(1,j)))  sigma3(1,j)];
        DeltaAX(:,j)=det(temp1)+det(temp2);
        
        temp1=[1 lambda1*sigma1(1,j)-2*lambda1*lambda1 sigma1(1,j);1 lambda2*sigma2(1,j)-2*lambda2*lambda2 sigma2(1,j);1 lambda3*sigma3(1,j)-2*lambda3*lambda3 sigma3(1,j)];
        DeltaA=det(temp1);
        
        if UWFlag==1
            temp1=[1 sigma1(1,j) lambda1*sigma1(1,j);1 sigma2(1,j) lambda3*sigma2(1,j);1 sigma3(1,j) lambda3*sigma3(1,j)];
            DeltaT(:,j)=det(temp1);
        
            temp1=[1 sigma1(1,j)*(2*lambda1-sigma1(1,j)) lambda1*sigma1(1,j);1 sigma2(1,j)*(2*lambda2-sigma2(1,j)) lambda2*sigma2(1,j);1 sigma3(1,j)*(2*lambda3-sigma3(1,j)) lambda3*sigma3(1,j)];
            temp2=[1 sigma1(1,j) lambda1*(1+sigma1(1,j)*(2*lambda1-sigma1(1,j)));1 sigma2(1,j) lambda2*(1+sigma2(1,j)*(2*lambda2-sigma2(1,j)));1 sigma3(1,j) lambda3*(1+sigma3(1,j)*(2*lambda3-sigma3(1,j)))];
            DeltaTX(:,j)=det(temp1)+det(temp2);
            
            temp1=[1 lambda1*(sigma1(1,j)-2*lambda1) 2*(lambda1-sigma1(1,j))*(1+sigma1(1,j)*(2*lambda1-sigma1(1,j)));1 lambda2*(sigma2(1,j)-2*lambda2) 2*(lambda2-sigma2(1,j))*(1+sigma2(1,j)*(2*lambda2-sigma2(1,j)));1 lambda3*(sigma3(1,j)-2*lambda3) 2*(lambda3-sigma3(1,j))*(1+sigma3(1,j)*(2*lambda3-sigma3(1,j)))];
            temp2=[1 2*lambda1*(lambda1-sigma1(1,j))*(1+sigma1(1,j)*(2*lambda1-sigma1(1,j))) sigma1(1,j);1 2*lambda2*(lambda2-sigma2(1,j))*(1+sigma2(1,j)*(2*lambda2-sigma2(1,j))) sigma2(1,j);1 2*lambda3*(lambda3-sigma3(1,j))*(1+sigma3(1,j)*(2*lambda3-sigma3(1,j))) sigma3(1,j)];
            DeltaAXX(:,j)=det(temp1)+det(temp2);
        
            temp1(:,2)=[lambda1 lambda2 lambda3]';
            DeltaXX(:,j)=det(temp1);
        end
    end
     
    if UWFlag~=1    
        w=(Delta.*(Delta-DeltaAX)+DeltaA.*DeltaX)./(Delta.*Delta);
        %surf(X,Y,w);
        plot(w(1,:));
    else
        u=0.5*(DeltaX./Delta-DeltaTX./DeltaT)-1/6*(DeltaAXX.*Delta-2*DeltaAX.*DeltaX-DeltaA.*DeltaXX+2*DeltaA.*DeltaX.*DeltaX./Delta)./(Delta.*Delta-DeltaAX.*Delta+DeltaA.*DeltaX);
        %surf(X,Y,u);
        plot(u(1,:));
    end
    
    title(t);
	xlabel('X');
	ylabel('Y');
	zlabel('Z');   
    pause(2);
    
    xi1=xi1+deltapt1;   %增量
    xi2=xi2+deltapt2;   %增量
    xi3=xi3+deltapt3;   %增量
end