Python/SciPyで積分（1重/2重）した結果を可視化する方法

Pythonのscipy.integrate.quadやscipy.integrate.dblquadを使えば簡単に数値積分ができますが、当然結果はシンプルに積分値が返ってくるだけです。ここでは他者へ説明する目的で基本的な1重積分をはじめ、2重積分までのmatplotlibによる可視化方法を紹介します。

こんにちは。wat（＠watlablog）です。今回は積分の可視化をテーマにscipy.integrateを使ってみます！

目次（項目クリックでジャンプできます）

本記事の目標

積分結果を可視化する

曲線を積分すると面積を、曲面を2重積分すると体積を求めることができます。
一般的に積分がどのような処理を行っているかを他者へ説明する時、下図に示すような、関数と求積した部分を共にプロットするのではないでしょうか。

動作環境

この記事のコードは以下のPython環境で動作を確認しました。

Python	Python 3.9.6
PyCharm (IDE)	PyCharm CE 2020.1
Numpy	1.23.0
Scipy	1.4.1
matplotlib	3.7.1

1重積分の可視化を行うPythonコード

まずは最も簡単な式（1）を使って積分と可視化を行ってみましょう。

\[ f(x)=x \tag{1} \]

scipyのintegrateからquadをimportすることで積分を行います。関数fを定義しておき、積分区間a, bを使ってquad(f, a, b)と書くことで求積が可能です。

積分結果の可視化は.fill_betweenで行っています。1重積分の可視化はこれ一択だと思いました。

scipy.integrate.quadの公式ドキュメント：https://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.quad.html

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import quad

def f(x):
    """ 被積分関数 """

    y = x

    return y

def plot(x, y, a, b, result):
    """ グラフ描画する関数 """

    # フォントの種類とサイズを設定する。
    plt.rcParams['font.size'] = 14
    plt.rcParams['font.family'] = 'Times New Roman'

    # 目盛を内側にする。
    plt.rcParams['xtick.direction'] = 'in'
    plt.rcParams['ytick.direction'] = 'in'

    # グラフの上下左右に目盛線を付ける。
    fig = plt.figure()
    ax1 = fig.add_subplot(111)
    ax1.yaxis.set_ticks_position('both')
    ax1.xaxis.set_ticks_position('both')

    # 軸のラベルを設定する。
    ax1.set_xlabel('x')
    ax1.set_ylabel('y')

    # データをプロットする。
    label_text = "$\int_{{a}}^{{b}} f(x) \, dx = $" + str(result)
    ax1.plot(x, y, label=label_text)
    ax1.fill_between(x, y, where=[a <= xi <= b for xi in x], color='red', alpha=0.5)
    ax1.legend(loc='upper left')

    # グラフを表示する。
    fig.tight_layout()
    plt.show()
    plt.close()

    return


if __name__ == '__main__':
    """ Main """

    # 積分：x ∈ [a, b]
    a, b = 1, 2
    result, _ = quad(f, a, b)
    print('Quadrature result=', result)

    # グラフの描画
    x = np.linspace(0, 4, 100)
    y = f(x)
    plot(x, y, a, b, result)

import numpy as np

import matplotlib.pyplot as plt

from scipy.integrate import quad

def f(x):

""" 被積分関数 """

y = x

return y

def plot(x, y, a, b, result):

""" グラフ描画する関数 """

# フォントの種類とサイズを設定する。

plt.rcParams['font.size'] = 14

plt.rcParams['font.family'] = 'Times New Roman'

# 目盛を内側にする。

plt.rcParams['xtick.direction'] = 'in'

plt.rcParams['ytick.direction'] = 'in'

# グラフの上下左右に目盛線を付ける。

fig = plt.figure()

ax1 = fig.add_subplot(111)

ax1.yaxis.set_ticks_position('both')

ax1.xaxis.set_ticks_position('both')

# 軸のラベルを設定する。

ax1.set_xlabel('x')

ax1.set_ylabel('y')

# データをプロットする。

label_text = "$\int_{{a}}^{{b}} f(x) \, dx = $" + str(result)

ax1.plot(x, y, label=label_text)

ax1.fill_between(x, y, where=[a <= xi <= b for xi in x], color='red', alpha=0.5)

ax1.legend(loc='upper left')

# グラフを表示する。

fig.tight_layout()

plt.show()

plt.close()

return

if __name__ == '__main__':

""" Main """

# 積分：x ∈ [a, b]

a, b = 1, 2

result, _ = quad(f, a, b)

print('Quadrature result=', result)

# グラフの描画

x = np.linspace(0, 4, 100)

y = f(x)

plot(x, y, a, b, result)

結果はこちら。出力ウィンドウにもquadの結果が書かれていますが、labelにも載せるようにしてみました。手計算してみるとわかりますが、quadの数値1.5は正解です（この積分は暗算でもできますね）。

次はもう少し変曲点が多くて面白い関数として、式（2）で試してみましょう。

\[ f(x)= x^{4} - 10 x^{3} + 25 x^{2} \tag{2} \]

関数fの部分と、積分区間a, bを変更します。

def f(x):
    """ 被積分関数 """

    y = x ** 4 - 10 * x ** 3 + 25 * x ** 2

    return y

# 積分：x ∈ [a, b]
a, b = 2, 3

def f(x):

""" 被積分関数 """

y = x ** 4 - 10 * x ** 3 + 25 * x ** 2

return y

# 積分：x ∈ [a, b]

a, b = 2, 3

こちらが結果です。式を変更するだけで正常に機能していることがわかりました。

これで1重積分の計算方法や可視化方法はわかったので、次に進みましょう。

2重積分の可視化を行うPythonコード

次は式（3）の2重積分を行います。式（3）は平面の方程式です。ここで$w_{0}=1, w_{1}=1, w_{2}=1$とします。

\[ f(x, y) = w_{0} + w_{1} x + w_{2} y \tag{3} \]

コードはこちら。2重積分はscipy.integrate.dblquadを使います（double quadratureという意味らしい）。2重積分の可視化は色々試行錯誤しましたが、積分区間を黒いエッジで囲い、面にカラーマップを貼るという見せ方が最も美しいとの結論に達しました（個人の感想です）。

scipy.integrate.dblquadの公式ドキュメント：https://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.dblquad.html

積分区間をエッジで囲うために.plot_wireframeを使っていますが、これだけだと細かい線が多数表示されるためrstrideとcstrideを使ってメッシュの間隔を大きくする工夫をしました。

また本筋とは関係ないですが、*x_range_integrateとアスタリスクをつけて引数を渡すことで一つずつ記載しなくて良いという小ワザも使ってみました。

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import dblquad

def f(x, y):
    """ 被積分関数 """

    w0 = 1
    w1 = 1
    w2 = 1

    z = w0 + w1 * x + w2 * y

    return z


def plot(x_vals, y_vals, z_vals, x_vals_integrate, y_vals_integrate, z_vals_integrate, result):
    """ グラフ描画する関数 """

    # フォントの種類とサイズを設定する。
    plt.rcParams['font.size'] = 14
    plt.rcParams['font.family'] = 'Times New Roman'

    # 目盛を内側にする。
    plt.rcParams['xtick.direction'] = 'in'
    plt.rcParams['ytick.direction'] = 'in'

    # グラフの上下左右に目盛線を付ける。
    fig = plt.figure()
    ax1 = fig.add_subplot(111, projection='3d')
    ax1.yaxis.set_ticks_position('both')
    ax1.xaxis.set_ticks_position('both')

    # 軸のラベルを設定する。
    ax1.set_xlabel('x')
    ax1.set_ylabel('y')
    ax1.set_zlabel('z')

    # データをプロットする。
    # 被積分関数
    label_text = "$\int_{{a}}^{{b}}\int_{{c}}^{{d}} f(x, y) \, dydx = $" + str(result)
    ax1.plot_wireframe(x_vals, y_vals, z_vals, alpha=0.2, edgecolor='gray')

    # 積分区間の面
    ax1.plot_wireframe(x_vals_integrate, y_vals_integrate, z_vals_integrate, rstride=99, cstride=99, color='black', label=label_text)
    ax1.plot_surface(x_vals_integrate, y_vals_integrate, z_vals_integrate, cmap='jet', alpha=0.5)

    # グラフを表示する。
    ax1.legend()
    fig.tight_layout()
    plt.show()
    plt.close()

    return


if __name__ == '__main__':
    """ Main """

    # 描画の範囲を定義
    x_range = [-10, 10]
    y_range = [-10, 10]

    # 積分の範囲を定義
    a, b = 0, 5
    c, d = -5, 5
    x_range_integrate = [a, b]
    y_range_integrate = [c, d]

    # 積分：x ∈ [a, b], y ∈ [c, d]
    result, _ = dblquad(f, a, b, c, d)
    print('Quadrature result=', result)

    # グリッドを生成
    x_vals, y_vals = np.meshgrid(np.linspace(*x_range, 100), np.linspace(*y_range, 100))
    x_vals_integrate, y_vals_integrate = np.meshgrid(np.linspace(*x_range_integrate, 100), np.linspace(*y_range_integrate, 100))

    # zの値を計算
    z_vals = f(x_vals, y_vals)
    z_vals_integrate = f(x_vals_integrate, y_vals_integrate)

    # 可視化
    plot(x_vals, y_vals, z_vals, x_vals_integrate, y_vals_integrate, z_vals_integrate, result)

import numpy as np

import matplotlib.pyplot as plt

from scipy.integrate import dblquad

def f(x, y):

""" 被積分関数 """

w0 = 1

w1 = 1

w2 = 1

z = w0 + w1 * x + w2 * y

return z

def plot(x_vals, y_vals, z_vals, x_vals_integrate, y_vals_integrate, z_vals_integrate, result):

""" グラフ描画する関数 """

# フォントの種類とサイズを設定する。

plt.rcParams['font.size'] = 14

plt.rcParams['font.family'] = 'Times New Roman'

# 目盛を内側にする。

plt.rcParams['xtick.direction'] = 'in'

plt.rcParams['ytick.direction'] = 'in'

# グラフの上下左右に目盛線を付ける。

fig = plt.figure()

ax1 = fig.add_subplot(111, projection='3d')

ax1.yaxis.set_ticks_position('both')

ax1.xaxis.set_ticks_position('both')

# 軸のラベルを設定する。

ax1.set_xlabel('x')

ax1.set_ylabel('y')

ax1.set_zlabel('z')

# データをプロットする。

# 被積分関数

label_text = "$\int_{{a}}^{{b}}\int_{{c}}^{{d}} f(x, y) \, dydx = $" + str(result)

ax1.plot_wireframe(x_vals, y_vals, z_vals, alpha=0.2, edgecolor='gray')

# 積分区間の面

ax1.plot_wireframe(x_vals_integrate, y_vals_integrate, z_vals_integrate, rstride=99, cstride=99, color='black', label=label_text)

ax1.plot_surface(x_vals_integrate, y_vals_integrate, z_vals_integrate, cmap='jet', alpha=0.5)

# グラフを表示する。

ax1.legend()

fig.tight_layout()

plt.show()

plt.close()

return

if __name__ == '__main__':

""" Main """

# 描画の範囲を定義

x_range = [-10, 10]

y_range = [-10, 10]

# 積分の範囲を定義

a, b = 0, 5

c, d = -5, 5

x_range_integrate = [a, b]

y_range_integrate = [c, d]

# 積分：x ∈ [a, b], y ∈ [c, d]

result, _ = dblquad(f, a, b, c, d)

print('Quadrature result=', result)

# グリッドを生成

x_vals, y_vals = np.meshgrid(np.linspace(*x_range, 100), np.linspace(*y_range, 100))

x_vals_integrate, y_vals_integrate = np.meshgrid(np.linspace(*x_range_integrate, 100), np.linspace(*y_range_integrate, 100))

# zの値を計算

z_vals = f(x_vals, y_vals)

z_vals_integrate = f(x_vals_integrate, y_vals_integrate)

# 可視化

plot(x_vals, y_vals, z_vals, x_vals_integrate, y_vals_integrate, z_vals_integrate, result)

結果がこちら。3次元プロットはやはり映えますね。積分結果は175（数値積分なのでピッタリいかない）となりました。果たしてこの結果は合っているのでしょうか？2重積分の暗算は少々やっかいなので計算してみましょう。

積分は加算に関して線形であるため、3つの積分に分解できます。3つ目の項は奇関数であるため0になり簡単になるため、式（4）のように展開して答えは175を得ました。プログラムは正しく動作しているようですね。

\[ \begin{eqnarray} \int_{0}^{5}\int_{-5}^{5} (1 + x + y) dydx &=&　\int_{0}^{5} \int_{-5}^{5} 1 \, dy \, dx + \int_{0}^{5} \int_{-5}^{5} x \, dy \, dx + \int_{0}^{5} \int_{-5}^{5} y \, dy \, dx\\ &=&5 \times 10 + \frac{5^2}{2} \times 10 + 0\\ &=&175 \end{eqnarray} \tag{4} \]

自分でプログラムを書いた時は、手計算でも簡単に答えが出る問題を解いて検証することが重要です！今回のケースだとxとyの範囲を逆にしてしまうバグもあるあるです。

次は式（5）で示す2次元のガウス関数です。2次元ガウス関数は様々な用途で用いられるので、これも2重積分してみましょう。

\[ f(x, y) = e^{-\frac{(x - \mu_x)^2}{2 \sigma_x^2}} e^{-\frac{(y - \mu_y)^2}{2 \sigma_y^2}} \tag{5} \]

Pythonコードで変更するのは関数部分のみです。積分範囲は先ほどと同様にしました。

def f(x, y):
    """ 被積分関数 """

    mu_x = 0
    mu_y = 2
    sigma_x = 3
    sigma_y = 1

    z = np.exp(-(x - mu_x) ** 2 / (2 * sigma_x ** 2)) * np.exp(-(y - mu_y) ** 2 / (2 * sigma_y ** 2))

    return z

def f(x, y):

""" 被積分関数 """

mu_x = 0

mu_y = 2

sigma_x = 3

sigma_y = 1

z = np.exp(-(x - mu_x) ** 2 / (2 * sigma_x ** 2)) * np.exp(-(y - mu_y) ** 2 / (2 * sigma_y ** 2))

return z

結果がこちら。綺麗ですね。ガウス関数は天文学であれば星の光強度、工学であればベクトルデータのヒストグラムといった用途で2重積分の需要がありそうです。

まとめ

この記事ではscipy.integrate.quadやdblquadを使った積分方法を紹介しました。それだけだと類似の記事がたくさんあったので、このブログでは主に可視化に特化してみました。

積分結果を得るだけなら可視化は全く必要ありませんが、結果を他者へ説明したり理解を深めたりする時には有用だと思います。

今回は積分区間が定数の問題のみを扱いましたが、世の中には積分区間にxやyが含まれる問題もあります。次回はそのような一歩進んだ問題の可視化もトライしてみたいと思いました。

基本的な積分を求めるPythonコードが書けました！
X（旧Twitter）でも関連情報をつぶやいているので、wat(@watlablog)のフォローお待ちしています！

Python/SciPyで積分（1重/2重）した結果を可視化する方法

本記事の目標

積分結果を可視化する

関連記事

動作環境

1重積分の可視化を行うPythonコード

2重積分の可視化を行うPythonコード

まとめ

コメントを残すコメントをキャンセル

最近の投稿

アーカイブ

カテゴリー

本記事の目標

積分結果を可視化する

関連記事

動作環境

1重積分の可視化を行うPythonコード

2重積分の可視化を行うPythonコード

まとめ

SNSでもご購読できます。

コメントを残す コメントをキャンセル

最近の投稿

アーカイブ

カテゴリー

コメントを残すコメントをキャンセル