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

微分で極値を求める

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

from sympy import *

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

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

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

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

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

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

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

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

• 極値を求める

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)
        ))

積分で面積を求める

• 不定積分を求める

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

• 定積分を求める

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

計算例

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()

参考記事