計算数理A演習第3・4回

今週の目標

matplotlibによるグラフ描画

先週の演習では基本的なsin関数の描画を行い,まず描画するということになれてもらいました。
今週はmatplotlibで描画する際に必要になる基礎的な知識について学んでいきます。
少しハードルが高いですが,一番大事な情報はmatplotlibのホームページ(matplotlib.org)に書いてありますので,より高度なことをしたくなったら見てみてください。
また,実例を見てみることに勝る教科書はありませんので,どんなグラフがかけるのか興味があれば見てみてください。(https://matplotlib.org/gallery/index.html)サンプルコードもついているので,そのまま実行すれば手元で遊ぶことが出来ます。

matplotlibの構造

それでは,まず初めにmatplotlibの描画システムを理解するための基本構造を学びます。
下に示すのは公式ガイドに載っているmatplotlibの描画システムを説明するためのグラフです。

引用:https://matplotlib.org/1.5.1/faq/usage_faq.html

実際にはより複雑なシステムですが,ここでは最も頻度が高く重要な要素について学びます。
本当は全てオブジェクトと呼ばれていますが,本演習ではオブジェクト指向という考え方について触れませんので,適当な言葉で説明します。
また,この部分については対話的にプログラムを実行すると具体的にわかりやすいので,エディタやjupyter notebookではなく,ipython自体を利用しますので,以下のコマンドをターミナルで実行してください。

ipython --pylab

すると次の様な画面になり,ipythonが起動します。


まず,一番外側にあるFigureで,これをプログラムの中で作ることで皆さんがグラフを描画するためのキャンバスが作られます。
早速,ターミナル上で以下のようにコマンドを入力してください。

fig = plt.figure()

すると以下のような画面になったはずです。

位置は調整していますが,上手の右にあるような何も描かれていない真っ白なウィンドウが現れたと思います。これがまさにキャンバスを作るということであり,Figureの実体になります。
本物のキャンバスは布と木で作られますが,それらがあったとして,いきなりキャンバスをもし作れと言われたら皆さんすぐに作れるでしょうか?絵を描くためのキャンバスには大小様々なものがあり,もし作れと言われたら当然サイズを指定する必要があります。
このFigureという構造も同様です。もちろん何の指示を与えなくてもデフォルトの大きさは与えられていますが。
ここでは示しませんが手元で,次のコマンドを実行してみてください。先のコマンドに比べるとかなり大きなウィンドウが現れたと思います。

fig = plt.figure(figsize=(10,8))

次に,その内側にある構造がAxesについてです。これは,データをグラフにするための座標系に当たります。早速以下のコマンドを実行してみてください。

ax = fig.add_subplot(111)

見ての通り真っ白だったウィンドウの中に四角い枠と縦軸,横軸に[0,1]区間の目盛が加えられました。準備されたFigureに座標系を与えることで,数値を図示することが出来るわけです。実はFigureというキャンバスの中には複数の座標系を用意したり,他の字を描くための構造を書くことが出来ます。

最後に緑色の円で囲まれているのが,axisという構造です。これは,Axesのそれぞれの軸を示しています。AxesはAxisの複数形で,先にも座標系と書きましたが,二次元か三次元のグラフしか描画しないmatplotlibでは,2つか3つの軸が必要になるので座標系に当たる構造がAxesとなるわけです。
グラフを描く際には,軸の幅や範囲など様々な情報を指定することが必要になります。そのときに最もよく利用して制御する必要がある構造が,このAxisになるわけです。

次にサンプルのための三角関数の数値を準備して,グラフを描くために以下のコマンドを入力してください。

x = np.linspace(0, 2.0*np.pi, 101)
y = np.sin(x)
ax.plot(x, y)

準備されたFigureに座標系を与えることで,数値を図示することが出来るわけです。実はFigureというキャンバスの中には複数の座標系を用意したり,他の字を描くための構造を書くことが出来ます。

最後に緑色の円で囲まれているのが,axisという構造です。これは,Axesのそれぞれの軸を示しています。AxesはAxisの複数形で,先にも座標系と書きましたが,二次元か三次元のグラフしか描画しないmatplotlibでは,2つか3つの軸が必要になるので座標系に当たる構造がAxesとなるわけです。
グラフを描く際には,軸の幅や範囲など様々な情報を指定することが必要になります。そのときに最もよく利用して制御する必要がある構造が,このAxisになるわけです。

具体的な描画方法

簡易的なグラフ描画(pyplot方式)

ここでは簡易的にグラフを描画する方法を学びます。
基本的な描き方は,先週学んだsin関数の描画プログラムが基本になります。
このやり方でも十二分に有用なグラフを描くことが出来ますが,少し高度なことをしようとするとすぐに次の方式を部分的に用いることになります。
そのため,本演習では次のObject-Oriented API(オブジェクト指向API)による方法を基本として進めていき,ここでは簡単に必要最低限のグラフを示す方法を学びます。
早速次のコマンドを実行してみましょう。

%matplotlib inline
import matplotlib.pyplot as plt
import numpy as np

N = 100
dx = 2.0*np.pi / N

x = np.zeros(N)
y = np.zeros(N)
for i in range(N):
    x[i] = i*dx
    y[i] = np.sin(x[i])

plt.plot(x, y, label="sin(x)")
plt.xlabel("x", fontsize=20)
plt.ylabel("y", fontsize=20)
plt.legend(loc="best", fontsize=20)

先週のグラフに比べて,情報が増えたと思います。
この方法で行えば,自然科学の実験やデータ解析で必要な必要最低限の情報を載せることが出来ます。また,plt.figureとplt.show()というコマンドを消しましたが,問題なくグラフが出てきたと思います。
これはplt.plotという関数が,呼び出してグラフを描けば,細かいことを気にしない代わりに基本設定のまま,空気を読んでグラフを描いてくれる関数だからです。
これによって,データをすぐに可視化してくれるわけです。

課題1
上のサンプルコードを書き換え,媒介変数表示(x(θ), y(θ)) = (cosθ, sinθ),0<θ<2πとして,円を描くプログラムをかけ。

高度な描画方式(Object-Oriented API・オブジェクト指向API)

このやり方は,先のmatplotlibの構造を説明を行うときに用いた方法で,このやり方がmatplotlibの基本的なグラフ描画方法だと考えてください。
大まかな流れは,先ほど説明したようにFigureを作り,座標系となるAxesを準備して,データを描くことになります。

fig = plt.figure(figsize=(10,8))
ax = fig.add_subplot(111)

x = np.linspace(0, 2.0*np.pi, 101)
y = np.sin(x)

ax.plot(x, y, label="sin(x)")
ax.set_xlabel("x", fontsize=20)
ax.set_ylabel("y", fontsize=20)
ax.legend(loc="best", fontsize=20)
fig.show()

線の色・線のパターンの変更

matplotlibではデフォルトのグラフの色として,青色が選択され,線で結ばれてグラフが描画されます。
ここでは,任意の色に変更したり,点線にしたりしてみます。
変更方法はとても簡単で,上のサンプルコードのax.plotに以下のように追加するだけで変更できます。

ax.plot(x, y, "r--", label="sin(x)")

この例にあるようにしてもう一度グラフを描くと,赤い点線でグラフが描かれたと思います。matplotlibのplot関数では,データとなる二つの配列を与えた後に,文字列を書き込むことでグラフの色と線の種類を選択することが出来ます。
色は次に示す省略した文字で指定し,線の種類は記号の組み合わせで指定します。
基本的には,”[色][線の種類]”という構造になっています。
一文字のアルファベットで指定できる色には,

‘b’(青), ‘g’(緑), ‘r’(赤), ‘c’(シアン), ‘m’(マゼンダ), ‘y’(黄), ‘k’(黒), ‘w’(白)

があります。

次に線の種類を変更するには,一文字のアルファベットのあとに記号の組み合わせで指定できるのですが,簡単に指定できるものが4種類があります。
https://matplotlib.org/gallery/lines_bars_and_markers/line_styles_reference.html

‘-‘, ‘–‘, ‘-.’, ‘:’

があります。それぞれ試してみて線の種類を見てみてください。

色も種類もより高度な方法を用いれば細かな設定が可能ですが,ここではあまり高度な内容は必要ないので,興味がある人はマニュアルページを見てみてください。

Specifying Colors (https://matplotlib.org/tutorials/colors/colors.html#sphx-glr-tutorials-colors-colors-py)
Linestyles (https://matplotlib.org/gallery/lines_bars_and_markers/linestyles.html)

マーカーのプロット

線の種類はデータ間をつないで可視化しました。ここでは,そのデータ点自体を描画する方法を学びます。数値計算ではあまり表だって使うことはないかもしれませんが。

fig = plt.figure(figsize=(10,8))
ax = fig.add_subplot(111)

x = np.linspace(0, 2.0*np.pi, 101)
y = np.sin(x)

ax.plot(x, y, "o", label="sin(x)")
ax.set_xlabel("x", fontsize=20)
ax.set_ylabel("y", fontsize=20)
ax.legend(loc="best", fontsize=20)
fig.show()

サンプルコードを実行すると先ほどまでつながれていた線が消えて,データとしてxとyにあるデータ自体が点として描かれています。x軸方向は等間隔でもy軸のデータが山の頂点に近づくにつれて詰まっているのがよくわかります。
ここでは,プロットするタイプとして”o”(アルファベットのオー)を指定しました。線の種類と同様に他にもプロットの種類があります。
種類がかなり多いので,リンク先のmatplotlibのマニュアルページを参照してください。(https://matplotlib.org/gallery/lines_bars_and_markers/marker_reference.html

また,これはデータ点自体のある場所にプロットしているだけなので,データ間をつなぐ線は独立に描くことが出来ます。

fig = plt.figure(figsize=(10,8))
ax = fig.add_subplot(111)

x = np.linspace(0, 2.0*np.pi, 101)
y = np.sin(x)

ax.plot(x, y, "o-", label="sin(x)")
ax.set_xlabel("x", fontsize=20)
ax.set_ylabel("y", fontsize=20)
ax.legend(loc="best", fontsize=20)
fig.show()

描画範囲の指定

それでは次にグラフを描く範囲の指定方法を学びます。
特に指定しなければ,サンプルコードにあるようにデータを全て描ける範囲を自動で取ってくれます。
ここでは,あえて描画する範囲を絞って見ます。

fig = plt.figure(figsize=(10,8))
ax = fig.add_subplot(111)

x = np.linspace(0, 2.0*np.pi, 101)
y = np.sin(x)

ax.plot(x, y, label="sin(x)")
ax.set_xlim([0,2.0*np.pi])
ax.set_ylim([-1.0,1.0])
ax.set_xlabel("x", fontsize=20)
ax.set_ylabel("y", fontsize=20)
ax.legend(loc="best", fontsize=20)
fig.show()

このサンプルコードでは,描画するサイン関数がぴったりはまるように範囲を指定しました。

ax.set_xlim([0,2.0*np.pi])
ax.set_ylim([-1.0,1.0])

上ののコードがそれぞれax.set_xlim,ax.set_ylimが,x軸,y軸の範囲を制限しています。x軸y軸どちらを指定するかは見ての通りxlimかylimかの違いで使い方は全く同じです。
関数の引数に与えるリスト[a,b]でaが軸の下限値,bが上限値を表しています。
サンプルコードでは,元のコードにあった余白がなくなっていることがわかると思います。

課題2
先のサンプルコードにおいて,(0,π)近傍を拡大して直線的に見える様子を確認せよ。

円を描く・図形を描く

図形を描くことはいろいろな場面で必要になるかと思いますが,matplotlibには基本的な図形を描けるようにモジュールが用意されています。これを使えば色々と図形を描くことが簡単にできるのですが,簡単な円を描くだけであればplot関数を用いて出来ますので,その方法を学びます。また多角形の枠の作り方も学びます。

fig = plt.figure(figsize=(10,8))
ax = fig.add_subplot(111)

cx = 0.0
cy = 0.0

ax.plot(cx, cy, 'o', markersize=200.0)
ax.set_xlim([-0.5, 0.5])
ax.set_ylim([-0.5, 0.5])

ax.set_xlabel("x", fontsize=20)
ax.set_ylabel("y")
ax.legend(loc="best")
#fig.show()

これまでのplot関数に与えるデータはそれぞれいくつかの値を持つ配列だけでしたが,上のサンプルコードにあるようにそれぞれ一つのデータでもグラフを描くことができます。この場合データ間をつなぐ線の種類を指定しても一点しかデータがないのでつなぐことが出来ず,何も描かれません。そこで,データ点自体を描画するマーカーを指定してその座標点のデータを描画しています。ただし,折れ線グラフのようにデータをプロットしながら線を結ぶことを想定しているので,マーカーは線の太さよりも少し太いくらいです。なので,何かしら円を描くつもりでマーカー機能を使う場合には,サイズを大きくする必要があります。

課題3
ビリヤードの球の初期配置のように9個の円で三角形を作れ。(向きは好きにしてください。)

多角形は線をつなげばいいので,四角形は以下のようにすれば出来ます。

fig = plt.figure(figsize=(10,8))
ax = fig.add_subplot(111)

x = np.zeros(5)
y = np.zeros(5)

for i in range(5):
    x[i] = np.cos((i*2.0*np.pi/4.0)+np.pi/4.0)
    y[i] = np.sin((i*2.0*np.pi/4.0)+np.pi/4.0)

ax.plot(x, y, "-", label="sin(x)")
ax.set_xlabel("x", fontsize=20)
ax.set_ylabel("y", fontsize=20)

ax.axis('scaled')
ax.set_aspect('equal')
ax.legend(loc="best", fontsize=20)
fig.show()

matplotlibでは以下のようにいくつかの図形を簡単に描けるようになっています。
描くためのモジュールがpatchesで以下のサンプルコードにあるように必要な場合にインポートしてください。

import matplotlib.patches as patches

fig = plt.figure(figsize=(10,8))
ax = fig.add_subplot(111)

# 円を作る
c = patches.Circle(xy=(0, 0), radius=0.5, fc='g', ec='r')
# 楕円を作る
e = patches.Ellipse(xy=(-0.25, 0), width=0.5, height=0.25, angle=10, fc='b', ec='y')
# 長方形を作る
r = patches.Rectangle(xy=(0.3, 0), width=0.25, height=0.5, angle=-10, fc='#00FFFF', ec='#00FFFF', fill=False)

# 作った図形を座標系の中へ描く
ax.add_patch(c)
ax.add_patch(e)
ax.add_patch(r)
# 二点の座標データを直接指定して直線を引く
ax.plot([0, 0.8], [0, -0.1])

ax.axis('scaled')
ax.set_aspect('equal')

課題4
多角形として,六角形と八角形を描け。可能であればn角形を描けるプログラムを作れ。

任意の場所に文字を書く

次にグラフをよりよく説明するために必要な文章をグラフ内に書く方法を学びます。
ax.plot関数がデータを書き出すための関数でしたが,文字を描くにはax.text関数を使います。基本的なルールは,

ax.text(左下のx座標, 左下のy座標, [文字列])

のように文字を書き出す地点を指定して,次に書きたい文字列を書きます。

import matplotlib
import matplotlib.pyplot as plt

fig = plt.figure()
ax = fig.add_subplot(111)

ax.text(3, 9, '1 Simplest way')
ax.plot(3, 9, 'o')

ax.text(3, 8, r'2 an equation: $E=mc^2$', fontsize=15)
ax.plot(3, 8, 'o')

ax.text(2, 6, '3 boxed italics text in data coords', style='italic')
ax.plot(2, 6, 'o')

ax.text(7, 3, r'4 unicode: Institut für Festkörperphysik',
        verticalalignment='bottom', horizontalalignment='right')
ax.plot(7, 3, 'o')

ax.text(0.55, 0.01, '5 colored text in axes coords',
        transform=ax.transAxes,
        color='green', fontsize=15)
ax.plot(0.55, 0.01, 'o')

ax.axis([0, 10, 0, 10])

plt.show()

当然axは座標系を作るために呼び出したので,ax.textにまず初めに書き込む座標はこの座標系上での値になります。set_xlimなどを使い描画する範囲を指定するためには,text関数に入力する座標がその範囲に収まっている必要があります。

サンプルコードの一番最初にあるものが最も基本的な使い方です。

二つ目の書き方は,LaTeXコマンドを使う方法です。文字列は””か”で包めば文字列と認識されるのですが,その前にrとつけておくとLaTeXコマンドが利用できます。加えて,ラベルの指定と同様にfontsizeで文字の大きさを指定できます。

三つ目の書き方は,文字のフォントスタイルをstyle=’italic’とすることで斜体へ変更しています。
加えて,これまで座標位置が必ず左下だったものが,右下に変わっています。
これはverticalalignment=’bottom’, horizontalalignment=’right’によるもので,verticalalignmentには{‘center’, ‘top’, ‘bottom’, ‘baseline’, ‘center_baseline’}が指定でき,horizontalalignmentには{‘center’, ‘right’, ‘left’}が指定できます。

最後の書き方は,これまでのものと少し違います。上記の三つのtextでは,初めの座標が左下の位置と対応していることを見るためにax.plotで点を同じ座標に打っていて,ちゃんと左角にあることが見て取れていると思います。しかし,最後のものだけ一致していません。
これは,後ろの方に書いてあるtransform=ax.transAxesというコマンドのためです。
これを行うとaxによって作られた座標系の中の位置ではなく,figure全体に[0,1]×[0,1]の座標系があるような,絶対的な位置により指定されることになります。

また,本演習を通して数値を表せる必要があるので,数値を文字列として書き出す方法を述べておきます。
いくつかやり方はあるのですが,format形式が私の知る限り一番便利なのでその方法を示します。
format形式では,以下のようにint型やfloat型の変数を文字列変換します。これを行うと表示だけでなく,ファイル名として保存するときにもパラメータの数字が入ったファイル名を連続して呼び出すことが出来ます。

import matplotlib
import matplotlib.pyplot as plt

fig = plt.figure()
ax = fig.add_subplot(111)

a = 2.3
b = 5

print(f'Format Style: a = {a}, b = {b}')
ax.text(3, 2, f'Format Style: a = {a}, b = {b}')

a = 2
b = 1
ax.text(a, b, f'Format Style: a = {a}, b = {b}')

ax.axis([0, 5, 0, 5])

plt.show()

これまでと違うのは文字列の”もしくは””の前にfを付け加えることと,文字列に変更したい変数を{}で囲うことです。
print関数で文字列を書き出す場合でも,ax.text関数でグラフに書き出す場合でも文字列として数値が出ているのがわかると思います。
format形式だと最後の例のように変数の値を変更すると数値が変更されていることがわかると思います。

複数のグラフを一つの画面に描画する

複数のデータを一つのグラフ上にプロットするには以下のように二つの座標系を準備します。
fig.add_subplot(111)で一つの座標系をFigure全体に描画していました。ここでは書き込む座標系を二つにします。そのためにadd_subplotの引数について,以下のサンプルコードを読んでもらった後説明します。

fig = plt.figure(figsize=(10,8))
ax1 = fig.add_subplot(211)
ax2 = fig.add_subplot(212)

x = np.linspace(0, 2.0*np.pi, 101)
y = np.sin(x)
y2 = np.cos(x)

ax1.plot(x, y, label="sin(x)")
ax1.plot(x, y2, '--', label="cos(x)")
ax1.set_xlabel("x", fontsize=20)
ax1.set_ylabel("y", fontsize=20)
ax1.legend(loc="best", fontsize=20)

y2 = np.cos(x)
ax2.plot(x, y2, '--', label="cos(x)")
ax2.set_xlabel("x", fontsize=20)
ax2.set_ylabel("y", fontsize=20)
ax2.legend(loc="best", fontsize=20)

fig.show()

fig.add_subplot(211)にある謎の211という数字の組み合わせですが,実はこの部分は以下のように書き換えることが出来ます。

fig.add_subplot(2,1,1)

意味としては(行数,列数,場所)となりつながった三桁の数字も同じ意味です。
行数、列数によってfigureを縦横に必要な分だけ分割して,その中でどの領域に座標を設定するかを最後の数字で指定しています。分割された領域は行優先で数が割り振られています。
例えば,fig.add_subplot(3,4,10)という指示を出すと、ちょうど10の当りにだけ座標系が設定されます。
注意すべきは,分割した全ての領域を利用する必要はなく,あくまで作りたい座標系のサイズと位置を決めるだけです。なので、適当に与えるとグラフが重なってしまうので注意してください。

1
234
5678
9101112

課題5
上下にグラフを分割して、上に任意の二次以上の多項式関数のグラフを描き,下段にその導関数のグラフを描け。定義域は0<x<3とする。

matplotlibでのアニメーション

アニメーションの描画を行う方法には,先の二つの方式に対応する大きく分けて二つの方式が存在します。

jupyter notebookでアニメーションを見るための設定

これまでグラフを描画してjupyter上で見るために,初めに

%matplotlib inline

と書いていました。アニメーションを見るためには,この方法では見られないため,以下のように変更してください。

%matplotlib tk

この部分を変更して実行しても,すでに先のコマンドを実行してjupyter notebookを使っている場合には,変更がききません。
新しく二つ目のコマンドを有効にするためには,いったんnotebookを閉じるか,ツールバーのkernelを選択してRestartをクリックして再起動する必要があります。

データを直接変更するアニメーション

初めに行う方法は,一度描いたグラフのデータを更新することでアニメーションを作る方法です。題材として進行波を描くプログラムを用意したので実行してみてください。

import numpy as np
import matplotlib.pyplot as plt

fig = plt.figure()
ax = fig.add_subplot(111)
x = np.arange(0, 10, 0.1)
y = np.sin(x)
dt = 0.05
v = 1.0
lines, = ax.plot(x, y, "r")

for i in range(50):
    t = dt * i
    y = np.sin(x - v*t)
    lines.set_data(x,y)
    plt.pause(0.01)

plt.show()

どこでグラフのデータを更新しているかというと,lines.set_data(x,y)の部分になります。
グラフ自体の実態としてlinesがまず,lines, = ax.plot(x, y, “r”)によって作られています。ax.plot関数はグラフを描くと同時に,関数として描かれたグラフに関する様々な情報をreturnしてくれます。ここでは詳細に触れませんが,その中でlinesとして受け取ったものをこの例題の様に制御するとグラフが変わることだけ覚えておいてください。

このデータを更新するだけではグラフの描画は変更されないので,その次の行にあるplt.pause(0.01)によりグラフは更新されます。0.01は更新された後にそのまま何秒待つかを指定しています。変えてみてアニメーションの進む速度がどう変わるのか確かめてみてください。

パラパラ漫画的手法

次は一度描かれたグラフのデータを変更するのではなく,少しずつずらしたグラフを用意して,そのたくさんのグラフを一つの関数に渡すことでまさに「パラパラ漫画」の要領でアニメーションを作る方法です。

import numpy as np
import matplotlib.pyplot as plt
import matplotlib.animation as animation

fig = plt.figure()
ax = fig.add_subplot(111)
x = np.arange(0, 10, 0.1)
y = np.sin(x)
dt = 0.05
v = 1.0

ims = []
for i in range(50):
    t = dt * i
    y = np.sin(x - v*t)
    im = ax.plot(x, y, "r")
    ims.append(im)

ani = animation.ArtistAnimation(fig, ims)
plt.show()

先ほどのやり方では,ax.plot関数はfor文の外側で一度使っただけですが,今回はfor文の中で毎回呼び出します。そして,そのグラフに関する情報をimとして受け取ります。先ほどlinesを作るときにlinesの後ろにカンマがついていたのですが,今回はついていないことに注意してください。
そして,for文の直前に書いてあるims = []というコードは,そのパラパラ漫画を格納するためのリストとなり,ims.append(im)によって実際に格納していき,for文が終了するとimsの中に50枚の絵が用意されたと思ってください。
最後に用意されたimsとアニメーションを書き出したいキャンバスであるfigを引数として持つanimation.ArtistAnimation関数に二つの変数を渡すことでアニメーションを作ることが出来ます。

変化を定義する関数を与える方法

二番目の方法にかなり近いやり方ですが,少し違う方法を紹介します。

import numpy as np
import matplotlib.pyplot as plt
import matplotlib.animation as animation

fig = plt.figure()
ax = fig.add_subplot(111)

x = np.arange(0, 2*np.pi, 0.01)
line, = ax.plot(x, np.sin(x))

def init():  # only required for blitting to give a clean slate.
    line.set_ydata([np.nan] * len(x))
    return line,

def animate(i):
    line.set_ydata(np.sin(x + i / 100))  # update the data.
    return line,

ani = animation.FuncAnimation(fig, animate, init_func=init, interval=2, blit=True, save_count=50)

plt.show()

この方法ではfor文を直接使わずに,animation.FuncAnimation関数にfor文で使われるコードをanimate関数として与えることでアニメーションを描いています。
この方法は余裕のある人が勉強して見てください。

どのやり方も一長一短なのですが,たぶん一番上の方法が楽かと思います。
一応お勧めのやり方は2番目か3番目です。

練習
円がx軸上を0から10まで移動するアニメーションを作れ。

課題6
いずれかの方法を用いて,中心座標に円を描きその周りを円が公転するアニメーションを作成せよ。

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


今日は、計算数理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)をグラフとして表示するプログラムを作成することです。パラメータ\lambdaおよびN_0を与えれば、(S)のグラフをかけるはずです。
実は前回の内容でy=\sin (x)のグラフを描いてもらいましたが、前回のy=\sin (x)を描く方法を記述した部分にしたがって、(S)のグラフを描いてください。

課題7
N(t) = N_0 e ^{-\lambda t}のグラフを描け。ただし、N_0=1, \lambda=1,  0\le t \le 10とする。

課題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)と課題6のグラフから明らかなように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の近似値を求め、課題8で求めた値との差(誤差と呼ぶ)

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

を求めよ。

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

参考:

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

課題10

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

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

なるiを求め、2点(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の近似値を求め、課題8で求めた値との差

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

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

注意:

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


課題11

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