計算数理A演習第4回 – Masashi Shiraishi Homepage

第４回（2016年5月1日）

今日の演習

前回の復習
- $y=\sin(x)$ グラフの描画
- $y=\sin(x+t)$ のアニメーションの作成
簡単な常微分方程式とその解のグラフ化
- 常微分方程式の数値解法：オイラー法

y = sin(x) のグラフの描画

これから、GLSCを使って $y=\sin(x)$ のグラフを描きます。
実際には $y=\sin(x)$ のグラフは滑らかな曲線ですが、コンピュータではそのような滑らかな曲線を直接は扱えませんので折れ線で近似します。
つまり、変数 $x$ は、実際には連続的に変化するものですが、とびとびの値を使います。
後の例で見るように、とびとびとはいえ、十分に細かく取れば、折れ線で曲線をうまく近似することができます。

図のように $[0,L]$ 区間を $N-1$ 等分して、その等分した小区間の幅を $dx$ とします。
$N-1$ 等分すると図のように $N$ 個の区切りが現れますが、それらに $0, 1, \dots, i, \dots, N-1$ と番号をつけることにします。
また、 $N$ を分割数、 $dx$ を分割幅と呼びます。このとき、

の部分についてそれぞれ $y=\sin(X_i)$ を計算し、それらを折れ線で結ぶことを考えます。
C言語風に書けば、

X[i] = dx*i;
Y[i] = sin(X[i]);

をそれぞれ求め、(X[0], Y[0]), (X[1],Y[1]), …., (X[N-1], Y[N-1]) を結ぶ折れ線を描くことになります。
これをプログラムにすると次のようになります。（上記の説明にあわせるためにかなり無駄なことをしています。）
ファイル名を 0501-1.c として打ち込み、実行してみてください。

注意：　数学関数(sin, cos, tan 等)を使う場合には #include <math.h> が必要となります。

#include <glsc.h>
#include <stdio.h>
/* 次の文は数学関数を使う場合に必要 */
#include <math.h>

/* 定数の定義 */
#define N (5)
#define PI (3.1415926)
#define L (2*PI)

int main()
{
    int i;
    double X[N], Y[N], dx;

    /* dx を求める */
    dx = L/(N - 1);

    /* GLSCの初期化および仮想座標系の定義 */
    g_init("GRAPH", 200.0, 100.0);
    g_device(G_DISP);
    g_def_scale(1, 0.0, L, -1.0, 1.0, 10.0, 10.0, 180.0, 80.0);

    /* 外枠の描画 */
    g_sel_scale(1);
    g_area_color(G_WHITE);
    g_line_color(G_BLACK);
    g_line_width(2);
    g_box(0.0, L, -1.0, 1.0, G_YES, G_YES);
    g_move(0.0, 0.0);
    g_plot(L, 0.0);

    /* X[i] を求める */
    for(i = 0; i < N; i++)
    {
        X[i] = i*dx;
    }

    /* Y[i] を求める */
    for(i = 0; i < N; i++)
    {
        Y[i] = sin(X[i]);
    }

    /* 折れ線の種類を設定 */
    g_line_color(G_RED);
    g_line_type(G_LINE_SOLID);
    g_line_width(2);

    /* 折れ線を描く */
    g_move(X[0], Y[0]);
    for(i = 1; i < N; i++)
    {
        g_plot(X[i], Y[i]);
    }

    g_sleep(G_STOP);
    g_term();

    return 0;
}

上記プログラムに、グラフの左端右端等の情報を書き入れたプログラムの出力を使って、 $N$ とグラフの関係を見てみましょう。
青色の破線で書かれたグラフが描きたいグラフ $y=\sin(x)$ )です（とはいっても、もちろんこれも $N=100$ の折れ線です）。
赤色の実線が描こうとしている折れ線です。

$N=5$ の場合は次のようになります（赤線）。
これでは、まったく $y=\sin(x)$ のグラフには見えません。

つづいて、 $N=10$ の場合（赤線）。曲線っぽくみえてきました。

$N=50$ の場合（赤線）。見た目はほぼ $y=\sin(x)$ のグラフに見えます。

課題1

上のサンプルプログラム（0421-1.c）を表示例のように表示するプログラムに変更せよ。
（グラフの左端等の表示と $N=100$ のグラフを重ねて表示するように変更）

ヒント：　Y[N] とは別に Y0[100] という配列を作って、それぞれを異なる線種で描く。また、文字列の描画（その２）で説明した手法を用いてタイトル部分を書くと良い。(注意：g_text 関数の最初の２つの引数（座標は）標準座標系での値になります.)

課題2

$y=3\cos(0.5x)$ のグラフを $[-2\pi, 4\pi]$ の範囲で描け。滑らかに見えるように、適度に $N$ を調整すること.　また、適度な大きさの「丸」を使った点描画でのグラフも描け.

y = sin(x + t) のグラフのアニメーション

すでに、アニメーションの作り方を知っているので簡単です。ただし、いま $t$ は時間として扱いますが、当然連続量として扱うことができないので、 $x$ 同様とびとびの値として扱います。先ほどの $x$ と同様に考え、 $t$ の範囲を [0, T] として、それを K-1 等分します。その小区間の幅を $dt$ として（つまり、 dt = T/(K-1)）、

Tk = dt*k

とすることにより、 Tk を止めるごとに Y[i] = sin(Xi + Tk) のグラフを描くことを繰り返せば、アニメーションとなります。次のサンプルプログラムを 0428-2.c として打ち込み実行してみてください。前々回のアニメーションプログラムと仕組みが同じであることに気づくことが重要です。（先のプログラムにはあったX[N]の配列は不要なのでなくしています.（先のプログラムでも実際にはX[N] の配列は不要でした））。

#include <glsc.h>
#include <stdio.h>
#include <math.h>

/* 定数の定義 */
#define N (50)
#define K (100)
#define PI (3.1415926)
#define L (2*PI)
#define T (10.0)

int main()
{
    int i, k;
    double Y[N], dx;
    double dt;

    /* dx, dt を求める */
    dx = L/(N - 1);
    dt = T/(K - 1);

    g_init("GRAPH", 200.0, 100.0);
    g_device(G_DISP);

    /* 仮想座標系の定義 */
    g_def_scale(1, 0.0, L, -1.0, 1.0, 10.0, 10.0, 180.0, 80.0);

    g_sel_scale(1);

    /* 毎時間繰り返し */
    for(k = 0; k < K; k++)
    {
        /* グラフ等を消去 */
        g_cls();

        /* 外枠の描画 （毎時刻グラフの消去してから描き直す） */
        g_area_color(G_WHITE);
        g_line_color(G_BLACK);
        g_line_width(2);
        g_box(0.0, L, -1.0, 1.0, G_YES, G_YES);
        g_move(0.0, 0.0);
        g_plot(L, 0.0);

        /* Y[i] を求める */
        for(i = 0; i < N; i++)
        {
            Y[i] = sin(i*dx + k*dt);
        }

        /* 折れ線の種類*/
        g_line_color(G_RED);
        g_line_type(G_LINE_SOLID);
        g_line_width(1);

        /* 描画の開始点 */
        g_move(0.0, Y[0]);

        /* 描画（曲線をつなぐ） */
        for(i = 0; i < N; i++)
        {
            g_plot(i*dx, Y[i]);
        }

        g_sleep(0.05);
    }

    g_sleep(G_STOP);
    g_term();

    return 0;
}

上記プログラムでは, i について Y[i] を一度全部求めてから, 再び i について Y[i] を描く, というように２度手間になってますが, このように計算しながら描く事も出来ます. また, 点描画もできます.

課題3

上のアニメーションのプログラムで $N$ を 5 にするとどうなるか？やってみよ。
（好きな方法で構いません.）

課題4

課題2のプログラムにグラフの左端等の情報を表示する部分を追加せよ。

課題5

$0.5(\sin (x+t) + \sin (x-t))$ のアニメーションを作成せよ。

課題6

$\sin (x+t) + \sin (x-t)$ のアニメーションを作成せよ。

簡単な常微分方程式とその解のグラフ化

今日は、計算数理Aの講義に出てきた次の常微分方程式を用いて演習を進めます。この方程式は、講義で詳しく説明されたと思いますが、例えば放射性同位体の崩壊を記述するモデル方程式として導出されます。講義で紹介があったように、例えば放射性炭素による年代測定にはこういった微分方程式が役立ちます。 $N(t)$ を時刻 $t$ における放射性同位体の原子数とすると下記のような微分方程式が得られることは講義で見たとおりです。

(P) $\begin{equation*} \left\{\begin{array}{lll} \frac{d}{dt} N(t) &=& - \lambda N, \\ N(0) &=& N_0 \end{array}\right. \end{equation*}$

さて、この微分方程式は簡単に解けて、解は次のようになることも講義で見ました。

(S) $\begin{equation*} N(t) = N_0 e ^{-\lambda t} \end{equation*}$

今日の目標は、解(S)をグラフとして表示するGLSCプログラムを作成することです。パラメータ $\lambda$ および $N_0$ を与えれば、(S)のグラフをかけるはずです。実は前回の内容で $y=\sin (x)$ のグラフを描いてもらいましたが、そこで注意したようにコンピュータを用いてこのようなグラフを描く場合には折れ線近似で描く必要があります。前回の $y=\sin (x)$ を描く方法を記述した部分にしたがって、(S)のグラフを描いてください。

課題7

GLSCを用いて、 $N(t) = N_0 e ^{-\lambda t}$ のグラフを描け。(0428-1.c, 0428-2.c 等を参考にせよ.)
ただし、 $N_0=1, \lambda=1, 0\le t \le 10$ とする。
例えば次のようになる。（この例では、0≦t≦10 を 100等分している。つまり時間刻み幅を dt と書くことにすると、 dt = 10/100 = 0.1 である。曲線を100本の直線のつながりである折れ線で近似するのである。）

課題8

$N(\tau) = 0.5$ となる $\tau$ を求めよ。（（S）を用いて具体的に求めよ）

答え： $\log 2$ (約 0.693)

課題9

課題7で作成したプログラムを用いて課題8のτを近似的に次のような手法で求めたい。そのようなプログラムを作成せよ。

方法

$N(t)$ の値を（S)を用いて飛び飛びの $t$ の値 $t_i = dt\ast i$ （ただし $dt$ は時間刻み幅）について求めるのであるが、(S)と課題６のグラフから明らかなように $N(t)$ は単調減少関数である。であるから、

$N(t_i) > 0.5 \ge N(t_{i+1})$

となる $i$ があるであろう。このとき、 $t_i$ を $\tau$ の近似値 $\tau^\ast$ とする。

$0\le t \le 10$ を100等分して上記方法にて $\tau$ の近似値を求め、課題７で求めた値との差（誤差と呼ぶ）

$err^\ast = |\tau - \tau^\ast|$

を求めよ。

注意：当然ながら刻み数を100から1000に増やせばより正確な値が求まる。

参考：

課題８、課題９（後で出てきます）の方法は、 $N(t)$ が具体的に与えられていない場合（例えば実測データなど）にも有効なので、実用的です。

課題10

課題9の方法ではあまり良い精度で近似値を求めることができない。
次の方法でより精度良く近似値をもとめることを考える。

$N(t_i) > 0.5 \ge N(t_{i+1})$

なる $i$ を求め、２点 $(t_i, N(t_i)), (t_{i+1}, N(t_{i+1}))$ を結ぶ直線を考える。
その直線と $N = 0.5$ の直線の交点の $t$ の値 $\tau^{\ast\ast}$ をもって、 $\tau$ の近似値をする。
（このような手法を線形補間と呼ぶ。）

$0\le t \le 10$ を100等分して上記方法にて $\tau$ の近似値を求め、課題２で求めた値との差

$err^{\ast\ast} = |\tau - \tau^{\ast\ast}|$

を求めよ。課題9の結果よりも良いだろうか？

注意：

課題9の方法でも刻み数を増やせばより正確な値が求まるのであるが、課題10の補間を使えば、あまり多くない刻み数でも正確な近似値が得られることを実感してほしい。

課題11

課題10において刻み数100で近似値を求めた場合の誤差と同程度に小さい誤差とするためには課題9において刻み数をどの程度にしなければならないか。実験してみよ。

常微分方程式の数値解法

次のような常微分方程式の初期値問題を、オイラー法と呼ばれる方法を用いて数値的に解き、GLSCを用いて結果をグラフにしてもらいます。次の問題を考えます。

(1) $\begin{equation*} \left\{\begin{array}{lll} \frac{d}{dt}y(t) &=& f(t, y(t)), \\ y(0) &=& a \end{array}\right. \end{equation*}$

求めたいのは、 $t>0$ の範囲の $y(t)$ であり、関数 $f(t, y(t))$ および定数 $a$ は与えられているとします。また、 $f(t, y(t))$ の値は $t, y(t)$ が与えられればいつでも計算できるものとします。（関数として定義するとよい。）

コンピューターでは数値を離散的に（とびとびに）しか扱うことができません。そこで、（１）式をコンピューターで扱えるようにする必要があります。（１）をコンピュータで解く方法には色々とありますが、この演習では代表的なオイラー法を用います。その前に、準備として差分商について説明します。

差分商

$f(x)$ を $x$ の関数としたとき、 $f(x)$ の $x=a$ における微分係数は

(2) $\begin{equation*} \lim_{h\rightarrow 0}\frac{1}{h}\{f(a+h)-f(a)\} \end{equation*}$

で定義されます。ここで、定義中の極限操作を取り払い、 $h$ を有限にとどめた

(3) $\begin{equation*} D=\frac{1}{h}\{f(a+h)-f(a)\} \end{equation*}$

を考えると、 $h$ を十分小さな正の実数ととれば、（３）式は（２）式の近似となっていると考えられます。
この $D$ のように、関数のいくつかの点における値の差を用いてその関数の微分係数を近似することを差分近似といい、（３）式の右辺のような量を差分商といいます。いまの場合は１階微分係数を近似する１階差分商です。

では、このような差分商を用いて（１）式を離散化してみましょう。まず、（１）式の第一式の微分係数を上記の差分商で置き換えてみます。

(4) $\begin{equation*} \frac{1}{h}\{y(t+h)-y(t)\}=f(t, y(t)) \end{equation*}$

ただし、 $h$ は正の数とします。（４）式と（１）式の第一式は異なる方程式なので、（１）式の第一式の解 $y(t)$ は一般に（４）式を満たしません。
そこで、混乱を防ぐために（４）式の $y(t)$ を $Y(t)$ と書き換えましょう。

(5) $\begin{equation*} \frac{1}{h}\{Y(t+h)-Y(t)\}=f(t, Y(t)) \end{equation*}$

（５）式のように差分を含む方程式を差分方程式といいます。

　次に、（１）式の $y$ を $Y$ に置き換えた初期条件

(6) $\begin{equation*} Y(t) = a \end{equation*}$

の下で（５）式を解くことを考えます。（５）式は、

(7) $\begin{equation*} Y(t+h)=Y(t)+hf(t, Y(t)) \end{equation*}$

と変形できるので、t = 0 を代入すると、

(8) $\begin{equation*} Y(h)=Y(0)+hf(0, Y(0))=a+hf(0,a) \end{equation*}$

となり、Y(0) から直ちに Y(h) が求まります。同様にして（７）式を繰り返し用いると、

(9) $\begin{equation*} \left\{ \begin{array}{lll} Y(2h) &=& Y(h)+hf(h, Y(h)), \\ Y(3h) &=& Y(2h)+hf(2h, Y(2h)), \\ Y(3h) &=& Y(2h)+hf(2h, Y(2h)), \\ &\vdots& \\ Y((j+1)h) &=& Y(jh)+hf(jh, Y(jh)), \\ &\vdots& \\ \end{array} \right. \end{equation*}$

というように、 $j=1, 2, 3, \dots$ とすると $Y((j+1)h)$ を $Y(jh)$ から計算できることがわかります。

$h$ をいったん決めると、 $t=jh$ 以外の時刻の $Y$ の値は（７）式からは求めることができません。
このように、差分方程式を解くと従属変数はとびとびの時刻で値が定まります。そのようなとびとびの時刻を格子点と呼びます。
$Y$ は格子点でのみ意味があるので、そのことを明示するために $Y(jh)$ を $Y_j$ と書き換え、 $T_j = jh$ とすると、（５）式と初期条件は、

(10) $\begin{equation*} \left\{ \begin{array}{lll} \frac{1}{h}\left(Y_{j+1}-Y_{j}\right) &=& f(t_j, Y_j), \\ Y_0 &=& a \end{array} \right. \end{equation*}$

となり、結局（１）式の常微分方程式の初期値問題が、 $Y_j$ に関する漸化式の問題（１０）式に置き換えられました。この方法をオイラー法と呼びます。

オイラー法のアルゴリズム

オイラー法のアルゴリズムを以下に擬似コードで表す。

解を求める時刻の上限をTとし、適当な自然数Nを定めてh=T/Nとする。

T, h, a, Y, t 及び関数 f は実数型(double)

N, j は整数型(int)

T, Nを設定

h = T/N

aを設定（入力する：printf関数＋scanf関数）

Y = a

t = 0でのグラフの描画（g_move関数）

j = 0, 1, 2, ..., N-1の順に（for文）

    t = j * h

    Y = Y + h * f(t, Y)

   一定時間毎にグラフを描画（g_plot関数)

以上をj=N-1まで繰り返す。

戻る