C ====================================================================
C  ODE (Ordinary Differential Equation = 常微分方程式) 
C  応用数学演習(I) 夏休み課題 1997
C --------------------------------------------------------------------
C  鳥取大学 (工・応用数理) 大信田 丈志 <ooshida@damp.tottori-u.ac.jp>
C --------------------------------------------------------------------
C  (d/dt)**2 y + OMEGA**2 * y = f(t)
C  初期条件:  y[t=0] = 0,  (d/dt)y[t=0] = 0
C --------------------------------------------------------------------
C  高階の常微分方程式は 1階に直して解くのが普通である。
C  今の場合 v = (d/dt)y を導入し
C
C  (d/dt) [ y, v ]  = [ v,  -OMEGA**2*y + f(t) ]
C
C  とすればよいだろう。
C  (なお v = (1.0/OMEGA)*(d/dt)y とするほうがいい場合もある。)
C  初期条件は [ y, v ](t=0) = [ 0, 0 ] と なる。
C --------------------------------------------------------------------
C  以下は数値計算および制御ルーチンである。
C  ただし 関数 f だけは 別のファイルに書いてある。
C  (差し替えを容易にするため)
C ====================================================================
        subroutine rk2(y0,v0,ynew,vnew,t0,dt)
C --------------------------------------------------------------------
C
C  『t=t0 において (y,v)=(y0,v0) である』という情報を受け取り、
C  これに基づき 2階の Runge-Kutta法を用いて
C  t=t0+dt での (y,v)の値 (ynew,vnew) を求める。
C
C  解くべき微分方程式は
C
C         +-   -+    +-                     -+
C         |  y  |    |    v                  |
C  (d/dt) |     |  = |                       |
C         |  v  |    |  - OMEGA**2*y + f(t)  |
C         +-   -+    +-                     -+
C
C  であり、これを有限のdtに対して近似しなければならない。
C  ひとつの方法は t = t0+(dt/2)での (y,v)の値 (ym,vm) を用いて
C
C  +-    -+    +-    -+        +-                            -+
C  | ynew |    |  y0  |        |    vm                        |
C  |      |  = |      | + dt * |                              |
C  | vnew |    |  v0  |        |  - OMEGA**2*ym + f(t0+dt/2)  |
C  +-    -+    +-    -+        +-                            -+
C
C  のように書くことであろう。 ここで (ym,vm) は
C
C  +-  -+    +-    -+            +-                       -+
C  | ym |    |  y0  |            |    v0                   |
C  |    |  = |      | + (dt/2) * |                         |
C  | vm |    |  v0  |            |  - OMEGA**2*y0 + f(t0)  |
C  +-  -+    +-    -+            +-                       -+
C
C  として求めればよい。
C
C  こうして求めた (ynew,vnew) と 厳密な (y,v)[t=t0+dt] の値との
C  差は dt**3 に比例することが示される。 くわしくは数値計算の
C  教科書を参照されたい。
C --------------------------------------------------------------------
        implicit double precision (a-h,o-z)

        common /par/ omega2,gamma

        ym = y0 + 0.5*dt *    v0
        vm = v0 + 0.5*dt * ( -omega2*y0 + f(t0) )

        ynew = y0 + dt *    vm
        vnew = v0 + dt * ( -omega2*ym + f(t0+0.5*dt) )

        return
        end
C ====================================================================
        double precision function typein(vname)
        character*(*) vname
C --------------------------------------------------------------------
C  変数の名前などを示す vname を端末に表示し、
C  値をキーボード入力するよう求める。 
C  入力された値は 関数 typein の値として返される。 つまり
C    gamma = typein('gamma')
C  のような使い方を想定している。
C --------------------------------------------------------------------
        parameter(MAXTRY=5)
        ic = 0

 1200   ic = ic + 1
        if(ic.gt.MAXTRY)  goto 4000
        write(6,fmt='(1x,a,'' = ?'')') vname
        read(5,*,err=1200,end=4000) typein
        return

C       入力失敗の場合はゼロを返す
 4000   continue
        typein = 0.d0

        return
        end
C ====================================================================
C  以下がメインルーチンである。
C  計算結果は RESULT で示されるファイルに 文字列形式で記録される。
C --------------------------------------------------------------------
        implicit double precision (a-h,o-z)

        character RESULT*(*)
        parameter(OMEGA=1.0)
        parameter(YINIT=0.0, DYINIT=0.0)
        parameter(NOTE=7,RESULT='ode.out')
C --------------------------------------------------------------------
C  変数 omega2 (=OMEGA**2) や gamma は ほとんど定数のようなものなので
C  引数として渡すような面倒なことは したくない。
C  こういうときこそ common 文を用いる。
C --------------------------------------------------------------------
        common /par/ omega2,gamma
C --------------------------------------------------------------------
C  方程式を解くための時間刻みが dt で、出力の間隔が sample で
C  与えられる。 計算は 0 <= t <= tmax の範囲で行なわれる。
C  当然 0 < dt < sample < tmax で あろう。  
C --------------------------------------------------------------------
	call whatf(6)
        omega2 = OMEGA*OMEGA
        gamma  = typein('gamma')
        dt     = typein('dt')
        sample = typein('sample')
        tmax   = typein('tmax')

        print *,'omega = ', OMEGA
        print *,'gamma = ', gamma
        print *,'dt    = ', dt

        if( gamma.lt.0.d0 .or. dt.le.0.d0
     *      .or. sample.lt.(0.999*dt) .or. tmax.le.sample ) then
          write(6,*) '*** illegal data ? ***'
          stop
        endif
C --------------------------------------------------------------------
        open(NOTE,file=RESULT)
	write(6,*) 'result written in : ', RESULT
C --------------------------------------------------------------------
C  以下、2重ループになっている。 
C  外側のループの変数 nrec は『n番目の記録』の意味で、
C  ループ一回につき一度 write文によって結果を記録する。
C  内側のループは nstep回だけ rk2()を呼び出している。
C  なお ここではdt=[一定値]としているが、本格的な数値計算では 計算の
C  誤差を自動的にチェックしてdtの値を調節するような仕掛を組み込むのが
C  普通である。
C  今回のようにdtを固定したプログラムでは、たとえば dt=0.01 の場合と
C  その半分(dt=0.005)の場合の結果を比較し 実質的な差がないことを
C  確認すると良い。 (もし差があったら dtを小さくして計算しなおす。)
C --------------------------------------------------------------------
        t = 0.0
        y = YINIT
        v = DYINIT

        maxrec = int( 0.5 + (tmax/sample) )
        do 1200 nrec=0,maxrec

          if(t.gt.tmax) goto 3000

          write(NOTE,*) t,y
          tnext = (nrec+1)*sample
          nstep = int( 0.5 + (tnext-t)/dt )

          do 1500 nc=1,nstep
            call rk2(y,v,y1,v1,t,dt)
            y = y1
            v = v1
            t = t + dt
 1500     continue

 1200   continue
C --------------------------------------------------------------------
 3000   continue
        close(NOTE)
C --------------------------------------------------------------------
        stop
        end
C ====================================================================
