doublevk1k2.asv

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

clear all;
clc;

k1=1.0;
k2=-1.0;
w1=k1^2/2;
w2=k2^2/2;
p1=w1*k1/2;
p2=w2*k2/2;
A=-(k1-k2)^2/(k1*k2);

deltat=1;
ST=-10.0; 	%开始时间
ET=10.0;		%终止时间

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

y=x;	%y序列,和x相同

z10=-5.0;	%z1初值
z20=5.0;	%z2初值

[X,Y]=meshgrid(x,y);

q1=exp(k1*X+w1*Y+z10+p1*ST);
r1=exp(k2*X+w2*Y+z20+p2*ST);
deltapt1=exp(p1*deltat);
deltapt2=exp(p2*deltat);

%qa=A*ones(size(X));
qa=A;
q2a=qa*qa;

j=1;
figure(1);
pause(5);

for t=ST:deltat:ET
    q2=q1.*q1;
    q3=q2.*q1;
    q4=q3.*q1;
    
    r2=r1.*r1;
    r3=r2.*r1;
    r4=r3.*r1;
    
    s=(q1+r1+q1.*r1*qa).*(1+q1+r1);
    s=1./(s.*s);
    v =-0.5*s.*(3*qa*(q2.*r3)*k2^2+4*qa*(q2.*r2)*k1^2-4*k1*(q2.*r1)*k2-4*k1*(q1.*r2)*k2+2*qa*(q3.*r1)*k2^2-4*k1*(q3.*r1)*k2+4*qa*(q2.*r3)*k1^2-4*k1*(q2.*r3)*k2*qa-2*q2a*(q3.*r3)*k1*k2+k1^2*(q3)+qa*(q1.*r4)*k1^2+q2a*(q2.*r3)*k2^2+q2a*(q3.*r3)*k1^2+q2a*(q3.*r3)*k2^2+2*qa*(q1.*r3)*k2^2+2*qa*(q1.*r3)*k1^2-4*k1*(q1.*r3)*k2+qa*(q1.*r2)*k1^2+q2a*(q3.*r2)*k1^2-8*k1*(q2.*r2)*k2+3*qa*(q3.*r2)*k1^2+4*qa*(q3.*r2)*k2^2+4*k2^2*(q2.*r2)*qa+k2^2*(r3)+2*k1^2*(q1.*r3)+2*k2^2*(q1.*r3)-4*k1*(q3.*r2)*k2*qa+4*k2^2*(q1.*r2)+3*k1^2*(q1.*r2)+4*k2^2*(q2.*r2)+4*k1^2*(q2.*r2)+k1^2*(q1.*r1)+qa*(q4.*r1)*k2^2+qa*(q2.*r1)*k2^2-2*k1*(q1.*r1)*k2+2*qa*(q3.*r1)*k1^2+k2^2*(q1.*r1)+2*k1^2*(q3.*r1)+3*k2^2*(q2.*r1)+4*k1^2*(q2.*r1)+2*k2^2*(q3.*r1));
    
    surf(X,Y,v);
    
	xlabel('X');
	ylabel('Y');
	zlabel('Z');  
	view([30 50]);
	axis([-50 50 -50 50 -0.01 0]);
	pause(0.5);
	%F(j)=getframe;
	j=j+1;
    
    p1=p1*deltapt1;
    r1=r1*deltapt2;
end
%pause;
%movie(F,2);