「数学の言葉で世界を見たら」をPython で見たら… Vol.11 積分のいろいろ

第７話「微積は積分から」

書籍の著者　大栗博司 先生

書籍「数学の言葉で世界を見たら」第７話「微積は積分から」の Python写経活動記録 です。 

書籍にならって いろいろな関数を積分 します。

では書籍を開いて数学の旅に出発です🚀

はじめに

このブログシリーズは、書籍「数学の言葉で世界を見たら」（幻冬舎）で学んだ「数学の楽しさ」を「Python 写経の形式」でご紹介いたします。

【引用表記】
この記事は、出典に記載の書籍に掲載された文章とデータを引用し、適宜、掲載文章・データを改変して書いています。
【出典】
「数学の言葉で世界を見たら」 第3刷、著者 大栗博司、幻冬舎

５ いろいろな関数を積分してみる

学びポイント

いろいろな関数の積分について、書籍で「解析的に求める方法」を学び、ブログでは「Python の積分」を学びます。

この記事で用いるライブラリをインポートします。

## インポート

# 数値計算
import numpy as np
import pandas as pd

# 積分
from scipy.integrate import quad  # 数値積分
import sympy as sp

# TeX表示
from IPython.display import Math

# 可視化
import matplotlib.pyplot as plt
plt.rcParams['font.family'] = 'Meiryo'  # または import japanize_matplotlib

１次関数の積分

書籍 p.225～の

１次関数 $${f(x)=x}$$ を 区間 $${[0, a]}$$ で積分したらどうなるか

を $${a=2}$$ にして Python で実践します。

1️⃣ 関数の可視化
$${f(x)=x}$$ を可視化します。

### p.225~ １次関数 y=x の場合 

## 関数の定義 f=x
f1 = lambda x: x

## 可視化
# 関数fの曲線の描画
plt.plot(x_val:=np.linspace(0, 3, 1001), f1(x_val))
# 区間[1, 2]の面の塗りつぶし
plt.fill_between(x_val2:=np.linspace(0, 2, 101), f1(x_val2), alpha=0.1)
# 区間[1, 2]の両端の垂直線の描画
plt.vlines(0, 0, f1(0), color='tab:blue', lw=1)
plt.vlines(2, 0, f1(2), color='tab:blue', lw=1)
# 積分数式の表示
plt.text(x=1.5, y=0.5, s=r'$\int_0^2 f\ (x)\ dx$', ha='center', fontsize=14)
# 修飾
plt.xlabel('$x$', fontsize=12)
plt.ylabel('$f\ (x) = x$', fontsize=12)
plt.xticks(range(0, 4))
plt.ylim(0, 3);

【実行結果】
$${f(x)=x}$$ は傾き１の直線です。

青く塗られた領域について、面積算出：積分をします。
ちなみに答えは $${2}$$ です！

2️⃣ 「アルキメデスのはさみうち」による面積算出の形式で積分
前回までの記事で大活躍した「アルキメデスのはさみうち」で青く塗られた領域の面積を求めます。
横に並べた細い長方形の面積を計算するアレです。

%%time
## アルキメデスのはさみうちで関数の面積を推定

for n in range(3999995, 4000005):

    # 長方形の横の長さεの算出
    eps = 2 / n

    # n+1個の分割点ごとのf1(x)の値の算出
    fx = f1(np.linspace(0, 2, n + 1))
    
    # Cnの面積の算出
    Cn = sum([fx[i] for i in range(1, n + 1)]) * eps
    
    # Bnの面積の算出
    Bn = sum([fx[i] for i in range(n)]) * eps

    # 差が一定値未満になったら推定を終了
    if abs(Cn - Bn) < 1e-6:
        break

## 面積の推定値の表示
print(f'n = {n}\nCn = {Cn:.7f}\nBn = {Bn:.7f}\nCn-Bn = {Cn - Bn:.7f}\n'
      f'(Cn+Bn)/2 = {(Cn + Bn) / 2:.7f}')
print('-'*25)

【実行結果】
横に 4000000 個に分割した長方形の面積は $${2}$$ です。

三角形の面積「底辺 $${2 \times}$$ 高さ $${2 \div 2=2}$$」で求めることもできます。

3️⃣ scipy で数値積分
scipy.integrate の関数 quad() で積分します。

## scipyで数値積分
result = quad(f1, a=0, b=2)
print(f'数値積分による面積: {result[0]:.5f}')

【実行結果】

4️⃣ 書籍の極限で求めた積分
書籍 p.227 の「$${n}$$ が無限大の極限では面積は $${a^2/2}$$ になる」を利用する積分の解法を利用します。

## p.227 極限と積分で求めた式を用いて面積を算出
a = 2
print(f'極限と積分で求めた式による面積: {a**2 / 2:.6f}')

【実行結果】

２次関数の積分

書籍 p.227～の

２次関数 $${f(x)=x^2}$$ の 区間 $${[0, a]}$$ の面積も同じようにして…

を $${a=2}$$ にして Python で実践します。

1️⃣ 関数の可視化
$${f(x)=x^2}$$ を可視化します。

### p.227~ ２次関数 y=x² の場合

## 関数の定義 f=x²
f2 = lambda x: x**2

## 可視化
# 関数fの曲線の描画
plt.plot(x_val:=np.linspace(0, 3, 1001), f2(x_val))
# 区間[1, 2]の面の塗りつぶし
plt.fill_between(x_val2:=np.linspace(0, 2, 101), f2(x_val2), alpha=0.1)
# 区間[1, 2]の両端の垂直線の描画
plt.vlines(0, 0, f2(0), color='tab:blue', lw=1)
plt.vlines(2, 0, f2(2), color='tab:blue', lw=1)
# 積分数式の表示
plt.text(x=1.5, y=0.5, s=r'$\int_0^2 f\ (x)\ dx$', ha='center', fontsize=14)
# 修飾
plt.xlabel('$x$', fontsize=12)
plt.ylabel('$f\ (x) = x^2$', fontsize=12)
plt.xticks(range(0, 4))
plt.ylim(0, 5);

【実行結果】
$${f(x)=x^2}$$ は放物線です。
$${0}$$ 以上の区間なので、放物線の右半分が描かれています。

青く塗られた領域について、面積算出：積分をします。
ちなみに答えは $${2.666 \cdots}$$ です！

2️⃣ 「アルキメデスのはさみうち」による面積算出の形式で積分
前回までの記事で大活躍した「アルキメデスのはさみうち」で青く塗られた領域の面積を求めます。
横に並べた細い長方形の面積を計算するアレです。

%%time
## アルキメデスのはさみうちで関数の面積を推定

for n in range(7999995, 8000005):

    # 長方形の横の長さεの算出
    eps = 2 / n

    # n+1個の分割点ごとのf1(x)の値の算出
    fx = f2(np.linspace(0, 2, n + 1))
    
    # Cnの面積の算出
    Cn = sum([fx[i] for i in range(1, n + 1)]) * eps
    
    # Bnの面積の算出
    Bn = sum([fx[i] for i in range(n)]) * eps

    # 差が一定値未満になったら推定を終了
    if abs(Cn - Bn) < 1e-6:
        break

## 面積の推定値の表示
print(f'n = {n}\nCn = {Cn:.7f}\nBn = {Bn:.7f}\nCn-Bn = {Cn - Bn:.7f}\n'
      f'(Cn+Bn)/2 = {(Cn + Bn) / 2:.7f}')
print('-'*25)

【実行結果】
横に 8000001 個に分割した長方形の面積は $${2.666 \cdots}$$ です。

3️⃣ scipy で数値積分
scipy.integrate の関数 quad() で積分します。

## scipyで数値積分
result = quad(f2, a=0, b=2)
print(f'数値積分による面積: {result[0]:.6f}')

【実行結果】

4️⃣ 書籍の極限で求めた積分
書籍 p.228 の積分の解法を利用します。

# p.227 極限と積分で求めた式を用いて面積を算出
a = 2
print(f'極限と積分で求めた式による面積: {a**3 / 3:.6f}')

【実行結果】

高次の関数の積分

書籍 p.228 で言及されている $${f(x)=x^k}$$、区間 $${[0, a]}$$ の積分です！
もちろん $${a=2}$$ にして Python で実践します。

1️⃣ 書籍の解法
書籍 p.228 の積分の解法を利用します。
$${k=1, 2, \cdots, 20}$$ でやってみます。

### p.228 高次の関数 y=xᵏ の定積分 ※書籍の数式利用

## 積分関数の定義
integrate_xk_0_to_a = lambda a, k: (a**(k + 1) / (k + 1))

## 設定
# 区間 [0, a]
a = 2
# 指数 k
ks = np.arange(1, 21)

## 積分を実行してデータフレーム化
result_text = pd.DataFrame({'積分': integrate_xk_0_to_a(a, ks)}, index=ks)
result_text.index.name = 'k'
result_text

【実行結果】
次数 $${k}$$ が大きくなるにつれて、積分の結果はさらに大きくなっているようです。

可視化しましょう。

# 可視化
result_text.plot();

【実行結果】
（数学の言葉ではないかもですが）逓増してます！

3️⃣ scipy で数値積分
scipy.integrate の関数 quad() で積分します。

### p.228 高次の関数 y=xᵏ の定積分 ※scipy.integrate 利用

## 設定
# 区間 [0, a]
a = 2
# 指数 k
ks = np.arange(1, 21)
# 関数 f(x) = xᵏ
fx1 = lambda x, k: x**k

## 数値積分の計算
result = [quad(fx1, a=0, b=a, args=(k, ))[0] for k in ks]

## データフレーム化
result_scipy = pd.DataFrame({'積分': result}, index=ks)
result_scipy.index.name = 'k'
result_scipy

【実行結果】
書籍の結果と同じです！（当然！？）

もちろん可視化します。

# 可視化
result_scipy.plot();

【実行結果】

sympyで積分

高次の関数 $${f(x)=x^k}$$ の変数 $${x, k, a}$$ をそのまま活かして積分したいときは、sympy の登場です！

1️⃣ 定積分
関数 $${f(x)=x^k}$$、区間 $${[0, a]}$$ の定積分です。
sympy の integrate で積分できます。
引数は、関数、（積分区間の対象変数・下端・上端）です。
積分の数式を integ_func1 に格納しておきましょう。

## sympy で積分

# 変数の定義
x, a = sp.symbols('x a')
k = sp.symbols('k', positive=True)  # 正値

# 関数 f(x) = xᵏ
fx2 = x**k

# 定積分の場合
integ_func1 = sp.integrate(fx2, (x, 0, a))
display(Math(f'\int_0^a x^k dx = {sp.latex(integ_func1)}'))

【実行結果】
書籍 p.228 の数式を導くことができました！
変数 $${x}$$ は消えて、$${k, a}$$ が変数の形式で残っています！

可視化です！
sympy 用の可視化ライブラリ SymPy Plotting Backends（spb）を使えば、sympy 様式の関数 integ_func1 自体を可視化できます！
$${a=2}$$ で固定して $${k}$$ の関数とし、$${0 \leq k \leq 10}$$ の範囲で描画します。

# 可視化（a=2の場合）
spb.plot(integ_func1.subs(a, 2), (k, 0, 10));

【実行結果】

こちらは、$${k=2}$$ で固定して $${a}$$ の関数にしています。

# 可視化（k=2の場合）
spb.plot(integ_func1.subs(k, 2));

【実行結果】

2️⃣ 不定積分
（書籍に掲載されていないですが）関数 $${f(x)=x^k}$$ の不定積分です。
sympy の integrate で、積分区間を設定しない感じです。
積分の数式を integ_func2 に格納しておきましょう。

# 不定積分の場合 ※積分定数の表示は省略
integ_func2 = sp.integrate(fx2, x)
display(Math(f'\int x^k dx = {sp.latex(integ_func2)}'))

【実行結果】
$${x}$$ が変数の形式で残っています！

可視化です！
sympy 様式の関数 integ_func2 を活用します。
$${x=2}$$ で固定して $${k}$$ の関数とし、$${0 \leq k \leq 10}$$ の範囲で描画します。

# 可視化（x=2の場合）
spb.plot(integ_func2.subs(x, 2), (k, 0, 10));

【実行結果】

こちらは、$${k=2}$$ で固定して $${x}$$ の関数にしています。

# 可視化（k=2の場合）
spb.plot(integ_func2.subs(k, 2));

【実行結果】

ぜひ、$${a, k, x}$$ の設定値を変えて、Python を動かしてみましょう！
特別な何かを発見できるかもしれません🍀

おわり

この記事がシリーズ最終話となります。

著者の大栗博司先生の親切な文章のおかげで、数学を身近に感じることができました。
数学への恐怖心も少し和らぎました。
先生、誠にありがとうございました✨️

◆ ◆ ◆

このブログは、書籍の豊富なテーマから、個人的に Python で書きたくなったものを（気分次第で）選び取って始めたシリーズです。
書籍の実践と同時並行で綴った Python コードの全てを出し切りました。

将来、まだ手つかずのテーマを Python 化したくなったら、シリーズが再開するかもしれません。
そのときまで、しばらくお別れです。

最後までお読みいただき、ありがとうございました🍀

シリーズの記事

前の記事

目次

ブログの紹介

note で７つのシリーズ記事を書いています。
ぜひ覗いていってくださいね！

１．のんびり統計

統計検定２級の問題集を手がかりにして、確率・統計をざっくり掘り下げるブログです。
雑談感覚で大丈夫です。ぜひ覗いていってくださいね。
統計検定２級公式問題集CBT対応版に対応しています。
Python、EXCELのサンプルコードの配布もあります。

２．実験！たのしいベイズモデリング１＆２をPyMC Ver.5で

書籍「たのしいベイズモデリング」・「たのしいベイズモデリング２」の心理学研究に用いられたベイズモデルを PyMC Ver.5で描いて分析します。
この書籍をはじめ、多くのベイズモデルは Ｒ言語＋Stanで書かれています。
PyMCの可能性を探り出し、手軽にベイズモデリングを実践できるように努めます。
身近なテーマ、イメージしやすいテーマですので、ぜひぜひPyMCで動かして、一緒に楽しみましょう！

３．実験！岩波データサイエンス1のベイズモデリングをPyMC Ver.5で

書籍「実験！岩波データサイエンスvol.1」の４人のベイジアンによるベイズモデルを PyMC Ver.5で描いて分析します。
この書籍はベイズプログラミングのイロハをざっくりと学ぶことができる良書です。
楽しくPyMCモデルを動かして、ベイズと仲良しになれた気がします。
みなさんもぜひぜひPyMCで動かして、一緒に遊んで学びましょう！

４．楽しい写経 ベイズ・Python等

ベイズ、Python、その他の「書籍の写経活動」の成果をブログにします。
主にPythonへの翻訳に取り組んでいます。
写経に取り組むお仲間さんのサンプルコードになれば幸いです🍀

５．RとStanではじめる心理学のための時系列分析入門 を PythonとPyMC Ver.5 で

書籍「RとStanではじめる心理学のための時系列分析入門」の時系列分析をPythonとPyMC Ver.5 で実践します。
この書籍には時系列分析のテーマが盛りだくさん！
時系列分析の懐の深さを実感いたしました。
大好きなPythonで楽しく時系列分析を学びます。

６．データサイエンスっぽいことを綴る

統計、データ分析、AI、機械学習、Pythonのコラムを不定期に綴っています。
統計・データサイエンス書籍にまつわる記事が多いです。
「統計」「Python」「数学とPython」「R」のシリーズが生まれています。

７．Python機械学習プログラミング実践記

書籍「Python機械学習プログラミング PyTorch & scikit-learn編」を学んだときのさまざまな思いを記事にしました。
この書籍は、scikit-learnとPyTorchの教科書です。
よかったらぜひ、お試しくださいませ。

最後までお読みいただきまして、ありがとうございました。