![](http://oshiete.xgoo.jp/images/v2/pc/qa/question_title.png?e8efa67)
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
!サブルーチン:4次ルンゲ
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
!4次ルンゲクッタ法による計算
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
このプログラムでh(分割幅)を0.0000001以下にすると正確な値が出なくなってしまいました。
考えても原因が分からなかったのでどなたか教えていただけると本当に助かります!
お願いします!!
No.2ベストアンサー
- 回答日時:
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だけ、倍精度で計算して見れば、正確な値になってくれるのではないでしょうか。
お試しください。頑張ってね。
にんたまくんさんありがとうございます!!
hが小さすぎるとtが更新してくれないということだったのですね!
てっきりhは小さければ小さいほど良いもんだと勘違いしてました。
とってもスッキリしました。ありがとうございました!
No.1
- 回答日時:
こんばんは。
分割数が小さすぎると、有効桁数、誤差等の影響で、正確な値が出なくなってきた
のではないでしょうか???。
効果がどの程度あるか分かりませんが、変数をrealから倍精度へ変更して見て
比較するのもやり方かと思いますが。
4次ルンゲクッタ法がn次ルンゲクッタ法の中で有効的、刻み幅のことについて、
引用ですが、SIKINOさんという方が
(http://slpr.sakura.ne.jp/qp/runge-kutta-ex/#rk4e …
が解説されていますが、具体的な幅の数値等は記載がありませんでした。
倍精度の話は記載がありましたけど。
ご参考まで。m(_ _)m
お探しのQ&Aが見つからない時は、教えて!gooで質問しましょう!
似たような質問が見つかりました
- C言語・C++・C# C言語で再起関数とポインタを用いて文字列反転をする方法がわかりません。 4 2023/04/29 20:32
- Visual Basic(VBA) エクセルVBAで教えて頂きたいのですが? 2 2022/12/31 20:28
- その他(プログラミング・Web制作) 物理の斜方投射のシミュレーションにおける位置や速度の単位について 4 2023/05/31 09:50
- 統計学 t値の計算方法 1 2022/11/29 18:37
- 化学 化学の溶解度積Kspの計算に関しての質問です。MgF2の溶解度積を求める問題で Ksp=[Mg^2+ 1 2022/08/13 18:07
- 数学 写真について質問なのですが、 ①の図の面積Sを求めるとき、②と③の図の面積、つまりS=S2+S3で求 4 2023/04/27 17:20
- C言語・C++・C# c言語でユーザ関数を利用して入力された文字列を反転させるプログラムを作りたいです。 3 2023/01/29 19:47
- 数学 高校時代電離平衡の計算に関しての質問です。 問題集で、 酢酸は水溶液中で一部が電離し、次のような電離 2 2022/10/22 18:59
- ドライブ・ストレージ PCのSSD換装 2 2022/06/11 08:26
- C言語・C++・C# C 言語の Gauss Jordan 法について 2 2022/12/28 11:16
関連するカテゴリからQ&Aを探す
おすすめ情報
デイリーランキングこのカテゴリの人気デイリーQ&Aランキング
-
1、Rstudioで回帰直線を求める...
-
ビーリアルのユーザー名を変え...
-
教えてください
-
再起動後必ず2つのエラーが出...
-
英数字を含む文字列(0-9,A-Z)...
-
C言語の入力した文字を反転させ...
-
プログラミングの課題で1万円か...
-
Ruby on railsをrails sで立ち...
-
実行時エラー450:引数の数が一...
-
(再質問)エクセルのマクロボ...
-
エクセルvbaでチェックボックス...
-
100万件越えCSVから条件を満た...
-
これらは書誌情報だと思うので...
-
情報の表現。()内がどうしても...
-
pythonのerrorコード
-
三項でたとえば交換って
-
WinSCPで画像のように puttyを...
-
パソコンのスクリーンセーバー...
-
こういう問題分をよんだとき
-
バーチャルボックスが使えなく...
マンスリーランキングこのカテゴリの人気マンスリーQ&Aランキング
-
バーコードのチェックデジット...
-
rubyのエラー out of float ra...
-
ruby の Σ計算
-
fortran のプログラムが分かり...
-
【訂正】Ruby(Perlでも)で60進...
-
小文字wと大文字Wの区別
-
システムエンジニアの適正について
-
web上のhtmlファイルから文字デ...
-
VB.NETで階乗を求めるプログラ...
-
COBOLのIFの入れ子について
-
Ruby / passenger のインストー...
-
Ruby interpreter (CUI) 2.2.3 ...
-
rubyでバイナリファイルを直接...
-
Passengerがインストールできな...
-
RubyでNo such file or directo...
-
プログラミング言語で大文字と...
-
式?文?節?
-
CかC++どちらを覚えるべきですか?
-
Rubyで画像処理
-
Rudyを覚えたい
おすすめ情報