********************************************* * etpro.f 自励系微分方程式の解曲線群 * * (c) Kenrou Adachi, 2004/11/18 * * FreeBSD eggx version * ********************************************* ****** 主ルーチン ***** program etude real a, b, eps, h, sh integer n, l, i, nwins character key * *** データの入力 *** call inputdata(a, b, eps, n, l) * examples * a = 12.8, b = 8.0, eps = 0.2, n = 5, l = 100 * *** グラフィックスシステムの初期化 *** call gopen(640, 400, nwins) call gsetbgcolor(nwins, 'black\0') call gclr(nwins) * *** 座標軸の描画 *** call drawaxes(nwins) * *** 時間刻み幅の設定 *** sh = (a / 320.0) * (a /320.0) + (b / 200.0) * (b / 200.0) h = 2.0 * sqrt(sh) * *** 解曲線群描画 *** do 10 i = 1, n call drawing(a, b, eps, h, l, nwins) 10 continue * *** 処理終了 *** write(*,*) '終了しました。' write(*,*) '任意のキーに続いて Enter キーを押して下さい。' read(*,*) key call gclose(nwins) stop end ***** 解曲線群描画ルーチン ***** subroutine drawing(a, b, eps, h, l, nwins) real a, b, eps, h integer l, nwins integer j, count real e1, e2, e, x(1:40), y(1:40), xx(1:40), yy(1:40) real rnd 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 * *** 40 本の積分曲線を描く *** * * 時間の向きを正に * 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 do 60 count = 1, l call calcul(x, y, a, b, eps, h, e, nwins) 60 continue * * 時間の向きを逆転してもう一度 * e = -e * * 逆転は一度だけ * if (e .lt. 0.0) goto 40 return end ***** 微小積分曲線素描画ルーチン ***** subroutine calcul(x, y, a, b, eps, h, e, nwins) real x(1:40), y(1:40), a, b, eps, h, e integer nwins integer j real x0, y0, r0, x1, y1, r1, wx, wy, wx1, wy1 real u, v, u1, v1, q, qh do 70 j = 1, 40 x0 = x(j) y0 = y(j) if ((abs(x0).gt.(2.0 * a)).and.(abs(y0).gt.(2.0 * b))) goto 70 r0 = sqrt(x0 * x0 + y0 * y0) if (r0 .lt. h) goto 70 u = fx(x0, y0, r0) v = gy(x0, y0, r0) q = sqrt(u * u + v * v) if (q .lt. eps) goto 70 qh = h / q x1 = x0 + 0.5 * e * u * qh y1 = y0 + 0.5 * e * v * qh r1 = sqrt(x1 * x1 + y1 * y1) u1 = fx(x1, y1, r1) v1 = gy(x1, y1, r1) x(j) = x0 + e * u1 * qh y(j) = y0 + e * v1 * qh wx = 320.0 * x0 / a wy = 200.0 * y0 / b wx1 = 320.0 * x(j) / a wy1 = 200.0 * y(j) / b call line(nwins, wx + 320.0, 200.0 + wy, 3) call line(nwins, wx1 + 320.0, 200.0 + wy1, 2) 70 continue return end ***** 座標軸描画ルーチン ***** subroutine drawaxes(nwins) integer nwins, BLUE parameter (BLUE = 4) call newpencolor(nwins, BLUE) call line(nwins, 0.0, 200.0, 3) call line(nwins, 639.0, 200.0, 2) call line(nwins, 320.0, 0.0, 3) call line(nwins, 320.0, 399.0, 2) return end ***** 自励系微分方程式の定義 ***** ***** dx/dt = fx(x, y, r), dy/dt = gy(x, y, r) ***** * real function fx(x, y, r) * real x, y, r * fx = cos(x) - x * sin(y) * return * end * real function gy(x, y, r) * real x, y, r * gy = x * sin(y) + y * sin(x) * return * end ***** 描画データの入力 ***** subroutine inputdata(a, b, eps, n, l) real a, b, eps integer n, l write(*,*) '描画領域の範囲を入力して下さい。' write(*,*) 'x = [-a, a], y = [-b,b] です。' write(*,*) 'Input a = ?' read(*,*) a write(*,*) 'Input b = ?' read(*,*) b write(*,*) '精度限界を入力して下さい。' write(*,*) 'Input eps = ?' read(*,*) eps write(*,*) '描画回数を入力して下さい。' write(*,*) 'Input n = ?' read(*,*) n write(*,*) '計算回数を入力して下さい。' write(*,*) 'Input l = ?' read(*,*) l return end * g77 has rand() ***** 疑似乱数の発生 ***** ***** 大駒誠一 著、改訂 FORTRAN 77、サイエンス社 より ***** * real function rand() * integer c, k, r, m * real MAXR * parameter (c = 656329, k = 163, m = 12518383, MAXR = 12518383.0) * save r * data r /12345678/ * r = mod(k * r + c, m) * rand = real(r) / MAXR * return * end