RK4 のプログラムについて---追加説明

資料のプログラムが、間違ってはいないが難しいので、分かりやすいと 思われるプログラムにしてみる。まず、素直に書いてみる。 RK4 のサブルーチン自体はこんな感じ。なお、以下の話では微分方程式の 右辺が従属変数だけの関数であるものを考えている。独立変数も入るように 書き直すのは課題の一部ということで。
      subroutine rk4(x,h,df)
      implicit none
      integer N
      parameter (N=2)
      real*8 x(N), h
      real*8 kx1(N), kx2(N), kx3(N), kx4(N), xtmp(N), ftmp(N)
      integer i
      external df
      call df(x,ftmp)
      do i=1,N
         kx1(i) = ftmp(i)*h
         xtmp(i) = x(i)+kx1(i)/2
      enddo
      call df(xtmp,ftmp)
      do i=1,N
         kx2(i) = ftmp(i)*h
         xtmp(i) = x(i)+kx2(i)/2
      enddo
      call df(xtmp,ftmp)
      do i=1,N
         kx3(i) = ftmp(i)*h
         xtmp(i) = x(i)+kx3(i)
      enddo
      call df(xtmp,ftmp)
      do i=1,N
         kx4(i) = ftmp(i)*h
      enddo
      do i=1,N
         x(i) = x(i)+(kx1(i)+2*kx2(i)+2*kx3(i)+kx4(i))/6
      enddo
      end
つまり、 k1, k2, k3, k4 というのを使うので、その名前で素直に変数を宣言 し、式の通りに順番に計算するだけ。一 応念のために導関数を計算するプログラムとメインプログラムの例もつけると、 こんな感じ。
      subroutine calcdf(x,f)
      implicit none
      real*8 x(*),f(*)
      f(1) = x(2)
      f(2) =  -x(1)
      end

      program main
      implicit none
      integer N
      real*8 pi
      parameter (N=2)
      real*8 x(N),t,h,phase
      integer i, nsteps, norbits, print_interval
      external calcdf
      x(1)=1
      x(2)=0
      pi=4.0d0*atan(1d0)
      write(*,*)'PI=', pi
      write(*,*) 'Enter steps,  orbits, print_interval:'
      read(*,*) nsteps,norbits, print_interval
      h = 2*pi/nsteps
      write(*,*) 'X=', x(1), ' V=', x(1), ' h=', h,
     &      ' norbits= ', norbits
      do i=1,nsteps*norbits
         phase= h*(mod(i, nsteps))
        t = h* i
	call rk4(x,h,calcdf)
        if ( (mod(i, print_interval)) == 0) then
           write(6,*) t, x(1), x(2),  x(1)-cos(phase),x(2)+sin(phase)
        endif
      enddo
      end

こちらは説明は省略。

さて、上のRK4 はなにもいけないところはないが、なんとなくセンスがない気 がする。というのは、k1-k4 の計算はほとんど同じであるのに、それを繰り返 して書いているためにプログラムが必要以上に長いからである。もうちょっと なんとかしてみよう。

まず、 計算した導関数に h を掛けるのは必ずするので、これをするサブルー チンを作ってみる。

      subroutine dfbyh(x,h,df,k,n)
      implicit none
      integer n,i
      real*8 x(n),h,k(n)
      external df
      call df(x,k)
      do i=1,n
         k(i)=k(i)*h
      enddo
      end
つまり、元々は rk4 が df を直接呼ぶようになっていたが、それをやめて rk4 は上の dfbyh を呼ぶ。と、これが df*h を計算して k にしまう。これを 使うと rk4 は下のように少し簡単になる。
      subroutine rk4(x,h,df)
      implicit none
      integer N
      parameter (N=2)
      real*8 x(N), h
      real*8 kx1(N), kx2(N), kx3(N), kx4(N), xtmp(N), ftmp(N)
      integer i
      external df
      call dfbyh(x,h,df,kx1,N)
      do i=1,N
         xtmp(i) = x(i)+kx1(i)/2
      enddo
      call dfbyh(xtmp,h,df,kx2,N)
      do i=1,N
         xtmp(i) = x(i)+kx2(i)/2
      enddo
      call dfbyh(xtmp,h,df,kx3,N)
      do i=1,N
         xtmp(i) = x(i)+kx3(i)
      enddo
      call dfbyh(xtmp,h,df,kx4,N)
      do i=1,N
         x(i) = x(i)+(kx1(i)+2*kx2(i)+2*kx3(i)+kx4(i))/6
      enddo
      end
さて、 xtemp の計算も、同じようなことの繰り返しである。これもまとめて サブルーチンにできないだろうか? 2 で割るかどうかとか、最初は x をその ままつかうとかあるが、これは掛ける数も引数にすればよい。つまり、以下の ようなサブルーチンを作ってみる

      subroutine dfbyh2(x,dx,dxcoef,xwork,h,df,k,n)
      implicit none
      integer n,i
      real*8 x(n),dx(n),dxcoef,xwork(n),h,k(n)
      external df
      do i=1,n
         xwork(i)=x(i)+dx(i)*dxcoef
      enddo
      call df(xwork,k)
      do i=1,n
         k(i)=k(i)*h
      enddo
      end
k2 を計算する時には dx にk1 を、dxcoef には 0.5 を入れる。 k1 を計算する時には dxcoef を 0 にする必要がある。さらに、 dx になにか まともな数字がはいったもの(0を掛けるので何がはいっていてもいいが、「数 ではないもの」では困る)を渡すためにこれは x にしておく。とすると rk4 は以下のようになる。
      subroutine rk4(x,h,df)
      implicit none
      integer N
      parameter (N=2)
      real*8 x(N), h
      real*8 kx1(N), kx2(N), kx3(N), kx4(N), xtmp(N), ftmp(N)
      integer i
      external df
      call dfbyh2(x,x,  0d0,  xtmp,h,df,kx1,N)
      call dfbyh2(x,kx1,0.5d0,xtmp,h,df,kx2,N)
      call dfbyh2(x,kx2,0.5d0,xtmp,h,df,kx3,N)
      call dfbyh2(x,kx3,1.0d0,xtmp,h,df,kx4,N)
      do i=1,N
         x(i) = x(i)+(kx1(i)+2*kx2(i)+2*kx3(i)+kx4(i))/6
      enddo
      end

さて、上の rk4 は dfbyh2 をほとんど同じやり方で4回呼んでいる。これも繰 り返しにできないだろうか?ちょっと技を使えばできないことはない。やった 例が以下である。
      subroutine rk4(x,h,df)
      implicit none
      integer N
      parameter (N=2)
      real*8 x(N), h
      real*8 kx(N), xtmp(N), dxtmp(N),kcoef1(4),kcoef2(4)
      data kcoef1/0d0,0.5d0,0.5d0,1d0/
      data kcoef2/1d0,2d0,2d0,1d0/
      integer i,j
      external df
      do i=1,N
         kx(i)=0
         dxtmp(i)=0
      enddo
      do j=1,4
         call dfbyh2(x,kx, kcoef1(j), xtmp,h,df,kx,N)
         do i=1,N
            dxtmp(i)=dxtmp(i)+kx(i)*kcoef2(j)
         enddo
      enddo
      do i=1,N
         x(i) = x(i)+dxtmp(i)/6
      enddo
      end
確かに繰り返しで表現できてはいるが、プログラムはかえって長くなってしま う。ここまではしないほうがいいかもしれない。なお、このプログラムには少 しいいところがある。 kxが 1-4 の4個ではなく 1 つだけになっていることで ある。 kx を使い回すために、 k1+2k2+2k3+k4 を最後に計算するのではなく、 求まった順に dxtmp に積算していく。自由度が2とかなら別になんでもいいが、 100万変数とか 1億変数とかになると余計なメモリをつかわないですむので嬉 しいことがある。

さて、上の形にしてしまうと、 dfbyh2 は1度しか使われていない。それなら、 サブルーチンにするよりも中身をそのまま書いてしまったほうがわかりやすい。 そうしてみたのが以下である。

      subroutine rk4(x,h,df)
      implicit none
      integer N
      parameter (N=2)
      real*8 x(N), h
      real*8 kx(N),  xtmp(N), dxtmp(N),kcoef1(4),kcoef2(4)
      data kcoef1/0d0,0.5d0,0.5d0,1d0/
      data kcoef2/1d0,2d0,2d0,1d0/
      integer i,j
      external df
      do i=1,N
         kx(i)=0
         dxtmp(i)=0
      enddo
      do j=1,4
         do i=1,N
            xtmp(i)=x(i)+kx(i)*kcoef1(j)
         enddo
         call df(xtmp,kx)
         do i=1,N
            kx(i)=kx(i)*h
            dxtmp(i)=dxtmp(i)+kx(i)*kcoef2(j)
         enddo
      enddo
      do i=1,N
         x(i) = x(i)+dxtmp(i)/6
      enddo
      end