「数学の言葉で世界を見たら」をPython で見たら… Vol.10 積分の気分
第７話「微積は積分から」
書籍の著者 大栗博司 先生
書籍「数学の言葉で世界を見たら」第７話「微積は積分から」の Python写経活動記録 です。
関数 $${f(x)}$$ の区間 $${[a, b]}$$ の面積を求める過程を通じて 積分の気持ち に近づきます。
では書籍を開いて数学の旅に出発です🚀
はじめに
このブログシリーズは、書籍「数学の言葉で世界を見たら」（幻冬舎）で学んだ「数学の楽しさ」を「Python 写経の形式」でご紹介いたします。
【引用表記】
この記事は、出典に記載の書籍に掲載された文章とデータを引用し、適宜、掲載文章・データを改変して書いています。
【出典】
「数学の言葉で世界を見たら」 第3刷、著者 大栗博司、幻冬舎
４ 「積分」では何を計算しているのか
学びポイント
面積計算を通じて積分を体感シリーズの２作目です。
関数などで囲った形状の面積（＝積分）について、微小な幅の長方形をたくさん並べて計算します。
この記事で用いるライブラリをインポートします。
## インポート
# 数値計算
import numpy as np
# 数値積分
from scipy.integrate import quad
# 可視化
import matplotlib.pyplot as plt
from matplotlib import patches
plt.rcParams['font.family'] = 'Meiryo' # または import japanize_matplotlib主人公「図形Ａ」
「単調に増加する関数 $${f(x)}$$ の曲線」、「区間 $${[a, b]}$$ の両端からに伸びる垂直線」、「$${y=0}$$ の水平線 」で囲まれた図形が今回主役の 図形 $${A}$$ です。
と言葉の説明は混乱の原因なので…
可視化で一目瞭然！といきましょう。
図形 $${A}$$ を描画します。
書籍 p.222 図 7-9 に似せてみました。
### p.222 図 7-9 グラフの下の面積～図形A
## 関数の定義 f(x) = 1.1x³ - 4x² + 6.5x
f = lambda x: 1.1 * x**3 -4 * x**2 + 6.5 * x
## 可視化
# 描画領域の設定
fig, ax = plt.subplots(figsize=(4.2, 3.9))
# 区間の設定
a, b = 0.5, 2
# 関数fの曲線の描画
ax.plot(x_val:=np.linspace(a, b, 1001), f(x_val), lw=3)
# 区間[a, b]の面の塗りつぶし
ax.fill_between(x_val2:=np.linspace(a, b, 101), f(x_val2), alpha=0.1)
# 区間[a, b]の両端の垂直線・下端の水平線の描画
ax.vlines(a, 0, f(a), color='tab:blue')
ax.vlines(b, 0, f(b), color='tab:blue')
ax.hlines(0, a, b, color='tab:blue')
# 図形名の表示
ax.text(x=1.25, y=1.5, s='図形Ａ', ha='center', fontsize=16)
# 修飾
ax.set_title('関数 $f\ (x) = 1.1x^3 - 4x^2 +6.5x$ のグラフ')
ax.set_xlabel('$x$', fontsize=12)
ax.set_ylabel('$f\ (x)$', fontsize=12)
ax.set(xlim=(0, 2.5), ylim=(0, 7));【実行結果】
太い曲線は 関数 $${f(x)=1.1x^3 - 4x^2+6.5x}$$ です。
区間 $${[0.5, 2.0]}$$ と $${y=0}$$ で囲まれた青く塗られた部分が 図形 $${A}$$ です。
図形 $${A}$$ の面積を「アルキメデスのはさみうち」で計算します。
「アルキメデスのはさみうち」は前回記事でご紹介した方法です。
図形Ｂと図形Ｃではさみうち
はさみうちの計算イメージを図示しましょう。
書籍 p.223 図 7-10 に相当します。
最初に描画関数を定義します。
関数 $${f(x)}$$ を屋根にする $${n}$$ 個の長方形を横並びに描画します。
### p.223 図 7-10 アルキメデスのはさみうちでグラフの下の面積を計算する
# 図形 Bn, Cn 描画関数の定義：引数 分割数 n
def plot_integrate(n):
## 設定と準備
# 区間の設定
a, b = 0.5, 2
# 関数の定義 f(x) = 1.1x³ - 3x² + 6.5x
f = lambda x: 1.1 * x**3 -4 * x**2 + 6.5 * x
# 長方形作成ヘルパー関数の定義
def plot_rectangle(start_x, height, width, fc='none', ec='none', alpha=1):
return patches.Rectangle(
xy=(start_x, 0), # 左下の座標
width=width, # 幅
height=height, # 高さ
facecolor=fc, # 塗りつぶし色
edgecolor=ec, # 枠線色
alpha=alpha, # 透明度
)
## 描画
# 描画領域の設定
fig, axes = plt.subplots(1, 2, figsize=(8, 4), tight_layout=True)
# 図形B_nとC_nの描画
for i, (fig_name, ax) in enumerate(zip(['$B_n$', '$C_n$'], axes.flat)):
# 関数fの曲線の描画
ax.plot(x_val:=np.linspace(a, b, 1001), f(x_val), lw=3)
# n 個の長方形の描画
for j in range(n):
# 幅の算出
width = (b - a) / n
# 始点の算出
start_x = a + width * j
# 長方形の高さの算出
if i == 0:
height = f(start_x)
else:
height = f(start_x + width)
# 長方形の塗りつぶしの描画
rectangle = plot_rectangle(
start_x, height, width, 'tab:blue', 'none', 0.1)
ax.add_patch(rectangle)
# 長方形の辺の描画
rectangle = plot_rectangle(
start_x, height, width, 'none', 'tab:blue', 1)
ax.add_patch(rectangle)
# 区間[a, b]の両端の垂直線の描画
ax.vlines(a, 0, f(a), color='tab:blue', lw=1)
ax.vlines(b, 0, f(b), color='tab:blue', lw=1)
# 修飾
ax.set_title(f'図形 {fig_name}, $n$={n}', fontsize=14)
ax.set_xlabel('$x$', fontsize=12)
ax.set_ylabel('$f\ (x)$', fontsize=12)
ax.set(xlim=(0, 2.5), ylim=(0, 7))
plt.show();【実行結果】なし
1️⃣ $${n=4}$$ のケース
$${n}$$ は長方形を横方向に分割する数です。
最初に４分割のときの 図形 $${B}$$、図形 $${C}$$ を描きます。
４個の長方形が現れます。
# n=4 の場合
plot_integrate(4)【実行結果】
図形 $${B}$$ では、長方形は関数 $${f(x)}$$ の下側で接しています。
図形 $${C}$$ では、長方形は関数 $${f(x)}$$ の上側で接しています。
長方形たちの面積の合計が図形 $${A}$$ に近似する「はず」なのですが…
図形 $${A}$$ の面積は、関数 $${f(x)}$$ の「滑らかな」曲線に沿って計算する必要があります。
分割の数 $${n}$$ をどんどん大きくして、長方形の幅を小さくすれば、長方形群が形成する「上側のギザギザ」が滑らかになるはずです。
2️⃣ $${n=10}$$ のケース
# n=10 の場合
plot_integrate(10)【実行結果】
上部の少しギザギザが取れてきている印象です。
3️⃣ $${n=20}$$ のケース
# n=20 の場合
plot_integrate(20)【実行結果】
上部がかなり滑らかになり、長方形群が図形 $${A}$$ にとても似てきている印象です！
このまま $${n}$$ を大きくすれば（極限！）、図形 $${B}$$ も図形 $${C}$$ も、図形 $${A}$$ にぴったり当てはまるでしょう。
「はさみうち」は $${n}$$ を増やして、小さい方の図形 $${B}$$ の面積と、大きい方の図形 $${C}$$ の面積の「差」をキリキリと小さく詰めることで、図形 $${A}$$ の面積に近づけていく考え方です。
面積算出シミュレーション
分割数 $${n}$$ を大きくして、図形 $${C}$$ と図形 $${B}$$ の長方形面積の差が僅少（$${10^{-6}}$$ 未満）になったところで面積算出完了にするシミュレーションを行います。
どのくらい分割すれば図形 $${A}$$ の面積と言えるのでしょう？？？
%%time
## 面積算出シミュレーションの実行
## 設定
# 区間
a, b = 0.5, 2
# 関数の定義 f(x) = 1.1x³ - 3x² + 6.5x
f = lambda x: 1.1 * x**3 -4 * x**2 + 6.5 * x
## 探索区間を設定の上、面積算出シミュレーションを実行
for n in range(5118750, 5118755):
# 長方形の横の長さεの算出
eps = (b - a) / n
# n+1個の分割点ごとのf(x)の値の算出
fx = f(np.linspace(a, b, n + 1))
# Bnの面積の算出
Bn = sum([fx[i] for i in range(n)]) * eps
# Cnの面積の算出
Cn = sum([fx[i] for i in range(1, n + 1)]) * 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)【実行結果】
分割数 $${n=5118751}$$ です！
とても小さな幅になっています。
このときの面積は
・図形 $${B}$$：$${6.0703120}$$
・図形 $${C}$$：$${6.0703130}$$
になりました。
平均して $${6.0703125}$$ を図形 $${A}$$ の面積といたしましょう。
テキスト p.224 のまとめを引用いたします。
こうして計算した図形 $${A}$$ の面積は、「関数 $${f(x)}$$ の、区間 $${a \leq x \leq b}$$ での積分」と呼ばれ、
$$
\int_a^b f(x) dx
$$
と表される。
「積分」の誕生です。
Python のライブラリで面積を計算
Python のライブラリを使えば、サクッと積分計算ができます。
scipy.integrate の数値積分関数 quad() を利用して、積分で面積を計算しましょう。
引数は、関数、区間 $${a, b}$$ です。
## scipyで数値積分
## 設定
# 区間
a, b = 0.5, 2
# 関数の定義 f = 1.1x³ - 3x² + 6.5x
f = lambda x: 1.1 * x**3 -4 * x**2 + 6.5 * x
## 数値積分の実行
result = quad(f, a=0.5, b=2)
print(f'数値積分による面積: {result[0]:.7f}')【実行結果】
シミュレーションで求めた面積と一致しました！
おわり
シリーズの記事
次の記事
前の記事
目次
ブログの紹介
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の教科書です。
よかったらぜひ、お試しくださいませ。
最後までお読みいただきまして、ありがとうございました。
いいなと思ったら応援しよう！
応援ありがとうございます。これからもがんばって記事を作成します！