Jupyter notebook の記事を投稿してみた
2018-04-09

せっかく Jupyter の記事を投稿できるので試してみました。

Miyadaiku と Jupyter でこういう記事がかけるんだーくらいに思ってもらえれば御の字です。

Note

ありきたりな例では面白くないので、今回はサイクロイド曲線を題材にしていろいろやってみます。なつかしい。

図形描画 figure figure

描画ライブラリは色々ありますが、割と一般的な matplotlib を使います。

In [1]:
%autosave 0
%matplotlib inline
Autosave disabled
In [2]:
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
In [3]:
# x = r(θ - sinθ)
# y = r(1 - cosθ)
# となるような直交座標を計算する
In [4]:
data = np.linspace(0, np.pi * 2, 100) 

xs = []
ys = []

r = 1.0

for i in data:
    xs.append(r * (i - np.sin(i)))
    ys.append(r * (1 - np.cos(i)))

xs = np.array(xs)
ys = np.array(ys)

plt.gca().set_aspect('equal', adjustable='box')
plt.plot(xs, ys)
plt.fill_between(xs, ys, where=xs <= 3)
plt.show()
In [5]:
# 続いて X 軸を軸として回転させたもの
fig = plt.figure()
ax = fig.add_subplot(1,1,1,projection='3d')

r = 1
φ, θ = np.mgrid[0:2*np.pi:20j, 0:2*np.pi:20j]

x = r * (θ - np.sin(θ))
y = r * (1 - np.cos(θ)) * np.sin(φ)
z = r * (1 - np.cos(θ)) * np.cos(φ)

ax.plot_wireframe(x, y, z, color=('red', 'green', 'blue'))
plt.show()

面積と体積 surface_volume surface_volume

SymPy も使ってみたいので、(上記で描画した図形の)面積と体積を求めます。

In [3]:
%autosave 0
Autosave disabled
In [4]:
from sympy import sin, cos, symbols, integrate, Integral, diff, Eq, pi, init_printing
init_printing()
S, V, r, x, y, θ, dx, , a = symbols('S V r x y θ dx dθ a')
In [5]:
# 計算方法のおさらい
Eq(x, r(θ - sin(θ))), Eq(y, r(1 - cos(θ)))
Out[5]:
$$\left ( x = r{\left (θ - \sin{\left (θ \right )} \right )}, \quad y = r{\left (- \cos{\left (θ \right )} + 1 \right )}\right )$$
In [6]:
# まずは 面積から
Eq(S, Integral(y, x))
Out[6]:
$$S = \int y\, dx$$
In [7]:
# 上記の条件から dx と dθ の関係を求める
Eq(dx/, diff(r * (θ - sin(θ)), θ))
Out[7]:
$$\frac{dx}{dθ} = r \left(- \cos{\left (θ \right )} + 1\right)$$
In [8]:
# 変形
Eq(dx, r(1 - cos(θ)) * )
Out[8]:
$$dx = dθ r{\left (- \cos{\left (θ \right )} + 1 \right )}$$
In [9]:
# y と等しいので以下のようにまとめられる
Eq(S, r ** 2 * Integral((1 - cos(θ))**2, (θ, 0, a)))
Out[9]:
$$S = r^{2} \int_{0}^{a} \left(- \cos{\left (θ \right )} + 1\right)^{2}\, dθ$$
In [10]:
# 定積分するには Integral ではなく integrate 関数を使う (あとはaがわかれば面積がわかる)
result = r ** 2 * integrate((1 - cos(θ))**2, (θ, 0, a))
result
Out[10]:
$$r^{2} \left(\frac{a}{2} \sin^{2}{\left (a \right )} + \frac{a}{2} \cos^{2}{\left (a \right )} + a + \frac{1}{2} \sin{\left (a \right )} \cos{\left (a \right )} - 2 \sin{\left (a \right )}\right)$$
In [11]:
# 仮に 求める面積が 0 -> 2π までなら以下のようになる 
result.subs({a: 2 * pi})
Out[11]:
$$3 \pi r^{2}$$
In [ ]:
# 続いて体積 (X軸を軸に図形を回転させたときの体積)
In [12]:
# (X軸の)回転体の体積は以下のように求められる
pi * Integral(y ** 2, x)
Out[12]:
$$\pi \int y^{2}\, dx$$
In [13]:
# 先程の情報 (dx dθ の関係, y の媒介変数表示) を適用すると以下のように変形できる
r ** 3 * pi * Integral((1 - cos(θ)) ** 3, (θ, 0, 2 * pi))
Out[13]:
$$\pi r^{3} \int_{0}^{2 \pi} \left(- \cos{\left (θ \right )} + 1\right)^{3}\, dθ$$
In [14]:
V = r ** 3 * pi * integrate((1 - cos(θ)) ** 3, (θ, 0, 2 * pi))
V
Out[14]:
$$5 \pi^{2} r^{3}$$

ニュートン法 newton_method newton_method

最後に先程の塗りつぶした部分(一番最初の図形)の面積の近似値を求めます。

Note

近似解を求めるためにニュートン法というものを使います。

簡単に言うと 接線の X 切片が解に近づいていく性質を利用するということだそうです。

汎用的なコードが見つからなかったので自作しました。

pandas も使いたいので DataFrame を使って値の遷移を表してみます。

In [1]:
%autosave 0
%matplotlib inline
Autosave disabled
In [2]:
from itertools import count
from functools import wraps

from IPython.display import display
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from sympy import sin, cos, Symbol, symbols, integrate, Integral, diff, Eq, pi, init_printing
r, x, y, θ, a = symbols('r x y θ a')
init_printing()
In [3]:
# 最終的に求めたい値 は以下を満たす θ
expr = Eq(θ - sin(θ) - 3, 0)
expr
Out[3]:
$$θ - \sin{\left (θ \right )} - 3 = 0$$
In [10]:
X = Symbol('x')


def get_variable(expr: Symbol):
    """式で使われている変数の抽出(1変数関数のみを想定)"""
    return list(expr.free_symbols)[0]


class NewtonMethod:
    def __init__(
        self,
        plotting: bool=False,  # True にするとプロットする
        framing: bool=False,  # True にすると DataFrame を表示する
        xlim=[-10, 10], 
        ylim=[-10, 10],
        interval=0.1
    ):
        self.plotting = plotting
        self.framing = framing
        
        if plotting:
            self.interval = interval
            self.xs = np.arange(xlim[0], xlim[1], interval)
            plt.xlim(xlim)
            plt.ylim(ylim)
            plt.gca().set_aspect('equal', adjustable='box')

    def plot(self, expr: Symbol, label='', plotting=True):
        """グラフの描画"""
        if not (self.plotting and plotting):
            return
        v = get_variable(expr)
        ys = [expr.subs({v: x}) for x in self.xs]
        plt.plot(self.xs, ys, label=label)
        
    def frame(self, x, x_intercept, y, y_intercept, framing=True):
        """レコード追加"""
        if self.framing and framing:
            self.rows.append([x, x_intercept, y, y_intercept])
        
    def pre_process(self):
        """前処理"""
        self.plot(expr, label='Target')
        self.rows = []
        
    def post_process(self):
        """後処理"""
        if self.plotting:
            plt.legend()
            plt.axhline(0, color='lightgray')
            plt.axvline(0, color='lightgray')
        
        if self.framing:
            display(pd.DataFrame(self.rows, columns=['X', 'X切片', 'Y', 'Y切片']))
    
    def miscellaneous_process(f):
        @wraps(f)
        def _deco(self, *args, **kwargs):
            self.pre_process()
            result = f(self, *args, **kwargs)
            self.post_process()
            return result
        return _deco
        
    @miscellaneous_process
    def solve(
        self,
        expr: Symbol,
        tolerance: float=0.001,
        initial_value: float=100.0,
        # 収まりきらないときに表示条件(index)を指定する
        index_condition=lambda index: True
    ):
        """解く"""
        v = get_variable(expr)
        x = initial_value
        
        for i in count(1):
            gradient = diff(expr).subs({v: x})  # 傾き
            assert gradient, '傾きが 0 になってしまうと収束しないので initial_value をずらしてね'
            y = expr.subs({v: x})  # Y 座標
            y_intercept = y - (gradient * x)  # Y 切片
            x_intercept = - y_intercept / gradient  # X 切片
            
            self.plot(gradient * X + y_intercept, label=f'Line {i}', plotting=index_condition(i))
            self.frame(x, x_intercept, y, y_intercept, framing=index_condition(i))
            
            if abs(x - x_intercept) <= tolerance:
                return x_intercept
            x = x_intercept
In [11]:
# 試しに2次関数で 
expr = x ** 2  - 1
expr
Out[11]:
$$x^{2} - 1$$
In [12]:
# だんだんと答え(1)に近づいていくのがわかる
NewtonMethod(plotting=True, framing=True).solve(expr)
X X切片 Y Y切片
0 100 50.0050000000000 9999.00000000000 -10001.0000000000
1 50.0050000000000 25.0124990001000 2499.50002500000 -2501.50002500000
2 25.0124990001000 12.5262395058466 624.625106230003 -626.625106230003
3 12.5262395058466 6.30303596239437 155.906676157833 -157.906676157833
4 6.30303596239437 3.23084483304812 38.7282623432367 -40.7282623432367
5 3.23084483304812 1.77018069983297 9.43835833523375 -11.4383583352337
6 1.77018069983297 1.16754738949848 2.13353971006116 -4.13353971006116
7 1.16754738949848 1.01202183653539 0.363166906724712 -2.36316690672471
8 1.01202183653539 1.00007140387117 0.0241881976244729 -2.02418819762447
9 1.00007140387117 1.00000000254907 0.000142812840862039 -2.00014281284086
Out[12]:
$$1.00000000254907$$
In [13]:
# 逆向きから (今度は 別解(-1) に近づく)
NewtonMethod(plotting=True, framing=True).solve(expr, initial_value=-100.0)
X X切片 Y Y切片
0 -100 -50.0050000000000 9999.00000000000 -10001.0000000000
1 -50.0050000000000 -25.0124990001000 2499.50002500000 -2501.50002500000
2 -25.0124990001000 -12.5262395058466 624.625106230003 -626.625106230003
3 -12.5262395058466 -6.30303596239437 155.906676157833 -157.906676157833
4 -6.30303596239437 -3.23084483304812 38.7282623432367 -40.7282623432367
5 -3.23084483304812 -1.77018069983297 9.43835833523375 -11.4383583352337
6 -1.77018069983297 -1.16754738949848 2.13353971006116 -4.13353971006116
7 -1.16754738949848 -1.01202183653539 0.363166906724712 -2.36316690672471
8 -1.01202183653539 -1.00007140387117 0.0241881976244729 -2.02418819762447
9 -1.00007140387117 -1.00000000254907 0.000142812840862039 -2.00014281284086
Out[13]:
$$-1.00000000254907$$
In [14]:
# よさそうなので 本命を解いてみよう
# (サイクロイド曲線の中央に近いので π に近い値になるはず)
expr = θ - sin(θ) - 3
a_value = NewtonMethod(plotting=True, framing=True).solve(expr, initial_value=0.1, tolerance=0.001)
a_value
X X切片 Y Y切片
0 0.1 600.566905650820 -2.99983341664683 -3.00033300011903
1 600.566905650820 280.140581309612 598.066134724575 -522.874002277437
2 280.140581309612 130.728391351840 277.653717842245 -242.933484186483
3 130.728391351840 -65.7259080721581 128.666974217383 43.0469261509143
4 -65.7259080721581 -30.9555393418859 -68.4808906588905 60.9675129246774
5 -30.9555393418859 299.432585449740 -34.3998343617261 -31.1767602073194
6 299.432585449740 108.383643910971 297.263796995853 -168.640209475768
7 108.383643910971 3.86384774469680 104.383644759404 -3.85881455169729
8 3.86384774469680 2.99261912556238 1.52492612901676 -5.23803166991841
9 2.99261912556238 3.07095494251531 -0.155803982563299 -6.10789583778577
10 3.07095494251531 3.07076672776844 0.000375960122687966 -6.13387555931756
Out[14]:
$$3.07076672776844$$
In [16]:
# 期待通りの出力です
# これが 塗りつぶした面積の近似解ですね :)
S = r ** 2 * integrate((1 - cos(θ)) **2, (θ, 0, a))
Eq(S, S.subs({a: a_value, r: 1}))
Out[16]:
$$r^{2} \left(\frac{a}{2} \sin^{2}{\left (a \right )} + \frac{a}{2} \cos^{2}{\left (a \right )} + a + \frac{1}{2} \sin{\left (a \right )} \cos{\left (a \right )} - 2 \sin{\left (a \right )}\right) = 4.42932198525805$$

なんとか求めることができました。

Jupyter と Miyadaiku 、 相性がいいのでおすすめですよ。

興味ある方は 勉強会 へ!

参考