「SymPy Plotting Backends でグラフ作成」と同等の内容を Matplotlib を使って。
ライブラリの import
必要なライブラリを import します。
import numpy as np
import matplotlib.pyplot as plt
# 以下はグラフを SVG で Notebook にインライン表示させる設定
%config InlineBackend.figure_formats = ['svg']
plt.plot() の基本形
Matplotlib の plt.plot() でグラフを描く際の基本形は，以下の通りです。
点 (x0, y0), (x1, y1), … (xn, yn) を結んで折れ線グラフを描くには，
- $x$ 座標だけのリスト
[x0, x1, ..., xn]， - $y$ 座標だけのリスト
[y0, y1, ..., yn]をつくり，以下のように…plt.plot([x0, x1, ..., xn], [y0, y1, ..., yn]);
陽関数 $y = f(x)$ や媒介変数表示 $x(t), y(t)$ のグラフを描く場合も，描画範囲の $x$ のリストや $t$ のリストを作って，上記のパターンにする。
等間隔の数値配列の作成法
描画範囲を細かく分ければ，滑らかな曲線に見えてきます。しかし，座標値のリストを [0, 0.1, 0.2, ...] などと手で一つ一つ入力するのは大変なので，以下のように等間隔の数値配列を作成する関数を使います。
np.arange() は間隔を指定
# 0以上5未満，間隔1ごとの数値要素を作成
x = np.arange(0, 5, 1)
print(x)
print(x**2)
np.linspace() は要素数を指定
# 0から2までを4等分し，端点を含めて5個の数値要素を作成
x = np.linspace(0, 2, 5)
print(x)
# NumPy の関数を使う
print(np.sin(x))
陽関数のグラフ
例として，$y = \sin x$ のグラフを $-2\pi \leq x \leq 2\pi$ の範囲で描きます。三角関数を含む全ての関数は，NumPy の関数を使います。基本的な定数の一つである円周率 $\pi$ は NumPy の np.pi を使います。
まずは，$x$ の範囲 $-2\pi \leq x \leq 2\pi$ を 10 等分し，端点も含めて 11 個の要素のリスト x を作成します。
# x 座標のリストを作成
x = np.linspace(-2*np.pi, 2*np.pi, 11)
print(x)
次に，このリスト x を NumPy の関数 np.sin(x) に入れると，x の各要素 $x_i$ についての $\sin x_i$ の値を要素にしたリストが出来上がります。
# y 座標のリストを作成
y = np.sin(x)
print(y)
# x 座標，y 座標を入れて plt.plot()
plt.plot(x, y);
$x$ の範囲をもっと細かく分割して，滑らかな曲線に見えるグラフにしてみます。
# x 座標のリストを作成
x = np.linspace(-2*np.pi, 2*np.pi, 101) # 100 等分
# y 座標のリストを作成
y = np.sin(x)
plt.plot(x, y);
いくつかオプションを設定する例。
plt.plot(x, y,
# 凡例
label = "正弦関数",
# 線の太さと色
linewidth = 2, color = "red")
# 横軸縦軸の表示範囲
plt.xlim(-10, 10)
plt.ylim(-1.1, 1.1)
# グリッド（格子線）の表示
plt.grid()
# 凡例の表示
plt.legend();
陰関数のグラフ
例として，$x^2 + y^2 = 1$ のグラフを描きます。$y$ について解いて陽関数の形にしたり，以下で述べるような媒介変数表示にしなくても，陰関数のまま，グラフに描くことができます。
$z \equiv x^2 + y^2$ として，$z = 1$ の等高線をプロットします。
delta = 0.02
xrange = np.arange(-2, 2, delta)
yrange = np.arange(-2, 2, delta)
x, y = np.meshgrid(xrange, yrange)
z = x**2 + y**2
# z = 1 の等高線を描く
plt.contour(x, y, z, [1]);
縦横比・サイズの設定
縦横の比を設定して円らしく見えるようにします。グラフのサイズ figsize も設定してみます。
# グラフのサイズの設定：縦横を同じ長さに
fig = plt.figure(figsize=[6,6])
# 縦横のアスペクト比を等しくして円が円らしく見えるように
plt.axes().set_aspect('equal')
# 軸の設定 [xmin, xmax, ymin, ymax]
plt.axis([-1.1, 1.1, -1.1, 1.1])
delta = 0.02
xrange = np.arange(-2, 2, delta)
yrange = np.arange(-2, 2, delta)
x, y = np.meshgrid(xrange, yrange)
z = x**2 + y**2
# z = 1 の等高線を描く
plt.contour(x, y, z, [1])
# グリッド（格子線）の表示
plt.grid();
媒介変数表示のグラフ
半径 $1$ の円の方程式は $x^2 +y^2 = 1$ です。この円を以下のような媒介変数表示にして描きます。
$$ x=\cos t, \quad y=\sin t \quad (0 \leq t \leq 2\pi) $$
# 縦横のアスペクト比を等しくして円が円らしく見えるように
plt.axes().set_aspect('equal')
# 軸の設定
plt.axis([-1.1, 1.1, -1.1, 1.1])
# 媒介変数 t の数値配列を作成
t = np.linspace(0, 2*np.pi, 100)
# グラフを描く
plt.plot(np.cos(t), np.sin(t))
# x 軸の目盛設定
plt.xticks(np.linspace(-1, 1, 5))
# y 軸の目盛設定
plt.yticks(np.linspace(-1, 1, 5))
# グリッド（格子線）の表示
plt.grid();
点・数値データのグラフ
以下の例では，配列 X に 座標の値，配列 Y に 座標の値を入れて，6個の点をつないだ折れ線グラフを描きます。
X = np.linspace(0, 5, 6)
Y = X**2
plt.plot(X, Y)
# グリッド（格子線）
plt.grid();
各点を線でつながずに赤色でプロットする例。点の種類を赤いr 円 o に設定。
X = np.linspace(0, 5, 6)
Y = X**2
plt.plot(X, Y, "ro", label = "データ")
# グリッド（格子線）
plt.grid()
# 凡例の表示
plt.legend();
plt.scatter() でも，各点を線でつながずにプロットします。
plt.scatter(X, Y);
ベクトルを描く
原点を始点とし，成分が $v_x = 0.5, v_y = 1$ のベクトルを描く例。
# 始点の x 座標
X = [0]
# 始点の y 座標
Y = [0]
# ベクトルの x 成分
Vx = [0.5]
# ベクトルの y 成分
Vy = [1]
# 縦横のアスペクト比を等しくして円が円らしく見えるように
plt.axes().set_aspect('equal')
# ベクトルを描く
plt.quiver(X, Y, Vx, Vy, color="blue",
# 以下の3点セットを書かないと自動スケーリングされる
angles='xy', scale_units='xy', scale=1
)
# x 軸の目盛設定
plt.xticks(np.linspace(-2, 2, 9))
# y 軸の目盛設定
plt.yticks(np.linspace(-2, 2, 9))
plt.grid();
複数のベクトルを描く例。斜方投射の速度ベクトルを描いてみます。
t = np.linspace(0, 4, 5)
X = t
Y = t - 0.25*t**2
Vx = 0.5 + t*0
Vy = 0.5 - 0.25*t
print(X)
print(Y)
print(Vx)
print(Vy)
# 縦横のアスペクト比を等しく
plt.axes().set_aspect('equal')
# 軸の設定
plt.axis([-0.5, 4.5, -0.5, 2])
plt.quiver(X, Y, Vx, Vy, color="blue",
# 以下の3点セットを書かないと自動スケーリングされる
angles='xy', scale_units='xy', scale=1
)
plt.grid();
複数のグラフを重ねて表示
複数のグラフを重ねて表示する例を示します。
$y = x^2 – 1$ と $y = 4 x – 5$ の 2 つのグラフを重ねて描きます。
$x$ の範囲は $-5 \leq x \leq 5$ です。
x = np.linspace(-5, 5, 101)
y1 = x**2 - 1
y2 = 4*x - 5
plt.plot(x, y1, x, y2)
plt.grid();
それぞれの曲線ごとに，いくつかオプションを設定して描く例です。
x = np.linspace(-5, 5, 101)
y1 = x**2 - 1
y2 = 4*x - 5
plt.plot(x, y1, linewidth=2, color="red", label="$x^2-1$")
plt.plot(x, y2, linewidth=1, color="black", label="$4x-5$")
plt.grid()
plt.legend();
数値データファイルを読み込む
あらかじめ作成された数値データファイルを読み込んでグラフを描くこともできます。
まず，数値データファイル data.txt をあらかじめ作成しておきます。数値データの書き込み・読み込みには，NumPy の savetxt() や loadtxt() を使うことにします。
# ファイルに書き込むデータを作成
dat = []
for i in range(6):
dat.append([i, i**2])
# できた dat を表示
dat
# ファイル data.txt に dat の内容を書き込む
np.savetxt('data.txt', dat, fmt='%i') # 整数として書き出し
確かにファイルが作成されているか，確認します。
（弘大 JupyterHub は Ubuntu サーバですから，! に続けてテキストファイルの中身を見るコマンド cat data.txt を打ち込んで，内容を表示させます。）
# ファイルの内容表示
!cat data.txt
# loadtxt() で読み込み，変数 data に代入。
data = np.loadtxt('data.txt', dtype = 'int') # 整数として読み込み
# Python で作成した dat はリスト
type(dat)
# NumPy の loadtxt() で読み込んだ data は配列になる
type(data)
# numpy.ndarray 配列であれば，以下のように一挙に列データを指定可能
# data の1列目 data[:,0] を X, 2列目 data[:,1] を Y
X = data[:,0]
Y = data[:,1]
print(X)
print(Y)
plt.plot(X, Y, "ro", label="データ")
plt.legend();
数値データと理論曲線を重ねて表示
前節の数値データファイル data.txt と理論曲線 $y = x^2$ の 2 つのグラフを重ねて表示してみます。
plt.plot(X, Y, "ro", label="データ")
x = np.linspace(0, 6)
plt.plot(x, x**2, label="理論曲線 $y=x^2$")
plt.grid()
plt.legend();
海王星と冥王星の軌道
焦点を原点とした楕円のグラフ
太陽からの万有引力を受けて運動する惑星は，太陽（二体問題では質量中心）を焦点の一つとした楕円運動を描いて運動します。
焦点の一つを原点とし，長半径 $a$，離心率 $e$ の楕円の方程式は，極座標 $r$，$\varphi$ を使って以下のように表すことができます。
$$ r= \frac{a (1−e^2)}{ 1 + e\cos\varphi}$$
さて，かつては第9惑星，現在では準惑星の一つである冥王星も楕円軌道を描きます。冥王星の軌道長半径 $𝑎_P=39.445 \mbox{au}$，離心率 $𝑒_P=0.250$
を使って冥王星の軌道を描きます。
まず，極座標表示の楕円の式を関数として定義します。
def r(a, e, phi):
return a*(1-e**2)/(1+e*np.cos(phi))
def x(a, e, phi):
return r(a, e, phi)*np.cos(phi)
def y(a, e, phi):
return r(a, e, phi)*np.sin(phi)
aP = 39.445
eP = 0.250
# グラフのサイズの設定
#fig = plt.figure(figsize=[6,6])
# 縦横のアスペクト比を等しくして
plt.axes().set_aspect('equal')
# 軸の設定
#plt.axis([-60, 60, -60, 60])
phi = np.linspace(0, 2*np.pi, 200)
plt.plot(x(aP, eP, phi), y(aP, eP, phi), label="冥王星")
plt.legend();
では，この冥王星の軌道と，その内側をまわる海王星の軌道を重ねて描いてみましょう。
海王星の軌道長半径は $a_N=30.181 \mbox{au}$，離心率は $e_N=0.0097$ と小さいので簡単のために $e_N=0$ として扱います。
実際の2つの天体の軌道は同じ平面上にありませんが，太陽からの距離のみをみるために，ここでは同一平面上に描きます。
aN = 30.181
eN = 0
# 縦横のアスペクト比を等しくして
plt.axes().set_aspect('equal')
# 軸の設定
#plt.axis([-60, 60, -60, 60])
phi = np.linspace(0, 2*np.pi, 200)
plt.plot(x(aP, eP, phi), y(aP, eP, phi), label="冥王星")
plt.plot(x(aN, eN, phi), y(aN, eN, phi), label="海王星")
plt.legend();
x 軸 y 軸を引く
Matplotlib 業界では $x$ 軸や $y$ 軸を引くという慣習がないのか，そのようなオプションが見当たりません。仕方がないので，axhline() と axvline() を使ってそれらしい線を引いてみます。
# グラフのサイズの設定
fig = plt.figure(figsize=[6,6])
# 縦横のアスペクト比を等しくして
plt.axes().set_aspect('equal')
# 軸の設定
plt.axis([-60, 40, -50, 50])
phi = np.linspace(0, 2*np.pi, 200)
plt.plot(x(aP, eP, phi), y(aP, eP, phi), label="冥王星")
plt.plot(x(aN, eN, phi), y(aN, eN, phi), label="海王星")
# x軸 y軸
plt.axhline(0, color='black', dashes=(3, 3), linewidth=0.6)
plt.axvline(0, color='black', dashes=(3, 3), linewidth=0.6)
plt.legend();
月別平年気温の関数フィット
弘前市の月別平年気温のデータを使って配列 HiroT を作ります。
Month = [i for i in range(1,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]
plt.plot(Month, HiroT, "o", label="月別平年気温")
plt.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)$ のグラフを重ねてプロットしてみます。
plt.plot(Month, HiroT, "o", label="月別平年気温")
x = np.linspace(1, 12, 100)
a, b, c = popt
plt.plot(x, f(x, a, b, c), label="関数フィット")
# x 軸の目盛設定
plt.xticks(np.arange(1, 13, 1))
plt.xlabel("月")
# y 軸の目盛設定
plt.yticks(np.arange(-5, 26, 5))
plt.ylabel("°C")
plt.grid()
plt.legend();
練習 3-1
弘前市の月別平年気温のデータを使って
- フィットさせる関数を（周期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))
領域の塗りつぶし
2本の陽関数で挟まれた領域を塗りつぶす
$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)
# f(x) の表示範囲
x1 = np.linspace(0.25, 2.25, 100)
# 塗りつぶしの範囲
x2 = np.linspace(0.5, 2, 100)
plt.plot(x1, f(x1))
plt.fill_between(x2, f(x2), 0, fc="yellow")
plt.axis([-0.2, 2.5, -0.2, 1.5])
# グリッドを細めの灰色点線で
plt.grid(color="lightgray", dashes=(3, 5), linewidth=0.5)
# x軸 y軸
plt.axhline(0, color='black', dashes=(3, 5), linewidth=0.6)
plt.axvline(0, color='black', dashes=(3, 5), linewidth=0.6);
多角形の内部を塗りつぶす
# 多角形の頂点の x 座標
Xpoly = [0, 1, 2, 1, 0]
# 多角形の頂点の y 座標
Ypoly = [1, 2, 1, 0, 1]
# 縦横のアスペクト比を等しくして
plt.axes().set_aspect('equal')
# 多角形を描く
plt.plot(Xpoly, Ypoly);
# 縦横のアスペクト比を等しくして
plt.axes().set_aspect('equal')
# 多角形の内部を塗りつぶす
plt.fill(Xpoly, Ypoly);
円の内部を塗りつぶす
円周上の多数の頂点からなる多角形と考えて…
def enx(t):
return np.cos(t)
def eny(t):
return np.sin(t)
t = np.linspace(0, 2*np.pi, 200)
Xen = enx(t)
Yen = eny(t)
# 縦横のアスペクト比を等しくして
plt.axes().set_aspect('equal')
# 多角形として円を描く
plt.plot(Xen, Yen);
# 縦横のアスペクト比を等しくして
plt.axes().set_aspect('equal')
# 多角形として円を塗りつぶす
plt.fill(Xen, Yen, fc = "yellow", edgecolor = "tab:blue");
t = np.linspace(np.pi/6, np.pi/3, 100)
# list にしてから結合
Xen = [0] + enx(t).tolist() + [0]
Yen = [0] + eny(t).tolist() + [0]
# 縦横のアスペクト比を等しくして
plt.axes().set_aspect('equal')
plt.xlim(-1.1, 1.1)
plt.ylim(-1.1, 1.1)
# 多角形として扇形を描く
plt.plot(Xen, Yen);
# 縦横のアスペクト比を等しくして
plt.axes().set_aspect('equal')
plt.xlim(-1.1, 1.1)
plt.ylim(-1.1, 1.1)
# 扇形の内部を塗りつぶす
plt.fill(Xen, Yen, fc = "yellow")
# edgecolor を使わずに縁を黒線にする
# 多角形として扇形を塗りつぶす
plt.plot(Xen, Yen, c = "black");
練習 3-2：メイカイにカイメイせよ
冥王星と海王星の軌道のグラフをみると，軌道が交差している。これは，ある期間冥王星のほうが海王星より太陽に近いことを意味する。そこで，以下の指示に従ってグラフを完成させてください。
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):
global aP, eP, aN
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])
phieq = answer.root
# 参考までにラジアンから度に変換してみると...
np.degrees(phieq);
# f = 0 の数値解は -π/2 < x < 0 にもある
answer2 = root_scalar(f, bracket = [-np.pi/2, 0])
answer2.root;
軌道の対称性から，$- \phi_{\rm eq} < \phi < + \phi_{\rm eq}$ の間は冥王星のほうが海王星より太陽に近いことになる。
2. Matplotlib で海王星と冥王星の軌道のグラフを描く
# 念のため，すでに定義しているがここでもう一度コピペ
# 楕円の極座標表示
def r(a, e, phi):
return a*(1-e**2)/(1+e*np.cos(phi))
# x 座標
def x(a, e, phi):
return r(a, e, phi)*np.cos(phi)
# y 座標
def y(a, e, phi):
return r(a, e, phi)*np.sin(phi)
# グラフのサイズの設定
fig = plt.figure(figsize=[6,6])
# 縦横のアスペクト比を等しくして
plt.axes().set_aspect('equal')
# 軸の設定
plt.axis([-60, 40, -50, 50])
phi = np.linspace(0, 2*np.pi, 200)
plt.plot(x(aP, eP, phi), y(aP, eP, phi), label="冥王星")
plt.plot(x(aN, eN, phi), y(aN, eN, phi), label="海王星")
# x軸 y軸
plt.axhline(0, color='black', dashes=(3, 3), linewidth=0.6)
plt.axvline(0, color='black', dashes=(3, 3), linewidth=0.6)
plt.legend();
3. 原点と2つの軌道の交点を結ぶ直線を2本描く
原点 $(0, 0)$ と $(x1, y1)$ を結ぶ直線を描くのは
plt.plot([0, x1], [0, y1])
であった。
交点の1つは
- $(x(a_N, e_N, \phi_{\rm eq}), y(a_N, e_N, \phi_{\rm eq}))$，
もう一つは
- $(x(a_N, e_N, -\phi_{\rm eq}), y(a_N, e_N, -\phi_{\rm eq}))$
であるから，1本ずつ直線を引く。同じ色で描くには，
plt.plot([0, x1], [0, y1], color='black')
plt.plot([0, x2], [0, y2], color='black')
4. 海王星が冥王星より太陽から遠い期間を求める
海王星の公転周期（1周 $2\pi$ ラジアンで）$T_N$ のうち，$\displaystyle \frac{2\phi_{\rm eq}}{2 \pi}$ の間は海王星が冥王星より遠いから
$$t_1 = T_N \times \frac{\phi_{\rm eq}}{\pi}$$
で与えられる時間 $t_1$ の間だけ，海王星が冥王星より太陽から遠いことになる。以下のように…
# 海王星の公転周期は何年？
TN = ...
# round() で数値を丸める
t1 = round(TN * phieq/(np.pi))
# グラフの x=7, y=0 からはじまるテキストを書く
plt.text(7, 0, "この期間は"+str(t1)+"年")
# グラフのタイトルもつけてみよう
plt.title("冥王星が海王星より太陽に近い期間")
5. 扇形の部分を塗りつぶす
扇形を，円周の一部と原点からなる多角形だと考えて…
完成イメージ
だいたい以下のようなグラフができると思います。
