************************************************************************** * selflinear.f * * 定数係数線形自励系常微分方程式の解のグラフ * * (C) Kenrou Adachi 2004/11/22 * * * * dv/dt = Av * * v = (x, y), A = [[s, t], [u, v]] * * solution v(t) = v(0) + exp(tA)v(0) * * * * You needs to compile etude_x.c with eggx graphics library * ************************************************************************** program selflin real a, b, h, sh real s, t, u, v character key integer l, m, n integer nwins * 挨拶 call message * 係数行列の入力 * d [x] [ s t ][x] * --- [ ] = [ ][ ] * dt [y] [ u v ][y] call inputcoeff(s, t, u, v) * データの設定 a = 6.4 b = 4.0 m = 3 l = 50 * グラフィックの初期化 call gopen(640, 400, nwins) call gsetbgcolor(nwins, 'black\0') call gclr(nwins) * きざみ幅の設定 sh = (a / 320.0) * (a / 320.0) + (b / 200.0) * (b / 200.0) h = 2.0 * sqrt(sh) * 座標軸の描画 call axes(nwins) * 繰り返し処理 do 10 n = 1, m call drawing(a, b, h, l, s, t, u, v, nwins) 10 continue * 終了処理 write(*,*) '終了しました。' write(*,*) '任意の key に続いてEnter key を押してください。' read(*,*) key call gclose(nwins) stop end subroutine drawing(a, b, h, l, s, t, u, v, nwins) real a, b, h, s, t, u, v integer l, nwins integer j, count real rnd, e1, e2, e, x(1:40), y(1:40), xx(1:40), yy(1:40) real eh, AA, BB, CC, DD integer GREEN, RED parameter (GREEN = 3, RED = 2) * 描画範囲に４０個の初期点をランダムに生成する j = 1 do 20 e1 = -1.0, 0.9, 0.25 do 30 e2 = -1.0, 0.9, 0.4 rnd = rand() xx(j) = a * (e1 + 0.25 * rnd) rnd = rand() yy(j) = b * (e2 + 0.4 * rnd) j = j + 1 30 continue 20 continue * 軌道計算 routine を呼出す。 e = 1.0 40 continue if (e .gt. 0.0) then call newpencolor(nwins, GREEN) else call newpencolor(nwins, RED) endif do 50 j = 1, 40 x(j) = xx(j) y(j) = yy(j) 50 continue * Ah + (Ah)^2/2 + (Ah)^3/6 + (Ah)^4/24 + (Ah)^5/120 eh = e * h call detectcoeff(s, t, u, v, AA, BB, CC, DD, eh) * 軌道計算 do 60 count = 1, l call culcul(x, y, a, b, h, AA, BB, CC, DD, nwins) 60 continue * 時間を逆転して過去を計算 e = -e if (e.lt.0.0) goto 40 return end subroutine culcul(x, y, a, b, h, AA, BB, CC, DD, nwins) real x(1:40), y(1:40), a, b, h, AA, BB, CC, DD integer nwins integer j, PENUP, PENDOWN real x0, y0, r0, wx, wy, wx1, wy1 parameter (PENUP = 3, PENDOWN = 2) do 70 j = 1, 40 x0 = x(j) y0 = y(j) if ((abs(x0) .gt. 2.0 * a) .or. (abs(y0) .gt. 2.0 * b)) goto 70 r0 = sqrt(x0 * x0 + y0 * y0) if (r0 .lt. h) goto 70 * x(t) = x(0) + exp(tA)x(0) wx1 = x0 + (AA * x0 + BB * y0) wy1 = y0 + (CC * x0 + DD * y0) x(j) = wx1 y(j) = wy1 * Window に描画範囲を割当るように座標変換する wx = 320.0 * x0 / a wy = 200.0 * y0 / b wx1 = 320.0 * wx1 / a wy1 = 200.0 * wy1 / b * 線素の描画 call line(nwins, wx + 320.0, wy + 200.0, PENUP) call line(nwins, wx1 + 320.0, wy1 + 200.0, PENDOWN) 70 continue return end * 座標軸の描画 subroutine axes(nwins) integer nwins integer BLUE parameter (BLUE = 4) call newpencolor(nwins, BLUE); call line(nwins, 0.0, 200.0, PENUP) call line(nwins, 639.0, 200.0, PENDOWN) call line(nwins, 320.0, 0.0, PENUP) call line(nwins, 320.0, 399.0, PENDOWN) return end * 係数行列の入力 subroutine inputcoeff(s, t, u, v) real s, t, u, v write(*,*) ' 係数行列 A = [[s, t], [u, v]] を入力して下さい。' write(*,*) ' s = ? ' read(*,*) s write(*,*) ' t = ? ' read(*,*) t write(*,*) ' u = ? ' read(*,*) u write(*,*) ' v = ? ' read(*,*) v return end * 挨拶 subroutine message write(*,*) ' このプログラムは2次元定数係数線形自励系微分方程式: ' write(*,*) ' dx/dt = A * x の解軌道群を描画するものです。 ' write(*,*) ' 計算には解公式: ' write(*,*) ' x(t) = x(0) + exp(t * A) * x(0)' write(*,*) ' の近似式を使います。' return end subroutine detectcoeff(s, t, u, v, AA, BB, CC, DD, eh) real s, t, u, v, AA, BB, CC, DD, eh real s2, t2, u2, v2, s3, t3, u3, v3 real s4, t4, u4, v4, s5, t5, u5, v5 * dx/dt = x + hAx + ((hA)^2/2)x * + ((hA)^3/6)x + ((hA^4)/24)x + ((hA)^5/120)x * A^2 s2 = s * s + t * u t2 = s * t + t * v u2 = u * s + v * u v2 = u * t + v * v * A^3 s3 = s2 * s + t2 * u t3 = s2 * t + t2 * v u3 = u2 * s + v2 * u v3 = u2 * t + v2 * v * A^4 s4 = s2 * s2 + t2 * u2 t4 = s2 * t2 + t2 * v2 u4 = u2 * s2 + v2 * u2 v4 = u2 * t2 + v2 * v2 * A^5 s5 = s2 * s3 + t2 * u3 t5 = s2 * t3 + t2 * v3 u5 = u2 * s3 + v2 * u3 v5 = u2 * t3 + v2 * v3 * exp(hA) - I の Maclaurin 近似 AA = eh * (s + eh * (s2 / 2.0 + eh * (s3 / 6.0 + eh $ * (s4 / 24.0 + eh * s5 / 120.0)))) BB = eh * (t + eh * (t2 / 2.0 + eh * (t3 / 6.0 + eh $ * (t4 / 24.0 + eh * t5 / 120.0)))) CC = eh * (u + eh * (u2 / 2.0 + eh * (u3 / 6.0 + eh $ * (u4 / 24.0 + eh * u5 / 120.0)))) DD = eh * (v + eh * (v2 / 2.0 + eh * (v3 / 6.0 + eh $ * (v4 / 24.0 + eh * v5 / 120.0)))) return end ***** end of selflinear.f *****