double precision qx,qy,qz,qx0,qy0,qz0,qx1,qy1,qz1,a,b,c,h,t,d
double precision qkx1,qkx2,qkx3,qkx4,qky1,qky2,qky3,qky4,qkz1,qkz2
&,qkz3,qkz4
double precision px,py,pz,px0,py0,pz0,px1,py1,pz1
double precision pkx1,pkx2,pkx3,pkx4,pky1,pky2,pky3,pky4,pkz1,pkz2
&,pkz3,pkz4
double precision k1,k2 !耦合强度
double precision delta, eps ! 误差精度
a=0.66d0
b=0.201d0
c=0.165d0
d=1d0/0.407d0
h=1.0d-3
px0=0.3d0
py0=0.5d0
pz0=1.0d0
qx0=0.8d0
qy0=1.5d0
qz0=0.2d0
eps=1d-4
open(10,file='delta.dat')
k1=0.4
k2=0.2
! do k1=0.001d0,1.0d0,0.05d0
! do k2=0.001d0,1.0d0,0.05d0
do 30 t=h,2.0d2,h
qkx1=h*qx(qy0,qz0,a)
qky1=h*qy(qx0,qy0,qz0,py0,b,c,k1,k2)
qkz1=h*qz(qx0,qz0,d)
pkx1=h*px(py0,pz0,a)
pky1=h*py(px0,py0,pz0,qy0,b,c,k1,k2)
pkz1=h*pz(px0,pz0,d)
qx1=qx0+qkx1/2.0d0
qy1=qy0+qky1/2.0d0
qz1=qz0+qkz1/2.0d0
px1=px0+pkx1/2.0d0
py1=py0+pky1/2.0d0
pz1=pz0+pkz1/2.0d0
qkx2=h*qx(qy1,qz1,a)
qky2=h*qy(qx1,qy1,qz1,py1,b,c,k1,k2)
qkz2=h*qz(qx1,qz1,d)
pkx2=h*px(py1,pz1,a)
pky2=h*py(px1,py1,pz1,qy1,b,c,k1,k2)
pkz2=h*pz(px1,pz1,d)
qx1=qx0+qkx2/2.0d0
qy1=qy0+qky2/2.0d0
qz1=qz0+qkz2/2.0d0
px1=px0+pkx2/2.0d0
py1=py0+pky2/2.0d0
pz1=pz0+pkz2/2.0d0
qkx3=h*qx(qy1,qz1,a)
qky3=h*qy(qx1,qy1,qz1,py1,b,c,k1,k2)
qkz3=h*qz(qx1,qz1,d)
pkx3=h*px(py1,pz1,a)
pky3=h*py(px1,py1,pz1,qy1,b,c,k1,k2)
pkz3=h*pz(px1,pz1,d)
qx1=qx0+qkx3
qy1=qy0+qky3
qz1=qz0+qkz3
px1=px0+pkx3
py1=py0+pky3
pz1=pz0+pkz3
qkx4=h*qx(qy1,qz1,a)
qky4=h*qy(qx1,qy1,qz1,py1,b,c,k1,k2)
qkz4=h*fz(qx1,qz1,d)
pkx4=h*px(py1,pz1,a)
pky4=h*py(px1,py1,pz1,qy1,b,c,k1,k2)
pkz4=h*pz(px1,pz1,d)
qx1=qx0+(qkx1+2.0*qkx2+2.0*qkx3+qkx4)/6.0d0
qy1=qy0+(qky1+2.0*qky2+2.0*qky3+qky4)/6.0d0
qz1=qz0+(qkz1+2.0*qkz2+2.0*qkz3+qkz4)/6.0d0
px1=px0+(pkx1+2.0*pkx2+2.0*pkx3+pkx4)/6.0d0
py1=py0+(pky1+2.0*pky2+2.0*pky3+pky4)/6.0d0
pz1=pz0+(pkz1+2.0*pkz2+2.0*pkz3+pkz4)/6.0d0
qx0=qx1
qy0=qy1
qz0=qz1
px0=px1
py0=py1
pz0=pz1
大家看看我 我要用四阶龙格-库塔方法解一个六变量的方程组(六个方程) 如何去解 这里我写了一点了 
这个是方程组!求大侠帮忙 |