Python sympy 微分・積分で曲線の極値や面積を求める

# 微分で極値を求める

- 関数 $f(x)$ を定義

```py
from sympy import *

x = Symbol('x')  # 変数を設定
f = x ** 3 - 2 * x
f
```

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/621414/2b2c8acf-ea42-7edd-acef-4ddba88fbc3f.png)

- $f'(x), f''(x)$ を求める

```py
diff(f, x)  # 第二引数は変数を渡す
```

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/621414/d1572785-99de-c410-8aac-8b08f4cc4782.png)

```py
diff(f, x, 2)  # 第三引数は微分の階数を渡す
```

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/621414/70860f1a-efab-560c-7187-99e4dc0c45a3.png)


- $f'(x) = 0$ の解 $x_1, x_2$ を求める

```py
critical_points = solve(diff(f, x))  # solveメソッドには関数を渡す
x_1, x_2 = critical_points
```

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/621414/0a296dbf-a6c2-0f45-62a6-c11408d16419.png)

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/621414/ece15c3c-4113-410f-707a-b04382f5d371.png)

- $f''(x_1), f''(x_2)$ を求める

```py
f2 = diff(f, x, 2)
f2.subs(x, x_1)  # 関数に値を代入するにはsubsメソッドを使う
f2.subs(x, x_2)
```

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/621414/b8e9e975-eda0-99c6-6b9d-577885cf8c46.png)

- 極値を求める

```py
for p in critical_points:
    if f2.subs(x, p) < 0:
        print('f(x)はx={}で極大値{}をとる'.format(
            p, f2.subs(x, p)
        ))
    elif f2.subs(x, p) > 0:
        print('f(x)はx={}で極小値{}をとる'.format(
            p, f2.subs(x, p)
        ))
```

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/621414/9ef274ff-f469-56da-26b1-79c514a38a95.png)

# 積分で面積を求める

- 不定積分を求める

```py
integrate(f, x)  # 第二引数は変数を渡す
```

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/621414/f387608f-6bdf-9670-41a1-bb2c2e8b3257.png)

- 定積分を求める

```py
# 第二引数はタプルを渡す。タプルは(変数, 左端の区間, 右端の区間)
integrate(f, (x, -1, 0))
```

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/621414/5bd75adb-1297-21b4-354c-7a4e2cb6fe44.png)

# 計算例

```py
import sympy as *
import numpy as np
import matplotlib.pyplot as plt

EPSILON = 1e-10
plt.figure(figsize=(12.0, 10.0))

def cal_integral(curve1, curve2):
    """
    積分の区間を定義するメソッド
    """
    x = Symbol('x')
    x_interval = np.linspace(-2, 2, 100)
    # 2曲線で囲まれた領域の面積を求める
    new_function = curve1 - curve2
    # 2曲線の交点のx座標を求める
    x_intersections = sorted(solve(new_function))
    print('交点のx座標: {}\n'.format(x_intersections))
    # 2曲線の交点のy座標を求める
    y_intersections = np.empty(0)
    for x_point in x_intersections:
        y_intersections = np.append(y_intersections, curve1.subs(x, x_point))
    # 区間に応じて積分を計算する
    for n in range(len(x_intersections)):
        if n == 0:
            _cal_integral(new_function, x_interval[n], x_intersections[n])
        if n == len(x_intersections) - 1:
            _cal_integral(new_function, x_intersections[n], x_interval[-1])
        else:
            _cal_integral(
                new_function, x_intersections[n], x_intersections[n + 1])
    # 交点をプロットする
    plt.scatter(x_intersections, y_intersections, marker='o', color='black', s=200)
            
            
def _cal_integral(curve, left_edge, right_edge):
    """
    積分を計算するメソッド
    """
    x = Symbol('x')
    area = integrate(
        curve, (x, left_edge, right_edge)
    )
    print('積分区間: {} to {}\n面積: {}\n'.format(
        left_edge, right_edge, area
    ))


def cal_extremes(curve, length):
    """
    極値を計算するメソッド
    """
    x = Symbol('x')
    # f'(x) = 0の解を求める
    critical_points = solve(diff(curve, x))
    # 解の虚部が十分小さければ実数とみなす
    critical_points = [re(point) for point in critical_points
                       if abs(im(point)) < EPSILON]
    # 実数とみなせるf'(x) = 0の解について極値を求める
    for critical_point in critical_points:
        if im(critical_point) != 0\
           or length < critical_point or critical_point < -length:
            continue
        extreme_value = curve.subs(x, critical_point)
        diff1 = diff(curve, x)
        diff2 = diff(diff1, x)
        if diff2.subs(x, critical_point) < 0:
            print('x = {}のとき極大値 f(x) = {}\n'.format(
                critical_point, extreme_value))
        elif diff2.subs(x, critical_point) > 0:
            print('x = {}のとき極小値 f(x) = {}\n'.format(
                critical_point, extreme_value))
        else:
            # do samething
            pass
    plot_extremes(curve, critical_points, length)

def plot_extremes(curve, critical_points, length):
    """
    極値を描画するメソッド
    """
    x = Symbol('x')
    critical_images = np.empty(0)
    real_critical_points = np.empty(0)
    # 実数とみなせるf'(x) = 0の解xについて、極値をプロットする
    for point in critical_points:
        if im(point) != 0\
           or length < point:
            continue
        real_critical_points = np.append(real_critical_points, point)
        critical_images = np.append(critical_images, curve.subs(x, point))
    plt.scatter(real_critical_points, critical_images, s=200, color='red')


def plot_curve():
    """
    関数を描画するメソッド
    """
    x = Symbol('x')
    # 2つの曲線（または直線）を定義する
    # 関数の係数を定義する
    coeff = [0, -2, 0, 1]
    curve = 0
    # 係数をもとに関数を定義する（曲線1つ目）
    for i in range(len(coeff)):
        curve += coeff[i] * x ** i 
    # 定数関数を定義する（曲線2つ目）
    const_line = 0

    curve_images = np.empty(0)
    const_images = np.empty(0)

    # xの区間を定義する
    x_interval = np.linspace(-2, 2, 100)
    # 関数に対して、xに対応するf(x)の値を計算する
    for x_point in x_interval:
        curve_images = np.append(curve_images, curve.subs(x, x_point))
        const_images = np.append(const_images, const_line)
    # (x, f(x))をもとにグラフをプロットする（曲線1つ目）
    plt.plot(x_interval, curve_images, 'blue')
    # (x, f(x))をもとにグラフをプロットする（曲線2つ目）
    plt.plot(x_interval, const_images, 'green')
    plt.bar(x_interval, curve_images, width=0.01)
    # 積分を計算
    cal_integral(curve, const_line)
    # 極値を計算
    cal_extremes(curve, x_interval[-1] - x_interval[0])


plot_curve()
```

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/621414/5dc0a738-9c83-4557-5a32-5730081dfb92.png)

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/621414/868f4871-3f9f-2032-b330-cb5dbf7a72b5.png)


# 参考記事

https://docs.sympy.org/latest/tutorial/calculus.html

https://home.hirosaki-u.ac.jp/jupyter/python-sk/

https://thr3a.hatenablog.com/entry/20190916/1568571135