グラフのプロット

以下に述べる方法でグラフのプロットができる. 必要な関数はライブラリファイル drill に入っているので, まず最初に

[1] load("drill");

を行う. 描画に必要な関数は次の通り.

1. 座標系を設定してウィンドウを生成する

   
   [2] glib_window(X0,Y0,X1,Y1);
   

   (X0,Y0) で左上, (X1,Y1) で右下の座標を与える. 座標系は, x 座標は左から右, y 座標は上から下に向かって増える. y 方向に関して通常と逆転していることに注意. したがって, 例えば glib_window(0,-5,50,5) は正しい呼び出し方.

2. 点を打つ

   
   [2] glib_putpixel(X,Y);
   

   座標 (X,Y) に点を打つ.

3. 線を引く

   
   [2] glib_line(X0,Y0,X1,Y1);
   

   座標 (X0,Y0) と (X1,X1) を直線で結ぶ.

4. ウィンドウのクリア

   現在書かれている点,線を全て消去する.

例 13.3   単振動の方程式を解いて, グラフを書くプログラム.

def osci0(K) {
    glib_window(0,-2,100,2);
    glib_clear();
    glib_line(0,0,100,0);

    X1 = 1.0; X2 = 1.0; T = 0; H = 0.01;
    while (T<100) {
        X3 = (2-K*H^2)*X2-X1;
        glib_putpixel(T,X1);
        T+=H;
        X1=X2; X2=X3;
    }
}

K を変えると何が起こるか?

外力が加わる場合でも同様である. $f(t)$ を時刻 $t$ に依存する外力とする とき, 方程式は

$$y'' + y = f(t)$$

この場合も

$${{y_{k+1}-2y_k+y_{k-1}} \over h^2} + y_k = f(kh)$$

より

$$y_{k+1}=(2-h^2)y_k - y_{k-1}+h^2f(kh)$$

からプログラムを書けばよい. もっと一般に

$$y'' + a(t)y'+b(t)y = f(t)$$

でも同じ. $y'$ を $(y_{k+1}-y_k)/h$ で 置き換えて

$$y_{k+1}={{(2+ha(t)-h^2b(t))y_k-y_{k-1}+h^2f(t)}\over {1+ha(t)}}$$

演習 13.3   強制振動の方程式

$$y''+y =\sin(at), \quad y(0) = 1, \quad y'(0) = 0 \quad (a > 0)$$

の数値近似解を求めるプログラムを書け. $a$ を引数とし, $a$ をいろいろ変えて 結果のグラフを鑑賞する. 特に, $a=1$ のとき何が起こるか? $a \neq 1$ での解

$$y = \cos(t) + {1 \over {1-a^2}} (\sin(at)-a\sin(t))$$

から $a=1$ の解を導き納得せよ.

演習 13.4   減衰振動の方程式

$$y''+ay'+y =0, \quad y(0) = 1, \quad y'(0) = 0$$

の数値近似解を求めるプログラムを書け. $a$ を引数とし, $a$ をいろいろ変えて 様子を見よ. ($a < 0$, $a = 0$, $a > 0$ を比較する.)