06/07/18 15:03:58
>>680
うまく行かないのは、subroutine rk4 に誤りがあるから。
以下のループ部分は、i = 1 について積分し、次にi = 2について積分しているが、
これはおかしい。i=1,2は独立な量ではなく、相互に依存している。
i=1,2を同時に積分して、ちょっとづつ前進させなければならない。
do i = 1, 2
yt(i) = y(i) + hh*dydx(i)
call derive(xh, yt, dyt)
yt(i) = y(i) + hh*dyt(i)
call derive(xh, yt, dym)
yt(i) = y(i) + h*dym(i)
dym(i) = dyt(i) + dym(i)
call derive(x+h, yt, dyt)
yout(i) = y(i) + h6*(dydx(i) + dyt(i) + 2.0d0*dym(i))
end do
こう直せばおｋ
yt = y + hh*dydx
call derive(xh, yt, dyt)
yt = y + hh*dyt
call derive(xh, yt, dym)
yt = y + h*dym
dym = dyt + dym
call derive(x+h, yt, dyt)
yout = y + h6*(dydx + dyt + 2.0d0*dym)
数値誤差は１０＾－１０程度になるが、これは数値微分を使っているのだからまぁおｋ