戻る
ルンゲクッタ法のサンプル
ルンゲクッタ法(4次) とオイラー法との比較をしてみる計算です。
(2009.6)
プログラムのダウンロード
解いている問題は、dy/dx = f(x,y) で初期値 x0 y0 になっています。
具体的には、x0=0; y0=2; x_end=2; f(x,y) = 2*x*(1-y*y) で解いています。
答えは、(3*exp(2*x*x)+1)/(3*exp(2*x*x)-1) のはずです。
比較してみると、既に１６分割くらいで良い答えになっていますが、
その後、Runge-Kutta で 128分割、Euler 法で 8192 あたりで最高の精度を
示し、それ以上分割すると、数値誤差のために却って精度が落ちて行きます。
どんどん分割数を増やせば、精度がどんどん上るということを、実際の
プログラム上で実装するには、数値誤差に対する考察が必要なことが良く分かります。
gfortran43 rungekutta.f90
./a.out
# n Euler Runge-Kutta Answer
2 -4.0000000 -0.56600535 1.0002236
4 0.40625000 1.1356248 1.0002236
8 1.0000620 1.0007145 1.0002236
16 1.0000006 1.0002341 1.0002236
32 1.0000477 1.0002241 1.0002236
64 1.0001163 1.0002238 1.0002236
128 1.0001647 1.0002236 1.0002236
256 1.0001930 1.0002238 1.0002236
512 1.0002081 1.0002236 1.0002236
1024 1.0002156 1.0002239 1.0002236
2048 1.0002197 1.0002235 1.0002236
4096 1.0002218 1.0002239 1.0002236
8192 1.0002229 1.0002241 1.0002236
16384 1.0002277 1.0002282 1.0002236
32768 1.0002435 1.0002439 1.0002236
65536 1.0002440 1.0002440 1.0002236
131072 1.0004882 1.0004882 1.0002236
262144 1.0009761 1.0009761 1.0002236
524288 1.0019512 1.0019512 1.0002236
1048576 1.0038986 1.0038986 1.0002236
プログラム上の注意点は、関数副プログラムを
ちゃんと external 指定してやって、関数を引数にうける
サブルーチンに渡していることです。call する側、受ける側双方で
external 指定が必要なようです。
program main
implicit none
real :: x0, y0, x_end, y_end_r, y_end_e
real, external :: func, answer
integer :: i, n
x0=0; y0=2; x_end=2
write(*,*) '# n Euler Runge-Kutta Answer'
n=1
do i=1, 20
n = n*2
call euler(x0, y0, n, x_end, func, y_end_e)
call rungekutta(x0, y0, n, x_end, func, y_end_r)
write(*,*) n, y_end_e, y_end_r, answer(x_end)
end do
end
subroutine euler( x0, y0, n, x_end, f, y_end)
implicit none
real :: x0, y0, x_end, y_end
real, external :: f
real :: k1, k2, k3, k4
real :: x, y, dx
integer :: i, n
dx = (x_end - x0)/n
x=x0
y=y0
do i=1, n
y = y + dx*f(x,y)
x = x + dx
! write(*,*) 'Euler = ',i, x, y
end do
y_end=y
! write(*,*) 'Euler ', x0, dx, n, x_end, y_end
end subroutine
subroutine rungekutta( x0, y0, n, x_end, f, y_end)
implicit none
real :: x0, y0, x_end, y_end
real, external :: f
real :: k1, k2, k3, k4
real :: x, y, dx
integer :: i, n
dx = (x_end - x0)/n
x=x0
y=y0
do i=1, n
k1 = dx * f(x, y)
k2 = dx * f( x+dx/2, y+k1/2 )
k3 = dx * f( x+dx/2, y+k2/2 )
k4 = dx * f( x+dx , y+k3 )
y = y + (1.0/6.0)*(k1+2*k2+2*k3+k4)
x = x + dx
! write(*,*) x, y
end do
y_end=y
! write(*,*) 'Rungekutta ', x0, dx, n, x_end, y_end
end subroutine
real function func(x, y)
real :: x, y
func = 2*x*(1-y*y)
end function
real function answer( x )
real :: x
answer = (3*exp(2*x*x)+1)/(3*exp(2*x*x)-1)
end function
玉川一郎