コンテンツ

第３・４回（第２週）

今週の目標

グラフを描画する
簡単な常微分方程式を解いて描画する
アニメーションを描画する

【重要】グラフを別ウィンドウで描画するためのSpyderの設定

特に何も設定していない状態でグラフなどを描画すると、iPythonコンソールの中に描画されます。静止画の場合はこれで問題ないのですが、アニメーションを使う際には正しく動作してくれません。これでは少々（演習の進行上）勝手が悪いので、別ウィンドウで開くように設定を変えます。

上部にある「ツール」から「設定」を選ぶ
左側の一覧から「iPythonコンソール」を選ぶ
右側上部のタブから「グラフィクス」を選ぶ
右側の「グラフィックのバックエンド」の中の「バックエンド」の項目を「インライン」から「自動」に変更
すぐ下の「OK」
spyderを再起動

Jupyter Notebookを使う人は、このページの最後にある付録を参照してください。

matplotlibによるグラフ描画

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

まずは次のコードを実行してみてください。

import numpy as np
import matplotlib.pyplot as plt
x = np.linspace(0, 2.0*np.pi, 101)
y = np.sin(x)
plt.plot(x, y)

1行目は前回学んだ数学関数・定数を使うためのモジュールであるnumpyの読み込みです。
2行目が今回学ぶ、グラフ描画のためのモジュールであるmatplotlibの読み込みですが、その中でもpyplotというものを使いますので、matplotlib.pyplotを直接pltと略すようにしています。
4行目でxのリストを、5行目でそれに対応するyのリストを作っています。
5行目では、得られたリストxとリストyの組を折線で繋げて描画しています。

グラフの横軸や縦軸のラベルや、凡例を追加するには、次のようにします。

import matplotlib.pyplot as plt
import numpy as np

x = np.linspace(0, 2*np.pi, 101)
y = np.sin(x)
plt.plot(x, y, label="sin(x)")
plt.xlabel("x", fontsize=20)
plt.ylabel("y", fontsize=20)
plt.legend(loc="best", fontsize=20)

4行目までは同じです。5行目にlabel="sin(x)"が追加されています。ここで凡例での表示を指定しておきます。6行目にx軸のラベル、7行目にy軸のラベル、8行目に凡例の表示位置が指定されています。

計算数学演習で学んだgnuplotとは異なり、先にplotで描画をした後にラベルや凡例の指定をしていることに注意！

単純に必要最低限を描画するだけであればこれだけで済みますが、今後のためにもう少し踏み込んで、「オブジェクト指向」な方法を学びましょう。

小ネタ：オブジェクト指向って？

オブジェクト指向とは、プログラミングの一種の手法であり、プログラム内のデータとそれに関連する操作を、オブジェクトという単位で組織化する考え方です。

クラスに関する完全な表現ではないのですが、具体例を通して考えてみましょう。あるサークルに10人の会員がいるとします。それぞれの会員の情報として、名前（name）・年齢（age）・所属学部（faculty）・競技歴（career）・昨年度の成績（score）などの項目があるとしましょう。これまでの知識をベースにこの情報を変数のセットにしようとすると、name・age・…について、それぞれ人数分のリストを用意する、という方法をとると思います（上図左）。
これに対してオブジェクト指向では、これらの変数などをまとめたmemberなるものを定義し、会員ごとにこのmemberという変数の組（クラス）を割り当てるということをします（上図右）。仮に会員が変化しても、いちいち個別の配列の個数を増減する必要はなく、memberさえ増やしてしまえば良いわけです。

実はこれだけではオブジェクト指向の説明としては不十分で、例えばPythonにおけるオブジェクト指向プログラミングを理解するには、「クラス」や「メソッド」、「インスタンス化」や階層構造などについても学ぶ必要があります。
この授業ではそこまで学ばなくても、まぁ、複数の変数や関数がまとまって階層構造を作っているよ、くらいに理解しておけば十分です。

それではmatplotlibのオブジェクト指向な方法を学んでいきます。

matplotlibの構造

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

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

先程の例との対応で言うと、AxisやAxesがnameやageに対応し、Figureがmemberに対応していると思ってください。
実際にはAxisの中にもさらに色々な変数やメソッドが組み込まれて、より複雑なシステムですが、ここでは最も頻度が高く重要な要素に絞って学びます。

早速、コンソール上で以下のようにコマンドを入力してください。

import matplotlib.pyplot as plt
fig = plt.figure()

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

※ Spyderにおいて、図をインライン表示に設定している場合、この段階ではウィンドウは表示されません。この後の軸の設定をしないと表示されません。このページの最初に記載している方法で、グラフを別ウィンドウで開くように設定するか、コンソールではなくエディタからファイルを直接編集して実行してください。

左上にあるような何も描かれていない真っ白なウィンドウが現れたと思います。これのウインドウがFigureの実体になります。ウインドウの大きさを変えることもできます。次のコマンドを実行してみてください。さきほどのコマンドに比べるとかなり大きなウィンドウが現れたと思います。

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

一見するとただの何も書かれていないウインドウですが、変数エクスプローラーでfigの中身を見てみると、膨大な変数やメソッド（オブジェクト内の関数のこと）で嫌になると思います。
（要 「非サポートのデータ型を除外」のチェックを外す）

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

ax = fig.add_subplot(111)

先ほど定義したfigのメソッドであるadd_subplotにより、ウィンドウの中に四角い枠と縦軸、横軸に[0,1]区間の目盛が加えられ、それらを含んだオブジェクトをaxと置きました。準備されたFigureに座標系を与えることで、数値を図示することが出来るわけです。実はFigureの中には複数の座標系を用意したり、他の字を描くための構造を書くことが出来ます。

fig.add_subplot(111)の引数の意味は、100の位が縦方向のグラフの分割数、10の位が横方向のグラフの分割数で、1の位がその何番目かを表します。たとえば、fig.add_subplot(211)とすれば、上下に2分割されたうちの上半分にグラフが現れます。詳細はこの下の「複数のグラフを一つの画面に描画する」にあります。

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

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

import numpy as np
x = np.linspace(0, 2.0*np.pi, 101)
y = np.sin(x)
ax.plot(x, y)
plt.show()

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

グラフ描画の2つの流儀：オブジェクト指向インターフェース vs. pyplotインターフェース

これまで作成したの2つのプログラムを実行して比較してみましょう。

プログラム1（pyplotインターフェース）

import matplotlib.pyplot as plt
import numpy as np

x = np.linspace(0, 2*np.pi, 101)
y = np.sin(x)
plt.plot(x, y, label="sin(x)")
plt.xlabel("x", fontsize=20)
plt.ylabel("y", fontsize=20)
plt.legend(loc="best", fontsize=20)

プログラム2（オブジェクト指向インターフェース）

import matplotlib.pyplot as plt
import numpy as np

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

x = np.linspace(0, 2*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)
plt.show()

実行してみると、どちらも同じsin波が現れたと思います。 2つのプログラムを比べると、大きな違いは2点です。

まずプログラム2では、

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

のようにfigやaxというインスタンスを定義し、それらを利用してplotや設定を行います。一方プログラム1では、すべてpltを利用しており、上で説明したようなオブジェクトを有効に利用していません。

また、プログラム1にはplt.figureやax.add_subplotに相当するものが明示的に書かれていないですが、問題なくグラフが出てきたと思います。これはplt.plotという関数が、ひとたび呼び出されれば、細かいことを気にしない代わりに基本設定のまま、空気を読んでグラフを描いてくれる関数だからです。これによって、plt.figureやax.add_subplotなしにデータをすぐに可視化してくれるわけです。

今後の授業では、どちらの方法についても紹介する予定ではいますが、できるだけオブジェクト指向インターフェイスに慣れてください。

これらのプログラムで描かれたグラフには、自然科学の実験やデータ解析で必要な必要最低限の情報が記載されています。今後レポート課題等はこのように必要な情報が載ったグラフの提出を心がけてください。

小ネタ：matplotlib？pyplot？

matplotlibとは、Matlabという科学技術計算ソフトの描画システムを真似て作られたもの、と言われています。今回の例の場合、Matlabでは、plt.〇〇の「plt.」を取り除けばそのまま使えるほど、使い方は近いです。ところで最近では、オブジェクト指向が主流になりつつあり、Matlabでもオブジェクト指向での記述が標準仕様になりました（現状は、従来式でも可）。グラフの詳細設定に関しては、オブジェクト指向インターフェースの方が、従来式より有用であることは、MatlabでもPythonでも同様です。おそらく今後は、様々なプログラミングの場面において、オブジェクト指向主義が完全な主流になります。早い段階からオブジェクト指向プログラミングを身につけていくと良いでしょう。

課題1

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

オブジェクト指向インターフェースによる描画の詳細設定

ここからはオブジェクト指向のこのやり方は、先の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)
plt.show()

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

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

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

この例にあるようにしてもう一度グラフを描くと、赤い点線でグラフが描かれたと思います。matplotlibのplot関数では、データとなる二つの配列を与えた後に、文字列を書き込むことでグラフの色と線の種類を選択することが出来ます。色は次に示す省略した文字で指定し、線の種類は記号の組み合わせで指定します。基本的には、”[色][線の種類]”という構造になっています。一文字のアルファベットで指定できる色には、 ‘b’（青）, ‘g’（緑）, ‘r’（赤）, ‘c’（シアン）, ‘m’（マゼンダ）, ‘y’（黄）, ‘k’（黒）, ‘w’（白）があります。

次に線の種類を変更するには、一文字のアルファベットのあとに記号の組み合わせで指定できるのですが、簡単に指定できるものとして　’-‘, ‘–‘, ‘-.’, ‘:’ があります。それぞれ試してみて線の種類を見てみてください。（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)

マーカーのプロット

線の種類はデータ間をつないで可視化しました。ここでは、そのデータ点自体を描画する方法を学びます。散布図などを描画する際には非常に重宝します。
先ほどと同様に、ベースコード7行目のax.plotの部分を次のように変更してみてください。

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

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

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

描画範囲の指定

それでは次にグラフを描く範囲の指定方法を学びます。特に指定しなければ、ベースコードで実行したようにデータを全て描ける範囲を自動で取ってくれます。ここでは、あえて描画する範囲を絞って見ます。次の2行を8行目以降（ax.plot(~)を実行したあと）に追加してください。

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

ベースコードを改良し、線の色を赤くかつマーカーを打ち、x軸の範囲を[0, 1]に拡大して直線的に見える様子を確認せよ。

複数の線を一つのグラフの中に描く

複数のグラフを描くには、2つの方法があります。

方法1

import matplotlib.pyplot as plt
import numpy as np

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)
y = np.cos(x)
ax.plot(x, y)
plt.show()

1つ目の方法は、ax.plotを何回も書く方法です。for文などで複数の線を繰り返し描く場合に重宝します。

方法2

import matplotlib.pyplot as plt
import numpy as np

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

x = np.linspace(0, 2.0 * np.pi, 101)
y1 = np.sin(x)
y2 = np.cos(x)
ax.plot(x, y1, x, y2)
plt.show()

2つ目の方法は、ax.plotの中に繰り返して描く方法です。
直接書く分、あまりにも線の数が多い場合には向いていませんが、線が少ない場合は簡潔に描けます。

どちらの使い方も利点・欠点はありますので、どちらも使えるようになっておくと良いと思います。

図形を描く

円を描く

図形を描くことはいろいろな場面で必要になるかと思いますが、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")
plt.show()

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

課題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, "-")
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)
plt.show()

ただし、この場合多角形の内部の塗りつぶしができません。塗りつぶしできる図形描くためのモジュールが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', fill=False)

# 長方形を作る
r = patches.Rectangle(xy=(0.3, 0), width=0.25, height=0.5, angle=-10, fc='#00FFFF', ec='#00FFFF', fill=true)

# 作った図形を座標系の中へ描く
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.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関数に入力する座標がその範囲に収まっている必要があります。サンプルコードを順に見ていきましょう 1つ目の書き方が基本形です。特に指定することもないのであれば、これで十分です。 2つ目の書き方は、LaTeXコマンドを使う方法です。文字列は ” ” か ‘ ‘ で包めば文字列と認識されるのですが、その前にrとつけておくとLaTeXコマンドが利用できます。加えて、ラベルの指定と同様にfontsizeで文字の大きさを指定できます。 3つ目の書き方は、文字のフォントスタイルをstyle=’italic’とすることで斜体へ変更しています。 4つ目の書き方は、これまで座標位置が必ずテキストの左下だったものが、テキストの右下に変わっています。これはverticalalignment='bottom', horizontalalignment='right'によるもので、verticalalignmentには{‘center’, ‘top’, ‘bottom’, ‘baseline’, ‘center_baseline’}が指定でき、horizontalalignmentには{‘center’, ‘right’, ‘left’}が指定できます。 5つ目の書き方は、これまでのものと少し違います。上記の4つのtextでは、初めの座標がテキストの左下の位置と対応していることを見るためにax.plotで点を同じ座標に打っていて、ちゃんと左下角にあることが見て取れていると思います。しかし、最後のものだけ一致していません。これは、後ろの方に書いてあるtransform=ax.transAxesというコマンドのためです。これを行うとaxによって作られた座標系の中の位置ではなく、figure全体に[0,1]×[0,1]の座標系があるような、絶対的な位置により指定されることになります。

グラフへの数値の描画

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

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形式だと最後の例のように変数の値を変更すると数値が変更されていることがわかると思います。

小ネタ：print関数再訪

上の例で、format形式を学びました。その上で、以下の例を見てみましょう。

a = 1
b = 2
print('a = ', a, ', b = ', b, sep='')
print(f'a = {a}, b = {b}')

実行してみるとどうでしょうか？a = 1, b = 2が2回出てきたと思います。3行目は第2回で学んだやり方、4回目は今回学んだやり方ですが、結果に違いはありません。では、この形式、つまり、

ax.text(3, 2, ['a = ', a, ', b = ', b])

と置き換えたらどうなったのでしょうか？（※ sep=''を含めるとエラーになるので、あらかじめ消しています。） print文で同じ結果だったので…？それとも…？実際にどうなるか自分で確かめてみてください。

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

複数のデータを一つのグラフ上にプロットするには以下のように二つの座標系を準備します。 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)

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)
plt.show()

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

fig.add_subplot(2,1,1)

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

1	2	3	4
5	6	7	8
9	10	11	12

課題5

上下にグラフを分割して、上に任意の二次以上の多項式関数のグラフを描き、下段にその導関数のグラフを描け。定義域は0<x<3とする。ただし、導関数は数値的（近似的）に求めよ。

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

今日は、計算数理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)のグラフを描いてください。

課題6

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

※ pythonではlambdaが既に違う用途で定義されており、lambda = 1.0とすることはできない。lambなど、別の名前を使う。

課題7

課題6の式において、 $N(\tau) = 0.5N_0$ となる $\tau$ （半減期といいます）を手計算で求めよ。答え： $\log 2$ (約 0.693)

課題8

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

方法

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

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

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

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

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

を求めよ。

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

参考：
課題8、課題9（次に出てきます）の方法は、が具体的に与えられていない場合 （例えば実測データなど）にも有効なので、実用的です。

課題9

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

$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$ の近似値を求め、課題7で求めた値との差

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

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

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

課題10

課題9において刻み数100で近似値を求めた場合の誤差と同程度に小さい誤差とするためには課題8において刻み数をどの程度にしなければならないか。実験してみよ。（実は、刻み数を1ずつ増やしていくと、刻み幅によってたまたま上手くいってしまうことがあるので、何回か連続で課題9のときの値を下回るまで傾向をみると良い（少なくとも刻み数300000程度））

課題11（余力のある人はトライ！）

課題8および課題9の方法において、横軸刻み数・縦軸err（半減期の真の解と数値解との誤差の絶対値）とした次のような図を描け。

横軸・縦軸を対数軸にするには、それぞれ

ax.set_xscale('log')
ax.set_yscale('log')

matplotlibでのアニメーション

この授業の範囲で手軽にアニメーションの描画を行う方法には、以下の方式で十分です。ただし、matplotlibのanimationの関数を使った、より高速な方法もあります。下の方に付録として記載しておくので、余力がある人は見てみてください。

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

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(100):
    t = dt * i
    y = np.sin(x - v*t)
    lines[0].set_data(x,y)
    y = np.sin(x - v*t)
    lines[0].set_data(x,y)
    plt.pause(0.01)

plt.show()

グラフ自体の実態としてlinesがまず、lines = ax.plot(x, y, “r”)によって作られています。ax.plot関数はグラフを描くと同時に、関数として描かれたグラフに関する様々な情報をreturnしてくれおり、この値をlinesに代入しています。このlines、実はlistになっています。そこで、linesの1個目の要素について、データを更新しています。ここでは詳細に触れませんが、その中でlineとして受け取ったものを変更するとグラフが変わることだけ覚えておいてください。

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

小ネタ：linesがリストの理由

linesはなぜリストなのでしょうか？
実は複数のグラフを描く方法について、上で紹介した「複数の線を一つのグラフの中に描く」に関係します。
動かす線が2本の場合について考えましょう。
「複数の線を一つのグラフの中に描く」の方法1に対応するアニメーションは、次のようなプログラムになります。

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
line1 = ax.plot(x, y, "r")
line2 = ax.plot(x, y, "b")
for i in range(100):
    t = dt * i
    y = np.sin(x - v*t)
    line1[0].set_data(x,y)
    y = np.sin(x + v*t)
    line2[0].set_data(x,y)
    plt.pause(0.01)

plt.show()

一方、「複数の線を一つのグラフの中に描く」の方法2の場合は、次のようになります。

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
line = ax.plot(x, y, "r", x, y, "b")
for i in range(100):
    t = dt * i
    y = np.sin(x - v*t)
    line[0].set_data(x,y)
    y = np.sin(x + v*t)
    line[1].set_data(x,y)
    plt.pause(0.01)

plt.show()

課題12

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

課題13

いずれかの方法を用いて、中心座標に円を描きその周りを円が公転するアニメーションを作成せよ（運動方程式を考える必要はない）。

付録：matplotlibのanimation関数を使ったアニメーションの描画

方法1：全ての図を変数に格納し、パラパラ漫画的に描画

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

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(500):
    t = dt * i
    y = np.sin(x - v*t)
    im = ax.plot(x, y, "r")
    ims.append(im)

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

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

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

方法1にかなり近いやり方ですが、少し違う方法を紹介します。

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
line = ax.plot(x, y)

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

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

ani = animation.FuncAnimation(fig, animate, init_func = init, interval=1, blit=True, save_count=100)
plt.show()

この方法ではfor文を直接使わずに、animation.FuncAnimation関数にfor文で使われるコードをanimate関数として与えることでアニメーションを描いています。

どのやり方も一長一短なのですが、たぶん方法0が一番楽かと思います。ただし、作ったアニメーションを動画形式で保存するには、方法1か方法2であると楽でしょう。

付録：Jupyter Notebookでアニメーションを見るための設定

※macOS Mojave以前ではこの方法は使えません！強制ログアウトされます！

 macOS Catalina以降では解決しているようです。

前回の最後に、Jupyter Notebookを使ってグラフを描画し、Jupyter Notebook上で見るためには、初めにおまじないとして

%matplotlib inline

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

%matplotlib tk

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