* dynamics-f.f
**********************************************************************
**********************************************************************
*
* 2 次元の定数係数線形常微分方程式
* x'(t) = a x + b y
* y'(t) = c x + d y
* に初期条件
* x(0)=x0
* y(0)=y0
* を課した常微分方程式の初期値問題を解いて、相図を描く。
*
* このプログラムは次の4つの部分から出来ている。
* program dyna
* 主プログラム
* 行列の係数の入力、ウィンドウのオープン等の初期化をした後に、
* ユーザーにメニュー形式でコマンドを入力してもらう。
* 実際の計算のほとんどは他の副プログラムに任せている。
* subroutine draworbit(x0,y0,h,tlimit)
* (x0,y0) を初期値とする解の軌道を描く。
* 刻み幅を h、追跡時間を tlimit とする。
* 近似解の計算には Runge-Kutta 法を用いる。
* function fx(x,y)
* 微分方程式の右辺の x 成分
* function fy(x,y)
* 微分方程式の右辺の y 成分
**********************************************************************
**********************************************************************
program dyna
* 初期値
real x0,y0
* 時間の刻み幅、追跡時間
real h,tlimit,newh,newT
* 係数行列 A の成分（common 文により function fx,fy と共有する）
real a, b, c, d
common /coef/a,b,c,d
* ウィンドウに表示する範囲（common 文により draworbit と共有する）
real xleft, xright, ybottom, ytop
common /region/xleft, xright, ybottom, ytop
* メニューに対して入力されるコマンドの番号
integer cmd
* マウスのボタンの状態
integer lbut, mbut, rbut
* ウィンドウに表示する範囲の設定
xleft = -1.0
xright = 1.0
ybottom = -1.0
ytop = 1.0
* 時間刻み幅、追跡時間（とりあえず設定）
h=0.01
tlimit=10.0
* 行列の成分を入力
write(*,*) ' a,b,c,d='
read(*,*) a,b,c,d
* ウィンドウを開く
call openpl
call fspace(xleft, ybottom, xright, ytop)
* x 軸、y 軸を描く
call fline(xleft, 0.0, xright, 0.0)
call fline(0.0, ybottom, 0.0, ytop)
call xsync
* メイン・ループの入口
100 continue
* メニューを表示して、何をするか、番号で選択してもらう
write(*,*) ' したいことを番号で選んで下さい。'
* write(*,*) ' -1:finish, 0:keyboard input x0,y0, 1:mouse input,'
write(*,*) ' -1:メニュー終了, 0:初期値のキーボード入力,',
- ' 1:初期値のマウス入力,'
write(*,1000) h,tlimit
1000 format(' 2:change h,T(h=',f7.4,',T=',f7.4,')')
* 1000 format(' 2:刻み幅,追跡時間変更(現在 h=',f7.4,',T=',f7.4,')')
read(*,*) cmd
* 番号 cmd に応じて、指示された仕事をする
if (cmd .eq. 0) then
* 初期値の入力
write(*,*) ' 初期値 x0,y0='
read(*,*) x0,y0
call draworbit(x0,y0,h,tlimit)
else if (cmd .eq. 1) then
*
200 continue
write(*,*) 'マウスの左ボタンで初期値を指定して下さい',
- '（右ボタンで中止）。'
call fmouse(lbut, mbut, rbut, x0, y0)
if (lbut .eq. 1) then
write(*,*) '(x0,y0)=', x0,y0
call draworbit(x0,y0,h,tlimit)
goto 200
else if (mbut .eq. 1) then
write(*,*) '(x0,y0)=', x0,y0
call draworbit(x0,y0,-h,tlimit)
goto 200
endif
write(*,*) 'マウスによる初期値の入力を打ち切ります。'
else if (cmd .eq. 2) then
* 時間刻み幅、追跡時間の変更
write(*,*) ' 時間刻み幅 h, 追跡時間 T='
read(*,*) newh, newT
if (newh .ne. 0 .and. newT .ne. 0) then
h = newh
tlimit = newT
write(*,*) ' 新しい時間刻み幅 h = ', h,
- ' 新しい追跡時間 T = ', tlimit
else
write(*,*) ' h=', h, ' T=', newT, ' は不適当です。'
endif
else if (cmd .eq. -1) then
* 終了 -- メイン・ループを抜ける
goto 999
endif
* メイン・ループの先頭に戻る
goto 100
*
999 continue
call mkplot("dynamics-f.plot")
write(*,*) 'fplotウィンドウを左ボタンでクリックして下さい'
call closepl
end
**********************************************************************
**********************************************************************
* 指示された初期値に対する解軌道を描く
subroutine draworbit(x0,y0,h,tlimit)
real x0,y0,h,tlimit
real x,y,fx,fy
real k1x,k1y,k2x,k2y,k3x,k3y,k4x,k4y,t
real xleft, xright, ybottom, ytop
common /region/xleft, xright, ybottom, ytop
* 時刻を 0 にセットする
t = 0.0
* 初期値のセット
x = x0
y = y0
* 初期点を描く
call fpoint(x,y)
* ループの入口
100 continue
* Runge-Kutta 法による計算
* k1 の計算
k1x = h * fx(x, y)
k1y = h * fy(x, y)
* k2 の計算
k2x = h * fx(x + k1x / 2.0, y + k1y / 2.0)
k2y = h * fy(x + k1x / 2.0, y + k1y / 2.0)
* k3 の計算
k3x = h * fx(x + k2x / 2.0, y + k2y / 2.0)
k3y = h * fy(x + k2x / 2.0, y + k2y / 2.0)
* k4 の計算
k4x = h * fx(x + k3x, y + k3y)
k4y = h * fy(x + k3x, y + k3y)
* (Xn+1,Yn+1) の計算
x = x + (k1x + 2.0 * k2x + 2.0 * k3x + k4x) / 6.0
y = y + (k1y + 2.0 * k2y + 2.0 * k3y + k4y) / 6.0
* 解軌道を延ばす
call fcont(x,y)
* 時刻を 1 ステップ分進める
t = t + h
* まだ範囲内かどうかチェック
if (abs(t) .le. abs(tlimit)) goto 100
* サブルーチンを抜ける
999 continue
* 確実に描画が終るのを待つ
call xsync
end
*********************************************************************
*********************************************************************
* 微分方程式の右辺のベクトル値関数 f の x 成分
real function fx(x,y)
real x,y
real a, b, c, d
common /coef/a,b,c,d
fx = a * x + b * y
end
* 微分方程式の右辺のベクトル値関数 f の y 成分
real function fy(x,y)
real x,y
real a, b, c, d
common /coef/a,b,c,d
fy = c * x + d * y
end