fortran のプログラムが分かりません！！助けてください。ルンゲクッタ法についてです。

program runge4
!変数と関数の宣言
implicit none
integer i,N
real t,x,v,h,x0,v0,r,tmax
!初期値
x0=0.0
v0=0.0
tmax=10.0
!小区間の幅
h=0.1
N=int(tmax/h)
!開始時刻と初期値
t=0.0
x=x0
v=v0
write(*,*) t,x,v
!繰り返し計算
do i=0,N-1
t=t+h
call rgkt4(t,x,v,h)
write(*,*) t,x,v
end do
stop
end program runge4

!サブルーチン：４次ルンゲ
subroutine rgkt4(t,x,v,h)
implicit none
!サブルーチン内の変数
real t,x,v,h,f1,f2,m,g,r
real s1(2),s2(2),s3(2),s4(2)
!文関数（運動方程式）の設定
f1(t,x,v)=v
f2(t,x,v)=g-r*v*v/m
!定数
r=0.1
g=9.81
m=1.0
!４次ルンゲクッタ法による計算
s1(1)=h*f1(t,x,v)
s1(2)=h*f2(t,x,v)
s2(1)=h*f1(t+h/2.0,x+s1(1)/2.0,v+s1(2)/2.0)
s2(2)=h*f2(t+h/2.0,x+s1(1)/2.0,v+s1(2)/2.0)
s3(1)=h*f1(t+h/2.0,x+s2(1)/2.0,v+s2(2)/2.0)
s3(2)=h*f2(t+h/2.0,x+s2(1)/2.0,v+s2(2)/2.0)
s4(1)=h*f1(t+h,x+s3(1),v+s3(2))
s4(2)=h*f2(t+h,x+s3(1),v+s3(2))
x=x+(s1(1)+2.0*s2(1)+2.0*s3(1)+s4(1))/6.0
v=v+(s1(2)+2.0*s2(2)+2.0*s3(2)+s4(2))/6.0
return
end

このプログラムでｈ（分割幅）を0.0000001以下にすると正確な値が出なくなってしまいました。
考えても原因が分からなかったのでどなたか教えていただけると本当に助かります！
お願いします！！

No.2 ベストアンサー 回答者： にんたまくん 回答日時： 2018/12/15 09:37

No.1です。

補足です。久しぶりにFortranいじってみました。30年ぶり。相当忘れていました。。
（ルンゲクッタなんて、懐かしいな～）

原因は多分、hが非常に小さいと、tが更新してくれないんだと思います。
多分、現状の計算で、tが小さい方は合っていて、tが大きくなると合わなくなるのではないでしょうか。

t=0.0
x=x0
v=v0
の箇所を
t=10.0
x=92.114006
v=9.90454292
で動作させれ見れば、原因が分かると思います。
doループ掛けても、tが更新せずに、10のままになっていると思います。

例、
t=10.1　→　更新してくれる
t=10.0000001 桁数が足りず、計算時に、下位が削られ、t=10になってしまう。

なので、tだけ、倍精度で計算して見れば、正確な値になってくれるのではないでしょうか。
お試しください。頑張ってね。

この回答へのお礼

にんたまくんさんありがとうございます！！
ｈが小さすぎるとｔが更新してくれないということだったのですね！
てっきりｈは小さければ小さいほど良いもんだと勘違いしてました。
とってもスッキリしました。ありがとうございました！

お礼日時：2018/12/18 15:34

No.1 回答者： にんたまくん 回答日時： 2018/12/14 17:24

こんばんは。

分割数が小さすぎると、有効桁数、誤差等の影響で、正確な値が出なくなってきた
のではないでしょうか？？？。
効果がどの程度あるか分かりませんが、変数をrealから倍精度へ変更して見て
比較するのもやり方かと思いますが。
4次ルンゲクッタ法がn次ルンゲクッタ法の中で有効的、刻み幅のことについて、
引用ですが、SIKINOさんという方が
(http://slpr.sakura.ne.jp/qp/runge-kutta-ex/#rk4e …
が解説されていますが、具体的な幅の数値等は記載がありませんでした。
倍精度の話は記載がありましたけど。

ご参考まで。m(_ _)m