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