せっかく Jupyter の記事を投稿できるので試してみました。
Miyadaiku と Jupyter でこういう記事がかけるんだーくらいに思ってもらえれば御の字です。
Note
当サイトでは 独自にデザインをあてているため、 ipynb ファイルを置いてもまったく同じようにはなりません。
jupyter-themes onedork をベースにしています
ベースのテーマを継承したテンプレートさえ利用していれば、 ipynb ファイルを置くだけで通常の描画はできると思います。
当記事は ipynb ファイルだけで記事を構成しているわけではなく rst ファイルから ipynb ファイルを snippet として インクルードして描画しています。
やり方は Miyadaiku スニペット リファレンス を参照
ありきたりな例では面白くないので、今回はサイクロイド曲線を題材にしていろいろやってみます。なつかしい。
図形描画
描画ライブラリは色々ありますが、割と一般的な matplotlib を使います。
In [1]:
%autosave 0
%matplotlib inline
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()
面積と体積
SymPy も使ってみたいので、(上記で描画した図形の)面積と体積を求めます。
In [3]:
%autosave 0
In [4]:
from sympy import sin, cos, symbols, integrate, Integral, diff, Eq, pi, init_printing
init_printing()
S, V, r, x, y, θ, dx, dθ, a = symbols('S V r x y θ dx dθ a')
In [5]:
# 計算方法のおさらい
Eq(x, r(θ - sin(θ))), Eq(y, r(1 - cos(θ)))
Out[5]:
In [6]:
# まずは 面積から
Eq(S, Integral(y, x))
Out[6]:
In [7]:
# 上記の条件から dx と dθ の関係を求める
Eq(dx/dθ, diff(r * (θ - sin(θ)), θ))
Out[7]:
In [8]:
# 変形
Eq(dx, r(1 - cos(θ)) * dθ)
Out[8]:
In [9]:
# y と等しいので以下のようにまとめられる
Eq(S, r ** 2 * Integral((1 - cos(θ))**2, (θ, 0, a)))
Out[9]:
In [10]:
# 定積分するには Integral ではなく integrate 関数を使う (あとはaがわかれば面積がわかる)
result = r ** 2 * integrate((1 - cos(θ))**2, (θ, 0, a))
result
Out[10]:
In [11]:
# 仮に 求める面積が 0 -> 2π までなら以下のようになる
result.subs({a: 2 * pi})
Out[11]:
In [ ]:
# 続いて体積 (X軸を軸に図形を回転させたときの体積)
In [12]:
# (X軸の)回転体の体積は以下のように求められる
pi * Integral(y ** 2, x)
Out[12]:
In [13]:
# 先程の情報 (dx dθ の関係, y の媒介変数表示) を適用すると以下のように変形できる
r ** 3 * pi * Integral((1 - cos(θ)) ** 3, (θ, 0, 2 * pi))
Out[13]:
In [14]:
V = r ** 3 * pi * integrate((1 - cos(θ)) ** 3, (θ, 0, 2 * pi))
V
Out[14]:
ニュートン法
最後に先程の塗りつぶした部分(一番最初の図形)の面積の近似値を求めます。
Note
近似解を求めるためにニュートン法というものを使います。
簡単に言うと 接線の X 切片が解に近づいていく性質を利用するということだそうです。
汎用的なコードが見つからなかったので自作しました。
pandas も使いたいので DataFrame を使って値の遷移を表してみます。
In [1]:
%autosave 0
%matplotlib inline
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]:
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]:
In [12]:
# だんだんと答え(1)に近づいていくのがわかる
NewtonMethod(plotting=True, framing=True).solve(expr)
Out[12]:
In [13]:
# 逆向きから (今度は 別解(-1) に近づく)
NewtonMethod(plotting=True, framing=True).solve(expr, initial_value=-100.0)
Out[13]:
In [14]:
# よさそうなので 本命を解いてみよう
# (サイクロイド曲線の中央に近いので π に近い値になるはず)
expr = θ - sin(θ) - 3
a_value = NewtonMethod(plotting=True, framing=True).solve(expr, initial_value=0.1, tolerance=0.001)
a_value
Out[14]:
In [16]:
# 期待通りの出力です
# これが 塗りつぶした面積の近似解ですね :)
S = r ** 2 * integrate((1 - cos(θ)) **2, (θ, 0, a))
Eq(S, S.subs({a: a_value, r: 1}))
Out[16]:
なんとか求めることができました。
Jupyter と Miyadaiku 、 相性がいいのでおすすめですよ。
興味ある方は 勉強会 へ！