precisetimeintegration.txt

精细时程积分的FORTRAN源程序,可以精确进行积分运算

!精细时程积分FORTRAN源程序

主程序main.f90，搞了两个例子进行调用，例子是钟院士那篇文献里的，附后的文献里有的。一个程序是reciseTIM.f90，里面包含了精细时程积分的子程序。 


主程序Program Main
CODE:

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
DO i=2,7
k(i,i)=8
END DO
K(1,1)=4
K(8,8)=4
! WRITE (*,*) M
! WRITE (*,*) K
! WRITE (*,*) C
! INVK=INV(K,8)
! STOP

U=0
U(8,1)=10
V=0
A=0
call PreciseTIM (M, C, k, R, 8, 4, dt, U, V, A)
write (*,*) U
! write (*,*) V
! write (*,*) A

end


PreciseTIM.f90

!-----------------------------------------------------------------
    Module PreciseTimeIntegration

    Contains
CODE:

!//////////////////////////////////////////////////////////////
    Subroutine PreciseTIM (M, G, K, R, N, Nt, tao, X, XX, XXX)
    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