第4回(2009年5月20日)


今日の演習


report と work フォルダの作成

 計算数学演習を受講しなかった人は、ファイルマネージャーを用いて、 report および work という名前のフォルダをホームディレクトリ(ホームフォルダー)に作成してください。 

 以前、計算数学演習を受けていた人は、既に work および report フォルダがあると思います。それらをそれぞれ、work_old および report_old と名前を変えてください。その後、新規に report および work フォルダを作成してください。

 以下、より詳細な説明を書いておきます。よくわからない人は、我々に質問してください。

 端末室のパソコンにログインした直後には、ファイルマネージャー(Gnome Midnight Commander という名前だそうです)が起動しているはずです。これを用いて、フォルダーをつくり、これまで作ってきたファイルをフォルダー内にまとめましょう。ファイルマネージャーの使い方の詳細はマニュアルにまかせます。work というフォルダーを作成し、先々週作成した0422-1.c などのファイルを移動してください。また、自動回収によるレポート提出の為のフォルダー(reportという名前のフォルダー)を作成してください。

 以前、計算数学演習を受けていた人は、既に work および report フォルダがあると思います。それらをそれぞれ、work_old および report_old と名前を変えましょう。その後、新規に report および work フォルダを作成してください。

 実際に work というフォルダーを作成してもらいました。そして、そのwork フォルダーにこれまで作成したソースファイル(0422-1.c)とそれをコンパイルしたもの(0422-1)を移動してもらいました。こんなことをして、ターミナルでの作業(コンパイルおよび実行)には影響ないのでしょうか?(当然影響があります。)

 ターミナルを開きます。先週同様、

cglsc  0422-1.c

としてみましょう。フォルダー作成等の作業がうまくいっていれば、エラーがでます。なぜなら、コンパイルするべきファイル 0421-1.c を別の場所である work フォルダに移動してしまったからです。

ホームディレクトリ(ホームフォルダー)

 ターミナルを起動した時点では、自分はホームディレクトリと呼ばれる位置にいます(ディレクトリとフォルダーは同じものだと思ってください)。それは、起動時のファイルマネージャーも同様で、起動時にはホームディレクトリの内容を表示しています。先ほどの作業で、ホームディレクトリにあったファイルを workフォルダーにファイルを移動してしまったため、ターミナルを起動した後、

cglsc  0422-1.c

としても、0422-1.c というファイルが見つからないのでエラーとなったわけです(エラーメッセージもファイルが見つからないといったものだったはずです)。正しくコンパイルするためには、work フォルダー内に自分が移動する必要があります。それには、ターミナルで次のコマンドを実行します。

cd  work

この新たに出てきたコマンド cd は Change Directory の略で、文字通り自分の位置を変更します。上記コマンドで、work フォルダーの中に移動できたので、改めて

cglsc 0422-1.c

としてみましょう。今度はうまくいくはずです。ついでに、コンパイルされたプログラムが実行できるかどうかも確かめておきましょう。

./0422-1

正しく実行できたと思います。

 ところで、ターミナル上で、そのフォルダー内のファイルの一覧を見たい場合もあるでしょう。次のコマンドで見ることができます(ls は LiSt の略です)。

ls  -l

ファイル名と、そのサイズ変更日時などが一覧として表示されます。hello.c や、そのコンパイルされた結果の hello もでていると思います。

 さて、今 work フォルダーの中にいますが、先のホームディレクトリに戻るにはどうするのか?次のコマンドで、ホームディレクトリに戻れます。

cd

当面はこれだけのことができれば十分でしょう。実際ターミナル上でもフォルダーを作成することができるのですが、この演習ではフォルダーの作成ファイルの移動ファイルのコピーはファイルマネージャーで行うことにします。(Windows や Macintosh を使いなれた人は何も問題ないでしょう)

 ターミナルで使えるフォルダーやファイル操作のコマンドのより詳しい情報はこちらにあります。

 今後、この演習で作成するファイルは全て work フォルダーに保存するようにしてください。ただし、自動回収によるレポート提出をする場合には、対象となるファイルを report フォルダーにコピーしてもらう必要があります。(計算数学演習と同じ形式です。)


y = sin(x) のグラフの描画(前回の復習)

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

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

Xi = i*dx

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

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

をそれぞれ求め、(X[0], Y[0]), (X[1],Y[1]), ...., (X[N-1], Y[N-1]) を結ぶ折れ線を描くことになります。これをプログラムにすると次のようになります。(上記の説明にあわせるためにかなり無駄なことをしています。)ファイル名を 0513-2.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)

main()
{
  int i;
  float 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();
}

--------------------

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

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

 

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

 

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

 


課題2

 上のサンプルプログラム(0513-2.c, もしくは 0513-2.c)を表示例のように表示するプログラムに変更せよ。(グラフの左端等の表示とN=100のグラフを重ねて表示するように変更)

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


課題3

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



今回 (5/20) はここから




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) のグラフを描くことを繰り返せば、アニメーションとなります。次のサンプルプログラムを 0520-1.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)

main()
{
  int i, k;
  float Y[N], dx; 
  float 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();
}

--------------------

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


課題4

 課題4のプログラムにグラフの横軸、縦軸等の情報を表示する部分を追加せよ。 


課題5

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


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


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

 (P)

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

(S)

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


課題6

GLSCを用いて、

のグラフを描け。(0513-2.c, 0520-1.c 等を参考にせよ.) ただし、

N0=1, λ=1, 0≦t≦10

とする。例えば次のようになる。(この例では、0≦t≦10 を 100等分している。つまり時間刻み幅を dt と書くことにすると、 dt = 10/100 = 0.1 である。曲線を100本の直線のつながりである折れ線で近似するのである。)


課題7

N(τ) = 0.5 となる τ を求めよ。((S)を用いて具体的に求めよ)

答え: log2    (約 0.693)


課題8

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

方法

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

N(ti) > 0.5 ≧ N(ti+1

となる i があるであろう。このとき、 ti を τ の近似値 τ* とする。

0≦t≦10 を 100等分して上記方法にてτの近似値を求め、課題7で求めた値との差(誤差と呼ぶ)

err* = |τ- τ*|

を求めよ。

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


参考:

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


レポート課題

GLSCを用いて、

のグラフを描け。ただし、

N0=1, λ=1, 0≦t≦10

とする。また半減値τも表示する事。次に作成したプログラムの実行例を示す。(この例では、 0≦t≦10 を 1000等分している。つまり時間刻み幅を dt と書くことにすると、 dt = 10/1000 = 0.01 である。曲線を1000本の直線のつながりである折れ線で近似するのである。)この例の出力に近くなるよう作成すること。

 半減値の表示については、例えば課題7の結果を用いて上記例と似たグラフを表示するプログラムでも構いませんが、例えば課題8または課題9(この後にあります.)の内容を組み込んで、τの近似値を実際に数値的に計算して表示するプログラムの方を、より高く評価しようと思います。

ファイル名は rep0520.c として、report フォルダ内におくこと。なお、ファイルの先頭にコメント文として、自分の学籍番号と名前を記入し、プログラムの途中にもコメント文として処理の内容を書き込むこと

自動回収日時:2009年5月26日(火)午後4時


課題9

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

N(ti) > 0.5 ≧ N(ti+1

なる i を求め、2点 (ti, N(ti)), (ti+1, N(ti+1)) を結ぶ直線を考える。その直線と N = 0.5 の直線の交点の t の値 τ** をもって、τ の近似値をする。(このような手法を線形補間と呼ぶ。)

0≦t≦10 を 100等分して上記方法にて τ の近似値を求め、課題2で求めた値との差

err** = |τ- τ**|

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

注意:

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


課題10

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


常微分方程式の数値解法


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

           (1)

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

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


差分商

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

         (2)

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

         (3)

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

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

         (4)

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

         (5)

(5)式のように差分を含む方程式を差分方程式といいます。

 次に、(1)式の y を Y に置き換えた初期条件

         (6)

の下で(5)式を解くことを考えます。(5)式は、

         (7)

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

         (8)

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

         (9)

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

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

         (10)

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


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

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

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

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 関数)

以上を繰り返す

 


戻る