fouru1111.m

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

clear all;

clear all;

k1=1;
k2=1;
k3=1;
k4=1;


A12=-(k1-k2)^2/(k1*k2);
A13=-(k1-k3)^2/(k1*k3);
A14=-(k1-k4)^2/(k1*k4);
A23=-(k2-k3)^2/(k2*k3);
A24=-(k2-k4)^2/(k2*k4);
A34=-(k3-k4)^2/(k3*k4);

w1=k1^2/2;
w2=k2^2/2;
w3=k3^2/2;
w4=k4^2/2;

p1=k1^3/4;
p2=k2^3/4;
p3=k3^3/4;
p4=k4^3/4;

z10=-40.0;	%z1初值
z20=-20.0;	%z2初值
z30=20.0;	%z1初值
z40=40.0;	%z2初值

ST=-2.0; 	%开始时间
ET=2.0;		%终止时间

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

y=x;



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

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

z1_1=k1*X+w1*Y+z10;
z2_1=k2*X+w2*Y+z20;
z3_1=k3*X+w3*Y+z30;
z4_1=k4*X+w4*Y+z40;

one=ones(size(z1_1));

q_12a=A12*one;
q2_12a=q_12a.*q_12a;

q_13a=A13*one;
q2_13a=q_13a.*q_13a;

q_14a=A14*one;
q2_14a=q_14a.*q_14a;

q_23a=A23*one;
q2_23a=q_23a.*q_23a;

q_24a=A24*one;
q2_24a=q_24a.*q_24a;

q_34a=A34*one;
q2_34a=q_34a.*q_34a;

for t=ST:0.2:ET
	z1=z1_1+p1*t;
	z2=z2_1+p2*t;
	z3=z3_1+p3*t;
	z4=z4_1+p4*t;
    
	q1=exp(z1);
	q2=q1.*q1;
    
    r1=exp(z2);
	r2=r1.*r1;

    s1=exp(z3);
	s2=s1.*s1;
    
    t1=exp(z4);
	t2=t1.*t1;
    
    u=(1+(q1)+(r1)+(s1)+(t1));
    u=1./u;
    
    surf(X,Y,u);
    
	xlabel('X');
	ylabel('Y');
	zlabel('Z');
    
	%view([30 50]);
	%axis([-40 40 -40 40 -1 1]);
	pause(0.5);
	%F(j)=getframe;
	j=j+1;
end
%pause;
%movie(F,2);