「数学の言葉で世界を見たら」をPython で見たら… Vol.9 アルキメデスのはさみうち
第７話「微積は積分から」
書籍の著者 大栗博司 先生
書籍「数学の言葉で世界を見たら」第７話「微積は積分から」の Python写経活動記録 です。
アルキメデスのはさみうち で円の面積を求めます。
では書籍を開いて数学の旅に出発です🚀
はじめに
このブログシリーズは、書籍「数学の言葉で世界を見たら」（幻冬舎）で学んだ「数学の楽しさ」を「Python 写経の形式」でご紹介いたします。
【引用表記】
この記事は、出典に記載の書籍に掲載された文章とデータを引用し、適宜、掲載文章・データを改変して書いています。
【出典】
「数学の言葉で世界を見たら」 第3刷、著者 大栗博司、幻冬舎
３ 「アルキメデスのはさみうち」
学びポイント
円と正多角形と三角形と面積が混ざり合う合理的な世界を堪能できます。
この記事で用いるライブラリをインポートします。
## インポート
# 数値計算
import numpy as np
import pandas as pd
# 可視化
import matplotlib.pyplot as plt
from matplotlib import patches
plt.rcParams['font.family'] = 'Meiryo' # または import japanize_matplotlibアルキメデスのはさみうち
この節は「アルキメデスのはさみうち」で図形の面積を計算します。
アルキメデスのはさみうちは、だいたいこんな感じ…
面積をはかりたい図形をＡとします。
図形Ａに対して「内側で接する折れ線図形Ｂ」と「外側で接する折れ線図形Ｃ」で図形Ａを「はさみうち」します。
図形Ｂと図形Ｃは三角形の組み合わせで構成されるので、三角形の面積を求めれば、図形Ｂ・Ｃの面積が求まります。
図形Ｂと図形Ｃの頂点をどんどん増やして、図形Ａへのフィット率を高めます。
図形Ｂの面積と図形Ｃの面積の差が無くなったとき、図形Ａ＝図形Ｂ＝図形Ｃになります。
書籍 p.220 の例をお借りして、図示します。
図形Ａは円、図形Ｂは円に内接する正多角形、図形Ｃは円に外接する正多角形です。
図形Ａ、Ｂ、Ｃを Python で描いて、２つの正多角形が円を内側と外側から「はさみうち」する様子を確認しましょう！
描画関数を定義します。
## 円に内接する正多角形の描画関数の定義 p.220 図 7-7 対応
# 参考サイト
# 正多角形の1辺の長さ https://keisan.casio.jp/exec/system/1166416582
# 内接円の半径 https://manabitimes.jp/math/2854
def plot_polygon_and_circle(num_Vert, rad_circ=1):
## 描画領域の設定
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(8, 5),
sharex=True, sharey=True, tight_layout=True)
## 左のプロット：円に内接する正多角形 Bn の描画
# 正多角形の描画
poly1 = patches.RegularPolygon(
xy=(0, 0), # 中心座標
numVertices=num_Vert, # 頂点数
radius=rad_circ, # 半径
orientation=np.pi / 6, # 回転角度（ラジアン）
facecolor='none', # 塗りつぶし色
edgecolor='black', # 枠線色
)
ax1.add_patch(poly1)
# 円の描画
circ1 = patches.Circle(
xy=(0, 0), # 中心座標
radius=rad_circ, # 半径
facecolor='none', # 塗りつぶし色
edgecolor='tab:red', # 枠線色
linewidth=2, # 枠線の太さ
)
ax1.add_patch(circ1)
## 右のプロット：円に外接する正多角形 Cn の描画
# 正多角形の半径の算出
a = rad_circ * (2 * np.tan(np.pi / num_Vert)) # 正多角形の１辺の長さ
rad_poly = a / (2 * np.sin(np.pi / num_Vert)) # 正多角形の半径の長さ
# 正多角形の描画
poly2 = patches.RegularPolygon(
xy=(0, 0), # 中心座標
numVertices=num_Vert, # 頂点数
radius=rad_poly, # 半径
orientation=np.pi / 6, # 回転角度（ラジアン）
facecolor='none', # 塗りつぶし色
edgecolor='black', # 枠線色
)
ax2.add_patch(poly2)
# 円の描画
circ2 = patches.Circle(
xy=(0, 0), # 中心座標
radius=rad_circ, # 半径
facecolor='none', # 塗りつぶし色
edgecolor='tab:red', # 枠線色
linewidth=2, # 枠線の太さ
)
ax2.add_patch(circ2)
## 共通修飾
for ax_, tex1, tex2 in zip([ax1, ax2], ['内接', '外接'], ['$B_n$', '$C_n$']):
ax_.set(xticks=[], yticks=[],
title=f'円 $A$ に{tex1}する正{num_Vert}角形 {tex2}')
ax_.set_aspect('equal', adjustable='box')
ax_.autoscale()
ax_.axis('off')
plt.show()【実行結果】なし
正四角形から徐々に頂点の数を増やしていきます。
1️⃣ 正四角形
# 正四角形の場合
plot_polygon_and_circle(4)【実行結果】
内接する正四角形 $${B}$$ は円の内側に位置して、頂点が円と接します。
外接する正四角形 $${C}$$ は円の外側に位置して、辺に円が接します。
正四角形では、円の面積と $${B}$$・$${C}$$ の面積が、あまり近くない印象です。
2️⃣ 正八角形
# 正八角形の場合
plot_polygon_and_circle(8)【実行結果】
正八角形では、円と $${B, C}$$ の面積がグッと近づいた印象です。
3️⃣ 正十六角形
# 正十六角形の場合
plot_polygon_and_circle(16)【実行結果】
正十六角形になると、円と $${B, C}$$ の線がかなり被っていて、面積もいい感じに近似しているように見えます。
$${B, C}$$ の面積の大小は、外側から接する $${C}$$ の方が大きいです。
このようにして正 $${n}$$ 角形の $${n}$$ を増やしていけば、いつしか（極限で） $${B, C}$$ は一致するはずです。
一致するときの面積が円の面積になります。
このように、面積をはかりたい図形があって、内側と外側から直線を辺でもつ図形で「はさみうち」して、図形の面積を計算するのが「アルキメデスのはさみうち」なのです！
【謝辞】
図形の計算について、以下のWebサイトの情報を参考にいたしました。
ありがとうございます！
◆ 正多角形の１辺の長さ
◆ 内接円の半径の長さ
正多角形のはさみうちで円の面積を計算
円 $${A}$$ の面積を、内接する正多角形 $${B}$$ と外接する正多角形 $${C}$$ ではさみうちして、計算しましょう。
もちろん Python で。
正４角形 ⇒ 正８角形 ⇒ 正16角形 のように、$${2^{n+1}}$$ で頂点の数を増やします。
$${B}$$ と $${C}$$ の面積の差が僅少になったときの面積で、円の面積を推定します。
このコードでは、円の半径を１、「僅少」の基準を誤差 $${10^{-10}}$$ にしています。
また、正多角形の面積は、頂点と中心を直線で結んでできる「三角形」の面積を計算して、頂点の数＝三角形の数だけ足し合わせて求めています。
### CnとBnのはさみうちによる面積の算出
## 設定
# 半径の長さ
radius = 1
# Bn用：三角形の面積算出関数の定義 ※2辺の長さとその角度（ラジアン）から算出
triangle_area_Bn = lambda a, b, θ: 1/2 * a * b * np.sin(θ)
# Cn用：三角形の面積算出関数の定義 ※高さと円の中心に接する角度の数から算出
triangle_area_Cn = lambda num_corner, radius: radius * np.tan(np.pi / num_corner)
# 結果を格納するリストの初期化
num_corners, area_Bn, area_Cn = [], [], []
## 正多角形の面積を計算
for n in range(1, 28):
# 正多角形の角の数 2^{n+1}の算出
num_corner = 2**(n + 1)
# 円に内接する正多角形Bnの面積の算出
θ = 2*np.pi / num_corner # 円の中心と接する三角形の角の角度の算出
est_area_Bn = triangle_area_Bn(radius, radius, θ) * num_corner
# 円に外接する正多角形Cnの面積の算出
est_area_Cn = triangle_area_Cn(num_corner, radius) * num_corner
# 結果の格納
num_corners.append(num_corner)
area_Bn.append(est_area_Bn)
area_Cn.append(est_area_Cn)
# 差が一定値未満になったら推定を終了
if abs(est_area_Cn - est_area_Bn) < 1e-10:
break
## データフレーム化して表示
result = pd.DataFrame({
'多角形': num_corners, 'Cn面積': area_Cn, 'Bn面積': area_Bn,
'差' : np.array(area_Cn) - np.array(area_Bn)})
# 小数点10桁で表示
with pd.option_context('display.float_format', '{:.10f}'.format):
display(result)
## 円の面積の表示（正解値）
circle_area = radius**2 * np.pi
print(f'円: {circle_area:.10f}')【実行結果】
あっさり求まりました。
正 1048576 角形で、面積は 3.1415926536 です。
おわり
シリーズの記事
次の記事
前の記事
目次
ブログの紹介
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の教科書です。
よかったらぜひ、お試しくださいませ。
最後までお読みいただきまして、ありがとうございました。
いいなと思ったら応援しよう！
応援ありがとうございます。これからもがんばって記事を作成します！