Python を使って，関数のグラフを描くことができます。また，数値データもグラフにすることができます。
Python の2次元グラフ作成は matplotlib.pyplot.plot() でもできますが，ここでは，SymPy Plotting Backends (SPB) について説明します。
SPB でプロットできるオブジェクトは以下のとおりです。
- 陽関数 $y = f(x)$:
plot(f(x), (x, xmin, xmax))
- 陰関数 $f(x, y) = 0$:
plot_implicit(f(x, y), (x, xmin, xmax), (y, ymin, ymax))
- 媒介変数表示 $x(t), y(t)$:
plot_parametric(x(t), y(t), (t, tmin, tmax))
- 点，$x$ 座標 $y$ 座標の数値データ:
plot_list([x1, ..., xn], [y1, ..., yn])
- ベクトル
ax.quiver(X, Y, Vx, Vy)
また，2本の陽関数で挟まれた領域などを塗りつぶすこともできます。
from sympy import *
# 1文字変数の Symbol の定義が省略できる
from sympy.abc import *
# π
from sympy import pi
Pi = float(pi)
# SymPy Plotting Backends (SPB)
from spb import *
# グラフを SVG で Notebook にインライン表示
%config InlineBackend.figure_formats = ['svg']
陽関数のグラフ
例として，$y = \sin x$ のグラフを $−2\pi \leq x \leq 2\pi$ の範囲で描きます。基本的な定数の一つである円周率 $\pi$ は SPB では pi と書きます。
plot(sin(x), (x, -2*pi, 2*pi));
いくつかオプションを設定する例。
plot(sin(x), (x, -2*pi, 2*pi),
# 凡例, 線の太さと色
"正弦関数", {"linewidth":2, "color":"red"},
# 表示範囲
xlim = (-10, 10), ylim = (-1.1, 1.1), legend = True);
陰関数のグラフ
例として，$x^2 + y^2 = 1$ のグラフを描きます。$y$ について解いて陽関数の形にしたり，以下で述べるような媒介変数表示にしなくても，陰関数のまま，グラフに描くことができます。
陰関数 $f(x, y) = 0$ を (x, xmin, xmax), (y, ymin, ymax) の範囲で描く書式は…
plot_implicit(f(x, y), (x, xmin, xmax), (y, ymin, ymax))
$f(x, y) \equiv x^2 + y^2 – 1$ としてプロットします。
plot_implicit(x**2 + y**2 - 1,
(x, -1.2, 1.2), (y, -1.2, 1.2));
縦横比・サイズの設定
縦横の比 aspect を設定して円らしく見えるようにします。size も設定してみます。
plot_implicit(x**2 + y**2 - 1, (x, -1.2, 1.2), (y, -1.2, 1.2),
aspect = "equal", size = (6, 6));
媒介変数表示のグラフ
半径 $1$ の円の方程式は $x^2 +y^2 = 1$ です。この円を以下のような媒介変数表示にして描きます。
$$ x=\cos t, \quad y=\sin t \quad (0 \leq t \leq 2\pi) $$
書式は…
plot_parametric(x(t), y(t), (t, tmin, tmax))
デフォルトで use_cm = True になっているので，カラーマップが不要なら False に。
plot_parametric(cos(t), sin(t), (t, 0, 2*pi),
aspect = "equal", use_cm = False);
点・数値データのグラフ
以下の例では，配列 X に 座標の値，配列 Y に 座標の値を入れて，6個の点をつないだ折れ線グラフを描きます。
X = [i for i in range(6)]
Y = [i**2 for i in range(6)]
plot_list(X, Y);
オプション is_point = True を設定して，各点を線でつながずに赤色でプロットする例。
plot_list(X, Y, "データ", {"color":"red"}, is_point = True, legend = True);
ベクトルを描く
SPB には plot_vector() がありますが，これはベクトル場を描く関数です。以下のように一旦 plot() すれば matplotlib.pyplot.quiver() が使えるので，Matplotlib 流にベクトルを描いてみます。
原点を始点とし，成分が $v_x = 0.5, v_y = 1$ のベクトルを描く例。
# 始点の x 座標
X = [0]
# 始点の y 座標
Y = [0]
# ベクトルの x 成分
Vx = [0.5]
# ベクトルの y 成分
Vy = [1]
p=plot(xlim = (-2, 2), ylim = (-2, 2),
aspect = "equal", size = (5, 5), show=False)
ax=p.ax
# ax を定義したら，あとは Matplotlib 流に
# x の目盛を 0.5 刻みに
ax.set_xticks([0.5*i for i in range(-4,5)])
# y の目盛を 0.5 刻みに
ax.set_yticks([0.5*i for i in range(-4,5)])
# x軸 y軸
ax.axhline(0, color='black', dashes=(3, 3), linewidth=0.7)
ax.axvline(0, color='black', dashes=(3, 3), linewidth=0.7);
# ベクトル
ax.quiver(X, Y, Vx, Vy, color="blue",
# 以下の3点セットを書かないと自動スケーリングされる
angles='xy', scale_units='xy', scale=1);
複数のベクトルを描く例。斜方投射の速度ベクトルを描いてみます。
X = [t for t in range(5)]
Y = [(t-0.25*t**2) for t in range(5)]
Vx = [0.5 for t in range(5)]
Vy = [(0.5-0.25*t) for t in range(5)]
X, Y, Vx, Vy
p = plot(xlim = (-0.5, 4.5), ylim = (-0.5, 2),
size=(6,3), ylabel = "$y$", show=False)
ax = p.ax
ax.quiver(X, Y, Vx, Vy, color="blue",
# 以下の3点セットを書かないと自動スケーリングされる
angles='xy', scale_units='xy', scale=1);
複数のグラフを重ねて表示
複数のグラフを重ねて表示する例を示します。
$y = x^2 – 1$ と $y = 4 x – 5$ の 2 つのグラフを重ねて描きます。
$x$ の範囲は $-5 \leq x \leq 5$ です。
var('x')
plot(x**2 - 1, 4*x - 5, (x, -5, 5));
それぞれの曲線ごとに，いくつかオプションを設定して描く例です。
plot((x**2 - 1, {"linewidth":2, "color":"red"}),
(4*x - 5, {"linewidth":1, "color":"black"}), (x, -5, 5),
# 表示範囲
ylim = (-5, 10));
数値データファイルを読み込む
SPB では，あらかじめ作成された数値データファイルを読み込んでグラフを描くこともできます。
# ファイルを作成
import numpy as np
dat = []
for i in range(6):
dat.append([i, i**2])
np.savetxt('data.txt', dat, fmt='%i') # 整数として書き出し
# ファイルの内容表示
!cat data.txt
data = np.loadtxt('data.txt', dtype = 'int') # 整数として読み込み
# data の1列目 data[:,0] を x, 2列目 data[:,1] を y
plot_list(data[:,0], data[:,1], "データ", legend = True,
is_point = True, line_color = "red");
数値データと理論曲線を重ねて表示
前節の数値データファイル data.txt と理論曲線 $y = x^2$ の 2 つのグラフを重ねて表示してみます。
p1 = plot_list(data[:,0], data[:,1], "データ",
is_point = True, line_color = "red", legend = True,
show = False)
p2 = plot(x**2, (x, 0, 6), "理論曲線 $y=x^2$", show = False)
(p1 + p2).show();
海王星と冥王星の軌道
焦点を原点とした楕円のグラフ
太陽からの万有引力を受けて運動する惑星は，太陽（二体問題では質量中心）を焦点の一つとした楕円運動を描いて運動します。
焦点の一つを原点とし，長半径 $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*cos(phi))
def x(a, e, phi):
return r(a, e, phi)*cos(phi)
def y(a, e, phi):
return r(a, e, phi)*sin(phi)
aP = 39.445
eP = 0.250
plot_parametric(
x(aP, eP, phi), y(aP, eP, phi), (phi, 0, 2*pi), "冥王星",
aspect = "equal", legend = True, use_cm = False);
では，この冥王星の軌道と，その内側をまわる海王星の軌道を重ねて描いてみましょう。
海王星の軌道長半径は $a_N=30.181 \mbox{au}$，離心率は $e_N=0.0097$ と小さいので簡単のために $e_N=0$ として扱います。
実際の2つの天体の軌道は同じ平面上にありませんが，太陽からの距離のみをみるために，ここでは同一平面上に描きます。
aN = 30.181
eN = 0
plot_parametric(
(x(aP, eP, phi), y(aP, eP, phi), "冥王星"),
(x(aN, eN, phi), y(aN, eN, phi), "海王星"),
(phi, 0, 2*pi),
aspect = "equal", legend = True, use_cm = False);
参考：極座標表示でプロット
極座標表示でプロットする場合は，plot_polar() を使います。（plot_polar() の場合はデフォルトで aspect = "equal" になっているようです。）
x 軸 y 軸を引く
また，Matplotlib 業界では $x$ 軸や $y$ 軸を引くという慣習がないのか，そのようなオプションが見当たりません。仕方がないので，axhline() と axvline() を使ってそれらしい線を引いてみます。
グリッドの非表示
grid = False でグリッドを非表示にしてみます。
p = plot_polar(
r(aP, eP, phi), "冥王星",
r(aN, eN, phi), "海王星", (phi, 0, 2*pi),
xlim = (-60, 40), ylim = (-50, 50), size = (6, 6),
grid = False, show = False);
ax = p.ax
# x軸 y軸
ax.axhline(0, color='black', dashes=(3, 3), linewidth=0.6)
ax.axvline(0, color='black', dashes=(3, 3), linewidth=0.6);
月別平年気温の関数フィット
弘前市の月別平年気温のデータを使って配列 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]
p = plot_list(Month, HiroT, "月別平年気温",
is_point = True, legend = True);
グラフをみると，月別平年基本は 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.# SymPy の symbolic な関数として定義
var('x a b c')
def f(x, a, b, c):
return a + b*sin(pi*x/6) + c*cos(pi*x/6)
# scipy.optimize.curve_fit 用に NumPy の 関数に変換
args = (x, a, b, c)
func = lambdify(args, f(x, a, b, c), 'numpy')
from scipy.optimize import curve_fit
popt, pcov = curve_fit(func, Month, HiroT)
popt
最小二乗法によって求めた $a, b, c$ の値をもった $f(x, a, b, c)$ のグラフを重ねてプロットしてみます。
p1 = plot_list(Month, HiroT, "月別平年気温",
ylim = (-5, 25),
xlabel = "月", ylabel = "°C",
is_point = True, show = False);
a, b, c = popt
p2 = plot(f(x, a, b, c), (x, 1, 12), "関数フィット",
show = False)
p = p1 + p2
ax = p.ax
# x の主目盛を 1 刻みに
ax.set_xticks([i for i in range(1, 13)]);
領域の塗りつぶし
2本の陽関数で挟まれた領域を塗りつぶす
ここでは，plot_implicit() を使って，$x$ 軸（$y = 0$）と $y = f(x)$ で挟まれた領域を塗りつぶしてみます。
$0.5 \leq x \leq 2$ で $0 \leq y \leq f(x)$ の領域を塗りつぶす例。
def f(x):
return 0.6*x + 0.4*cos(x)
# 0 <= y <= f(x) の領域を 0.5 <= x <= 2 の範囲で塗りつぶす
# plot_implicit() の場合は x の範囲と y の範囲が必要
var('x y')
plot_implicit(
(0 <= y) & (y <= f(x)),
# 塗りつぶす x の範囲
(x, 0.5, 2),
# 不等式が満たされているかを探す y の範囲
(y, -0.1, 1.2),
adaptive=True, # これをつけないと警告が出る
# グラフの表示範囲
xlim = (-0.2, 2.5), ylim = (-0.2, 1.5),
# 塗りつぶしの色
color='yellow'
);
いくつかオプションを指定して描く例。
$x$ 軸 $y$ 軸を表示します。
デフォルトのグリッド線が少し目障りな場合は，以下のように細めの灰色点線にして目立たないようにすることもできます。
# y = f(x) を 0.25 <= x <= 2.25 の範囲で描く
p1 = plot(
f(x), (x, 0.25, 2.25),
# 凡例。$ で囲まれた部分は latex 記法
"$f(x)$",
xlim = (-0.2, 2.5), ylim = (-0.2, 1.5),
show = False # 後でまとめて表示
)
# 0 <= y <= f(x) の領域を 0.5 <= x <= 2 の範囲で塗りつぶす
# plot_implicit() の場合は x の範囲と y の範囲が必要
p2 = plot_implicit(
(0 <= y) & (y <=f(x)),
(x, 0.5, 2), (y, -0.2, 1.5),
# 凡例。$ で囲まれた部分は latex 記法
"$0 \leq y \leq f(x)$",
adaptive=True, # これをつけないと警告が出る
xlim = (-0.2, 2.5), ylim = (-0.2, 1.5), color='yellow',
show = False # 後でまとめて表示
)
# 2つのグラフを重ねて
p = p1 + p2
ax = p.ax
# グリッドを細めの灰色点線で
ax.grid(which="major", color="lightgray", dashes=(3, 3), linewidth=0.5)
# x軸 y軸
ax.axhline(0, color='black', dashes=(3, 3), linewidth=0.6)
ax.axvline(0, color='black', dashes=(3, 3), linewidth=0.6);
多角形の内部を塗りつぶす
plot_geometry(Polygon()) を使って，$(0, 0)$ を中心とし，「半径」（中心から頂点までの距離）が $2$ の正四角形 $n=4$（正方形）を描く例。
plot_geometry(
# 中心， 半径，頂点の数
Polygon((0, 0), 2, n=4),
# グリッドは無しに
grid = False,
# 塗りつぶさない
is_filled=False
);
# デフォルトでは内部を塗りつぶす
plot_geometry(
# 中心， 半径，頂点の数
Polygon((0, 0), 2, n=4),
grid = False
);
# 頂点の x, y 座標を指定して描く例
p1, p2, p3, p4 = [(0, 0), (1, 0), (1, 1), (0, 1)]
plot_geometry(
Polygon(p1, p2, p3, p4),
# 塗りつぶしの色
{'color':'yellow'}
);
円の内部を塗りつぶす
plot_geometry(Circle()) を使って，$(0, 0)$ を中心とし，半径が $1$ の円を描く。
# 原点 (0, 0) を中心とした半径 1 の円
plot_geometry(
Circle((0, 0), 1),
grid = False,
# 塗りつぶさない
is_filled=False
);
# デフォルトでは内部を揺りつぶす
plot_geometry(
Circle((0, 0), 1),
# 縁の色を赤に
{'edgecolor':'red'},
grid = False
);
color と edgecolor の同時設定
本稿執筆時点では，以下のように edgecolor と color を両方設定しようとしてもうまくいかないようだ。（追記：{'edgecolor':'red', 'facecolor':'yellow'} でうまくいきますよ。）
plot_geometry(
Circle((0, 0), 1), {'edgecolor':'red', 'color':'yellow'}
);
どうしてもというのであれば，以下のように別々に作成して重ねて表示させるという方法があるだろう。
# 円のみを赤色の線で
en1 = plot_geometry(
Circle((0, 0), 1), {'color':'red'},
grid = False, legend = False, is_filled=False,
show=False
)
# 円の内部を黄色で塗りつぶす
en2 = plot_geometry(
Circle((0, 0), 1), {'color':'yellow'},
grid = False, legend = False,
show=False
)
# 重ねて表示
(en1 + en2).show();
円の内部が $ x^2 + y^2 \leq 1$ であることから，plot_implicit() を使って塗りつぶすこともできます。
plot_implicit(
x**2 + y**2 <=1, (x, -1, 1), (y, -1, 1),
# color は塗りつぶす色，border_color は縁の線の色
color="yellow", border_color="red",
xlim = (-1.1, 1.1), ylim = (-1.1, 1.1),
aspect = "equal", grid = False);
扇形の内部を塗りつぶす
まず扇形を，円周の一部と原点からの2本の直線で囲まれた領域と考えて縁のグラフを描く。
# 円の一部を媒介変数表示で
p1 = plot_parametric(
cos(t), sin(t), (t, pi/6, pi/3), {"color":"black"},
xlim = (-1.1, 1.1), ylim = (-1.1, 1.1),
grid = False, legend = False,
aspect = "equal", use_cm = False, show = False)
# 2本の直線
p2 = plot_list(
([0, cos(pi/6)], [0, sin(pi/6)], {"color":"black"}),
([0, cos(pi/3)], [0, sin(pi/3)], {"color":"black"}), show = False)
# 重ねて表示
(p1+p2).show();
扇形の内部を，
\begin{eqnarray}
x^2 + y^2 &\leq& 1 \\
y &\geq& \tan\frac{\pi}{6} \ x \\
y &\leq& \tan\frac{\pi}{3} \ x
\end{eqnarray}
を満たす領域として，plot_implicit() で塗りつぶします。
p3 = plot_implicit(
(x**2+y**2 <= 1) & (tan(pi/6)*x <= y) & (y <= tan(pi/3)*x),
(x, 0, 1), (y, 0, 1), legend = False,
color = "yellow",
adaptive=True, # これをつけないと警告が出る
# グラフの表示範囲
xlim = (-1.1, 1.1), ylim = (-1.1, 1.1),
aspect = "equal", grid = False
);
縁の線も重ねて表示。
(p1+p2+p3).show();
