home *** CD-ROM | disk | FTP | other *** search
/ The Fred Fish Collection 1.5 / ffcollection-1-5-1992-11.iso / ff_progs / prog_oth / prep_065.lzh / PREP / KTEST.F < prev    next >
Text File  |  1991-08-16  |  3KB  |  113 lines
  1.       implicit real ( a-h, o-z )
  2.       real y(2 )
  3.       external simple
  4.       xi = 10
  5.       xf = 0
  6.       y(1) = exp(xi)
  7.       y(2) = exp(xi)
  8.       n = 2
  9.       del = .1
  10.       accurc = 1.e-8
  11.       imax = 10
  12.       write(*,*) kutta( xi, xf, y, n, del, accurc, imax, simple ),
  13.      *              ' iterations'
  14.       do 12500 i = 1, 2
  15.          write(*,100) i, xf, y(i)-1
  16. 12500 continue
  17.       stop
  18. 100   format( ' error in y(', i2, ') at x=', g12.5, ' is ', g12.5 )
  19.       end
  20.       subroutine simple( x, y, f )
  21.       implicit real ( a-h, o-z )
  22.       real y(2 ), f(2 )
  23.       f(1) = y(2)
  24.       f(2) = 2*y(2) - y(1)
  25.       return
  26.       end
  27.       integer function kutta( xi, xf, y, n, del, accurc, imax, equa )
  28.       implicit real (a-h,o-z)
  29.       real   y(6 ), yi(6 ), yn(6 ),
  30.      *              k1(6 ), k2(6 ), k3(6 ), k4(6 ), k5(6 ),
  31.      *              f(6 ), e(6 ), f1(6 )
  32.       real   test, xi, xf, del, accurc, xn, h
  33.       real   amax1, abs
  34.       logical quit
  35.       if ( n .gt.  6  ) then
  36.          write(*,*) 'KUTTA: too many equations: ', n, '>', 6
  37.          kutta = -1
  38.          return
  39.       end if
  40.       kutta = 0
  41.       nhalves = 0
  42.       xn = xi
  43.       h = del
  44.       quit = .false.
  45.       do 10000 i000 = 1, ( n ), 1
  46.       yn(i000) = y(i000)
  47. 10000 continue
  48.       if (  ( xf .gt.  xi .and.  h .lt.  0 )  .or.   ( xf .lt.  xi .and.
  49.      *  h .gt.  0 )  ) h = -h
  50. 17500 continue
  51.          if ( ( ( (xn+h.ge. xf) .and.  (xf.ge. xi) )  .or.  ( (xn+h.le. 
  52.      *xf) .and.  (xf.le. xi) )  )  ) then
  53.             del = h
  54.             h = xf - xn
  55.             quit = .true.
  56.          end if
  57.          call equa(xn, yn, f1)
  58. 17501 continue
  59.       do 10001 i000 = 1, ( n ), 1
  60.             k1(i000) = h*f1(i000)/3.
  61.             yi(i000) = yn(i000) + k1(i000)
  62. 10001 continue
  63.             call equa(xn + h/3., yi, f)
  64.       do 10002 i000 = 1, ( n ), 1
  65.             k2(i000) = h*f(i000)/3.
  66.             yi(i000) = yn(i000) + k1(i000)/2. + k2(i000)/2.
  67. 10002 continue
  68.             call equa(xn + h/3., yi, f)
  69.       do 10003 i000 = 1, ( n ), 1
  70.             k3(i000) = h*f(i000)/3.
  71.             yi(i000) = yn(i000) + 3.*k1(i000)/8. + 9.*k3(i000)/8.
  72. 10003 continue
  73.             call equa(xn + h/2., yi, f)
  74.       do 10004 i000 = 1, ( n ), 1
  75.             k4(i000) = h*f(i000)/3.
  76.             yi(i000) = yn(i000) + 3.*k1(i000)/2. - 9.*k3(i000)/2. + 6.*k
  77.      *4(i000)
  78. 10004 continue
  79.             call equa(xn + h, yi, f)
  80.             test = 0.0
  81.       do 10005 i000 = 1, ( n ), 1
  82.             k5(i000) = h*f(i000)/3.
  83.             e(i000) = (k1(i000) - 9.*k3(i000)/2. + 4.*k4(i000) - k5(i000
  84.      *)/2.)/5.
  85.             test = amax1(test, abs(e(i000)))
  86. 10005 continue
  87.       if (.not.( test .ge.  accurc )) goto 15001
  88.             nhalves = nhalves + 1
  89.             if( nhalves .ge.  imax ) then
  90.                del = h
  91.                kutta = -1
  92.                write(*,*) 'KUTTA: step made too small:'
  93.                write(*,*) '       del = ', del, '    at x = ', xn
  94.                return
  95.             end if
  96.             h = h/2.
  97.             quit = .false.
  98.       goto 17501
  99. 15001 continue
  100.       do 10006 i000 = 1, ( n ), 1
  101.          yn(i000) = yn(i000) + (k1(i000) + 4.*k4(i000) + k5(i000))/2.
  102. 10006 continue
  103.          xn = xn + h
  104.          kutta = kutta + 1
  105.          if ( test .lt.  accurc/32. ) h = 2.*h
  106.       if (.not.( quit )) goto 17500
  107. 15000 continue
  108.       do 10007 i000 = 1, ( n ), 1
  109.       y(i000) = yn(i000)
  110. 10007 continue
  111.       return
  112.       end
  113.