euler.f90

在fortran语言中

!     Last change:  CU    1 Oct 2008   11:38 am
!*                                                                  *
!*                                                                  *
!* This program solves a SECOND-ORDER differential equations          *
!* by the Euler algorithm.             *
!*                                      *
!*------------------------------------------------------------------*
!*

  module external_sub_procs
  implicit none
  private
  public  ::  derivs

  contains

    subroutine derivs(x,y,f)
    implicit none
          REAL, DIMENSION(2), INTENT(OUT) :: f
          REAL, INTENT(IN) :: x
          REAL, DIMENSION(2), INTENT(in) :: y
   !
   ! in this example the two first order coupled eqns are
   !    dx/dt = y and dy/dt = - omega**2 x
   ! with omega = 2
   !
   ! the solution is ofcourse: x = constant * sin (omega * t)
   ! or cos(omega * t) ;
   ! with the constants determined by the boundary conditions.
   !
   ! Note that in this program, I have chosen to write
   !               t = x
   !               x = y(1)
   !           dx/dt = y(2)
   ! and so
   !          dy(1)/dt = f1(y(1),y(2),x)   ; f(1) = f1
   !          dy(2)/dt = f2(y(1),y(2),x)   ; f(2) = f2
   ! This is just to make the coding simpler

             f(1) = y(2) + 0.0D0*X
             f(2) = -4.d0*y(1)
    return
    END SUBROUTINE DERIVS
  end module external_sub_procs



   Program euler
   USE external_sub_procs
   Implicit None
       REAL, DIMENSION(2) :: y, f0, y0

       REAL ::   h, x0

       INTEGER :: i , n, NMAX
!
!* Initialize the integration:
!*
        Open (10,FILE='euler.dat',STATUS='UNKNOWN')
      h = 1.0d-2
      x0 = 0.0D0
      NMAX = 100
!   specify the boundary conditions
!   in this example, we put x = 1 and dx/dt = 0
!   this corresponds to the cosine form of the solution
!
!   The y0 variable denotes the position and velocity at the current time
!   step in the iteration
!
      y0(1) = 1.0D0
      y0(2) = 0.0D0
!*
!* The current point is specified as (x0,Y0) and the
!* step size to be used is H.
!  The subroutine DERS evaluates
!* all the derivatives, e.g., all the components of the
!* derivative vector, at (X,Y). The solution is moved
!* forward one step by the specified number to points, N.
!*
    do n = 1,NMAX

    ! up to some end point

       call derivs(x0,y0,f0)

       DO i = 1, 2
         y(i) = y0(i) + h*f0(i)
       END Do
      ! write out the answer at each time point
      ! and compare with the real answer
       WRITE(10, *) x0, y(1), cos(2*x0)
      ! update the time coordinate
      x0 = x0 +  h

      DO i = 1, 2
         y0(i) = y(i)
      END DO
   END do
        Close(10)
      STOP
END program euler