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 endk2 を計算する時には 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