Return to Python でグラフ作成

SymPy Plotting Backends でグラフ作成：The Graphic Module 編

ここでは，SymPy Plotting Backends (SPB) の The Graphic Module を使った 2次元グラフ作成について説明します。

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 による方法は今後はアップデートされないようです。

SPB の The Graphic Module を使ってプロットできるオブジェクトは以下のとおりです。

  • 陽関数 $y = f(x)$:
    • graphics(line(f(x), (x, xmin, xmax)))
  • 陰関数 $f(x, y) = 0$:
    • graphics(implicit_2d(f(x, y), (x, xmin, xmax), (y, ymin, ymax)))
  • 媒介変数表示 $x(t), y(t)$:
    • graphics(line_parametric_2d(x(t), y(t), (t, tmin, tmax)))
  • 点，$x$ 座標 $y$ 座標の数値データ:
    • graphics(list_2d([x1, ..., xn], [y1, ..., yn]))
  • ベクトル
    • graphics(arrow_2d((x0, y0), (vx, vy)))

また，2本の陽関数で挟まれた領域を塗りつぶすこともできます。

SymPy および SPB の import

In [1]:
from sympy.abc import *
from sympy import *

# SymPy Plotting Backends (SPB)
from spb import *

# 日本語がトーフになる場合
# import japanize_matplotlib

# グラフを SVG で Notebook にインライン表示
%config InlineBackend.figure_formats = ['svg']

# デフォルト設定のため
import matplotlib.pyplot as plt
# mathtext font の設定
plt.rcParams['mathtext.fontset'] = 'cm'

# デフォルトの figsize の設定変更は，
# 一度何か描いてからにするとよいようだ。
plt.rcParams['figure.figsize'] = [6.4, 4.8]

The Graphics Module graphics()

SymPy Plotting Backends (SPB) の Graphic Module でグラフを作成する場合は graphics() を使います。

陽関数のグラフ line()

$x$ の陽関数 $y = f(x)$ を $x_{\rm min} \leq x \leq x_{\rm max}$ の範囲でグラフを作成するには line() を使って…

graphics(line(f(x), (x, xmin, xmax)))

(x, xmin, xmax) を省略すると (x, -10, 10) とみなします。

例として，$f(x) = \sin x$ のグラフを描きます。

In [2]:
plt.rcParams['figure.figsize']
Out[2]:
[6.4, 4.8]
In [3]:
graphics(line(sin(x)));

In [4]:
plt.rcParams['figure.figsize']
Out[4]:
[6.0, 4.0]

figsize の設定変更 plt.rcParams[]

SPB でデフォルトのグラフのサイズを変更するには，バックエンドの Matplotlib の設定変更方法を使います。例えば以下のように。

import matplotlib.pyplot as plt
plt.rcParams['figure.figsize'] = [6.4, 4.8]

しかし，本稿執筆時点の環境では，SPB で最初にグラフを描くとなぜか，[6.0, 4.0] に設定されてしまうようなので，あらためて [6.4, 4.8] に設定しておきます。

In [5]:
# デフォルトの figsize の設定変更は，
# 一度何か描いてからにするとよいようだ。

plt.rcParams['figure.figsize'] = [6.4, 4.8]
In [6]:
graphics(line(sin(x)));

In [7]:
plt.rcParams['figure.figsize']
Out[7]:
[6.4, 4.8]

グラフ全体のオプションの設定例

グラフ全体のオプションは graphics() のオプションとして設定します。

  • タイトル title=' '
  • 表示範囲 xlim=( , ) ylim=( , )
  • サイズ　 size =( , ) デフォルトは (6, 4)
  • 縦横比　 aspect = 'equal' 等
  • 軸ラベル xlabel=' ' ylabel=' '
  • グリッドの非表示 grid = False

なお，複数のグラフの際には凡例が表示されますが，1本のグラフの際はデフォルトでは表示されません。どうしても表示させたいなら

  • 凡例表示 legend = True 1本のグラフのときにも凡例を表示させる

ラインごとのオプション設定例

ラインごとのオプションは line() 内に rendering_kw={} として設定します。

  • 線の太さ 'lw': ('linewidth':)
  • 線の色　 'color':
  • 線種　　 'ls': ('linestyle':) ':', '--' 等

ラインごとの凡例は自動でつくが，カスタマイズしたいときは

  • 凡例ラベル label = ' ' あるいは ' ' のみでも可の場合もあり

いくつかオプションを設定する例。

In [8]:
graphics(
    line(sin(x), (x, -2*pi, 2*pi), '正弦関数',
         rendering_kw={"lw":2, "color":"red", "ls":":"}),
    legend = True, 
    xlim = (-2*pi, 2*pi), ylim = (-1.1, 1.1)
);

初等関数のグラフ

陽関数 $y = f(x)$ のグラフの描き方を使って，（双曲線関数および逆双曲線関数を含む）初等関数のグラフを描いてみます。

陽関数のグラフを重ねて表示

複数の陽関数のグラフを重ねて表示するには，

graphics(
    line(f(x), (x, xmin, xmax)), 
    line(g(x), (x, xmin, xmax)), 
    ...
)

のように。

また複数のグラフを描くと凡例が表示されます。

べき関数

まず，$y = x^{-1}$ のグラフ。

In [9]:
graphics(line(1/x), xlim = (-5, 5), ylim = (-10, 10));

$y=x^{-1}$ のグラフの $x＝0$ 付近，$y \rightarrow -\infty$ と $y \rightarrow +\infty$ を縦の直線でつないでしまう。

不連続点をつなげない対策

detect_poles=True オプションをつけ，さらに n の値を大きめにすると，以下のようになる。

In [10]:
graphics(line(1/x, detect_poles=True, n = 1e4), 
         xlim = (-5, 5), ylim = (-10, 10)
);

$y = x^{-2}, \ x^{-1}, \ x^2, \ x^3$ のグラフ例。（凡例 label に $\LaTeX$ 記法を使っています。）

In [11]:
graphics(
    line(1/x**2, label = '$x^{-2}$'), 
    line(1/x, label = '$x^{-1}$', detect_poles = True, n = 1e4), 
    line(x**2),
    line(x**3), 
    xlim = (-5, 5), ylim = (-10, 10), 
    title = 'べき関数のグラフ'
);

$\displaystyle y = \sqrt{x}, \ \frac{1}{\sqrt{x}}$ のグラフ例。
（\ で始まる $\LaTeX$ コマンドを使う際は，念のため r から始めます。）

In [12]:
graphics(
    line(sqrt(x), (x, 0, 5), r'$\sqrt{x}$'),
    line(1/sqrt(x), (x, 0.0001, 5), r'$1/\sqrt{x}$'),
    xlim = (0, 5), ylim = (0, 5)
);

指数関数

$y = e^{-x}, \ e^x$ のグラフ例。

In [13]:
graphics(
    line(exp(-x), (x, -5, 5), '$e^{-x}$'),
    line(exp(x), (x, -5, 5), '$e^{x}$'),
    xlim = (-5, 5), ylim = (0, 5)
);

対数関数

$y = \log x$ のグラフ例。

In [14]:
graphics(
    line(log(x), (x, 0.0001, 10), r'$\log x$'),
    xlim = (-0.5, 10), ylim = (-5, 5), legend = True
);

三角関数

$ y = \sin x, \ \cos x, \ \tan x $ のグラフ例。

$\tan x$ の不連続点の処理は detect_poles = True で。

In [15]:
graphics(
    line(sin(x), label = '$\sin x$'),
    line(cos(x), label = '$\cos x$'),
    line(tan(x), label = r'$\tan x$', detect_poles = True),
    xlim = (-2*pi, 2*pi+1), ylim = (-5, 5)
);

主目盛・副目盛のカスタマイズ

グラフのスタイルの詳細な設定変更は SPB には備わっていないようですので，バックエンドの Matplotlib の ax で詳細設定を行ってみます。

ここでは，横軸主目盛を $\pi$ ごと，横軸副目盛を $\frac{\pi}{2}$ ごと，縦軸副目盛を $1$ ごとにしてみます。

In [16]:
# おまじない。これで ax が使える。
fig, ax = plt.subplots()

# 横軸主目盛
Pi = float(pi)
ax.set_xticks(
    [-2*Pi,     -Pi,      0,     Pi,      2*Pi], 
    ['$-2\pi$', '$-\pi$', '$0$', '$\pi$', '$2\pi$'])

from matplotlib.ticker import MultipleLocator
# 横軸副目盛を π/4 ごとに
ax.xaxis.set_minor_locator(MultipleLocator(Pi/2))
# 縦軸副目盛を 1 ごとに
ax.yaxis.set_minor_locator(MultipleLocator(1))

graphics(
    line(sin(x), label = '$\sin x$'),
    line(cos(x), label = '$\cos x$'),
    line(tan(x), label = r'$\tan x$', detect_poles = True),
    xlim = (-2*pi, 2*pi+1), ylim = (-5, 5), 
    ax = ax # ax の設定を反映
);

参考：グリッドのカスタマイズ ax.grid()

上の出力結果をみると，主目盛と副目盛ではグリッドのスタイルが異なっているようです。このへんのカスタマイズは

  • ax.grid()
参考：$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)
In [17]:
# おまじない。これで ax が使える。
fig, ax = plt.subplots()

# 横軸主目盛
Pi = float(pi)
ax.set_xticks(
    [-2*Pi,     -Pi,      0,     Pi,      2*Pi], 
    ['$-2\pi$', '$-\pi$', '$0$', '$\pi$', '$2\pi$'])

from matplotlib.ticker import MultipleLocator
# 横軸副目盛を π/4 ごとに
ax.xaxis.set_minor_locator(MultipleLocator(Pi/2))
# 縦軸副目盛を 1 ごとに
ax.yaxis.set_minor_locator(MultipleLocator(1))

# x軸 y軸は dashed に
ax.axhline(0, c='k', ls='--', lw=0.5)
ax.axvline(0, c='k', ls='--', lw=0.5)

# グリッドは主目盛副目盛とも同じスタイルで目立たなく
ax.grid(c='lightgray', ls='--', lw=0.5, which='both')

graphics(
    line(sin(x), label = '$\sin x$'),
    line(cos(x), label = '$\cos x$'),
    line(tan(x), label = r'$\tan x$', detect_poles = True),
    xlim = (-2*pi, 2*pi+1), ylim = (-5, 5),
    grid=False, # 一旦 False にして...
    ax = ax     # ax の設定を反映
);

逆三角関数

  • $y = \sin^{-1} x = \arcsin x =$ asin(x)
    • 定義域は $\displaystyle-1 \leq x \leq 1$
    • 値域は $\displaystyle -\frac{\pi}{2} \leq y \leq \frac{\pi}{2}$
  • $y = \cos^{-1} x = \arccos x =$ acos(x)
    • 定義域は $\displaystyle-1 \leq x \leq 1$
    • 値域は $\displaystyle0 \leq y \leq \pi$
  • $y = \tan^{-1} x = \arctan x =$ atan(x)
    • 定義域は $\displaystyle-\infty < x < \infty$
    • 値域は $\displaystyle-\frac{\pi}{2} \leq y \leq \frac{\pi}{2}$
In [18]:
fig, ax = plt.subplots()

# 縦軸主目盛
Pi = float(pi)
ax.set_yticks(
    [-Pi/2,      0,     Pi/2,      Pi], 
    ['$-\pi/2$', '$0$', '$\pi/2$', '$\pi$'])
# 縦軸副目盛
ax.yaxis.set_minor_locator(MultipleLocator(Pi/4))

# 横軸主目盛
ax.xaxis.set_major_locator(MultipleLocator(5))
# 横軸副目盛
ax.xaxis.set_minor_locator(MultipleLocator(1))

# x軸 y軸は dashed に
ax.axhline(0, c='k', ls='--', lw=0.5)
ax.axvline(0, c='k', ls='--', lw=0.5)

# グリッドは主目盛副目盛とも同じスタイルで目立たなく
ax.grid(c='lightgray', ls='--', lw=0.5, which='both')

graphics(
    line(asin(x), (x, -1, 1), r'$\arcsin x$'), 
    line(acos(x), (x, -1, 1), r'$\arccos x$'), 
    line(atan(x), (x, -10, 10), r'$\arctan x$'), 
    xlim = (-7, 7), ylim = (-pi/2, pi), 
    grid=False, # 一旦 False にして...
    ax = ax     # ax の設定を反映
);

双曲線関数

$y = \sinh x, \ \cosh x, \ \tanh x$ のグラフ例。

In [19]:
fig, ax = plt.subplots()

# 縦軸主目盛
ax.yaxis.set_major_locator(MultipleLocator(2))
# 縦軸副目盛
ax.yaxis.set_minor_locator(MultipleLocator(1))

# 横軸主目盛
ax.xaxis.set_major_locator(MultipleLocator(2))
# 横軸副目盛
ax.xaxis.set_minor_locator(MultipleLocator(1))

# x軸 y軸は dashed に
ax.axhline(0, c='k', ls='--', lw=0.5)
ax.axvline(0, c='k', ls='--', lw=0.5)

# グリッドは主目盛副目盛とも同じスタイルで目立たなく
ax.grid(c='lightgray', ls='--', lw=0.5, which='both')

graphics(
    line(sinh(x), (x, -5, 5), r'$\sinh x$'), 
    line(cosh(x), (x, -5, 5), r'$\cosh x$'),
    line(tanh(x), (x, -5, 5), r'$\tanh x$'), 
    xlim = (-5, 5), ylim = (-7, 7), 
    grid=False, # 一旦 False にして...
    ax = ax     # ax の設定を反映
);

逆双曲線関数

  • $ y = \sinh^{-1} x = \mbox{arsinh}\ x = $ asinh(x)
    • 定義域は $ -\infty < x < \infty$
    • 値域は $-\infty < y < \infty$
  • $ y = \cosh^{-1} x = \mbox{arcosh}\ x = $ acosh(x)
    • 定義域は $ 1 \leq x < \infty$
    • 値域は $0 \leq y < \infty$
  • $ y = \tanh^{-1} x = \mbox{artanh}\ x = $ atanh(x)
    • 定義域は $ -1 < x < 1$
    • 値域は $-\infty < y < \infty$
In [20]:
fig, ax = plt.subplots()

# 縦軸主目盛
ax.yaxis.set_major_locator(MultipleLocator(1))

# 横軸主目盛
ax.xaxis.set_major_locator(MultipleLocator(1))

# x軸 y軸は dashed に
ax.axhline(0, c='k', ls='--', lw=0.5)
ax.axvline(0, c='k', ls='--', lw=0.5)

# グリッドは主目盛副目盛とも同じスタイルで目立たなく
ax.grid(c='lightgray', ls='--', lw=0.5, which='both')

graphics(
    line(asinh(x), (x, -10, 10), r'$\sinh^{-1} x$'), 
    line(acosh(x), (x, 1, 10), r'$\cosh^{-1} x$'),
    line(atanh(x), (x, -0.99999, 0.99999), r'$\tanh^{-1} x$'), 
    xlim = (-7, 7), ylim = (-5, 5),
    grid=False, # 一旦 False にして...
    ax = ax     # ax の設定を反映
);

陰関数のグラフ implicit_2d()

例として，$x^2 + y^2 = 1$ のグラフを描きます。$y$ について解いて陽関数の形にしたり，以下で述べるような媒介変数表示にしなくても，陰関数のまま，グラフに描くことができます。

陰関数 $f(x, y) = 0$ を (x, xmin, xmax), (y, ymin, ymax) の範囲で描く書式は…

graphics(implicit_2d(f(x, y), (x, xmin, xmax), (y, ymin, ymax)))

$f(x, y) \equiv x^2 + y^2 – 1 = 0$ としてプロットします。

In [21]:
graphics(
    implicit_2d(x**2 + y**2 - 1, 
                (x, -1.2, 1.2), (y, -1.2, 1.2)
               )
);

縦横比・サイズの設定 aspect size

縦横の比 aspect を設定して円らしく見えるようにします。size も設定してみます。

In [22]:
graphics(
    implicit_2d(x**2 + y**2 - 1, 
                (x, -1.2, 1.2), (y, -1.2, 1.2)
               ), 
    aspect = 'equal', size = (6.4, 6.4)
);

媒介変数表示のグラフ line_parametric_2d()

半径 $1$ の円の方程式は $x^2 +y^2 = 1$ です。この円を以下のような媒介変数表示にして描きます。

$$ x=\cos t, \quad y=\sin t \quad (0 \leq t \leq 2\pi) $$

書式は…

graphics(line_parametric_2d((x(t), y(t), (t, tmin, tmax)))

デフォルトで use_cm = True になっているので，カラーマップが不要なら False に。

In [23]:
graphics(
    line_parametric_2d(
        cos(t), sin(t), (t, 0, 2*pi), use_cm = False), 
    aspect = "equal", size = (6.4, 6.4)
);

点・数値データのグラフ list_2d()

graphics(list_2d([x1, ..., xn], [y1, ..., yn]))

折れ線グラフ

以下の例では，配列 X に 座標の値，配列 Y に 座標の値を入れて，6個の点をつないだ折れ線グラフを描きます。

In [24]:
X1 = [i for i in range(6)]
Y1 = [i**2 for i in range(6)]

graphics(list_2d(X1, Y1));

散布図（線でつながない）scatter = True

オプション scatter = True を設定して各点を線でつながずに，赤色 {'color':'red'} でプロットする例。

In [25]:
graphics(
    list_2d(X1, Y1, 'データ', {'color':'red'}, scatter = True), 
    legend = True
);

ファイルから数値データを読み込んでプロット

あらかじめ作成された数値データファイルからデータを読み込んでグラフを描くこともできます。

以下のようなテキストファイル mydata2.txt が（このノートブックファイルと同じフォルダに）あるとします。（授業では事前に mydata2.txt も同じフォルダにダウンロードしておきます。）

In [26]:
! cat mydata2.txt

0  0
1  1
2  4
3  9
4 16
5 25
np.loadtxt() で読み込み

数値データの書き込み・読み込みには，NumPy の savetxt() や loadtxt() を使うことにします。

In [27]:
import numpy as np

# ファイル mydata2.txt から読み込み，data に代入
data = np.loadtxt('mydata2.txt') 
data
Out[27]:
array([[ 0.,  0.],
       [ 1.,  1.],
       [ 2.,  4.],
       [ 3.,  9.],
       [ 4., 16.],
       [ 5., 25.]])

読み込んだ data の1列目 data[:,0] を $x$ 座標に，2列目 data[:,1] を $y$ 座標にして散布図を描きます。

In [28]:
# numpy.ndarray 配列であれば，以下のように一挙に列データを指定可能
# data の1列目 data[:,0] を x, 2列目 data[:,1] を y

graphics(
    list_2d(data[:,0], data[:,1], "データ", 
            {"color":"red"}, scatter = True),
    legend = True
);

ベクトルを描く arrow_2d()

点 $(x, y) =$ (x0, y0) を始点とし，成分 $\boldsymbol{v} = (v_x, v_y) = $ (vx, vy) のベクトルを描くには，

graphics(arrow_2d((x0, y0), (vx, vy)))

原点を始点とし，成分が $v_x = 0.5, v_y = 1$ のベクトルを描く例。

In [29]:
graphics(
    arrow_2d((0, 0), (0.5, 1)), 
    xlim = (-2, 2), ylim = (-2, 2), 
    aspect = 'equal', size = (6.4, 6.4),
);

複数のベクトルを描く例。斜方投射の速度ベクトルを描いてみます。

適宜無次元化を行って，斜方投射の位置 $(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}

のように書けるとします。

In [30]:
def X(t):
    return t
def Y(t):
    return t - t**2/4

def Vx(t):
    return 1
def Vy(t):
    return 1 - t/2
In [31]:
graphics(
    arrow_2d((X(0), Y(0)), (Vx(0), Vy(0)), 
             rendering_kw={'color':'b'}), 
    arrow_2d((X(1), Y(1)), (Vx(1), Vy(1)), 
             rendering_kw={'color':'b'}),
    arrow_2d((X(2), Y(2)), (Vx(2), Vy(2)), 
             rendering_kw={'color':'b'}),
    arrow_2d((X(3), Y(3)), (Vx(3), Vy(3)), 
             rendering_kw={'color':'b'}),
    arrow_2d((X(4), Y(4)), (Vx(4), Vy(4)), 
             rendering_kw={'color':'b'}),
    xlim = (-0.4, 5), ylim = (-1, 2), 
    aspect = 'equal', legend = False
);

このように，arrow_2d() では複数のベクトルを描く際に少し手間がかかります。また，好みの問題ではありますが，ベクトルを太くすると何となく丸みを帯びた矢印になります。

デフォルトの凡例もデバッグ用以外には特に必要とされる機会は少ないかと思われるし，本稿執筆時点では，フォントの問題で日本語の凡例を設定できません。

In [32]:
graphics(
    arrow_2d((0, 0), (0.5, 1)), 
    arrow_2d((1, 0), (0.5, 1), rendering_kw={'lw':5}), 
    xlim = (-1, 2), ylim = (-0.5, 2)
);

参考：Matplotlib 流にベクトルを描く ax.quiver()

以下のように Matplotlib の ax.quiver() を使ってベクトルを描くこともできます。

In [33]:
# ax.quiver() 用にリスト作成
# 始点の x, y 座標
Xlist = [X(i) for i in range(5)]
Ylist = [Y(i) for i in range(5)]
# ベクトルの成分
Vxlist = [Vx(i) for i in range(5)]
Vylist = [Vy(i) for i in range(5)]
In [34]:
# おまじない。これで ax が使える。
fig, ax = plt.subplots()

# Matplotlib でベクトルを描く
ax.quiver(Xlist, Ylist, Vxlist, Vylist, label = r'$\vec{v}$',
          # 全ベクトルに色を一括指定
          color="blue", 
          # ベクトルの形状もその気になれば細かく設定可能
          # width=0.02, 
          # headaxislength=2.5, headlength=2.8,
          # 以下の3点セットを書かないと自動スケーリングされる
          angles='xy', scale_units='xy', scale=1)

graphics(ax = ax, # ax を反映
         xlim = (-0.5, 5.5), ylim = (-1.5, 2), 
         aspect='equal', legend=True);

複数のグラフを重ねて表示

複数のグラフを重ねて表示する例を示します。

複数の陽関数を重ねて表示

陽関数 $y = x^2 – 1$ と $y = 4 x – 5$ の 2 つのグラフを重ねて描きます。

$x$ の範囲は独立に指定可能です。

In [35]:
graphics(
    line(x**2 - 1, (x, -5, 5)), 
    line(4*x - 5,  (x, 1, 3)), 
    ylim = (-5, 10)
);

数値データと理論曲線を重ねて表示

前節の数値データファイル mydata2.txt と理論曲線 $y = x^2$ の 2 つのグラフ

  • 数値データ list_2d()
  • 陽関数 line()

を重ねて表示してみます。

描く順番を指定 'zorder'

複数のグラフを重ねて表示する際，オプション 'zorder' を使って描かれる順番を調整できます。{'zorder':0} だと最初（つまり最背面）に描かれます。

In [36]:
graphics(
    list_2d(data[:,0], data[:,1], "データ", 
            {"color":"blue"}, scatter = True),
    line(x**2, (x, 0, 5), "理論曲線", {'zorder':0})
);

さまざまなグラフを重ねて表示

graphics() の強力なところは，陽関数 line() だけでなく，媒介変数表示 line_parametric_2d()，数値データ list_2d()，ベクトル arrow_2d() なども全て重ねて表示できるところです。

例えば，前掲の斜方投射

\begin{eqnarray}
x &=& t \\
y &=& t – \frac{1}{4} t^2 \\
\therefore\ \ y &=& x – \frac{1}{4} x^2
\end{eqnarray}

を使って，

  • 陽関数 line()
  • 媒介変数表示 line_parametric_2d()
  • 数値データ list_2d()
  • ベクトル arrow_2d()

を重ねて表示してみます。

arrow_2d() の凡例には日本語が使えないようです。また，label='' でも凡例は非表示にならず，show_in_legend=False とする必要があります。

In [37]:
graphics(
    # 陽関数を少し上にずらして
    line(x - x**2/4 + 0.5, (x, 0, 4), '陽関数'), 
    # 媒介変数表示
    line_parametric_2d(
        X(t), Y(t), (t, 0, 4), '媒介変数表示', use_cm = False), 
    # 数値データ
    list_2d([X(i) for i in range(5)], 
            [Y(i) for i in range(5)], '数値データ', scatter = True), 
    # ベクトル
    arrow_2d((X(0), Y(0)), (Vx(0), Vy(0)), 'vector',
             rendering_kw={'color':'b'}), 
    arrow_2d((X(1), Y(1)), (Vx(1), Vy(1)), show_in_legend=False,
             rendering_kw={'color':'b'}),
    xlim = (-0.4, 5), ylim = (-0.5, 2), 
    aspect = 'equal'
);

極座標表示のグラフ line_polar()

極座標表示 $r = r(\varphi)$ のグラフは

graphics(line_polar(r(phi), (phi, 2, 2*pi))

焦点を原点とした楕円のグラフ

太陽からの万有引力を受けて運動する惑星は，太陽（二体問題では質量中心）を焦点の一つとした楕円運動を描いて運動します。
焦点の一つを原点とし，長半径 $a$，離心率 $e$ の楕円の方程式は，極座標 $r$，$\varphi$ を使って以下のように表すことができます。

$$ r= \frac{a (1−e^2)}{ 1 + e\cos\phi}$$

冥王星，海王星の軌道

さて，かつては第9惑星，現在では準惑星の一つである冥王星も楕円軌道を描きます。

まず，焦点を原点とした極座標表示の楕円の式を関数として定義します。

In [38]:
var('a e phi')

# elliptical（楕円）の e をつけて...
def re(a, e, phi):
    return a*(1-e**2)/(1+e*cos(phi))

re(a, e, phi)
Out[38]:
$\displaystyle \frac{a \left(1 – e^{2}\right)}{e \cos{\left(\phi \right)} + 1}$

冥王星の軌道長半径 $𝑎_P=39.445 \mbox{au}$，離心率 $𝑒_P=0.250$
を使って冥王星の軌道を描きます。

In [39]:
# 冥王星の軌道長半径と離心率
aP = 39.445
eP = 0.250

Pluto = line_polar(re(aP, eP, phi), (phi, 0, 2*pi), 
               '冥王星', {'color':'green'})

graphics(
    Pluto, 
    aspect = 'equal', size = (6.4, 6.4)
);

では，この冥王星の軌道と，その内側をまわる海王星の軌道を重ねて描いてみましょう。

海王星の軌道長半径は $a_N=30.181 \mbox{au}$，離心率は $e_N=0.0097$ と小さいので簡単のために $e_N=0$ として扱います。

実際の2つの天体の軌道は同じ平面上にありませんが，太陽からの距離のみをみるために，ここでは同一平面上に描きます。

グリッドの非表示 grid = False

原点にある太陽の位置を表すために，$x$ 軸と $y$ 軸の線を引き，グリッドを非表示に。

In [40]:
# 海王星の軌道長半径と離心率
aN = 30.181
eN = 0
Neptune = line_polar(re(aN, eN, phi), (phi, 0, 2*pi), 
               '海王星', {'color':'blue'})

Sun = list_2d([0],[0], '太陽', 
              {'color':'red', 'zorder':10}, scatter = True)

# おまじない。これで ax が使える。
fig, ax = plt.subplots(figsize = (6.4, 6.4))

# x軸 y軸は dashed に
ax.axhline(0, c='k', ls='--', lw=0.5)
ax.axvline(0, c='k', ls='--', lw=0.5)

graphics(
    Pluto, Neptune, Sun,
    ax = ax, 
    aspect = 'equal', grid = False
);

参考：各種図形の表示と塗りつぶし geometry()

点 geometry(Point2D(x0, y0))

座標 $(x, y) = $ (x0, y0) に点を描く。

graphics(geometry(Point2D(x0, y0)))
In [41]:
graphics(
    geometry(Point2D(0,0), 
             label='デフォルトでは塗りつぶし'), 
    geometry(Point2D(1,0),
             label='fill=False で塗りつぶしなし', fill=False), 
    xlim = (-1, 2), ylim = (-1, 1)
);

線分 geometry(Segment())

点 (x0, y0) から (x1, y1) までの線分を描く。

graphics(geometry(Segment((x0, y0), (x1, y1))))
In [42]:
graphics(
    geometry(Segment((0, 0), (2, 2)), ''),  
    geometry(Segment((0, 1), (2, 1)), '', 
             {'color':'red', 'lw':3, 'zorder':3}) # グリッドより上に
);

多角形 geometry(Polygon())

(x0, y0) を中心とし，「半径」R の正 N 角形

graphics(geometry(Polygon((x0, y0), R, n = N)))

p1 = (x1, y1), p2, ... を頂点とした多角形

graphics(geometry(Polygon(p1, p2, ...)))
In [43]:
p1, p2, p3, p4 = (-3, -3), (1, -3), (3, 3), (-2, 2)

graphics(
    geometry(Polygon((0,0), 2, n=8), '正8角形'), 
    geometry(Polygon(p1, p2, p3, p4), 
             '頂点の座標を指定', fill=False),
    geometry(Point2D(p1)),
    geometry(Point2D(p2)),
    geometry(Point2D(p3)),
    geometry(Point2D(p4)),
    aspect = 'equal', xlim = (-3.5, 5)
);

楕円 geometry(Ellipse())

点 (x0, y0) を中心とし，長半径（横半径）a，単半径（縦半径）b の楕円

graphics(geometry(Ellipse((x0, y0), hradius = a, vradius = b)))

点 (x0, y0) を中心とし，長半径（横半径）a，離心率 e の楕円

graphics(geometry(Ellipse((x0, y0), hradius = a, eccentricity = e)))
In [44]:
graphics(
    geometry(Ellipse((-2, -2), hradius = 2, vradius = 1), 
             '$a=2, b=1$'), 
    geometry(Ellipse((2, 0), hradius = 3, eccentricity = 0.8), 
             '$a = 3, e = 0.8$', fill = False), 
    aspect = 'equal'
);

円 geometry(Circle())

点 (x0, y0) を中心とし，半径 r の円

graphics(geometry(Circle((x0, y0), r)))
枠線色 'edgecolor' と背景色 'facecolor'

Circle() に限らず，全ての図形の枠線の色を 'edgecolor' で，塗り潰し（背景）の色を 'facecolor' でそれぞれ独立に設定できます。

In [45]:
graphics(
    geometry(Circle((1, 0), 2), 
             rendering_kw={'lw':2, 
                           'edgecolor':'black', 
                           'facecolor':'yellow'}), 
    aspect = 'equal'
);

領域の塗りつぶし implicit_2d()

2本の陽関数で挟まれた領域を塗りつぶす

ここでは，implicit_2d() を使って，$x$ 軸（$y = 0$）と $y = f(x)$ で挟まれた領域を塗りつぶしてみます。

$0.5 \leq x \leq 2$ で $0 \leq y \leq f(x)$ の領域を塗りつぶす例。

In [46]:
def f(x):
    return 0.6*x + 0.4*cos(x)
In [47]:
graphics(
    implicit_2d(
        # 領域を示す不等式
        (0 <= y) & (y <=f(x)), 
        # 塗りつぶす x の範囲
        (x, 0.5, 2), 
        # 不等式が満たされているかを探す y の範囲
        (y, -0.1, 1.2), 
        adaptive=True, # これをつけないと警告が出る
        # 塗りつぶしの色
        color='yellow'
    ),
    xlim = (-0.2, 2.5), ylim = (-0.2, 1.5)
);

扇形を描いて塗りつぶす

まず，扇形を，円周の一部と原点からの2本の直線で囲まれた領域と考えて縁のグラフを描く。

In [48]:
# 外枠
Arc1 = line_parametric_2d(
            cos(t), sin(t), (t, pi/6, pi/3), '',
            use_cm = False, 
            rendering_kw = {'color':'black'})
Line1 = geometry(Segment((0, 0), (cos(pi/6), sin(pi/6))), 
            label = '',
            rendering_kw = {'color':'black'})
Line2 = geometry(Segment((0, 0), (cos(pi/3), sin(pi/3))), 
            label = '',
            rendering_kw = {'color':'black'})

graphics(
    Arc1, Line1, Line2,
    aspect = 'equal', 
    xlim = (-0.1, 1), ylim = (-0.1, 1)
);

扇形の内部を，

\begin{eqnarray}
x^2 + y^2 &\leq& 1 \\
y &\geq& \tan\frac{\pi}{6} \ x \\
y &\leq& \tan\frac{\pi}{3} \ x
\end{eqnarray}

を満たす領域として，implicit_2d() で塗りつぶします。

In [49]:
FilledArea1 = implicit_2d(
        ((x**2+y**2 <= 1) 
          & (tan(pi/6)*x <= y)
          & (y <= tan(pi/3)*x)), 
        (x, 0, 1), (y, 0, 1), 
        adaptive = True, 
        color = 'yellow') 

graphics(
    FilledArea1, 
    aspect = 'equal', 
    xlim = (-0.1, 1), ylim = (-0.1, 1), legend = False
);

枠線も重ねて表示。

In [50]:
graphics(
    # 扇形塗りつぶし部分
    FilledArea1, 
    # 枠線部分
    Arc1, Line1, Line2,
    aspect = 'equal', 
    xlim = (-0.1, 1), ylim = (-0.1, 1) , legend = False
);

○練習：ドッテンメイカイを塗りつぶし

冥王星と海王星の軌道のグラフをみると，軌道が交差している。これは，ある期間冥王星のほうが海王星より太陽に近いことを意味する。そこで，以下の指示に従ってグラフを作成せよ。

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) をnsolve() で数値的に求める。SymPy による方程式の数値解の求め方については，以下を参照：

In [51]:
# 楕円の極座標表示
def re(a, e, phi):
    return a*(1-e**2)/(1+e*cos(phi))
In [52]:
# 冥王星 Pluto の軌道長半径，離心率
aP = 39.445
eP = 0.250
# 海王星 Neptune の軌道半径
aN = 30.181
eN = 0
In [53]:
# nsolve() を使って数値的に解を求める
phi_eq = nsolve(Eq(re(aP, eP, phi), aN), phi, (0, float(pi)/2))
phi_eq
Out[53]:
$\displaystyle 0.448599704341525$
In [54]:
# もう一つの解は...
nsolve(Eq(re(aP, eP, phi), aN), phi, (-float(pi)/2, 0))
Out[54]:
$\displaystyle -0.448599704341525$

軌道の対称性から，$- \phi_{\rm eq} < \phi < + \phi_{\rm eq}$ の間は冥王星のほうが海王星より太陽に近いことになる。

2. 原点と2つの軌道の交点を結ぶ直線を2本定義

交点の1つは

  • $(a_N \cos\phi_{\rm eq}, a_N \sin \phi_{\rm eq})$，

もう一つの交点は

  • $(a_N \cos\phi_{\rm eq}, -a_N \sin \phi_{\rm eq})$

であるから，原点と交点を結ぶ直線（線分）を1本ずつ，計2本定義する。

In [55]:
LineA = geometry(Segment((0, 0), 
                         (aN*cos(phi_eq), aN*sin(phi_eq))), '',
                 rendering_kw = {'color':'black'})
LineB = geometry(Segment((0, 0), 
                         (aN*cos(phi_eq), -aN*sin(phi_eq))), '',
                 rendering_kw = {'color':'black'})
3. 塗りつぶす扇形の部分を定義

最後に，海王星のほうが冥王星より太陽から遠い期間の扇形の部分を，塗りつぶす領域として定義します。

In [56]:
FilledArea = implicit_2d(
        ((x**2+y**2 <= aN**2) 
         & (tan(-phi_eq)*x <= y) 
         & (y <= tan(phi_eq)*x)), 
        (x, 0, aN), (y, -aN, aN), 
        adaptive = True, 
        color = 'yellow', label = '海王星が遠い期間')
In [57]:
# おまじない。これで ax が使える。
fig, ax = plt.subplots(figsize = (6.4, 6.4))

# x軸 y軸は dashed に
ax.axhline(0, c='k', ls='--', lw=0.5)
ax.axvline(0, c='k', ls='--', lw=0.5)

graphics(
    FilledArea, 
    Pluto, Neptune, Sun, LineA, LineB,
    ax = ax, 
    aspect = 'equal', grid = False
);

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 * phi_eq/(pi))

# グラフに表示させる準備
mytext = "この期間は"+str(t1)+"年"

# グラフの x=10, y=0.5 からはじまるテキストを書く
ax.text(10, 0.5, mytext)

あるいはグラフのタイトルに書いてもよいだろう。

title = '海王星が冥王星よりも太陽から遠い期間は??年'

In [58]:
In [59]:
 
完成イメージ
だいたい以下のようなグラフができると思います。
月別平年気温のカーブフィット
弘前市の月別平年気温のデータを使って配列 HiroT を作ります。
In [60]:
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]

In [61]:
graphics(
    list_2d(Month, HiroT, '月別平年気温', scatter = True), 
    legend = True
);

最小二乗法によるカーブフィット scipy.optimize.curve_fit()
グラフをみると，月別平年基本は 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.
In [62]:
# fit される関数を 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)

In [63]:
# scipy.optimize.curve_fit 用に NumPy の 関数に変換
args = (x, a, b, c)
func = lambdify(args, f(x, a, b, c), 'numpy')

In [64]:
from scipy.optimize import curve_fit

popt, pcov = curve_fit(func, Month, HiroT)
popt

Out[64]:
array([10.53333333, -8.34025527, -9.16106713])
最小二乗法によって求めた値を $a, b, c$ に代入して $f(x, a, b, c)$ のグラフを重ねてプロットしてみます。
In [65]:
# 最小二乗法によって求めた値を a, b, c に代入
a, b, c = popt

# せっかくだから横軸目盛を調整
fig, ax = plt.subplots()

from matplotlib.ticker import MultipleLocator
ax.xaxis.set_major_locator(MultipleLocator(1))

graphics(
    list_2d(Month, HiroT, '月別平年気温', scatter = True), 
    line(f(x, a, b, c), (x, 1, 12), '関数フィット'), 
    ax = ax, # ax の設定を反映
    title = '弘前市の月別平年気温', 
    xlabel = '月', ylabel = '°C'
);

○練習：カーブフィット
弘前市の月別平年気温のデータを使って
  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}
  2. 月別平年気温の数値データとフィット曲線 $f(x)$ および $g(x)$ のグラフを重ねて描け。
  3. 月平年気温の平均値を求めてみよ。（$g(x)$ の定数部分 $A_0$ が平均値。）
In [66]:
# 変数の宣言
var('x A0 a1 a2 b1 b2')

# SymPy の symbolic な関数として定義
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)

g(x, A0, a1, a2, b1, b2)

Out[66]:
$\displaystyle A_{0} + a_{1} \cos{\left(\frac{\pi x}{6} \right)} + a_{2} \cos{\left(\frac{\pi x}{3} \right)} + b_{1} \sin{\left(\frac{\pi x}{6} \right)} + b_{2} \sin{\left(\frac{\pi x}{3} \right)}$
In [ ]:
In [ ]:
In [ ]:
なお，3. の別解として，リストの総和 sum() と要素数 len() を使って以下のように平均値を求めることもできましたね。
In [67]:
sum(HiroT)/len(HiroT)

Out[67]:
10.533333333333333
また，月平年気温の最大値や最小値も簡単に求められますね。
In [ ]: