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, "データ",
is_point = True, line_color = "red", 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).save("spbplot15.svg");
参考：極座標表示でプロット
極座標表示でプロットする場合は，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)]);
領域の塗りつぶし
SPB には plot_parametric_region() というのがありまずが，開発中らしく，本稿執筆時点では fill ができないようです。ここでは，plot() の fill = dict オプションを使って塗りつぶしてみます。
$0.5 \leq x \leq 2$ で $y = f(x)$ と $x$ 軸（$y = 0$）に囲まれた領域を塗りつぶす例。$x$ 軸 $y$ 軸を表示します。
デフォルトのグリッド線が少し目障りな場合は，grid = False で非表示にすることもできますし，以下のように細めの灰色点線にして目立たないようにすることもできます。
def f(x):
return 0.6*x + 0.4*cos(x)
x_list = np.linspace(0.5, 2, 100)
y1_list = lambdify(x, f(x))(x_list)
p = plot(f(x), (x, 0.25, 2.25), "$f(x)$", grid = False,
xlim = (-0.2, 2.5), ylim = (-0.2, 1.5),
fill={'x': x_list,
'y1':y1_list,
'y2':0,
'color':'yellow'}, show = False);
ax = p.ax
# グリッドを細めの灰色点線で
ax.grid(which="major", color="lightgray", dashes=(3, 5), linewidth=0.5)
# x軸 y軸
ax.axhline(0, color='black', dashes=(3, 5), linewidth=0.6)
ax.axvline(0, color='black', dashes=(3, 5), linewidth=0.6);
