Pythonからはじめる数学入門 - Chapter 7: Solving Calculus Problems を読む

2021-12-29PythoncasSymPy

いよいよ最終章！Python からはじめる数学入門 (Doing Math with Python)（まとめ記事）の第 7 章を読む。第 4 章「Algebra and Symbolic Math with SymPy（ まとめ記事）」に続き、今度は 微分積分学 の初歩的な計算を SymPy で行う。SymPy の強みは数値解だけでなく、厳密解 も計算できる点にある。また応用として、関数の極値を求める 勾配降下法 を実装したりする。

内容

  • （数学的な）関数
    • 標準ライブラリと SymPy の関数を振り返る
  • Assumptions in SymPy: 不完全な型推論システムに、ヒントを与えて推論を手助けするような感じ
  • sympy.Limit: 極限
    • 関数の極限
      • 例: 関数の微分
    • 数列の極限
      • 例: 自然対数の底
  • sympy.Derivative: 微分
    • 応用: 2 階微分による極値の計算（厳密解）
    • 応用: Gradient Ascent（数値解）
  • sympy.Integral: 積分
    • 高校数学と同様の展開: 不定積分 → 定積分
    • Indefinite Integral（不定積分）
      • → 微分の逆演算
    • Definite Integral（定積分）
      • → x 軸と曲線に囲まれた領域の面積
  • Probability Density Functions（確率密度関数）
    • 正規分布 N(10,0)N(10, 0)（平均値 10）の確率計算

→ Jupyter ノートブック

ポイント

SymPy のクラス vs 関数

極限計算、微分、積分について、本では一貫して クラス を利用しているが、SymPy の 公式チュートリアル では 関数 を前面に押し出している。対応関係を表でまとめる。

Class Function
Limit limit
Derivative diff
Integral integrate

例えば微分する場合、sympy.diff() 関数では微分が即座に行われる。

x = sympy.Symbol("x")
f = sympy.sin(x)
sympy.diff(f, x)
cos(x)cos(x)

一方クラスを利用する場合、sympy.Derivative クラスを利用する。

x = sympy.Symbol("x")
f = sympy.sin(x)
d = Derivative(f, x)
d
ddxsin⁡(x)\frac{d}{d x} \sin{\left(x \right)}

結果は遅延評価されるので、実際に微分するには doit() メソッドを呼ぶ。

d.doit()
cos(x)cos(x)

また、クラス版は式の LaTeX\LaTeX コードを取得したい時に役立つ。

print(sympy.latex(d))

→ \frac{d}{d x} \sin{\left(x \right)}

2 階微分による極値の計算（厳密解）

2 階微分まで利用して極値を求める。そのポイントは以下の通り。

  • 極値 → 1 階微分が 00
  • 二階導関数が負 → 極大値
  • 二階導関数が 0 → 変曲点？
  • 二階導関数が正 → 極小値

5 次関数 f(x)=x5−30x3+50xf(x) = x^5 - 30 x^3 + 50x を例に計算してみる。まずは一階微分が 00 になるクリティカルポイント（極値の候補）を計算してみる。

from sympy import Derivative, solve, Symbol

x = Symbol("x")
f = x ** 5 - 30 * x ** 3 + 50 * x

d1 = diff(f, x)
critical_points = solve(d1)
critical_points.sort()
critical_points

→ [-sqrt(sqrt(71) + 9), -sqrt(9 - sqrt(71)), sqrt(9 - sqrt(71)), sqrt(sqrt(71) + 9)]

ちゃんと代数的に解けている！！

それぞれの解を昇順に A〜D と呼ぶことにして、二階微分係数を計算する。

A, B, C, D = critical_points

d2 = diff(f, x, 2)
print("d2(A): {0}".format(d2.subs({x: A}).evalf()))
print("d2(B): {0}".format(d2.subs({x: B}).evalf()))
print("d2(C): {0}".format(d2.subs({x: C}).evalf()))
print("d2(D): {0}".format(d2.subs({x: D}).evalf()))
  • d2(A): -703.493179468151
  • d2(B): 127.661060789073
  • d2(C): -127.661060789073
  • d2(D): 703.493179468151

したがって、A と C が極大値、B と D が極小値となる。お馴染みの増減表にまとめる。

x A B C D
f′f' + 0 - 0 + 0 - 0 +
f′′f'' - + - +
f(x)f(x) Max. Min. Max. Min.

Gradient Ascent（数値解）

関数の極小値を数値計算する有名な Gradient Descent（勾配降下法） の逆、Gradient Ascent で関数の極大値を求める。その心は単純で、次のように整理できる。

  • 微分係数が正 → 増加中 → 極大値は先方にある
  • 微分係数が負 → 減少中 → 極大値は後方にある
  • 微分係数が 0 に近い → 極大値に近い

従って初期値 x0x_0 から始めて、

xnew=xold+λdfdxx_{new} = x_{old} + \lambda \frac{df}{dx}

と、設定したステップ値 λ\lambda と微分係数に基づいて、小刻みに xx を繰り返し更新して行く。微分係数が 00 に近づき、殆ど動きがなくなったら、それを極値とする。

一階微分が数値的に 計算できれば良いので、適用範囲が広い事が、この手法の利点と言える。しかし以下の点は、気に留めておく必要がある。

  • 初期値の選択次第で到達する極値が変わる
  • Step Size: 解を探索する刻み幅
    • 大きすぎる → 解を飛び越えて振動してしまう
    • 小さすぎる → 計算に時間がかかる
  • Epsilon: 数値計算の誤差の許容範囲

SymPy で実装すると、次のようになる。

from sympy import Derivative, Symbol


def grad_ascent(x0, f1x, x):
    epsilon = 1e-6
    step_size = 1e-4

    x_old = x0
    x_new = x_old + step_size * f1x.subs({x: x_old}).evalf()

    while abs(x_old - x_new) > epsilon:
        x_old = x_new
        x_new = x_old + step_size * f1x.subs({x: x_old}).evalf()

    return x_new


def find_max(x0, f, x):
    d = Derivative(f, x).doit()

    x_max = grad_ascent(x0, d, x)

    print("x_max: {0:.5f}".format(x_max))
    print("Maximum: {0:.5f}".format(f.subs({x: x_max})))

プログラミング・チャレンジ（章末問題）

#1: Verify the Continuity of a Function at a Point

関数の値 f(a)f(a) を左極限 lim⁡x→a−f(x)\lim_{x \to a^-} f(x) と右極限 lim⁡x→a+f(x)\lim_{x \to a^+} f(x) と比較して、関数の連続性を検証する。

from sympy import limit


def is_continuous(f, x, at):
    lim_left = limit(f, x, at, "-")
    fa = f.subs({x: at})
    lim_right = limit(f, x, at, "+")

    if lim_left == fa and fa == lim_right:
        return True
    else:
        return False

#2: Implement the Gradient Descent

関数の極小値を数値計算する、お馴染みの 勾配降下法 を実装する。

from sympy import diff


def grad_descent(x0, f1x, x):
    epsilon = 1e-6
    step_size = 1e-4

    x_old = x0
    x_new = x_old - step_size * f1x.subs({x: x_old}).evalf()

    while abs(x_old - x_new) > epsilon:
        x_old = x_new
        x_new = x_old - step_size * f1x.subs({x: x_old}).evalf()

    return x_new


def find_min(x0, f, x):
    d = diff(f, x)

    x_min = grad_descent(x0, d, x)

    print("x_min: {0:.5f}".format(x_min))
    print("Minimum: {0:.5f}".format(f.subs({x: x_min})))

#3: Area Between Two Curves

2 つの関数 f(x)f(x), g(x)g(x) で定義される 2 つの曲線に囲まれた領域の面積は、積分

∫ab∣f(x)−g(x)∣dx\int_{a}^{b} | f(x) - g(x) | dx

で計算できる。SymPy では次のように実装できる。

from sympy import Abs, integrate


def area_2_curves(f, g, x, interval):
    return integrate(Abs(f - g), (x, interval[0], interval[1]))

1 次関数と 2 次関数に囲まれ領域の面積は 1/61/6 となるはずだが、積分に絶対値が入ると上手く計算できない場合があるようだ。

from sympy import Symbol

x = Symbol("x")
f = x
g = x ** 2

area_2_curves(f, g, x, (0, 1))

このとき結果は下のようなってしまう。

∫01{x2−xfor x2−x≥0−x2+xotherwise dx\int_{0}^{1} \begin{cases} x^{2} - x & \text{for}\: x^{2} - x \geq 0 \\- x^{2} + x & \text{otherwise} \end{cases}\, dx

絶対値を外すと計算は上手くいく。

#4: Finding the Length of a Curve

曲線の長さの定義通りに、次のように実装できる。

from sympy import diff, integrate, sqrt


def length_curve(f, x, interval):
    return integrate(sqrt(1 + diff(f, x) ** 2), (x, interval[0], interval[1]))

無事に年内に読破できた！！全体のまとめ記事「Python からはじめる数学入門 (Doing Math with Python) のまとめ」を公開しました！