Python の Matplotlib によるグラフ作成には，plt.*** という関数を使った plt (pyplot) 流（pyplot インターフェースとも）と，ax.*** という関数を使った ax 流（オブジェクト指向インターフェースとも）の2つの方法があります。
ここでは，ax 流で ax.*** のみを使ってグラフを作成する方法についてまとめておきます。
plt (pyplot) 流で plt.*** のみを使ってグラフを作成する方法については，以下のページにまとめています。
ライブラリの import と初期設定
必要なライブラリを import し，適宜いくつか初期設定をしておきます。
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.ticker import MultipleLocator
# 以下はグラフを SVG で Notebook にインライン表示させる設定
%config InlineBackend.figure_formats = ['svg']
# mathtext font の設定
plt.rcParams['mathtext.fontset'] = 'cm'
# デフォルトの figsize の設定変更は，
# 一度何か描いてからにするとよいようだ。
plt.rcParams['figure.figsize'] = [6.4, 4.8]
ax.plot() の基本形
以下，他のグラフ作成ページとの比較のため，便宜上，
- 陽関数のグラフ
- 媒介変数表示のグラフ
- 点・数値データのグラフ
などと分類して説明していますが，Matplotlib の ax.plot() でグラフを描く際の基本形は，（ax.contour() で陰関数のグラフを描く場合を除けば）以下の通り，基本的に「数値データのグラフを作成する方法」のみで対応します。
まず最初に
fig, ax = plt.subplots()
として ax を定義する。
点 (x0, y0), (x1, y1), … (xn, yn) を結んで折れ線グラフを描くには，
- $x$ 座標だけのリスト
[x0, x1, ..., xn]， - $y$ 座標だけのリスト
[y0, y1, ..., yn]をつくり，以下のように…
ax.plot([x0, x1, ..., xn], [y0, y1, ..., yn]);
陽関数 $y = f(x)$ や媒介変数表示 $x(t), y(t)$ のグラフを描く場合も，描画範囲の $x$ のリストや $t$ のリストを作って，上記のパターンにする。
等間隔の数値リストの作成法
描画範囲を細かく分ければ，滑らかな曲線に見えてきます。しかし，座標値のリストを [0, 0.1, 0.2, ...] などと手で一つ一つ入力するのは大変なので，以下のように等間隔の数値配列を作成する関数を使います。
np.arange() は間隔を指定
start 以上，stop 未満，一定間隔 step の数値リストを作成するには，
np.arange(start, stop, step)
start の省略値は 0，step の省略値は 1。
# 0以上5未満，間隔1ごとの数値要素を作成
x = np.arange(0, 5, 1)
print(x)
print(x**2)
# start と step を省略すると...
print(np.arange(5))
np.linspace() は要素数を指定
start から stop までを num-1 等分し，端点を含めて合計 num 個の数値リストを作成するには，
np.linspace(start, stop, num)
num の省略値は 50。
# 0から2までを4等分し，端点を含めて5個の数値要素を作成
x = np.linspace(0, 2, 5)
print(x)
# NumPy の関数を使う
print(np.sin(x))
陽関数のグラフ ax.plot()
$x$ の陽関数 $y = f(x)$ を $x_{\rm min} \leq x \leq x_{\rm max}$ の範囲でグラフを作成するには ax.plot() を使って…
# ax を使う際の最初のおまじない
fig, ax = plt.subplots()
x = np.linspace(x_min, x_start)
y = f(x) # numpy の関数であること
ax.plot(x, y)
例として，$y = \sin x$ のグラフを $-2\pi \leq x \leq 2\pi$ の範囲で描きます。三角関数を含む全ての関数は，NumPy の関数を使います。
今の場合は $\sin x = $ np.sin(x)。
基本的な定数の一つである円周率 $\pi$ は NumPy の np.pi を使います。
# ax を使うさいのおまじない
fig, ax = plt.subplots()
x = np.linspace(-2*np.pi, 2*np.pi)
y = np.sin(x)
ax.plot(x, y);
plt.rcParams['figure.figsize']
# デフォルトの figsize の設定変更は，
# 一度何か描いてからにするとよいようだ。
plt.rcParams['figure.figsize'] = [6.4, 4.8]
デフォルトの num = 50 では少しギクシャクしたグラフになるようです。$x$ の範囲を少し細かく分割して，滑らかな曲線に見えるグラフにしてみます。
fig, ax = plt.subplots()
x = np.linspace(-2*np.pi, 2*np.pi, 200)
y = np.sin(x)
ax.plot(x, y);
グラフ全体のオプションの設定例
グラフ全体のオプションは ax. で設定します。
- タイトル
ax.set_title() - 軸ラベル
ax.set_xlabel()ax.set_ylabel() - 表示範囲
ax.set_xlim(xmin,xmax)ax.set_ylim(ymin,ymax)
一挙に設定するにはax.axis([xmin,xmax,ymin,ymax]) - 縦横比
ax.set_aspect() - グリッドの表示・設定
ax.grid() - 凡例表示
ax.legend()
ラインごとのオプション設定例
ラインごとのオプションは ax.plot() 内に…
- 凡例ラベル
label = ' ' - 線の太さ
lw (linewidth) = - 線の色
c (color) = - 線種
ls (linestyle) = ':','--'等
いくつかオプションを設定する例。
fig, ax = plt.subplots()
x = np.linspace(-2*np.pi, 2*np.pi, 200)
y = np.sin(x)
ax.plot(x, y,
label = '正弦関数',
lw = 2, c = 'red', ls = ':')
ax.set_title('タイトル')
ax.set_xlabel('$x$')
ax.set_ylabel('$f(x)$')
ax.set_xlim(-10, 10)
ax.set_ylim(-1.1, 1.1)
ax.grid()
ax.legend();
初等関数のグラフ
陽関数 $y = f(x)$ のグラフの描き方を使って，（双曲線関数および逆双曲線関数を含む）初等関数のグラフを描いてみます。
べき関数
まず，$y = x^{-1}$ のグラフ。
fig, ax = plt.subplots()
# 滑らかな曲線にするには num = 200 程度は必要
x = np.linspace(-5, 5, 200)
y = x**(-1)
ax.plot(x, y)
ax.set_xlim(-5, 5)
ax.set_ylim(-10, 10)
ax.grid();
$y=x^{-1}$ のグラフの $x＝0$ 付近，$y \rightarrow -\infty$ と $y \rightarrow +\infty$ を縦の直線でつないでしまう。
不連続点をつなげない対策 np.nan
以下では，y の絶対値がある程度以上（以下の例では abs(y) > 20）であれば，その点はグラフに描かないように「欠損値」np.nan に置き換えている。
fig, ax = plt.subplots()
# 滑らかな曲線にするには num = 200 程度は必要
x = np.linspace(-5, 5, 200)
y = x**(-1)
y[abs(y) > 20] = np.nan
ax.plot(x, y)
ax.set_xlim(-5, 5)
ax.set_ylim(-10, 10)
ax.grid();
$y = x^{-2}, \ x^{-1}, \ x^2, \ x^3$ のグラフ例。$y = x^{-1}$ の不連続点をつなげないように。
fig, ax = plt.subplots()
# 滑らかな曲線にするには num = 200 程度は必要
x = np.linspace(-5, 5, 200)
ax.plot(x, x**(-2), label = r'$x^{-2}$')
y = x**(-1)
y[abs(y) > 20] = np.nan
ax.plot(x, y, label = r'$x^{-1}$')
ax.plot(x, x**2, label = r'$x^2$')
ax.plot(x, x**3, label = r'$x^3$')
ax.set_xlabel('$x$')
ax.set_ylabel('$f(x)$')
ax.set_xlim(-5, 5)
ax.set_ylim(-10, 10)
ax.grid()
ax.legend();
$\displaystyle y = \sqrt{x}, \ \frac{1}{\sqrt{x}}$ のグラフ例。
NumPy の関数 $\sqrt{x} = $ np.sqrt(x) を使います。
fig, ax = plt.subplots()
# x = 0 は除く。
x = np.linspace(0.001, 5, 200)
ax.plot(x, np.sqrt(x), label = r'$\sqrt{x}$')
ax.plot(x, 1/np.sqrt(x), label = r'$1/\sqrt{x}$')
ax.set_xlabel('$x$')
ax.set_ylabel('$f(x)$')
ax.set_xlim(0, 5)
ax.set_ylim(0, 5)
ax.grid()
ax.legend();
指数関数
$y = e^{-x}, \ e^x$ のグラフ例。
fig, ax = plt.subplots()
# x = 0 は除く。
x = np.linspace(-5, 5, 200)
ax.plot(x, np.exp(-x), label = '$e^{-x}$')
ax.plot(x, np.exp(x), label = '$e^{x}$')
ax.set_xlabel('$x$')
ax.set_ylabel('$f(x)$')
ax.set_xlim(-5, 5)
ax.set_ylim(-1, 30)
ax.grid()
ax.legend();
対数関数
$y = \log x$ のグラフ例。
fig, ax = plt.subplots()
# x = 0 は除く。
x = np.linspace(0.0001, 10, 200)
ax.plot(x, np.log(x), label = r'$\log x$')
ax.set_xlabel('$x$')
ax.set_ylabel('$f(x)$')
ax.set_xlim(-0.5, 10)
ax.set_ylim(-5, 5)
ax.grid()
ax.legend();
三角関数
$ y = \sin x, \ \cos x, \ \tan x $ のグラフ例。
$\tan x$ の不連続点をつなげない対策 np.nan
y[abs(y) > 10] = np.nan によって $y = \tan x$ の不連続点の処理を行っている。
fig, ax = plt.subplots()
x = np.linspace(-2*np.pi, 2*np.pi+1, 200)
ax.plot(x, np.sin(x), label = r'$\sin x$')
ax.plot(x, np.cos(x), label = r'$\cos x$')
y = np.tan(x)
y[abs(y) > 10] = np.nan # tan(x) の不連続点の処理
ax.plot(x, y, label = r'$\tan x$')
ax.set_xlabel('$x$')
ax.set_ylabel('$f(x)$')
ax.set_xlim(-2*np.pi, 2*np.pi+1)
ax.set_ylim(-5, 5)
ax.grid()
ax.legend();
主目盛・副目盛・グリッドのカスタマイズ
ここでは，横軸主目盛を $\pi$ ごと，横軸副目盛を $\frac{\pi}{2}$ ごと，縦軸副目盛を $1$ ごとにしてみます。
fig, ax = plt.subplots()
x = np.linspace(-2*np.pi, 2*np.pi+1, 200)
ax.plot(x, np.sin(x), label = r'$\sin x$')
ax.plot(x, np.cos(x), label = r'$\cos x$')
y = np.tan(x)
y[abs(y) > 10] = np.nan # tan(x) の不連続点の処理
ax.plot(x, y, label = r'$\tan x$')
ax.set_xlabel('$x$')
ax.set_ylabel('$f(x)$')
ax.set_xlim(-2*np.pi, 2*np.pi+1)
ax.set_ylim(-5, 5)
# 横軸主目盛
ax.set_xticks(
[-2*np.pi, -np.pi, 0, np.pi, 2*np.pi],
['$-2\pi$', '$-\pi$', '$0$', '$\pi$', '$2\pi$'])
# 横軸副目盛
ax.xaxis.set_minor_locator(MultipleLocator(np.pi/2))
# 縦軸副目盛
ax.yaxis.set_minor_locator(MultipleLocator(1))
# 主・副目盛両方に同じスタイルで目立たなくグリッド
ax.grid(c='lightgray', ls='--', lw=0.5, which='both')
ax.legend();
参考：$x$ 軸・$y$ 軸の表示 ax.axhline(0) ax.axvline(0)
また，Python のグラフ作成界隈では，$x$ 軸（$y = 0$）や $y$ 軸（$x = 0$）を（グリッドとは区別して）表示させる習慣はないようです。このあたりは以下のようにして対応できます。
- $x$ 軸（$y = 0$）の表示
ax.axhline(0) - $y$ 軸（$x = 0$）の表示
ax.axvline(0)
fig, ax = plt.subplots()
x = np.linspace(-2*np.pi, 2*np.pi+1, 200)
ax.plot(x, np.sin(x), label = r'$\sin x$')
ax.plot(x, np.cos(x), label = r'$\cos x$')
y = np.tan(x)
y[abs(y) > 10] = np.nan # tan(x) の不連続点の処理
ax.plot(x, y, label = r'$\tan x$')
ax.set_xlabel('$x$')
ax.set_ylabel('$f(x)$')
ax.set_xlim(-2*np.pi, 2*np.pi+1)
ax.set_ylim(-5, 5)
# 横軸主目盛
ax.set_xticks(
[-2*np.pi, -np.pi, 0, np.pi, 2*np.pi],
['$-2\pi$', '$-\pi$', '$0$', '$\pi$', '$2\pi$'])
# 横軸副目盛
ax.xaxis.set_minor_locator(MultipleLocator(np.pi/2))
# 縦軸副目盛
ax.yaxis.set_minor_locator(MultipleLocator(1))
# 主・副目盛両方に同じスタイルで目立たなくグリッド
ax.grid(c='lightgray', ls='--', lw=0.5, which='both')
# x軸 y軸は dashed に。k は black
ax.axhline(0, c='k', ls='--', lw=0.5)
ax.axvline(0, c='k', ls='--', lw=0.5)
ax.legend();
逆三角関数
- $y = \sin^{-1} x = \arcsin x =$
np.arcsin(x)- 定義域は $\displaystyle-1 \leq x \leq 1$
- 値域は $\displaystyle -\frac{\pi}{2} \leq y \leq \frac{\pi}{2}$
- $y = \cos^{-1} x = \arccos x =$
np.arccos(x)- 定義域は $\displaystyle-1 \leq x \leq 1$
- 値域は $\displaystyle0 \leq y \leq \pi$
- $y = \tan^{-1} x = \arctan x =$
np.arctan(x)- 定義域は $\displaystyle-\infty < x < \infty$
- 値域は $\displaystyle-\frac{\pi}{2} \leq y \leq \frac{\pi}{2}$
fig, ax = plt.subplots()
x = np.linspace(-1, 1, 100)
ax.plot(x, np.arcsin(x), label = r'$\arcsin x$')
ax.plot(x, np.arccos(x), label = r'$\arccos x$')
xt = np.linspace(-10, 10, 200)
ax.plot(xt, np.arctan(xt), label = r'$\arctan x$')
ax.set_xlabel('$x$')
ax.set_ylabel('$f(x)$')
ax.set_xlim(-7, 7)
ax.set_ylim(-np.pi/2, np.pi)
# 縦軸主目盛
ax.set_yticks(
[-np.pi/2, 0, np.pi/2, np.pi],
['$-\pi/2$', '$0$', '$\pi/2$', '$\pi$'])
# 横軸主目盛を 2 ごとに
ax.xaxis.set_major_locator(MultipleLocator(2))
# 横軸副目盛を 1 ごとに
ax.xaxis.set_minor_locator(MultipleLocator(1))
# 縦軸副目盛を π/4 ごとに
ax.yaxis.set_minor_locator(MultipleLocator(np.pi/4))
# 主・副目盛両方に同じスタイルで目立たなくグリッド
ax.grid(c='lightgray', ls='--', lw=0.5, which='both')
# x軸 y軸は dashed に。k は black
ax.axhline(0, c='k', ls='--', lw=0.5)
ax.axvline(0, c='k', ls='--', lw=0.5)
ax.legend();
双曲線関数
- $y = \sinh x =$
np.sinh(x) - $y = \cosh x =$
np.cosh(x) - $y = \tanh x =$
np.tanh(x)
fig, ax = plt.subplots()
x = np.linspace(-5, 5, 200)
ax.plot(x, np.sinh(x), label = r'$\sinh x$')
ax.plot(x, np.cosh(x), label = r'$\cosh x$')
ax.plot(x, np.tanh(x), label = r'$\tanh x$')
ax.set_xlabel('$x$')
ax.set_ylabel('$f(x)$')
ax.set_xlim(-5, 5)
ax.set_ylim(-10, 10)
# 横軸主目盛
ax.xaxis.set_major_locator(MultipleLocator(1))
# 縦軸主目盛
ax.yaxis.set_major_locator(MultipleLocator(2))
# 縦軸副目盛
ax.yaxis.set_minor_locator(MultipleLocator(1))
# 主・副目盛両方に同じスタイルで目立たなくグリッド
ax.grid(c='lightgray', ls='--', lw=0.5, which='both')
# x軸 y軸は dashed に。k は black
ax.axhline(0, c='k', ls='--', lw=0.5)
ax.axvline(0, c='k', ls='--', lw=0.5)
ax.legend();
逆双曲線関数
- $ y = \sinh^{-1} x = \mbox{arsinh}\ x = $
np.arcsinh(x)の定義域は $ -\infty < x < \infty$ - $ y = \cosh^{-1} x = \mbox{arcosh}\ x = $
np.arccosh(x)の定義域は $ 1 \leq x < \infty$ - $ y = \tanh^{-1} x = \mbox{artanh}\ x = $
np.arctanh(x)の定義域は $ -1 < x < 1$
（Mathematica や）NumPy の逆双曲線関数の表記は確信犯なのかなぁ。逆双曲線関数の表記については
fig, ax = plt.subplots()
x = np.linspace(-10, 10, 200)
ax.plot(x, np.arcsinh(x), label = r'$\sinh^{-1} x$')
x = np.linspace(1, 10, 200)
ax.plot(x, np.arccosh(x), label = r'$\cosh^{-1} x$')
x = np.linspace(-0.99999, 0.99999, 100)
ax.plot(x, np.arctanh(x), label = r'$\tanh^{-1} x$')
ax.set_xlabel('$x$')
ax.set_ylabel('$f(x)$')
ax.set_xlim(-10, 10)
ax.set_ylim(-5, 5)
# 横軸主目盛
ax.xaxis.set_major_locator(MultipleLocator(2))
# 横軸副目盛
ax.xaxis.set_minor_locator(MultipleLocator(1))
# 縦軸主目盛
ax.yaxis.set_major_locator(MultipleLocator(1))
# 主・副目盛両方に同じスタイルで目立たなくグリッド
ax.grid(c='lightgray', ls='--', lw=0.5, which='both')
# x軸 y軸は dashed に。k は black
ax.axhline(0, c='k', ls='--', lw=0.5)
ax.axvline(0, c='k', ls='--', lw=0.5)
ax.legend();
陰関数のグラフ ax.contour()
例として，$x^2 + y^2 = 1$ のグラフを描きます。$y$ について解いて陽関数の形にしたり，以下で述べるような媒介変数表示にしなくても，陰関数のまま，グラフに描くことができます。
$z \equiv x^2 + y^2$ として，$z = 1$ の等高線をプロットします。
fig, ax = plt.subplots()
xrange = np.linspace(-2, 2, 200)
yrange = np.linspace(-2, 2, 200)
x, y = np.meshgrid(xrange, yrange)
z = x**2 + y**2
# z = 1 の等高線を描く
ax.contour(x, y, z, [1]);
縦横比 ax.set_aspect()，サイズの設定 figsize
縦横の比を設定して円らしく見えるようにします。グラフのサイズ figsize も設定してみます。
# グラフのサイズの設定：縦横を同じ長さに
fig, ax = plt.subplots(figsize=[6.4, 6.4])
xrange = np.linspace(-2, 2, 200)
yrange = np.linspace(-2, 2, 200)
x, y = np.meshgrid(xrange, yrange)
z = x**2 + y**2
# z = 1 の等高線を描く
cs = ax.contour(x, y, z, [1])
# z = 1 のラベル
ax.clabel(cs)
# 縦横のアスペクト比を等しくして円らしく見えるように
ax.set_aspect('equal')
ax.set_xlabel('$x$')
ax.set_ylabel('$y$')
ax.set_title('$x^2 + y^2 = 1$ の等高線')
# 表示範囲 xmin, xmax, ymin, ymax を一挙に設定
ax.axis([-1.1, 1.1, -1.1, 1.1]);
# 横軸主目盛
ax.xaxis.set_major_locator(MultipleLocator(0.5))
# 縦軸主目盛
ax.yaxis.set_major_locator(MultipleLocator(0.5))
# 主・副目盛両方に同じスタイルで目立たなくグリッド
ax.grid(c='lightgray', ls='--', lw=0.5, which='both')
# x軸 y軸
ax.axhline(0, c='k', ls='--', lw=0.5)
ax.axvline(0, c='k', ls='--', lw=0.5);
媒介変数表示のグラフ ax.plot()
半径 $1$ の円の方程式は $x^2 +y^2 = 1$ です。この円を以下のような媒介変数表示にして描きます。
$$ x=\cos t, \quad y=\sin t \quad (0 \leq t \leq 2\pi) $$
# グラフのサイズの設定：縦横を同じ長さに
fig, ax = plt.subplots(figsize=[6.4, 6.4])
t = np.linspace(0, 2*np.pi, 100)
ax.plot(np.cos(t), np.sin(t), label = '$x^2 + y^2 = 1$')
# 縦横のアスペクト比を等しくして円らしく見えるように
ax.set_aspect('equal')
ax.set_xlabel('$x$')
ax.set_ylabel('$y$')
# 表示範囲 xmin, xmax, ymin, ymax を一挙に設定
ax.axis([-1.1, 1.1, -1.1, 1.1]);
# 横軸主目盛
ax.xaxis.set_major_locator(MultipleLocator(0.5))
# 縦軸主目盛
ax.yaxis.set_major_locator(MultipleLocator(0.5))
# 主・副目盛両方に同じスタイルで目立たなくグリッド
ax.grid(c='lightgray', ls='--', lw=0.5, which='both')
# x軸 y軸
ax.axhline(0, c='k', ls='--', lw=0.5)
ax.axvline(0, c='k', ls='--', lw=0.5)
ax.legend();
点・数値データのグラフ
折れ線グラフ ax.plot()
以下の例では，配列 Xp に 座標の値，配列 Yp に 座標の値を入れて，6個の点をつないだ折れ線グラフを描きます。
fig, ax = plt.subplots()
Xp = np.arange(0, 6)
Yp = Xp**2
ax.plot(Xp, Yp)
ax.grid();
散布図（線でつながない） ax.scatter()
各点を線でつながずに赤色でプロットする例。c = 'red' で点の色を赤に設定。
fig, ax = plt.subplots()
Xp = np.arange(0, 6)
Yp = Xp**2
ax.scatter(Xp, Yp, c = 'red', label = 'データ')
ax.grid()
ax.legend();
ax.plot() でも散布図
点種（'ro' で赤い丸）を指定することで，ax.plot() でも線でつながない散布図を描くことができます。
fig, ax = plt.subplots()
Xp = np.arange(0, 6)
Yp = Xp**2
ax.plot(Xp, Yp, 'ro', label = 'データ')
ax.grid()
ax.legend();
ファイルから数値データを読み込んでプロット
あらかじめ作成された数値データファイルからデータを読み込んでグラフを描くこともできます。
以下のようなテキストファイル mydata2.txt が（このノートブックファイルと同じフォルダに）あるとします。（授業では事前に mydata2.txt も同じフォルダにダウンロードしておきます。）
! cat mydata2.txt
np.loadtxt() で読み込み
数値データの書き込み・読み込みには，NumPy の savetxt() や loadtxt() を使うことにします。
# ファイル mydata2.txt から読み込み，data に代入
data = np.loadtxt('mydata2.txt')
print(data)
読み込んだ data の1列目 data[:,0] を $x$ 座標に，2列目 data[:,1] を $y$ 座標にして散布図を描きます。
fig, ax = plt.subplots()
ax.scatter(data[:,0], data[:,1], label = 'mydata2.txt のデータ')
ax.grid()
ax.legend();
ベクトルを描く ax.quiver()
原点を始点とし，成分が $v_x = 0.5, v_y = 1$ のベクトルを描く例。
fig, ax = plt.subplots(figsize=[6.4, 6.4])
# 縦横のアスペクト比を等しく
ax.set_aspect('equal')
# 表示範囲 xmin,xmax,ymin,ymax を一気に設定
ax.axis([-2, 2, -2, 2])
# 始点の x 座標
X = [0]
# 始点の y 座標
Y = [0]
# ベクトルの x 成分
Vx = [0.5]
# ベクトルの y 成分
Vy = [1]
# ベクトルを描く
ax.quiver(X, Y, Vx, Vy, color="blue",
# 以下の3点セットを書かないと自動スケーリングされる
angles='xy', scale_units='xy', scale=1
)
ax.grid();
複数のベクトルを描く例。斜方投射の速度ベクトルを描いてみます。
適宜無次元化を行って，斜方投射の位置 $(x(t), y(t))$ と速度 $(v_x(t), v_y(t))$ は
\begin{eqnarray}
x &=& t \\
y &=& t – \frac{1}{4} t^2 \\
v_x &=& 1 \\
v_y &=& 1 – \frac{1}{2} t
\end{eqnarray}
のように書けるとします。
t = np.arange(0, 5)
X = t
Y = t - 0.25*t**2
# Vx = 1 だけだとリストにならないので
Vx = 1 + t*0
Vy = (1 - 0.5*t)
fig, ax = plt.subplots(figsize=[6.4, 6.4])
# 縦横のアスペクト比を等しく
ax.set_aspect('equal')
# 表示範囲 xmin,xmax,ymin,ymax を一気に設定
ax.axis([-0.5, 5.5, -1.5, 2])
# ベクトルを描く
ax.quiver(X, Y, Vx, Vy, color="blue",
# 以下の3点セットを書かないと自動スケーリングされる
angles='xy', scale_units='xy', scale=1
)
ax.grid();
複数のグラフを重ねて表示
複数のグラフを重ねて表示する例を示します。
複数の陽関数を重ねて表示
陽関数 $y = x^2 – 1$ と $y = 4 x – 5$ の 2 つのグラフを重ねて描きます。
$x$ の範囲は独立に指定可能です。
fig, ax = plt.subplots()
x = np.linspace(-5, 5)
ax.plot(x, x**2 - 1, label = '$x^2 - 1$')
x = np.linspace(1, 3)
ax.plot(x, 4*x - 5, label = '$4 x - 5$')
ax.set_ylim(-5, 10)
ax.grid()
ax.legend();
数値データと理論曲線を重ねて表示
前節の数値データファイル mydata2.txt と理論曲線 $y = x^2$ の 2 つのグラフを重ねて表示してみます。
描く順番を指定 'zorder'
複数のグラフを重ねて表示する際，オプション zorder を使って描かれる順番を調整できます。zorder = 0 だと最初（つまり最背面）に描かれます。
fig, ax = plt.subplots()
ax.scatter(data[:,0], data[:,1], c = 'blue', label = 'データ')
x = np.linspace(0, 5)
ax.plot(x, x**2, c = 'red', label = '理論曲線', zorder = 0)
ax.grid()
ax.legend();
さまざまなグラフを重ねて表示
例えば，前掲の斜方投射
\begin{eqnarray}
x &=& t \\
y &=& t – \frac{1}{4} t^2 \\
\therefore\ \ y &=& x – \frac{1}{4} x^2
\end{eqnarray}
を使って，
- 陽関数 $\displaystyle y = f(x)$
- 媒介変数表示 $x = x(t), \ y = y(t) = $
- 数値データ
- ベクトル
を重ねて表示してみます。
fig, ax = plt.subplots()
# 陽関数を少し上にずらして
x = np.linspace(0, 4)
y = x - 0.25*x**2
ax.plot(x, y+0.5, label = '陽関数')
# 媒介変数表示，といっても陽関数と同じ書き方
t = np.linspace(0, 4)
x = t
y = t - 0.25*t**2
ax.plot(x, y, label = '媒介変数表示', zorder = 0)
# 数値データ
ax.scatter(X, Y, c = 'g', label = 'データ')
# ベクトル
ax.quiver(X, Y, Vx, Vy, color="blue", label = '速度ベクトル',
# 以下の3点セットを書かないと自動スケーリングされる
angles='xy', scale_units='xy', scale=1
)
ax.axis([-0.5, 5, -1, 2])
ax.set_aspect('equal')
ax.grid()
ax.legend();
極座標表示のグラフ
極座標表示 $r = r(\phi)$ のグラフは，以下のように媒介変数表示のグラフとして描きます。
\begin{eqnarray}
r &=& r(\phi) \\
x(\phi) &=& r(\phi) \cos\phi \\
y(\phi) &=& r(\phi) \sin\phi
\end{eqnarray}
焦点を原点とした楕円のグラフ
太陽からの万有引力を受けて運動する惑星は，太陽（二体問題では質量中心）を焦点の一つとした楕円運動を描いて運動します。
焦点の一つを原点とし，長半径 $a$，離心率 $e$ の楕円の方程式は，極座標 $r$，$\varphi$ を使って以下のように表すことができます。
$$ r= \frac{a (1−e^2)}{ 1 + e\cos\varphi}$$
海王星と冥王星の軌道
さて，さて，かつては第9惑星，現在では準惑星の一つである冥王星も楕円軌道を描きます。
まず，極座標表示の楕円の式を関数として定義します。
# elliptical（楕円）の e をつけて...
def re(a, e, phi):
return a*(1-e**2)/(1+e*np.cos(phi))
def xe(a, e, phi):
return re(a, e, phi)*np.cos(phi)
def ye(a, e, phi):
return re(a, e, phi)*np.sin(phi)
冥王星の軌道長半径 $𝑎_P=39.445 \mbox{au}$，離心率 $𝑒_P=0.250$
を使って冥王星の軌道を描きます。
fig, ax = plt.subplots(figsize = [6.4, 6.4])
aP = 39.445
eP = 0.250
phi = np.linspace(0, 2*np.pi, 200)
ax.plot(xe(aP, eP, phi), ye(aP, eP, phi),
c = 'g', label = '冥王星')
ax.set_aspect('equal');
では，この冥王星の軌道と，その内側をまわる海王星の軌道を重ねて描いてみましょう。
海王星の軌道長半径は $a_N=30.181 \mbox{au}$，離心率は $e_N=0.0097$ と小さいので簡単のために $e_N=0$ として扱います。
実際の2つの天体の軌道は同じ平面上にありませんが，太陽からの距離のみをみるために，ここでは同一平面上に描きます。
fig, ax = plt.subplots(figsize = [6.4, 6.4])
# 冥王星
aP = 39.445
eP = 0.250
phi = np.linspace(0, 2*np.pi, 200)
ax.plot(xe(aP, eP, phi), ye(aP, eP, phi),
c = 'g', label = '冥王星')
# 海王星
aN = 30.181
eN = 0
ax.plot(xe(aN, eN, phi), ye(aN, eN, phi),
c = 'b', label = '海王星')
# 太陽
ax.scatter([0], [0], c = 'r', label = '太陽')
# x軸 y軸は dashed に。k は black
ax.axhline(0, c='k', ls='--', lw=0.5)
ax.axvline(0, c='k', ls='--', lw=0.5)
ax.set_aspect('equal')
ax.legend();
参考：各種図形の表示と塗りつぶし
matplotlib.patches を以下のように import して使います。
from matplotlib.patches import *
正多角形 RegularPolygon()
(x0, y0) を中心とし，「半径」R の正 N 角形
p = RegularPolygon((x0, y0), N, radius = R)
ax.add_patch(p)
頂点の座標を指定した多角形 Polygon()
(x1, y1), (x2, y2), ... を頂点とした多角形
p = Polygon([(x1, y1), (x2, y2), ... ])
ax.add_patch(p)
デフォルトでは塗りつぶしなので，塗りつぶさないときは fill = False
fig, ax = plt.subplots(figsize = [6.4, 6.4])
# 正8角形
regPoly = RegularPolygon((0, 0), 8, radius = 2,
label = '正8角形')
ax.add_patch(regPoly)
# 頂点の座標を指定
ps = [(-3, -3), (1, -3), (3, 3), (-2, 2)]
poly = Polygon(ps, label = '頂点の座標を指定', fill = False,
color = 'tab:orange', lw = 1.5)
ax.add_patch(poly)
# 各頂点
ax.scatter(ps[0][0], ps[0][1], label = str(ps[0]))
ax.scatter(ps[1][0], ps[1][1], label = str(ps[1]))
ax.scatter(ps[2][0], ps[2][1], label = str(ps[2]))
ax.scatter(ps[3][0], ps[3][1], label = str(ps[3]))
ax.set_aspect('equal')
ax.axis([-3.5, 5, -3.5, 3.5])
ax.grid()
ax.legend();
長方形 Rectangle()
点 (x0, y0) を左下の頂点座標とし，幅（底辺）width，高さ height の長方形は，
p = Rectangle((x0, y0), width, height)
ax.add_patch(p)
fig, ax = plt.subplots(figsize = [6.4, 6.4])
# 底辺 2, 高さ 3 の長方形
p = Rectangle((-3, -3), 2, 3, label = '長方形')
ax.add_patch(p)
# 1辺の長さ 4 の正方形
# グリッドの上に描くように
p = Rectangle((0, -2), 4, 4, label = '正方形', zorder = 2,
fill = False, edgecolor = 'red', lw = 2)
ax.add_patch(p)
ax.set_aspect('equal')
ax.axis([-3.5, 5, -3.5, 3.5])
ax.grid()
ax.legend();
楕円 Ellipse()
点 (x0, y0) を中心とし，幅 width，高さ height の楕円は，
p = Ellipse((x0, y0), width, height)
ax.add_patch(p)
したがって，長半径（横半径）a，単半径（縦半径）b の楕円は，
p = Ellipse((x0, y0), 2*a, 2*b)
ax.add_patch(p)
塗りつぶさない場合は fill = False をつけて。
fig, ax = plt.subplots(figsize = [6.4, 6.4])
# 塗りつぶした楕円
a = 2
b = 1
e1 = Ellipse((-2, -2), 2*a, 2*b,
label = '$a = %d, b = %d$' % (a, b))
ax.add_patch(e1)
# 塗りつぶさない楕円
a = 3
e = 0.8
b = a*np.sqrt(1-e**2)
e2 = Ellipse((2, 0), 2*a, 2*b, fill = False,
color = 'tab:orange', lw = 1.5,
label = '$a = %d, e = %.1f$'%(a, e))
ax.add_patch(e2)
ax.set_aspect('equal')
ax.axis([-4.5, 5.5, -3.5, 2.5])
ax.grid()
ax.legend();
円 Circle()
点 (x0, y0) を中心とし，半径 r の円は，
p = Circle((x0, y0), r)
ax.add_patch(p)
枠線色 'edgecolor' と背景色 'facecolor'
Circle() に限らず，全ての図形の枠線の色を edgecolor (ec) で，塗り潰し（背景）の色を facecolor (fc) でそれぞれ独立に設定できます。
fig, ax = plt.subplots(figsize = [6.4, 6.4])
p = Circle((1, 0), 2, lw = 2,
ec = 'black',
fc = 'yellow')
ax.add_patch(p)
ax.axis([-1.3, 3.3, -2.3, 2.3])
ax.set_aspect('equal')
ax.grid();
弧 Arc()，扇形 Wedge()
点 (x0, y0) を中心とし，幅 width，高さ height の楕円の開始角度 t1 度から，終了角度 t2 度までの弧を描くには，
p = Arc((0, 0), width, height, theta1=t1, theta2=t2)
ax.add_patch(p)
長半径 a，単半径 b の楕円の場合は，width = 2*a，height = 2*b であるから
p = Arc((0, 0), 2*a, 2*b, theta1=t1, theta2=t2)
ax.add_patch(p)
グラフ内にテキスト ax.text()
グラフ内の座標 x = x1, y = y1 あたりにテキスト 'mytext' を表示させるには，
ax.text(x1, y1, 'mytext')
fig, ax = plt.subplots(figsize = [6.4, 6.4])
# 楕円の弧
a = 3
b = 2
width = 2*a
height = 2*b
p = Arc((0, 0), width, height,
theta1 = 180, theta2 = 270,
lw = 2, color = 'tab:blue', label = '楕円の一部')
ax.add_patch(p)
# 楕円の弧の残り
p = Arc((0, 0), width, height,
theta1 = -90, theta2 = 180,
lw = 2, color = 'tab:orange', label = '楕円の残り')
ax.add_patch(p)
# 扇形
p = Wedge((-1, -1), 3, theta1 = 30, theta2 = 60,
ec = 'black',
fc = 'yellow', label = '扇形')
ax.add_patch(p)
# 直線
ax.plot([-1, 4], [-1, -1], color = 'k')
# 円弧で角度
p = Arc((-1, -1), 3, 3, theta1 = 0, theta2 = 30,
color = 'blue', lw = 1.5)
ax.add_patch(p)
ax.text(0.1, -0.8, '30°', color = 'blue')
# 円弧で角度
p = Arc((-1, -1), 3.5, 3.5, theta1 = 0, theta2 = 60,
color = 'red', lw = 1.5)
ax.add_patch(p)
ax.text(0.7, -0.3, '60°', color = 'red')
# 円弧で角度
p1 = Arc((-1, -1), 1.7, 1.7, theta1 = 30, theta2 = 60)
p2 = Arc((-1, -1), 1.8, 1.8, theta1 = 30, theta2 = 60)
ax.add_patch(p1)
ax.add_patch(p2)
ax.text(-0.35, -0.3, '30°')
ax.axis([-4, 4, -3, 3])
ax.set_aspect('equal')
ax.grid()
ax.legend();
領域の塗りつぶし
2本の陽関数で挟まれた領域を塗りつぶす ax.fill_between()
$x$ の範囲 xrange で $y = f(x)$ と $y = g(x)$ で挟まれた領域を塗りつぶすには，
ax.fill_between(xrange, f(xrange), g(xrange))
$0.5 \leq x \leq 2$ で $y = f(x)$ と $x$ 軸（$y = 0$）に囲まれた領域を塗りつぶす例。$x$ 軸 $y$ 軸を表示します。
デフォルトのグリッド線が少し目障りな場合は，以下のように細めの灰色点線にして目立たないようにすることもできます。
def f(x):
return 0.6*x + 0.4*np.cos(x)
fig, ax = plt.subplots()
# 塗りつぶしの範囲
x2 = np.linspace(0.5, 2)
# y = f(x) と y = 0 の間を塗りつぶし
ax.fill_between(x2, f(x2), 0, fc="yellow")
ax.axis([-0.2, 2.5, -0.2, 1.5])
ax.grid();
図形の内部を塗りつぶす ax.fill()
一般的な図形の内部を塗りつぶすには，その図形を（滑らかな曲線で囲まれた図形であっても）多数の頂点からなる多角形と考えて塗りつぶします。
例えば，楕円の一部。楕円の中心を原点としたデカルト座標 $X, Y$ では，楕円の式は
$$\frac{X^2}{a^2} + \frac{Y^2}{b^2} = 1$$
これは以下のように媒介変数表示できます。
\begin{eqnarray}
X(\theta) &=& a \cos\theta \\
Y(\theta) &=& b \sin\theta
\end{eqnarray}
$\theta_1 = 10^{\circ}$ から $\theta_2 = 60^{\circ}$ までの楕円の弧と，原点からの直線で囲まれた領域を図示します。
ax.plot() は曲線のグラフを描くのに使いますが，始点と終点を同一にすると，閉曲線で囲まれた図形（の外枠）を描くことになります。
# 引数の角度は度で
def Xe(th):
return a*np.cos(np.radians(th))
def Ye(th):
return b*np.sin(np.radians(th))
th1 = 10
th2 = 60
a = 2
b = 1
fig, ax = plt.subplots(figsize = [6.4, 6.4])
# 楕円の弧と原点を結ぶ線で囲まれた図形
# theta の範囲
trange = np.linspace(th1, th2)
# 最初と最後に 0 を追加
# リストに変換してから前後に追加
Xlist = [0] + Xe(trange).tolist() + [0]
Ylist = [0] + Ye(trange).tolist() + [0]
ax.plot(Xlist, Ylist, color = 'black')
ax.axis([-1, 3, -1, 2])
ax.set_aspect('equal')
ax.grid();
この図形の内部を塗りつぶすには，ax.plot() を ax.fill() に変えるだけです。
fig, ax = plt.subplots(figsize = [6.4, 6.4])
# 楕円の弧と原点を結ぶ線で囲まれた図形
# theta の範囲
trange = np.linspace(th1, th2)
# 最初と最後に 0 を追加
# リストに変換してから前後に追加
Xlist = [0] + Xe(trange).tolist() + [0]
Ylist = [0] + Ye(trange).tolist() + [0]
# 内部を塗りつぶし
ax.fill(Xlist, Ylist, lw = 1.5,
ec = 'black', fc = 'yellow')
ax.axis([-1, 3, -1, 2])
ax.set_aspect('equal')
ax.grid();
○練習：ドッテンメイカイを塗りつぶし
冥王星と海王星の軌道のグラフをみると，軌道が交差している。これは，ある期間冥王星のほうが海王星より太陽に近いことを意味する。そこで，以下の指示に従ってグラフを作成せよ。
1. 海王星と冥王星の軌道が交差する点を求める
簡単のため，海王星は半径 $a_N$ の円軌道とし，実際の2つの天体の軌道は同じ平面上にないが，太陽からの距離のみをみるために，ここでは同一平面上に描くことにする。
$$r(a, e, \phi) = \frac{a (1-e^2)}{1 + e \cos\phi}$$
とし，
$$ r(a_P, e_P, \phi) = a_N$$
となる $\phi$ の値 $\phi_{\rm eq}$ (phieq) を数値的に求める。SymPy による方程式の数値解の求め方については，すでに解説している。以下を参照：
ここでは，Maplotlib で必要とされる NumPy と SymPy を混在させない（たとえば三角関数は両方で定義されているから）という方針で，後で別途説明する SciPy の root_scalar() を使って方程式の数値解を求める方法を紹介する。以下を参照：
$$f(x) \equiv r(a_P, e_P, x) – a_N = 0$$
となる $x$ の値 phieq を root_scalar() を使って数値的に求める。
from scipy.optimize import root_scalar
# 念のため，すでに定義しているがここでもう一度コピペ
# 冥王星 Pluto の軌道長半径，離心率
aP = 39.445
eP = 0.250
# 海王星 Neptune の軌道半径
aN = 30.181
eN = 0
# root_scalar に渡す関数の引数は x 決め打ち
def f(x):
return aP*(1-eP**2)/(1+eP*np.cos(x)) - aN
# 0 < x < π/2 の間で f(x) = 0 の数値解を探す
answer = root_scalar(f, bracket = [0, np.pi/2])
phi_eq = answer.root
phi_eq
# もう一つの解は...
root_scalar(f, bracket = [-np.pi/2, 0]).root
軌道の対称性から，$- \phi_{\rm eq} < \phi < + \phi_{\rm eq}$ の間は冥王星のほうが海王星より太陽に近いことになる。
2. あらためて海王星と冥王星の軌道のグラフを描く
# 念のため，すでに定義しているがここでもう一度コピペ
# elliptical（楕円）の e をつけて...
def re(a, e, phi):
return a*(1-e**2)/(1+e*np.cos(phi))
def xe(a, e, phi):
return re(a, e, phi)*np.cos(phi)
def ye(a, e, phi):
return re(a, e, phi)*np.sin(phi)
fig, ax = plt.subplots(figsize = [6.4, 6.4])
# 冥王星
aP = 39.445
eP = 0.250
phi = np.linspace(0, 2*np.pi, 200)
ax.plot(xe(aP, eP, phi), ye(aP, eP, phi),
c = 'g', label = '冥王星')
# 海王星
aN = 30.181
eN = 0
ax.plot(xe(aN, eN, phi), ye(aN, eN, phi),
c = 'b', label = '海王星')
# 太陽
ax.scatter([0], [0], c = 'r', label = '太陽')
# x軸 y軸は dashed に。k は black
ax.axhline(0, c='k', ls='--', lw=0.5)
ax.axvline(0, c='k', ls='--', lw=0.5)
ax.set_aspect('equal')
ax.legend();
3. 塗りつぶす扇形の部分を定義
最後に，冥王星のほうが海王星より太陽に近い期間の扇形の部分を，塗りつぶす領域として定義します。
# 楕円の弧と原点を結ぶ線で囲まれた図形
# phi の範囲
phirange = np.linspace(-phi_eq, phi_eq)
# 最初と最後に 0 を追加
# リストに変換してから前後に追加
Xlist = [0] + xe(aP, eP, phirange).tolist() + [0]
Ylist = [0] + ye(aP, eP, phirange).tolist() + [0]
4. 冥王星が海王星より太陽に近い期間を求める
$$t_1 = T_N \times \frac{\phi_{\rm eq}}{\pi}$$
で与えられる時間 $t_1$ の間だけ，冥王星が海王星より太陽に近いことになる。
以下のようにして $t_1$ (t1) を求める。
# 海王星の公転周期は何年？
TN = ...
# round() で数値を丸める
t1 = round(TN * phi_eq/(pi))
# グラフに表示させる準備
mytext = "この期間は"+str(t1)+"年"
# グラフの x=10, y=0.5 からはじまるテキストを書く
ax.text(10, 0.5, mytext)
あるいはグラフのタイトルに書いてもよいだろう。
ax.set_title('冥王星が海王星よりも太陽に近い期間は??年')
5. 扇形の部分を塗りつぶす
塗りつぶす部分を多角形だと考えて…
参考：凡例の位置 loc = 'upper left'
凡例の表示位置を左上にするには，
ax.legend(loc = 'upper left')
完成イメージ
だいたい以下のようなグラフができると思います。
月別平年気温のカーブフィット
弘前市の月別平年気温のデータを使って配列 HiroT を作ります。
Month = [i for i in range(1,13)] # 12 までのときは 13 と書く
HiroT = [-1.5, -1.0, 2.3, 8.6, 14.3, 18.3,
22.3, 23.5, 19.4, 12.9, 6.5, 0.8]
fig, ax = plt.subplots()
ax.scatter(Month, HiroT, label = '月別平年気温')
ax.grid()
ax.legend();
グラフをみると，月別平年基本は 12ヶ月周期の正弦関数あるいは余弦関数のようにみえます。では，以下のような関数でフィットしてみます。$a, b, c$ が最小二乗法で決定するパラメータです。
$$f(x, a, b, c) = a + b \sin\left(\frac{2\pi x}{12}\right)
+ c \cos\left(\frac{2\pi x}{12}\right)$$
scipy.optimize.curve_fit() で最小二乗法してみます。
help(curve_fit) から一部を引用：
curve_fit(f, xdata, ydata, ...)
Use non-linear least squares to fit a function, f, to data.
Returns
-------
popt : array
Optimal values for the parameters
pcov : 2-D array
The estimated covariance of popt.def f(x, a, b, c):
return a + b*np.sin(np.pi*x/6) + c*np.cos(np.pi*x/6)
from scipy.optimize import curve_fit
popt, pcov = curve_fit(f, Month, HiroT)
popt
最小二乗法によって求めた $a, b, c$ の値をもった $f(x, a, b, c)$ のグラフを重ねてプロットしてみます。
fig, ax = plt.subplots()
ax.scatter(Month, HiroT, label = '月別平年気温')
x = np.linspace(1, 12)
a, b, c = popt
ax.plot(x, f(x, a, b, c), label="関数フィット",
color = 'tab:orange')
# x 軸の目盛設定
ax.set_xticks(np.arange(1, 13))
ax.set_xlabel("月")
# y 軸の目盛設定
ax.set_yticks(np.arange(-5, 26, 5))
ax.set_ylabel("°C")
ax.set_title('弘前市の月別平年気温')
ax.grid()
ax.legend();
○練習：カーブフィット
弘前市の月別平年気温のデータを使って
- フィットさせる関数を（周期12ヶ月の周期関数をフーリエ級数に展開する要領で）以下のような
g(x)にし，最小二乗法で求めてみよ。
\begin{eqnarray}
g(x) = A_0 &+& a_1 \cos \frac{2\pi x}{12} + a_2 \cos \frac{4\pi x}{12} \\
&+& b_1 \sin \frac{2\pi x}{12} + b_2 \sin \frac{4\pi x}{12}
\end{eqnarray} - 月別平年気温の数値データとフィット曲線 $f(x)$ および $g(x)$ のグラフを重ねて描け。
- 月平年気温の平均値を求めてみよ。（$g(x)$ の定数部分 $A_0$ が平均値。）
def g(x, A0, a1, a2, b1, b2):
return (A0 + a1*np.cos(np.pi*x/6)
+ a2*np.cos(np.pi*x/3)
+ b1*np.sin(np.pi*x/6)
+ b2*np.sin(np.pi*x/3))
なお，3. の別解として，リストの総和 sum() と要素数 len() を使って以下のように平均値を求めることもできましたね。
sum(HiroT)/len(HiroT)
また，月平年気温の最大値や最小値も簡単に求められますね。