精细时程积分

    Module PreciseTimeIntegration
    Contains
!**************************************************************
    Subroutine PreciseTIM (M, G, K, R, N, Nt, tao, X, XX, XXX)
!**************************************************************
!    ?????±????·?·¨???ò??
!        ??ê????????ó?§?¨???¤??????2001.7.20
!!**************************************************************
    Real*8     M(N,N), G(N,N), K(N,N), R(N,Nt), X(N,Nt), &
               XX(N,Nt), XXX(N,Nt), H(2*N,2*N), T(2*N,2*N), &
               R0(2*N,Nt), R1(2*N,Nt), B(N,N), C(N,N), tao
    Call CalH (M, G, K, N, Nt, H, B, C)
    Call CalT (H, tao, N, T)
    Call CalR0R1 (R, R0, R1, N, Nt)
    Call CalX (T, H, R, R0, R1, tao, N, Nt, M, G, B, C, X, XX, XXX)
    End Subroutine PreciseTIM
!**************************************************************
    Subroutine CalH (M, G, K, N, Nt, H, B, C)
    Real*8     M(N,N), G(N,N), K(N,N), R(N,Nt), INVM(N,N), &
               A(N,N), B(N,N), C(N,N), D(N,N), H(2*N,2*N)
    INVM=INV(M,N)
    A=-0.5*Matmul(INVM,G)
    B=0.25*Matmul(Matmul(G,INVM),G)-K
    C=-0.5*Matmul(G,INVM)
    D=INVM
    H(1:N,1:N)=A
    H(N+1:2*N,1:N)=B
    H(1:N,N+1:2*N)=D
    H(N+1:2*N,N+1:2*N)=C
    End Subroutine CalH
!**************************************************************
    Subroutine CalT (H, tao, N, T)
    Real*8, Dimension(2*N,2*N) :: H, T, Ta, I
    Real*8 tao, dt
    Integer m
    m=2**20
    dt=tao/m
    I=0; T=0
    Do j=1,2*N
    I(j,j)=1
    End do
    Ta=Matmul(dt*H,(I+0.5*dt*H))
    Do j=1,20
    Ta=2*Ta+Matmul(Ta,Ta)
    End do
    T=I+Ta
    End Subroutine CalT
!**************************************************************
    Subroutine CalR0R1 (R, R0, R1, N, Nt)
    Real*8 R(N,Nt), R0(2*N,Nt), R1(2*N,Nt)
    R0=0; R1=0;
    R0(N+1:2*N,:)=R
    Do i=2,Nt
    R1(N+1:2*N,i-1)=R(:,i)-R(:,i-1)
    End do
    End Subroutine CalR0R1
!**************************************************************
    Subroutine CalX (T, H, R, R0, R1, tao, N, Nt, M, G, B, C, X, XX, XXX)
    Real*8     X(N,Nt), XX(N,Nt), XXX(N,Nt), H(2*N,2*N), &
               T(2*N,2*N), R0(2*N,Nt), R1(2*N,Nt), V(2*N,Nt), &
               p(N,Nt), q(N,Nt), R(N,Nt), B(N,N), C(N,N), M(N,N), &
               G(N,N), tao, INVH(2*N,2*N), index
    
    q=X
    p=Matmul(M,XX)+0.5*Matmul(G,X)
    V(1:N,:)=q
    V(N+1:2*N,:)=p
    index=0
    Do i=1,2*N
    Do j=1,Nt
    IF (ABS(R0(i,j)).GT.1E-8) THEN
    index=1; GOTO 10
    END IF
    End do
    End do
10    If (abs(index).gt.1e-8) then
    INVH=INV(H,2*N)
    Do i=2,Nt
    WRITE (*,'("***********   LOAD STEP: ",I5,"   ***********")') i
     V(:,i)=Matmul(T,(V(:,i-1)+Matmul(INVH,(R0(:,i-1)+ &
            Matmul(INVH,R1(:,i-1))))))-Matmul(INVH &
            ,(R0(:,i-1)+Matmul(INVH,R1(:,i-1))+ &
            tao*R1(:,i-1)))
    End do
    Else
    Do i=2,Nt
    WRITE (*,'("***********   LOAD STEP: ",I5,"   ***********")') i
    V(:,i)=Matmul(T,V(:,i-1))
    End do
    End if

    q=V(1:N,:)
    p=V(N+1:2*N,:)    
    Do i=1,Nt
    X(:,i)=q(:,i)
    XX(:,i)=Matmul(INV(M,N),p(:,i))-0.5*Matmul(Matmul(INV(M,N),G),X(:,i))
    XXX(:,i)=Matmul(Matmul(INV(M,N),B),X(:,i))-0.5*Matmul(Matmul(INV(M,N) &
             ,G),XX(:,i))+Matmul(Matmul(INV(M,N),C),p(:,i))+Matmul(INV(M,N),R(:,i))
    End do
    End Subroutine CalX
!**************************************************************
    Function INV (A, N)
    Use Numerical_Libraries
    REAL*8, DIMENSION(N,N) :: A, INV
    INV=0
    CALL DLINRG (N, A, N, INV, N)
    End Function INV
!**************************************************************
    End Module PreciseTimeIntegration
!##############################################################
       Program Main
    use PreciseTimeIntegration
    Real(8), Dimension(8, 8) :: K,M,C,INVK
    Real(8), Dimension(8, 4) :: R,U,V,A
       Real(8) dt
       dt=1
       R=0
      C=0
    C(7,7)=5
    C(7,8)=-5
    C(8,7)=-5
       C(8,8)=5
    DO i=1,8
     M(i,i)=8
       END DO
     DO i=1,7
       K(i,i+1)=-4
      k(i+1,i)=-4
       END DO