您所在位置：主题：我采用4阶的龙格库塔解方程组，计算结果绘制的图有很多小振动，可能是哪儿出问题了？？
function f1(t,x,a,b,c,d)
implicit none
real*8::v1,vb,x0,x,t
complex*16::f1,a,b,c,d
real*8::g1,g2,g,xh,pi
common/group/g1,g2,g,xh,pi
pi=3.141592653d0
vb=10.0d0
v1=x**2+vb*(cos(pi*x/5.2d0))**2
f1=(0.0d0,1.0d0)*((a+b)/(2.0d0*xh**2)-(1.0d0/xh**2+v1+g1/pi*c*conjg(c)+g/pi*d*conjg(d))*c)
end function
function f2(t,x,a,b,c,d)
implicit none
real*8::g1,g2,g,xh,pi
common/group/g1,g2,g,xh,pi
real*8::v2,x,t
complex*16::f2,a,b,c,d
pi=3.141592653d0
v2=2.0d0*x**2-3.0d0
f2=(0.0d0,1.0d0)*((a+b)/(2.0d0*xh**2)-(1.0d0/xh**2+v2+g2/pi*c*conjg(c)+g/pi*d*conjg(d))*c)
end function
program main
implicit none
integer,parameter::m=500
complex*16::y1(-m-1:m+1),y2(-m-1:m+1)
complex*16::K1(-m-1:m+1),L1(-m-1:m+1),M1(-m-1:m+1),N1(-m-1:m+1)
complex*16::K2(-m-1:m+1),L2(-m-1:m+1),M2(-m-1:m+1),N2(-m-1:m+1)
complex*16::a,b,c,d
real*8::x,t,th
complex*16,external::f1,f2
real*8::limitt
integer::l,j
real*8::p(-m-3:m+3),im(-m-3:m+3),re(-m-3:m+3)
real*8::ilm,rle,irm,rre,NL,NR,N,cta1
real*8::g1,g2,g,xh,pi
common/group/g1,g2,g,xh,pi
xh=0.01d0
th=0.0001d0
g1=0.1d0
g2=0.0d0
g=0.0d0
limitt=100.0d0
pi=3.141592653d0
l=0
open(unit=1,file='1.txt')
open(unit=2,file='2.txt')
do j=-m-1,m+1,1
y1(j)=(exp(-(j*xh+2.0d0)**2))*(0.55d0*pi)**(-0.25)+(exp(-(j*xh-2.0d0)**2))*(220.0d0*pi)**(-0.25)
y2(j)=exp(-0.5d0*(j*xh)**2)*pi**(-0.25)
enddo
do t=0.0d0,limitt,th
do j=-m,m,1
x=j*xh
K1(j)=f1(t,j*xh,y1(j+1),y1(j-1),y1(j),y2(j))
K2(j)=f2(t,j*xh,y2(j+1),y2(j-1),y2(j),y1(j))
enddo
do j=-m,m,1
x=j*xh
L1(j)=f1(t+0.5d0*th,j*xh,y1(j+1)+0.5d0*th*K1(j+1),y1(j-1)+0.5d0*th*K1(j-1),y1(j)+0.5d0*th*K1(j),y2(j)+0.5d0*th*K2(j))
L2(j)=f2(t+0.5d0*th,j*xh,y2(j+1)+0.5d0*th*K2(j+1),y2(j-1)+0.5d0*th*K2(j-1),y2(j)+0.5d0*th*K2(j),y1(j)+0.5d0*th*K1(j))
enddo
do j=-m,m,1
x=j*xh
M1(j)=f1(t+0.5d0*th,j*xh,y1(j+1)+L1(j+1)*0.5d0*th,y1(j-1)+L1(j-1)*0.5d0*th,y1(j)+L1(j)*0.5d0*th,y2(j)+L2(j)*0.5d0*th)
M2(j)=f2(t+0.5d0*th,j*xh,y2(j+1)+L2(j+1)*0.5d0*th,y2(j-1)+L2(j-1)*0.5d0*th,y2(j)+L2(j)*0.5d0*th,y1(j)+L1(j)*0.5d0*th)
enddo
do j=-m,m,1
x=j*xh
N1(j)=f1(t+th,j*xh,y1(j+1)+th*M1(j+1),y1(j-1)+th*M1(j-1),y1(j)+th*M1(j),y2(j)+th*M2(j))
N2(j)=f2(t+th,j*xh,y2(j+1)+th*M2(j+1),y2(j-1)+th*M2(j-1),y2(j)+th*M2(j),y1(j)+th*M1(j))
enddo
do j=-m,m,1
x=j*xh
y1(j)=y1(j)+th/6.0d0*(K1(j)+2.0d0*L1(j)+2.0d0*M1(j)+N1(j))
y2(j)=y2(j)+th/6.0d0*(K2(j)+2.0d0*L2(j)+2.0d0*M2(j)+N2(j))
p(j)=y1(j)*conjg(y1(j))
im(j)=p(j)*aimag(y1(j))
re(j)=p(j)*real(y1(j))
enddo
do j=-m,-1,1
ilm=ilm+0.5d0*xh*(im(j)+im(j+1))
rle=rle+0.5d0*xh*(re(j)+re(j+1))
NL=NL+0.5d0*xh*(p(j)+p(j+1))
enddo
do j=0,m,1
irm=irm+0.5d0*xh*(im(j)+im(j+1))
rre=rre+0.5d0*xh*(re(j)+re(j+1))
NR=NR+0.5d0*xh*(p(j)+p(j+1))
end do
N=NL-NR
cta1=(datan2(ilm,rle)-datan2(irm,rre))/pi
if (cta1>=0.9d0)then
cta1=cta1-2.0d0
endif
if (cta1<=-0.9d0)then
cta1=cta1+2.0d0
endif
if (mod(l,100)==0) then
write(1,*)t,cta1
write(2,*)t,N
endif
l=l+1
ilm=0.0d0
rle=0.0d0
irm=0.0d0
rre=0.0d0
NL=0.0d0
NR=0.0d0
end do
stop
end
implicit none
real*8::v1,vb,x0,x,t
complex*16::f1,a,b,c,d
real*8::g1,g2,g,xh,pi
common/group/g1,g2,g,xh,pi
pi=3.141592653d0
vb=10.0d0
v1=x**2+vb*(cos(pi*x/5.2d0))**2
f1=(0.0d0,1.0d0)*((a+b)/(2.0d0*xh**2)-(1.0d0/xh**2+v1+g1/pi*c*conjg(c)+g/pi*d*conjg(d))*c)
end function
function f2(t,x,a,b,c,d)
implicit none
real*8::g1,g2,g,xh,pi
common/group/g1,g2,g,xh,pi
real*8::v2,x,t
complex*16::f2,a,b,c,d
pi=3.141592653d0
v2=2.0d0*x**2-3.0d0
f2=(0.0d0,1.0d0)*((a+b)/(2.0d0*xh**2)-(1.0d0/xh**2+v2+g2/pi*c*conjg(c)+g/pi*d*conjg(d))*c)
end function
program main
implicit none
integer,parameter::m=500
complex*16::y1(-m-1:m+1),y2(-m-1:m+1)
complex*16::K1(-m-1:m+1),L1(-m-1:m+1),M1(-m-1:m+1),N1(-m-1:m+1)
complex*16::K2(-m-1:m+1),L2(-m-1:m+1),M2(-m-1:m+1),N2(-m-1:m+1)
complex*16::a,b,c,d
real*8::x,t,th
complex*16,external::f1,f2
real*8::limitt
integer::l,j
real*8::p(-m-3:m+3),im(-m-3:m+3),re(-m-3:m+3)
real*8::ilm,rle,irm,rre,NL,NR,N,cta1
real*8::g1,g2,g,xh,pi
common/group/g1,g2,g,xh,pi
xh=0.01d0
th=0.0001d0
g1=0.1d0
g2=0.0d0
g=0.0d0
limitt=100.0d0
pi=3.141592653d0
l=0
open(unit=1,file='1.txt')
open(unit=2,file='2.txt')
do j=-m-1,m+1,1
y1(j)=(exp(-(j*xh+2.0d0)**2))*(0.55d0*pi)**(-0.25)+(exp(-(j*xh-2.0d0)**2))*(220.0d0*pi)**(-0.25)
y2(j)=exp(-0.5d0*(j*xh)**2)*pi**(-0.25)
enddo
do t=0.0d0,limitt,th
do j=-m,m,1
x=j*xh
K1(j)=f1(t,j*xh,y1(j+1),y1(j-1),y1(j),y2(j))
K2(j)=f2(t,j*xh,y2(j+1),y2(j-1),y2(j),y1(j))
enddo
do j=-m,m,1
x=j*xh
L1(j)=f1(t+0.5d0*th,j*xh,y1(j+1)+0.5d0*th*K1(j+1),y1(j-1)+0.5d0*th*K1(j-1),y1(j)+0.5d0*th*K1(j),y2(j)+0.5d0*th*K2(j))
L2(j)=f2(t+0.5d0*th,j*xh,y2(j+1)+0.5d0*th*K2(j+1),y2(j-1)+0.5d0*th*K2(j-1),y2(j)+0.5d0*th*K2(j),y1(j)+0.5d0*th*K1(j))
enddo
do j=-m,m,1
x=j*xh
M1(j)=f1(t+0.5d0*th,j*xh,y1(j+1)+L1(j+1)*0.5d0*th,y1(j-1)+L1(j-1)*0.5d0*th,y1(j)+L1(j)*0.5d0*th,y2(j)+L2(j)*0.5d0*th)
M2(j)=f2(t+0.5d0*th,j*xh,y2(j+1)+L2(j+1)*0.5d0*th,y2(j-1)+L2(j-1)*0.5d0*th,y2(j)+L2(j)*0.5d0*th,y1(j)+L1(j)*0.5d0*th)
enddo
do j=-m,m,1
x=j*xh
N1(j)=f1(t+th,j*xh,y1(j+1)+th*M1(j+1),y1(j-1)+th*M1(j-1),y1(j)+th*M1(j),y2(j)+th*M2(j))
N2(j)=f2(t+th,j*xh,y2(j+1)+th*M2(j+1),y2(j-1)+th*M2(j-1),y2(j)+th*M2(j),y1(j)+th*M1(j))
enddo
do j=-m,m,1
x=j*xh
y1(j)=y1(j)+th/6.0d0*(K1(j)+2.0d0*L1(j)+2.0d0*M1(j)+N1(j))
y2(j)=y2(j)+th/6.0d0*(K2(j)+2.0d0*L2(j)+2.0d0*M2(j)+N2(j))
p(j)=y1(j)*conjg(y1(j))
im(j)=p(j)*aimag(y1(j))
re(j)=p(j)*real(y1(j))
enddo
do j=-m,-1,1
ilm=ilm+0.5d0*xh*(im(j)+im(j+1))
rle=rle+0.5d0*xh*(re(j)+re(j+1))
NL=NL+0.5d0*xh*(p(j)+p(j+1))
enddo
do j=0,m,1
irm=irm+0.5d0*xh*(im(j)+im(j+1))
rre=rre+0.5d0*xh*(re(j)+re(j+1))
NR=NR+0.5d0*xh*(p(j)+p(j+1))
end do
N=NL-NR
cta1=(datan2(ilm,rle)-datan2(irm,rre))/pi
if (cta1>=0.9d0)then
cta1=cta1-2.0d0
endif
if (cta1<=-0.9d0)then
cta1=cta1+2.0d0
endif
if (mod(l,100)==0) then
write(1,*)t,cta1
write(2,*)t,N
endif
l=l+1
ilm=0.0d0
rle=0.0d0
irm=0.0d0
rre=0.0d0
NL=0.0d0
NR=0.0d0
end do
stop
end