第5回(2009年5月27日)


今日の演習


GLSC で作成した静止画像のプリンターからの出力

 GLSC で作成した画像をプリンターから出力する方法です.

 (下の動画も参考にしてください.(TA 秋山氏 作))

1. プログラムを実行するターミナル(端末)と画像を取り込むためのターミナル(端末)を準備.

2. それぞれのターミナル上で、プログラムのあるフォルダ、画像を保存したいフォルダ(下の動画では Desktop フォルダ)に移動.

3. プログラムを実行(GLSC で作成された画像が現れる.)

4. 画像を取り込むためのターミナル上で

  ]$ import test.ps

 と入力して Enter キーを押す(ここではファイル名を「test」としましたが, 実際は何でもよいです. ただし 最後の「.ps」は必要です.)

5. すると「十字型」のマウスポインタが現れるので、それを実行の結果現れた window に移動して左クリック.
  (画像がプリンターから出力できるpostscript file として保存されます.)

6. ファイルマネージャーを開き, 保存した ps file を開き印刷.


 プリントしてみてはみ出してたりしてたら, 以下のようにして画像サイズ小さいファイルを作成します.

  6−1. ターミナル上で画像ファイル(この場合test.ps)のあるフォルダに移動

  6−2. ターミナル上で

   ]$ convert -resize 50% test.ps test1.ps

  と入力して Enter キーを押す. すると test.ps で表示されるものの 50% の大きさの画像ファイル test1.ps が作成されます.
  それをプリントしてください.



(参考) (ファイル名).ps とするとプリンターから出力するための Postscript ファイルが生成されますが, (ファイル名).jpg や (ファイル名).png とすると、それに応じて JPEG ファイル, PNG ファイルが生成されます.


常微分方程式の数値解法


今回は次のような常微分方程式の初期値問題を、オイラー法と呼ばれる方法を用いて数値的に解き、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)式に置き換えられました。この方法をオイラー法と呼びます。見やすいように(10)を書き直すと(11)になります.

   (11)


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

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

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

N, j は整数型(int)

Y[N+1], t[N+1] は実数型配列(float)


T, N, aを設定 (#define文で定数として定義すれば良い)

配列 Y[N+1], t[N+1] を用意する

h = T/N(TおよびN双方が整数型でないか注意)

Y[0] = a (初期値の設定)

t[0] = 0

j = 0,1,2,.......,N-1 の順に (for 文)

t[j+1] = (j+1)*h

Y[j+1] = Y[j] + h*f(t[j], Y[j])

以上を繰り返す

点(t[0], Y[0]), (t[1], Y[1]), ... , (t[N], Y[N]) を直線で結ぶことにより数値解をグラフとして描画 (g_move, g_plot 関数)

f(t,Y) は関数(上で復習してます)として定義するとよい.

計算数学演習での「関数の定義」の解説とサンプルプログラムはこちらにあります.


以下のアルゴリズム2では,数値解を求めながらグラフを描いています.(アルゴリズム1では,数値解を全て求めた後グラフを描いている.)


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

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

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

N, j は整数型(int)


T,N,aを設定

h = T/N

Y = a

t = 0 でのグラフの描画 (g_move 関数)

j = 0,1,2,.......,N-1 の順に (for 文)

t = j*h

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

一定間隔毎にグラフを描画 (g_plot 関数)

Y = Y_new ( Y の値を更新 )

以上を繰り返す

f(t,Y) は関数(上で復習してます)として定義するとよい.

計算数学演習での「関数の定義」の解説とサンプルプログラムはこちらにあります.


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


さて,実際にオイラー法をつかって常微分方程式を解いてみましょう.

今回も、計算数理Aの講義に出てきた次の常微分方程式を用いて演習を進めます。(ただし N を y と書き換えた.上記アルゴリズム等の分割数 N とまぎらわしいので.)

 (P)

前回は,この微分方程式を解析的に解き,得られた解をグラフ化してもらいましたが,今回は(P)をオイラー法を用いて数値計算で解いてもらい,得られた数値解をグラフ化してもらいます.


計算数学演習での「オイラー法の解説とサンプルプログラム」はこちらにあります. 参考にしてください. 以下の課題では,できればアルゴリズム2を用いて欲しいですが,分かりにくい人はアルゴリズム1を用いても良いです.


課題1

(P)をオイラー法で解き,得られた数値解をグラフとして表示するプログラムを作成せよ.ただし,

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

とする。例えば次のようになる。(この例では、0≦t≦10 を 100等分している。つまり時間刻み幅を h と書くことにすると、 h = 10.0/100 = 0.1 である。)(縦軸がN(t)となっているがy(t)と読み替えてください.)

先週,解析的に求めた解をグラフとして表示したものは以下であるが,殆ど違いがない事に注意.(課題2でみるが,微妙に異なっている)


課題2

前回解析的に求めた解(S)とオイラー法で求めた数値解を重ねて描くプログラムを作成せよ.例えば,次のようになる.この例では,解(S)を緑色で描き,数値解を赤色で描いている.(この例では、0≦t≦10 を 100等分している)

刻み幅を小さくすると両者の差は縮まる.(次の例では、0≦t≦10 を 300等分している)


レポート課題

課題2で作成したプログラムを rep0527.c というファイル名で report フォルダにおいてください.例えば次のような出力がえられます.0≦t≦10 を 100等分とすること.

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


課題3

次の常微分方程式の初期値問題をオイラー法を用いて解いて、t と y の関係を GLSC でグラフ化せよ.

 (1) dy/dt = -y + sin(t), y(0) = 1  (時刻 0≦ t ≦20 の範囲を20000等分割、100ステップ毎に値を表示)

 (2) dy/dt = y(5-y), y(0) = 1  (時刻 0≦ t ≦20 の範囲を20000等分割、50ステップ毎に値を表示)

 (3) dy/dt = y - 2y3, y(0) = 0.1 もしくは y(0) = -0.1 (時刻 0≦ t ≦20 の範囲を20000等分割、100ステップ毎に値を表示)


戻る