古典的4次Runge-Kutta法汎用ソルバー

　　古典的4次Runge-Kutta法汎用ソルバー クスクス 2006/12/13 02:06:08 
　　├RK4ソルバサンプル説明書き。 クスクス 2006/12/13 02:08:37 
　　│└RK4ソルバプログラム本体その1。 クスクス 2006/12/13 02:11:18 
　　│　└RK4ソルバプログラム本体その2。 クスクス 2006/12/13 02:12:25 
　　│　　├有難うございます。 SECOND 2006/12/13 05:24:55 
　　│　　└!これだけに、なりました。これと照合すると... SECOND 2006/12/16 00:32:26 
　　└RK4ソルバ動作確認(等速円運動で検査)。 クスクス 2006/12/13 15:26:27 
　　　└等速円運動シミレーションその１． クスクス 2006/12/13 15:28:20 
　　　　└その２。 クスクス 2006/12/13 15:29:37 
　　　　　└追加訂正分。 クスクス 2006/12/13 17:42:59 
　　　　　　└これで十分です。詳細に追跡します。有難う... SECOND 2006/12/14 02:34:16

　　古典的4次Runge-Kutta法汎用ソルバー クスクス 2006/12/13 02:06:08  ツリーへ

古典的4次Runge-Kutta法汎用ソルバー 返事を書く

クスクス 2006/12/13 02:06:08

4次のルンゲ・クッタ法プログラムは、Fortran,C,Java のコード
ならば、ネットのいたるところで手に入れることができますが、
FullBASIC のコードは探しても見当たりません。これは、RK4 の
原理の理解を目的とした、ｎ元連立１階常微分方程式対応のRK4
ソルバです。使用法を説明するためのサンプルは、単振子のシミ
レーションです。このソルバは
「パソコンで見る天体の動き」長沢工著　地人書刊
のN88BASICで書かれたコードCLARK.BASの十進BASICへの移植です。

　　├RK4ソルバサンプル説明書き。 クスクス 2006/12/13 02:08:37  ツリーへ

Re: 古典的4次Runge-Kutta法汎用ソルバー 返事を書く

クスクス 2006/12/13 02:08:37

RK4 ソルバ　サンプル説明書き。

! 単振子シミレーション
!
! 振子の諸元
! ・L = 長さ, M = 質量.
!
! 変数の表す系の状態
! ・TT = 現在時刻, H = 時間間隔(定数).
!
! ・XX(1) = X 軸との角度,
! XX(2) = 角速度.
!
! このシミレーションでは
! ・XX(1) = θ (振子の角度),
! ・XX(2) = ω = dθ/dt (振子の角速度).
!
! 系を定める微分方程式
! ・d(c(t))/dt = DFEQ(c(t))
! ・4次のルンゲ・クッタ法により積分曲線(微分方程式の解) c(t) を求める。

　　│└RK4ソルバプログラム本体その1。 クスクス 2006/12/13 02:11:18  ツリーへ

Re: RK4ソルバサンプル説明書き。 返事を書く

クスクス 2006/12/13 02:11:18

RK4 ソルバプログラム本体その1。

DECLARE NUMERIC TI,TF,TT,H ! 初期時刻,終了時刻,現在時刻,時間間隔
DECLARE NUMERIC G,L,M,R ! 重力加速度,振子長,錘質量,錘半径
DECLARE STRING FN$ ! ログファイル名
DECLARE NUMERIC MARGIN,W ! 余白率
DECLARE NUMERIC I,mx,my,mlb,mrb ! ループ変数,マウス用変数
DEF N = 2 ! 方程式の数
DIM XX(N), FM(N) ! 現在の相点, 次の相点への変位ベクトル

CALL ICOND ! 初期条件設定

LET XX(1) = (PI/180)*XX(1) ! ラジアンに変換
LET R = SQR(M) ! 錘の半径。
LET MARGIN = MAX(R/(2*L),0.1) ! 余白率
LET W = (1+MARGIN)*L ! ウィンドウのサイズ
SET BITMAP SIZE 501,501
SET WINDOW -W,W, -W,W
DRAW AXES(L/2,L/2) WITH ROTATE(-PI/2)

! OPEN #1: NAME FN$ ! ログを取るファイルを開く
! ERASE #1 ! 以前の計算内容を消去

LET TT = TI ! 初期時刻をセット
DO
! PRINT #1: XX(1); ","; XX(2) ! ログ書き込み
SET DRAW MODE HIDDEN
CLEAR
PLOT TEXT ,AT 0,L:"マウス右ボタンで、終了。"
SET LINE WIDTH 1
DRAW AXES(L/2,L/2) WITH ROTATE(-PI/2) ! 座標軸を描く
SET LINE WIDTH 2
DRAW PENDULUM(XX(1)) ! θの位置に振子を描く
CALL RUNGE(XX,TT,H,FM) ! XX = 現在の相点, TT = 現在時刻, H = 刻み幅,
! FM = (現在位置 XX から次の位置への変位)/H
!! 現在の位置を更新する
FOR I=1 TO N
LET XX(I) = XX(I) + FM(I)*H
NEXT I
!! 時刻を更新する
LET TT = TT + H
SET DRAW MODE EXPLICIT
MOUSE POLL mx,my,mlb,mrb !マウスの状態取得。
WAIT DELAY 1/(7*ABS(LOG(H)))
LOOP UNTIL mrb=1

! CLOSE #1 ! ログファイルパスを閉じる。

!! 振子を描く
PICTURE PENDULUM(θ)
DRAW CIRCLE WITH SCALE(W/80)
DRAW SINGLE_PENDULUM(L,R) WITH ROTATE(θ-PI/2)
END PICTURE

!! 単振子イメージ
PICTURE SINGLE_PENDULUM(L,R) ! L = 長さ, R = 錘の半径
PLOT LINES: 0,0; L,0
DRAW DISK WITH SCALE(R)*SHIFT(L,0)
END PICTURE

　　│　└RK4ソルバプログラム本体その2。 クスクス 2006/12/13 02:12:25  ツリーへ

Re: RK4ソルバプログラム本体その1。 返事を書く

クスクス 2006/12/13 02:12:25

RK4 ソルバプログラム本体その2。

! Runge Kutta 4
! ・古典的4次ルンゲ・クッタ法。
! ・与えられた相空間の1点 XX に対し, 現在の時刻 TT に点 XX = c(TT) を
! 通る積分曲線を c とする。 H 時間後 c(TT) は c(TT+H) に変位する。
! RUNGE は単位時間当たりの変位ベクトル (c(TT+H)-c(TT))/H の近似値 φ(TT,H)
! を配列 FM 入れて返すソルバーである。
! ・φ(TT,H) の近似の程度は
! (c(TT+H)-c(TT))/H = φ(TT,H) + O(H^4).
! よって, ソルバー RUNGE は4次の近似を与える。
! ・c(TT+H) ≒ XX() + H*FM()
!
! ・入力：XX = 現在の相点, TT = 現在時刻, H = 刻み幅
! ・返り値：FM = 返り値(変位ベクトル)
SUB RUNGE(XX(),TT,H,FM())
DIM FF(N,4)
DIM X(N),F(N)
DECLARE NUMERIC I,J,T,HH

FOR J=1 TO 4 !! J=1 のとき HH=0 だからループに入る前に F() を零初期化する必要がない。
!! 時間間隔 H を更に時間間隔 HH に刻む
LET HH = IP(J/2)*H/2 ! J=1,2,3,4 に応じて, HH=0,H/2,1H/2,H
!! TT より HH 時間後の時刻 T
LET T = TT + HH
!! HH 時間後の積分曲線上の点の1次近似 X
FOR I=1 TO N
LET X(I) = XX(I) + F(I)*HH
NEXT I
!! 位置 X にある接ベクトル F を計算し, FF に保存
CALL DFEQ(X,F)
FOR I=1 TO N
LET FF(I,J) = F(I)
NEXT I
NEXT J
FOR I=1 TO N ! 変位ベクトル FF(I,J) (J=1,2,3,4) の平均が, 求める変位ベクトル FM(I)
LET FM(I) = (FF(I,1)+2*FF(I,2)+2*FF(I,3)+FF(I,4))/6
NEXT I
END SUB

! 微分方程式 (ベクトル場)
! ・　入力： X = 相点
! ・返り値： F = X における接ベクトル
SUB DFEQ(X(),F())
LET F(1) = X(2)
LET F(2) = -(G/L)*SIN(X(1)) ! -(G/L)*X(1)
END SUB

! 初期条件、計算条件の設定
SUB ICOND
READ TI, TF, H ! 初期時刻, 終了時刻, 時間間隔
DATA 0, 30, 0.02
READ XX(1), XX(2), G, L, M ! 初期角度, 初速, 重力加速度, 振子の長さ, 錘の質量
DATA 30, 0, 9.8, 1, 0.001
READ FN$
DATA "single_pendulum.log" ! ログファイル名
END SUB

END

　　│　　├有難うございます。 SECOND 2006/12/13 05:24:55  ツリーへ

Re: RK4ソルバプログラム本体その2。 返事を書く

SECOND 2006/12/13 05:24:55

有難うございます。
話の腰を折ったあげくに、
こんなにやってもらって申訳なく思います。
しっかり勉強します。

　　│　　└!これだけに、なりました。これと照合すると... SECOND 2006/12/16 00:32:26  ツリーへ

Re: RK4ソルバプログラム本体その2。 返事を書く

SECOND 2006/12/16 00:32:26

!これだけに、なりました。これと照合すると、
!θの扱いが、私の方に間違いが有ったようです。
!t は、無視されており、取除きました。

! 単振子シミレーション

!　初速, 重力加速度, 振子の長さ, 錘の質量, 時間間隔
READ ω1, g, L1, m1, dt
DATA 0, 9.8, 1, 0.002, 0.02

LET θ1= PI/6 !初期角度

!dθ1/dt= ω1
!dω1/dt= -g/L1*SIN(θ)

DEF α1(θ)= -g/L1 *SIN(θ)

SUB Runge_Kutta_4
LET ω11= ω1
LET α11= α1(θ1)
!----
LET ω12= ω1+α11*dt/2
LET α12= α1(θ1+ω11*dt/2)
!----
LET ω13= ω1+α12*dt/2
LET α13= α1(θ1+ω12*dt/2)
!----
LET ω14= ω1+α13*dt
LET α14= α1(θ1+ω13*dt)
!----
LET w1= (ω11+2*ω12+2*ω13+ω14)/6 !Δ角速度
LET a1= (α11+2*α12+2*α13+α14)/6 !Δ角加速度
!----
LET θ1= θ1+w1*dt
LET ω1= ω1+a1*dt
END SUB

LET w= (L1+r1)*1.1
SET WINDOW -w,w, -w,w
SET COLOR MIX(15) 0.3,0.3,0.3
LET r1= SQR(m1) ! 錘の半径。
DO
SET DRAW mode hidden
CLEAR
PLOT TEXT ,AT w*0.1,w*0.9:"マウス右ボタンで、終了。"
DRAW grid(L1/2, L1/2) WITH ROTATE(-PI/2)
DRAW 振子12 WITH ROTATE(θ1) !θ1の位置に振子。
CALL Runge_Kutta_4 ! dt= 刻み幅,
SET DRAW mode explicit
MOUSE POLL mx,my,mlb,mrb !マウスの状態取得。
WAIT DELAY dt
LOOP UNTIL mrb=1

PICTURE 振子12
DRAW circle WITH SCALE(w/80)
DRAW 振子(L1,r1) !(m1)
!なし DRAW 振子(L2,r2) WITH ROTATE(θ2-θ1)*SHIFT(0,-L1) !(m2)
END PICTURE

PICTURE 振子(L,r)
PLOT LINES: 0,0;0,-L
DRAW disk WITH SCALE(r)*SHIFT(0,-L)
END PICTURE

END

　　└RK4ソルバ動作確認(等速円運動で検査)。 クスクス 2006/12/13 15:26:27  ツリーへ

Re: 古典的4次Runge-Kutta法汎用ソルバー 返事を書く

クスクス 2006/12/13 15:26:27

RK4 ソルバ動作確認(等速円運動で検査)。

RK4 ソルバが正しく作動しているか, 等速円運動の微分方程式を
解かせ, その厳密解(cos(t),sin(t))と比較することにより, 検査
します。RK4 ソルバの使用サンプルとしては, 単振子よりは, こちら
の方が妥当でした。

動かすと刻み幅Hを聞いてきます。最初はリーターンのみの入力で
暗黙値の0.02で動かしてみてください。一昼夜動かしても, 数値解
が, 厳密解に近いことが観測できるでしょう。テキスト画面では
数字による比較が観測できます。大きな刻み幅 H=1 等では急速に
厳密解から離れていく様子が見えます。

　　　└等速円運動シミレーションその１． クスクス 2006/12/13 15:28:20  ツリーへ

Re: RK4ソルバ動作確認(等速円運動で検査)。 返事を書く

クスクス 2006/12/13 15:28:20

等速円運動シミレーション　その１．

! 等速回転運動　--- RK4 ソルバの動作テスト ---
! ・等速円運動の微分方程式を RK4 ソルバで数値的に解いて
! 　その厳密解と比較する.

DECLARE NUMERIC TT,H ! 現在時刻,時間間隔
DECLARE NUMERIC R ! 質点半径
DECLARE NUMERIC L ! 原点と質点を結ぶワイヤー長
DECLARE STRING FN$ ! ログファイル名
DECLARE NUMERIC MARGIN,W ! 余白率,余白付きウィンドウ幅
DECLARE NUMERIC I,mx,my,mlb,mrb ! ループ変数,マウス用変数
DECLARE STRING L$ ! 行バッファ
DEF N = 4 ! 方程式の数
DIM XX(N), FM(N) ! 現在の相点,次の相点への変位ベクトル

CALL ICOND ! 初期条件設定

LET FN$ = "double_pendulum.log" ! ログファイル名
LET H = 0.02 ! 時間間隔
LET R = 0.02
LET L = 1
LET MARGIN = MAX(R,0.1) ! 余白率
LET W = 1+MARGIN
SET BITMAP SIZE 501,501
SET WINDOW -W,W, -W,W
DRAW AXES(W/2,W/2)

LINE INPUT PROMPT "時間間隔(リターンで0.02) H = ":L$
IF L$<>"" THEN LET H = VAL(L$)
PRINT
PRINT " ************** 数値解の厳密解との比較 **************"

!OPEN #1: NAME FN$ ! ログを取るファイルを開く
!ERASE #1 ! 以前の計算内容を消去

PRINT " TT COS(TT) SIN(TT) X座標 Y座標"
DO
!PRINT #1: XX(1); ","; XX(2); ",", XX(3); ","; XX(4) ! ログへ記録
PRINT USING "####.### #.####### #.####### ##.####### ##.#######": TT,COS(TT),SIN(TT), XX(1),XX(2)
SET DRAW MODE HIDDEN
CLEAR
PLOT TEXT ,AT -W,L:"マウス右ボタンで終了。黒: 厳密解/2, 赤: 数値解"
DRAW AXES(L/4,L/4) ! 座標軸を描く
SET LINE COLOR "RED"
SET AREA COLOR "RED"
DRAW CIRCLE WITH SCALE(L)
DRAW PENDULUM(XX(1),XX(2),R)
SET LINE COLOR "BLACK"
SET AREA COLOR "BLACK"
DRAW CIRCLE WITH SCALE(L/2)
DRAW PENDULUM(COS(TT)/2,SIN(TT)/2,R)
CALL RUNGE(XX,TT,H,FM) ! FM(I) = (現在位置 XX(I) から次の位置への変位)/H
!! 現在の位置を更新する
FOR I=1 TO N
LET XX(I) = XX(I) + FM(I)*H
NEXT I
!! 時刻を更新する
LET TT = TT + H
SET DRAW MODE EXPLICIT
MOUSE POLL mx,my,mlb,mrb !マウスの状態取得。
WAIT DELAY H/20
LOOP UNTIL mrb=1

!CLOSE #1
! *********************** 主プログラム終わり ***********************

　　　　└その２。 クスクス 2006/12/13 15:29:37  ツリーへ

Re: 等速円運動シミレーションその１． 返事を書く

クスクス 2006/12/13 15:29:37

その２。

! 振子を描く
PICTURE PENDULUM(X,Y,R) ! (X,Y) = 錘の座標, L = 長さ, R = 錘の半径
PLOT LINES: 0,0; X,Y
DRAW DISK WITH SCALE(R)*SHIFT(X,Y)
END PICTURE

! Runge Kutta 4
SUB RUNGE(XX(),TT,H,FM())
DIM FF(N,4)
DIM X(N),F(N)
DECLARE NUMERIC I,J,T,HH

FOR J=1 TO 4
LET HH = IP(J/2)*H/2
LET T = TT + HH
FOR I=1 TO N
LET X(I) = XX(I) + F(I)*HH
NEXT I
CALL DFEQ(X,F)
FOR I=1 TO N
LET FF(I,J) = F(I)
NEXT I
NEXT J
FOR I=1 TO N
LET FM(I) = (FF(I,1)+2*FF(I,2)+2*FF(I,3)+FF(I,4))/6
NEXT I
END SUB

! DFEQ(X(),F()): 等速円運動の微分方程式
! ・シンボルの対応関係
! 　(F(1), F(2), F(3), F(4)) = (x'(t), y'(t), ξ'(t), η'(t)),
! 　(X(1), X(2), X(3), X(4)) = (x(t), y(t), ξ(t), η(t)).
! ・x'(t) = ξ(t),
! 　y'(t) = η(t),
! ξ'(t) = -x(t),
! η'(t) = -y(t).
! ・向きは正, 速さは1の単位円上の等速回転運動を表す.
! ・これは、初期条件 x(0)=1, y(0)=0, ξ(0)=0, η(0)=1 のもと, 厳密解
! 　x(t)=cos(t), y(t)=sin(t), ξ(t)=-sin(t), η(t)=cos(t)
! 　を持つ.
! ・この微分方程式は初等解析で三角関数の定義に採用されることがある.
SUB DFEQ(X(),F())
LET F(1) = X(3)
LET F(2) = X(4)
LET F(3) = -X(1)
LET F(4) = -X(2)
END SUB

! 初期条件の設定
SUB ICOND
READ XX(1), XX(2), XX(3), XX(4)
DATA 1, 0, 0, 1
END SUB

END

　　　　　└追加訂正分。 クスクス 2006/12/13 17:42:59  ツリーへ

Re: その２。 返事を書く

クスクス 2006/12/13 17:42:59

追加訂正分。

プログラムの行
　PRINT " TT COS(TT) SIN(TT) X座標 Y座標"
および
　PRINT USING "####.### #.####### #.####### ##.####### ##.#######": TT,COS(TT),SIN(TT), XX(1),XX(2)
をコメント・アウトし、その部分に各々行

PRINT " TT X座標 - COS(TT) Y座標 - SIN(TT)"

PRINT USING "####.### ###.##### ###.#####": TT,XX(1)-COS(TT),XX(2)-SIN(TT)

を挿入してください。誤差の推移がより見易くなります。

　　　　　　└これで十分です。詳細に追跡します。有難う... SECOND 2006/12/14 02:34:16  ツリーへ

Re: 追加訂正分。 返事を書く

SECOND 2006/12/14 02:34:16

これで十分です。詳細に追跡します。有難うございました。

インデックスへ EXIT
新規発言を反映させるにはブラウザの更新ボタンを押してください。

古典的4次Runge-Kutta法汎用ソルバー	返事を書く
クスクス 2006/12/13 02:06:08
4次のルンゲ・クッタ法プログラムは、Fortran,C,Java のコードならば、ネットのいたるところで手に入れることができますが、 FullBASIC のコードは探しても見当たりません。これは、RK4 の原理の理解を目的とした、ｎ元連立１階常微分方程式対応のRK4 ソルバです。使用法を説明するためのサンプルは、単振子のシミレーションです。このソルバは「パソコンで見る天体の動き」長沢工著　地人書刊のN88BASICで書かれたコードCLARK.BASの十進BASICへの移植です。

Re: 古典的4次Runge-Kutta法汎用ソルバー	返事を書く
クスクス 2006/12/13 02:08:37
RK4 ソルバ　サンプル説明書き。 ! 単振子シミレーション ! ! 振子の諸元 ! ・L = 長さ, M = 質量. ! ! 変数の表す系の状態 ! ・TT = 現在時刻, H = 時間間隔(定数). ! ! ・XX(1) = X 軸との角度, ! XX(2) = 角速度. ! ! このシミレーションでは ! ・XX(1) = θ (振子の角度), ! ・XX(2) = ω = dθ/dt (振子の角速度). ! ! 系を定める微分方程式 ! ・d(c(t))/dt = DFEQ(c(t)) ! ・4次のルンゲ・クッタ法により積分曲線(微分方程式の解) c(t) を求める。

Re: RK4ソルバサンプル説明書き。	返事を書く
クスクス 2006/12/13 02:11:18
RK4 ソルバプログラム本体その1。 DECLARE NUMERIC TI,TF,TT,H ! 初期時刻,終了時刻,現在時刻,時間間隔 DECLARE NUMERIC G,L,M,R ! 重力加速度,振子長,錘質量,錘半径 DECLARE STRING FN$ ! ログファイル名 DECLARE NUMERIC MARGIN,W ! 余白率 DECLARE NUMERIC I,mx,my,mlb,mrb ! ループ変数,マウス用変数 DEF N = 2 ! 方程式の数 DIM XX(N), FM(N) ! 現在の相点, 次の相点への変位ベクトル CALL ICOND ! 初期条件設定 LET XX(1) = (PI/180)XX(1) ! ラジアンに変換 LET R = SQR(M) ! 錘の半径。 LET MARGIN = MAX(R/(2L),0.1) ! 余白率 LET W = (1+MARGIN)L ! ウィンドウのサイズ SET BITMAP SIZE 501,501 SET WINDOW -W,W, -W,W DRAW AXES(L/2,L/2) WITH ROTATE(-PI/2) ! OPEN #1: NAME FN$ ! ログを取るファイルを開く ! ERASE #1 ! 以前の計算内容を消去 LET TT = TI ! 初期時刻をセット DO ! PRINT #1: XX(1); ","; XX(2) ! ログ書き込み SET DRAW MODE HIDDEN CLEAR PLOT TEXT ,AT 0,L:"マウス右ボタンで、終了。" SET LINE WIDTH 1 DRAW AXES(L/2,L/2) WITH ROTATE(-PI/2) ! 座標軸を描く SET LINE WIDTH 2 DRAW PENDULUM(XX(1)) ! θの位置に振子を描く CALL RUNGE(XX,TT,H,FM) ! XX = 現在の相点, TT = 現在時刻, H = 刻み幅, ! FM = (現在位置 XX から次の位置への変位)/H !! 現在の位置を更新する FOR I=1 TO N LET XX(I) = XX(I) + FM(I)H NEXT I !! 時刻を更新する LET TT = TT + H SET DRAW MODE EXPLICIT MOUSE POLL mx,my,mlb,mrb !マウスの状態取得。 WAIT DELAY 1/(7ABS(LOG(H))) LOOP UNTIL mrb=1 ! CLOSE #1 ! ログファイルパスを閉じる。 !! 振子を描く PICTURE PENDULUM(θ) DRAW CIRCLE WITH SCALE(W/80) DRAW SINGLE_PENDULUM(L,R) WITH ROTATE(θ-PI/2) END PICTURE !! 単振子イメージ PICTURE SINGLE_PENDULUM(L,R) ! L = 長さ, R = 錘の半径 PLOT LINES: 0,0; L,0 DRAW DISK WITH SCALE(R)SHIFT(L,0) END PICTURE

Re: RK4ソルバプログラム本体その1。	返事を書く
クスクス 2006/12/13 02:12:25
RK4 ソルバプログラム本体その2。 ! Runge Kutta 4 ! ・古典的4次ルンゲ・クッタ法。 ! ・与えられた相空間の1点 XX に対し, 現在の時刻 TT に点 XX = c(TT) を ! 通る積分曲線を c とする。 H 時間後 c(TT) は c(TT+H) に変位する。 ! RUNGE は単位時間当たりの変位ベクトル (c(TT+H)-c(TT))/H の近似値 φ(TT,H) ! を配列 FM 入れて返すソルバーである。 ! ・φ(TT,H) の近似の程度は ! (c(TT+H)-c(TT))/H = φ(TT,H) + O(H^4). ! よって, ソルバー RUNGE は4次の近似を与える。 ! ・c(TT+H) ≒ XX() + HFM() ! ! ・入力：XX = 現在の相点, TT = 現在時刻, H = 刻み幅 ! ・返り値：FM = 返り値(変位ベクトル) SUB RUNGE(XX(),TT,H,FM()) DIM FF(N,4) DIM X(N),F(N) DECLARE NUMERIC I,J,T,HH FOR J=1 TO 4 !! J=1 のとき HH=0 だからループに入る前に F() を零初期化する必要がない。 !! 時間間隔 H を更に時間間隔 HH に刻む LET HH = IP(J/2)H/2 ! J=1,2,3,4 に応じて, HH=0,H/2,1H/2,H !! TT より HH 時間後の時刻 T LET T = TT + HH !! HH 時間後の積分曲線上の点の1次近似 X FOR I=1 TO N LET X(I) = XX(I) + F(I)HH NEXT I !! 位置 X にある接ベクトル F を計算し, FF に保存 CALL DFEQ(X,F) FOR I=1 TO N LET FF(I,J) = F(I) NEXT I NEXT J FOR I=1 TO N ! 変位ベクトル FF(I,J) (J=1,2,3,4) の平均が, 求める変位ベクトル FM(I) LET FM(I) = (FF(I,1)+2FF(I,2)+2FF(I,3)+FF(I,4))/6 NEXT I END SUB ! 微分方程式 (ベクトル場) ! ・　入力： X = 相点 ! ・返り値： F = X における接ベクトル SUB DFEQ(X(),F()) LET F(1) = X(2) LET F(2) = -(G/L)SIN(X(1)) ! -(G/L)*X(1) END SUB ! 初期条件、計算条件の設定 SUB ICOND READ TI, TF, H ! 初期時刻, 終了時刻, 時間間隔 DATA 0, 30, 0.02 READ XX(1), XX(2), G, L, M ! 初期角度, 初速, 重力加速度, 振子の長さ, 錘の質量 DATA 30, 0, 9.8, 1, 0.001 READ FN$ DATA "single_pendulum.log" ! ログファイル名 END SUB END

Re: RK4ソルバプログラム本体その2。	返事を書く
SECOND 2006/12/13 05:24:55
有難うございます。話の腰を折ったあげくに、こんなにやってもらって申訳なく思います。しっかり勉強します。

Re: RK4ソルバプログラム本体その2。	返事を書く
SECOND 2006/12/16 00:32:26
!これだけに、なりました。これと照合すると、 !θの扱いが、私の方に間違いが有ったようです。 !t は、無視されており、取除きました。 ! 単振子シミレーション !　初速, 重力加速度, 振子の長さ, 錘の質量, 時間間隔 READ ω1, g, L1, m1, dt DATA 0, 9.8, 1, 0.002, 0.02 LET θ1= PI/6 !初期角度 !dθ1/dt= ω1 !dω1/dt= -g/L1SIN(θ) DEF α1(θ)= -g/L1 SIN(θ) SUB Runge_Kutta_4 LET ω11= ω1 LET α11= α1(θ1) !---- LET ω12= ω1+α11dt/2 LET α12= α1(θ1+ω11dt/2) !---- LET ω13= ω1+α12dt/2 LET α13= α1(θ1+ω12dt/2) !---- LET ω14= ω1+α13dt LET α14= α1(θ1+ω13dt) !---- LET w1= (ω11+2ω12+2ω13+ω14)/6 !Δ角速度 LET a1= (α11+2α12+2α13+α14)/6 !Δ角加速度 !---- LET θ1= θ1+w1dt LET ω1= ω1+a1dt END SUB LET w= (L1+r1)1.1 SET WINDOW -w,w, -w,w SET COLOR MIX(15) 0.3,0.3,0.3 LET r1= SQR(m1) ! 錘の半径。 DO SET DRAW mode hidden CLEAR PLOT TEXT ,AT w0.1,w0.9:"マウス右ボタンで、終了。" DRAW grid(L1/2, L1/2) WITH ROTATE(-PI/2) DRAW 振子12 WITH ROTATE(θ1) !θ1の位置に振子。 CALL Runge_Kutta_4 ! dt= 刻み幅, SET DRAW mode explicit MOUSE POLL mx,my,mlb,mrb !マウスの状態取得。 WAIT DELAY dt LOOP UNTIL mrb=1 PICTURE 振子12 DRAW circle WITH SCALE(w/80) DRAW 振子(L1,r1) !(m1) !なし DRAW 振子(L2,r2) WITH ROTATE(θ2-θ1)SHIFT(0,-L1) !(m2) END PICTURE PICTURE 振子(L,r) PLOT LINES: 0,0;0,-L DRAW disk WITH SCALE(r)*SHIFT(0,-L) END PICTURE END

Re: 古典的4次Runge-Kutta法汎用ソルバー	返事を書く
クスクス 2006/12/13 15:26:27
RK4 ソルバ動作確認(等速円運動で検査)。 RK4 ソルバが正しく作動しているか, 等速円運動の微分方程式を解かせ, その厳密解(cos(t),sin(t))と比較することにより, 検査します。RK4 ソルバの使用サンプルとしては, 単振子よりは, こちらの方が妥当でした。動かすと刻み幅Hを聞いてきます。最初はリーターンのみの入力で暗黙値の0.02で動かしてみてください。一昼夜動かしても, 数値解が, 厳密解に近いことが観測できるでしょう。テキスト画面では数字による比較が観測できます。大きな刻み幅 H=1 等では急速に厳密解から離れていく様子が見えます。

Re: RK4ソルバ動作確認(等速円運動で検査)。	返事を書く
クスクス 2006/12/13 15:28:20
等速円運動シミレーション　その１． ! 等速回転運動　--- RK4 ソルバの動作テスト --- ! ・等速円運動の微分方程式を RK4 ソルバで数値的に解いて ! 　その厳密解と比較する. DECLARE NUMERIC TT,H ! 現在時刻,時間間隔 DECLARE NUMERIC R ! 質点半径 DECLARE NUMERIC L ! 原点と質点を結ぶワイヤー長 DECLARE STRING FN$ ! ログファイル名 DECLARE NUMERIC MARGIN,W ! 余白率,余白付きウィンドウ幅 DECLARE NUMERIC I,mx,my,mlb,mrb ! ループ変数,マウス用変数 DECLARE STRING L$ ! 行バッファ DEF N = 4 ! 方程式の数 DIM XX(N), FM(N) ! 現在の相点,次の相点への変位ベクトル CALL ICOND ! 初期条件設定 LET FN$ = "double_pendulum.log" ! ログファイル名 LET H = 0.02 ! 時間間隔 LET R = 0.02 LET L = 1 LET MARGIN = MAX(R,0.1) ! 余白率 LET W = 1+MARGIN SET BITMAP SIZE 501,501 SET WINDOW -W,W, -W,W DRAW AXES(W/2,W/2) LINE INPUT PROMPT "時間間隔(リターンで0.02) H = ":L$ IF L$<>"" THEN LET H = VAL(L$) PRINT PRINT " ************ 数値解の厳密解との比較 ***********" !OPEN #1: NAME FN$ ! ログを取るファイルを開く !ERASE #1 ! 以前の計算内容を消去 PRINT " TT COS(TT) SIN(TT) X座標 Y座標" DO !PRINT #1: XX(1); ","; XX(2); ",", XX(3); ","; XX(4) ! ログへ記録 PRINT USING "####.### #.####### #.####### ##.####### ##.#######": TT,COS(TT),SIN(TT), XX(1),XX(2) SET DRAW MODE HIDDEN CLEAR PLOT TEXT ,AT -W,L:"マウス右ボタンで終了。黒: 厳密解/2, 赤: 数値解" DRAW AXES(L/4,L/4) ! 座標軸を描く SET LINE COLOR "RED" SET AREA COLOR "RED" DRAW CIRCLE WITH SCALE(L) DRAW PENDULUM(XX(1),XX(2),R) SET LINE COLOR "BLACK" SET AREA COLOR "BLACK" DRAW CIRCLE WITH SCALE(L/2) DRAW PENDULUM(COS(TT)/2,SIN(TT)/2,R) CALL RUNGE(XX,TT,H,FM) ! FM(I) = (現在位置 XX(I) から次の位置への変位)/H !! 現在の位置を更新する FOR I=1 TO N LET XX(I) = XX(I) + FM(I)H NEXT I !! 時刻を更新する LET TT = TT + H SET DRAW MODE EXPLICIT MOUSE POLL mx,my,mlb,mrb !マウスの状態取得。 WAIT DELAY H/20 LOOP UNTIL mrb=1 !CLOSE #1 ! ********************* 主プログラム終わり *********************

Re: 等速円運動シミレーションその１．	返事を書く
クスクス 2006/12/13 15:29:37
その２。 ! 振子を描く PICTURE PENDULUM(X,Y,R) ! (X,Y) = 錘の座標, L = 長さ, R = 錘の半径 PLOT LINES: 0,0; X,Y DRAW DISK WITH SCALE(R)SHIFT(X,Y) END PICTURE ! Runge Kutta 4 SUB RUNGE(XX(),TT,H,FM()) DIM FF(N,4) DIM X(N),F(N) DECLARE NUMERIC I,J,T,HH FOR J=1 TO 4 LET HH = IP(J/2)H/2 LET T = TT + HH FOR I=1 TO N LET X(I) = XX(I) + F(I)HH NEXT I CALL DFEQ(X,F) FOR I=1 TO N LET FF(I,J) = F(I) NEXT I NEXT J FOR I=1 TO N LET FM(I) = (FF(I,1)+2FF(I,2)+2*FF(I,3)+FF(I,4))/6 NEXT I END SUB ! DFEQ(X(),F()): 等速円運動の微分方程式 ! ・シンボルの対応関係 ! 　(F(1), F(2), F(3), F(4)) = (x'(t), y'(t), ξ'(t), η'(t)), ! 　(X(1), X(2), X(3), X(4)) = (x(t), y(t), ξ(t), η(t)). ! ・x'(t) = ξ(t), ! 　y'(t) = η(t), ! ξ'(t) = -x(t), ! η'(t) = -y(t). ! ・向きは正, 速さは1の単位円上の等速回転運動を表す. ! ・これは、初期条件 x(0)=1, y(0)=0, ξ(0)=0, η(0)=1 のもと, 厳密解 ! 　x(t)=cos(t), y(t)=sin(t), ξ(t)=-sin(t), η(t)=cos(t) ! 　を持つ. ! ・この微分方程式は初等解析で三角関数の定義に採用されることがある. SUB DFEQ(X(),F()) LET F(1) = X(3) LET F(2) = X(4) LET F(3) = -X(1) LET F(4) = -X(2) END SUB ! 初期条件の設定 SUB ICOND READ XX(1), XX(2), XX(3), XX(4) DATA 1, 0, 0, 1 END SUB END

Re: その２。	返事を書く
クスクス 2006/12/13 17:42:59
追加訂正分。プログラムの行　PRINT " TT COS(TT) SIN(TT) X座標 Y座標" および　PRINT USING "####.### #.####### #.####### ##.####### ##.#######": TT,COS(TT),SIN(TT), XX(1),XX(2) をコメント・アウトし、その部分に各々行 PRINT " TT X座標 - COS(TT) Y座標 - SIN(TT)" PRINT USING "####.### ###.##### ###.#####": TT,XX(1)-COS(TT),XX(2)-SIN(TT) を挿入してください。誤差の推移がより見易くなります。

Re: 追加訂正分。	返事を書く
SECOND 2006/12/14 02:34:16
これで十分です。詳細に追跡します。有難うございました。