「数学の言葉で世界を見たら」をPython で見たら… Vol.7 ケプラーの第３法則

第３話「大きな数だって怖くない」

書籍の著者　大栗博司 先生

書籍「数学の言葉で世界を見たら」第３話「大きな数だって怖くない」の Python写経活動記録 です。 

「宇宙」にまつわる数学テーマです。
ケプラーの第３法則と対数の関係を Python で楽しく体感します。

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

はじめに

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

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

６ 自然の法則は「対数」で見抜ける

学びポイント

太陽系の惑星を一直線に並べる「対数グラフ」の威力を目の当たりにします。

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

## インポート

# 数値計算
import pandas as pd

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

ケプラーの第３法則

先生は、ケプラーが長い歳月を費やして発見した「ケプラーの第３法則」を優しく教えてくれます。
「太陽系惑星の軌道長半径と公転周期」の関係を次のように示す法則です。

① 両辺の正の平方根をとると

② 両辺の常用対数 $${\log_{10}}$$ をとると

③ 右辺の指数を前に出すと

③ の式が示す直線は傾きが $${3/2}$$ になるようです。
え？可視化しなきゃ！

惑星の軌道長半径と公転周期のデータ

可視化の前にデータを用意しなきゃ、とネットを調べたところ…
運良く 仙台市天文台 様 のサイトで 2007 年理科年表の惑星データを見つけました。
引用させていただきます。
ありがとうございます！

## p.98 惑星データの設定
# データ引用サイト:仙台市天文台 (理科年表2007年に基づく)
# https://www.sendai-astro.jp/nishikouen/photo/planet/planet_index.html

# データの登録 ※未使用項目はコメントアウトしてます
df = pd.DataFrame({
    '惑星名': ['水星', '金星', '地球', '火星', '木星', '土星', '天王星', '海王星'],
    # '赤道半径(km)':  [2440, 6052, 6378, 3396, 71492, 60268, 25559, 24764],
    # '体積（地球=1）':   [0.056, 0.857, 1.000, 0.151, 1321.0, 755.0, 63.0, 58.0],
    # '密度(g/cm3)': [5.43, 5.24, 5.52, 3.93, 1.33, 0.69, 1.27, 1.64],
    '軌道長半径(億km)': [0.579, 1.082, 1.496, 2.279, 7.783, 14.294, 28.75, 45.044],
    '公転周期(年)': [0.24, 0.62, 1, 1.88, 11.86, 29.46, 84.02, 164.77],
    # '自転周期(日)': [58.65, 243.02, 0.997, 1.026, 0.414, 0.444, 0.718, 0.671],
    # '衛星の数':  [0, 0, 1, 2, 63, 59, 27, 13],
})

# 軌道長半径を地球=1にして規格化（テキストに合わせる）
df['軌道長半径(地球=1)'] = df['軌道長半径(億km)'] / df.loc[2, '軌道長半径(億km)']
df = df.iloc[:, [0, 1, 3, 2]]

# データフレームの表示
df.round(3)

【実行結果】
軌道長半径は地球＝１を基準にして規格化しています。

いよいよ可視化

1️⃣ データそのままプロット
表の「軌道長半径（地球=1）」と「公転周期（年）」を「そのままのスケール」で可視化します。
テキスト p.98 図 3-1 に相当します。

## p.98 図3-1 の描画

# 描画領域の設定
plt.figure(figsize=(7, 7))
# 軌道長半径と公園周期の散布図の描画
plt.plot(df['軌道長半径(地球=1)'], df['公転周期(年)'], '--o', lw=1)
# 惑星名の表示
for _, (x, y, s) in df[['軌道長半径(地球=1)', '公転周期(年)', '惑星名']].iterrows():
    plt.text(x=x, y=y-5, s=s, ha='left', va='top')
# 修飾
plt.title('惑星の軌道半径と公転周期の関係')
plt.xlabel('軌道長半径 [地球=1]', fontsize=12)
plt.ylabel('公転周期 [年]', fontsize=12)
plt.ylim(-20, 180);

【実行結果】
欲張って海王星まで含めてしまいました！

全体的な雰囲気は「緩やかな曲線」でしょうか。
地球に近い惑星は団子のように偏っています。
天王星・海王星は軌道長半径が長すぎて、直線寄りに見えてしまいます。

2️⃣ 対数スケールでプロット
続いて、横軸・縦軸ともに常用対数 $${\log_{10}}$$ スケールで可視化します。
書籍 p.98 図 3-2 に相当します。

matplotlib の xscale、yscale で 軸のスケールを log, base=10 と指定するだけで、データを対数変換しなくても、対数スケールのチャートを描画できます。

## p.98 図3-2 の描画

# 描画領域の設定
plt.figure(figsize=(7, 7))
# 軌道長半径と公園周期の散布図の描画
plt.plot(df['軌道長半径(地球=1)'], df['公転周期(年)'], '--o', lw=1)
# 惑星名の表示
for _, (x, y, s) in df[['軌道長半径(地球=1)', '公転周期(年)', '惑星名']].iterrows():
    plt.text(x=x*1.03, y=y*0.9, s=s, ha='left', va='top')
# x=1, y=1の垂直線、水平線の描画
plt.axhline(1, color='black', lw=0.5)
plt.axvline(1, color='black', lw=0.5)
# x,y軸を対数log10スケールに変換
plt.xscale('log', base=10)
plt.yscale('log', base=10)
# 修飾
plt.title('惑星の軌道半径と公転周期の関係：対数スケール')
plt.xlabel('log(軌道長半径)', fontsize=12)
plt.ylabel('log(公転周期)', fontsize=12);

【実行結果】
おおおおお！惑星が一直線に並んでいます。しかもいい感じの間隔で。
美しい自然の法則です！
「③ の式の直線は傾きが $${3/2}$$ になる」です！

ケプラーの第３法則を計算する

①の式が成り立っているか、データで確認しましょう。

### p.99 惑星の公転周期＝起動長半径の3/2乗（p.99のケプラーの第３法則）
df['軌道長半径の3/2乗≒公転周期'] = df['軌道長半径(地球=1)']**(3/2)
df.round(3)

【実行結果】

右の２列に注目しましょう。
公転周期と軌道長半径の3/2乗はとても近い値になっています！
ケプラーの第３法則、すごいですね！

太陽系の天体の位置をWebサイトで

太陽系惑星が動く様子を見てみたい！とWebサイトを調べたところ…
STUDIO KAMADA 様の「太陽系の天体の位置」サイトを見つけました！

日時を指定して太陽系の天体の位置を見たり、「動く」天体シミュレーションを堪能できる素敵なサイトです！
じーっと見ていて、とても楽しい気持ちになりました。
ありがとうございます！
ぜひ皆さんもシミュレーションを体感して下さい！！！

ある日の天体の位置の画面ショットを引用いたします。

太陽系の惑星のとても大きな数も怖くなかったです！

おわり

シリーズの記事

次の記事

前の記事

目次

ブログの紹介

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の教科書です。
よかったらぜひ、お試しくださいませ。

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