doubleu.m

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

clear all;
clc;

deltat=0.01;
ST=-0.6; 	%开始时间
ET=0.4;	%终止时间

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

y=1:3;	%y序列,和x相同
[X,Y]=meshgrid(x,y);

lambda1=3;
r1=-10;
alpha1=-sqrt(0.25*lambda1*lambda1-1);
beta1=-4*(2+lambda1*lambda1);

lambda2=4;
r2=-20;
alpha2=-sqrt(0.25*lambda2*lambda2-1);
beta2=-4*(2+lambda2*lambda2);

lambda3=5;
r3=-31.5;
alpha3=-sqrt(0.25*lambda3*lambda3-1);
beta3=-4*(2+lambda3*lambda3);

xi1=-alpha1*(X+beta1*ST);
deltapt1=-alpha1*beta1*deltat;

xi2=-alpha2*(X+beta2*ST);
deltapt2=-alpha2*beta2*deltat;

xi3=-alpha3*(X+beta3*ST);
deltapt3=-alpha3*beta3*deltat;

Delta=zeros(size(X));
DeltaC=Delta;
DeltaD=Delta;
DeltaA=Delta;

figure(1);
pause(3);

for t=ST:deltat:ET
    sigma1=alpha1*(tanh(xi1)-r1)./(1-r1*tanh(xi1))+0.5*lambda1;
    sigma2=alpha2*(tanh(xi2)-r2)./(1-r2*tanh(xi2))+0.5*lambda2;    
    sigma3=alpha3*(tanh(xi3)-r3)./(1-r3*tanh(xi3))+0.5*lambda3;

    for j=1:size(X,2)
        temp=[1 sigma1(1,j) lambda1*sigma1(1,j);1 sigma2(1,j) lambda2*sigma2(1,j);1 sigma3(1,j) lambda3*sigma3(1,j)];
        Delta(:,j)=det(temp);
 
        temp=[-lambda1 sigma1(1,j) lambda1*sigma1(1,j);-lambda2 sigma2(1,j) lambda2*sigma2(1,j);-lambda3 sigma3(1,j) lambda3*sigma3(1,j)];
        DeltaC(:,j)=det(temp);
 
        temp=[1 sigma1(1,j) -lambda1;1 sigma2(1,j) -lambda2;1 sigma3(1,j) -lambda3];
        DeltaD(:,j)=det(temp);
 
        temp=[1 sigma1(1,j) -lambda1*lambda1;1 sigma2(1,j) -lambda2*lambda2;1 sigma3(1,j) -lambda3*lambda3];
        DeltaA(:,j)=det(temp);
    end
    
    DeltaC=DeltaC./Delta;
    DeltaA=DeltaA./DeltaD;
    DeltaD=DeltaD./Delta;
    
    u=DeltaC+DeltaD+DeltaA;
    
    %surf(X,Y,u);
    plot(u(1,:));
    %view(30,40);
    title(t);
    axis([0 30 -200 200]);
	xlabel('X');
	ylabel('Y');
	zlabel('Z');   
    pause;
    
    xi1=xi1+deltapt1;
    xi2=xi2+deltapt2;
    xi3=xi3+deltapt3;
end