第３・４回（第２週）

今週の目標

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

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

特に何も設定していない状態でグラフなどを描画すると、iPythonコンソールの中に描画されます。これでは少々（演習の進行上）勝手が悪いので、別ウィンドウで開くように設定を変えます。

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

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

matplotlibによるグラフ描画

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

matplotlibの構造

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

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

実際にはより複雑なシステムですが、ここでは最も頻度が高く重要な要素に絞って学びます。

※本当は全てオブジェクトと呼ばれていますが、本演習ではオブジェクト指向という考え方について深くは触れません。ひとまずグラフを描いたのちに、グラフ描画の2つの流儀について述べていますので、そちらを参考にしてください。何を言っているのか分からなくても、本演習の進め方には対して影響しませんので、構いません。

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

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

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

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

import matplotlib.pyplot as plt

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

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

import matplotlib.pyplot as plt
fig = plt.figure(figsize=(10,8))

ax = fig.add_subplot(111)

見ての通り真っ白だったウィンドウの中に四角い枠と縦軸、横軸に[0,1]区間の目盛が加えられました。準備された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
import matplotlib.pyplot as plt

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)

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

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

次の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*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)

プログラム2（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)

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

まずプログラム1では、

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

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

また、plt.figure（および第2回で使ったplt.show）というコマンドを消しましたが、問題なくグラフが出てきたと思います。これはplt.plotという関数が、ひとたび呼び出されれば、細かいことを気にしない代わりに基本設定のまま、空気を読んでグラフを描いてくれる関数だからです。これによって、plt.figureやax.add_subplotなしにデータをすぐに可視化してくれるわけです。

いずれの方法でも構いませんが、これらのプログラムで描かれたグラフには、自然科学の実験やデータ解析で必要な必要最低限の情報が記載されています。今後レポート課題等はこのように必要な情報が載ったグラフの提出を心がけてください。なお，上記プログラム中の”legend”関数についてはこちらのサイトに詳しく紹介しています．

小ネタ：オブジェクト指向？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ではデフォルトのグラフの色として、青色が選択され、線で結ばれてグラフが描画されます。ここでは、任意の色に変更したり、点線にしたりしてみます。変更方法はとても簡単で、上のサンプルコードの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)

マーカーのプロット

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

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

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

図形を描く

円を描く

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

課題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()

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.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とする。ただし、導関数は数値的（近似的）に求めよ。

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

グラフのデータが、line.set_data(x,y)の部分で更新されます。グラフ自体の実態としてlineがまず、line, = ax.plot(x, y, “r”)によって作られています。ax.plot関数はグラフを描くと同時に、関数として描かれたグラフに関する様々な情報をreturnしてくれます。ここでは詳細に触れませんが、その中でlineとして受け取ったものをこの例題の様に制御するとグラフが変わることだけ覚えておいてください（気になる人は、変数エクスプローラーの歯車をクリックし、「Exclude unsupported data types」のチェックを外すと、linesの中身を見ることができます）。

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

練習

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

（ラインの更新方法については述べたが，patchの更新方法については紹介していない．
patches.Circleを使った方法で，中心座標を更新して描画する場合には，
変数名.set_center(xy=(x, y))で中心座標の更新ができる．)

例えば
c=patches.Circle・・・・・
と置いたのであれば，
c.set_center(xy=(x, y))
で中心座標の更新が可能．

課題6

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

エクストラ課題1

課題6の公転する円を中心に，さらに一つ公転する円を描き（太陽，地球，月のような円を描き）アニメーションを作成せよ。その後，任意の円の数Nを入力すると，連鎖公転する円が描けるプログラムへ拡張せよ．ただし，各円の公転の角速度は（2π,4π,・・・2Nπ）となるものとする．（初期位相phaiはすべて0とし，公転軌道半径は任意のもので構わない．）

エクストラ課題2

エクストラ課題1で作成した”最後に生成される円”の軌道をアニメーション図中に併載させよ（軌道描画はアニメでなくても構わない）。この際，図の縦軸は”y”，図の横軸は”x”とする

エクストラ課題3

エクストラ課題2で描いた軌道および時間を基に，横軸が時間”t”，縦軸が”y”のアニメーション図(または，静止グラフ)を別途作成せよ。

エクストラ課題4

エクストラ課題3で，描いた軌道が近似的に矩形関数 – Wikipediaのような形にできるのか，またできるとしたら軌道半径・回転の角速度・初期位相をどのようにすれば良いか，考えよ。

エクストラ課題5

エクストラ課題4で矩形波を描けた原理を理解し，この動画を見て，連鎖公転円系を2軸分用意して任意の図形（ハートやミッキーなど）を2次元平面上に描画するプログラムを作成せよ(動画の2:40秒あたりに再生されているものを作る)。ただし，ハートを描くための関数はこちらを参照．（手順ヒント：①まずは2次元平面上に一筆書きの目標図形を描き，x座標用配列,・y座標用配列に値を格納する．②x座標用配列,・y座標用配列と同じサイズの時間t用配列を作成する．③得られた配列x(t), y(t)に対してそれぞれフーリエ級数展開を実装し，フーリエ係数を手に入れる．この際，同じ周波数成分のsinとcosのフーリエ係数から，位相差，および，ノルムを算出しておく．⑤得られた位相差・ノルム・周波数が公転系の何を意味するのかを考えてx軸用，y軸用の公転系をそれぞれ実装する．）

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

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

(P)

\[\left\{
\begin{array}{l}
\frac{d}{dt} N(t) &= – \lambda N, \\
N(0) &= N_0
\end{array}\right.\]

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

(S)
$$ N(t) = N_0e^{- \lambda t} $$

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

課題7

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

課題8

課題7の式において、$N( \tau ) \ \ = \ \ 0.5N_0 $となる$\tau$（半減期といいます）を手計算で求めよ．ただし，$\ \lambda = 1$とする．

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

課題9

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

方法

$N(t)$の値を(S)を用いて飛び飛びの$t$の値$t_i \ \ = \ \ dt * i$ （ただし$dt$は時間刻み幅）について求めるのであるが、(S)と課題６のグラフから明らかなように$N(t)$は単調減少関数である。であるから、$N(t_i) \ \ \geq 0.5N_0 \geq N(t_{i+1})$となる$i$があるであろう。このとき、$t_i$を$\tau$の近似値$\tau^*$とする。$0 \ \ \leq t \ \ \leq 10$を100等分して上記方法にて$\tau$の近似値$\tau^*$を求め、課題8で求めた値との差（誤差と呼ぶ）$$err^* \ \ = \ \ | \tau-\tau^* | $$を求めよ。

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

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

課題10

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

$N(t_i) \ \ \geq 0.5N_0 \geq N(t_{i+1})$となる$i$を求め、２点$(t_i, N(t_i)), \ \ (t_{i+1}, N(t_{i+1}))$ を結ぶ直線を考える。その直線と$N=0.5$の直線の交点の$t$の値$\tau^{**}$をもって、$\tau$の近似値とする。（このような手法を線形補間と呼ぶ。）$0 \ \ \leq t \leq 10$を100等分して上記方法にて$\tau$の近似値$\tau^{**}$を求め課題8で求めた値との差$$err^{**} \ \ = \ \ | \tau – \tau^{**} |$$

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

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

課題11

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

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

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

付録：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(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()

演習で行なったやり方（以下、方法0と呼びます）では、ax.plot関数はfor文の外側で一度使っただけですが、今回はfor文の中で毎回呼び出します。そして、そのグラフに関する情報をimとして受け取ります。方法0ではlinesを作るときにlinesの後ろに”, “がついていたのですが、今回は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.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, frames=range(0,2*314,1), interval=1, blit=True, save_count=100)

plt.show()

この方法ではfor文を直接使わずに、animation.FuncAnimation関数にfor文で使われるコードをanimate関数として与えることでアニメーションを描いています。init_func関数で初期化し，animate関数でlineのy座標を更新・描画、framesがiに渡される引数でfor loop文の代わりの役目をしています．

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

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

※macOS Mojaveではこの方法は使えません！強制ログアウトされます！ macOS Catalinaでは解決しているようですが…

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

%matplotlib inline

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

%matplotlib tk

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