[交流] 【求助】四阶龙格-库塔方法！！！！ 已有3人参与

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 大家看看我  我要用四阶龙格-库塔方法解一个六变量的方程组（六个方程） 如何去解  这里我写了一点了   这个是方程组！求大侠帮忙

回复此楼