fouru1234.m

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

clear all;

clear all;

k1=1;
k2=2;
k3=3;
k4=4;


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;

ST=0.0; 	%开始时间
ET=20.0;		%终止时间

SX=0.0;	%x开始
EX=20.0;	%x终止
x=SX:1:EX;	%x序列

SY=0.0;	%y开始
EY=20.0;	%y终止
y=SY:1:EY;	%y序列

z10=0.0;	%z1初值
z20=0.0;	%z2初值
z30=0.0;	%z1初值
z40=0.0;	%z2初值

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.5: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=((q1)+(r1)+(s1)+(t1)+(q1.*r1.*q_12a)+(q1.*s1.*q_13a)+(q1.*t1.*q_14a)+(r1.*s1.*q_23a)+(r1.*t1.*q_24a)+(s1.*t1.*q_34a)).*(1+(q1)+(r1)+(s1)+(t1));
    u=1./u;
    
    %u =-1/2*u.*(5*(q1.*t1.*q_14a)+3*(q1.*t1.*q_14a.*r1)+4*(q2.*t1.*q_14a)+(q1.*t2.*q_14a)+2*(q1.*t1.*q_14a.*s1)+4*(q1.*s1.*q_13a)+3*(q2.*s1.*q_13a)+(q1.*s2.*q_13a)+2*(q1.*s1.*q_13a.*r1)+(q1)+3*(s1)+2*(r1)+4*(t1)+5*(r1.*s1.*q_23a)+3*(r2.*s1.*q_23a)+(r1.*s1.*q_23a.*t1)+2*(r1.*s2.*q_23a)+4*(r1.*s1.*q_23a.*q1)+2*(r1.*t2.*q_24a)+3*(r1.*t1.*q_24a.*s1)+6*(r1.*t1.*q_24a)+4*(r2.*t1.*q_24a)+5*(r1.*t1.*q_24a.*q1)+7*(s1.*t1.*q_34a)+6*(s1.*t1.*q_34a.*q1)+5*(s1.*t1.*q_34a.*r1)+4*(s2.*t1.*q_34a)+3*(s1.*t2.*q_34a)+3*(q1.*r1.*q_12a)-(q1.*r1.*q_12a.*t1)+(q1.*r2.*q_12a)+2*(q2.*r1.*q_12a));
        
	surf(X,Y,u);
	xlabel('X');
	ylabel('Y');
	zlabel('Z');
    
	%view([30 50]);
	%axis([-40 40 -40 40 -50 50]);
	pause(0.5);
	%F(j)=getframe;
	j=j+1;
end
%pause;
%movie(F,2);