ここでは，SymPy Plotting Backends (SPB) の Plot functions を使った 2次元グラフ作成について説明します。
Python の SymPy 自身にもグラフ作成機能はありますが，SymPy Plotting Backends (SPB) は，言わばその機能拡張版に相当します。 SPB によるグラフ作成には，The Graphics Module による方法と，Plot functions による方法の2種類の方法があります。ドキュメントによると，The Graphics Module のほうがより直感的で簡単であり，異なったグラフを重ねて表示させることも簡単にできるのでオススメのようです。Plot functions による方法は今後はアップデートされないようです。
以下のページの内容とほぼ同等。
Python を使って，関数のグラフを描くことができます。また，数値データもグラフにすることができます。
ここでは，SymPy Plotting Backends (SPB) の Plot functions を使った 2次元グラフ作成について説明します。
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.abc import *
from sympy import *
# SymPy Plotting Backends (SPB)
from spb import *
# グラフを SVG で Notebook にインライン表示
%config InlineBackend.figure_formats = ['svg']
#
# init_printing()
陽関数のグラフ
例として，$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 では，あらかじめ作成された数値データファイルを読み込んでグラフを描くこともできます。
まず，数値データファイル data.txt をあらかじめ作成しておきます。数値データの書き込み・読み込みには，NumPy の savetxt() や loadtxt() を使うことにします。
# ファイルに書き込むデータを作成
dat = []
for i in range(6):
dat.append([i, i**2])
dat
import numpy as np
# ファイル data.txt に dat を書き込む
np.savetxt('data.txt', dat, fmt='%i') # 整数として書き出し
# ファイルの内容表示
!cat data.txt
# ファイル data.txt から読み込み，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
plot_list(data[:,0], data[:,1], "データ", {"color":"red"},
legend = True, is_point = True);
数値データと理論曲線を重ねて表示
前節の数値データファイル data.txt と理論曲線 $y = x^2$ の 2 つのグラフを重ねて表示してみます。
p1 = plot_list(data[:,0], data[:,1], "データ", {"color":"red"},
is_point = True, legend = True, show = False)
p2 = plot(x**2, (x, 0, 6), "理論曲線 $y=x^2$", {"color":"blue"}, 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)] # 12 までのときは 13 と書く
Month
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)]);
練習 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$ が平均値。）
# SymPy の symbolic な関数として定義
var('x A0 a1 a2 b1 b2')
def g(x, A0, a1, a2, b1, b2):
return A0+a1*cos(pi*x/6)+a2*cos(pi*x/3)+b1*sin(pi*x/6)+b2*sin(pi*x/3)
領域の塗りつぶし
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 を両方設定しようとしてもうまくいかないようだ。
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();
練習 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 による方程式の数値解の求め方については，以下を参照：
# 楕円の極座標表示
def r(a, e, phi):
return a*(1-e**2)/(1+e*cos(phi))
# x 座標
def x(a, e, phi):
return r(a, e, phi)*cos(phi)
# y 座標
def y(a, e, phi):
return r(a, e, phi)*sin(phi)
# 冥王星 Pluto の軌道長半径，離心率
aP = 39.767
eP = 0.254
# 海王星 Neptune の軌道半径
aN = 30.1104
eN = 0
# nsolve() を使って数値的に解を求める
phieq = nsolve(Eq(r(aP, eP, phi), aN), phi, (0, float(pi)/2))
# もう一つの解は...
# nsolve(Eq(r(aP, eP, phi), aN), phi, (-float(pi)/2, 0))
軌道の対称性から，$- \phi_{\rm eq} < \phi < + \phi_{\rm eq}$ の間は冥王星のほうが海王星より太陽に近いことになる。
2. 海王星と冥王星の軌道のグラフを描く
3. 原点と2つの軌道の交点を結ぶ直線を2本描く
原点 $(0, 0)$ と $(x1, y1)$ を結ぶ直線を描くのは
plot_list([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本ずつ直線を引く。原点 $(0, 0)$ と $(x2, y2)$ を結ぶ直線も一緒に同じ色で描くには，
plot_list(
([0, x1], [0, y1], {'color':'black'}),
([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$ の間だけ，海王星が冥王星より太陽から遠いことになる。
以下のようにして $t_1$ (t1) を求める。
# 海王星の公転周期は何年？
TN = ...
# round() で数値を丸める
t1 = round(TN * phieq/(pi))
# グラフに表示させる準備
mytext = "この期間は"+str(t1)+"年"
# グラフの x=7, y=0 からはじまるテキストを書く
ax.text(7, 0, mytext)
5. 扇形の部分を塗りつぶす
最後に，冥王星のほうが海王星より太陽に近い期間の扇形の部分を塗りつぶし，グラフを完成させてください。
完成イメージ
だいたい以下のようなグラフができると思います。
