まずはお絵描きしてみよう／Python で数学を

このサイトは学部では早稲田で物理を, 修士では東大で数学を専攻し, 今も非アカデミックの立場で数学や物理と向き合っている一市民の奮闘の記録です. 運営者情報および運営理念についてはこちらをご覧ください.

中高の数学の復習から専門的な数学・物理までいろいろな情報を発信しています.
中高数学に関しては自然を再現しようや役に立つ中高数学中高数学お散歩コース
大学数学に関しては現代数学観光ツアーなどの無料の通信講座があります.
その他にも無料の通信講座はこちらのページにまとまっています.
ご興味のある方はぜひお気軽にご登録ください!

これは Python, 特に numpy, scipy, matplotlib を使った,
数学や物理ネタと Python プログラミングに関する記事だ.
これから勉強の記録をつけること,
それを共有することに目的に記事を書いていく.
GitHub にも関係するコードを置いていくし,
このサイトでもタグで記事を一覧できるようにしていく.

ご興味ある方はぜひ継続的にチェックしてほしい.

readme を書いた時点では次のバージョンで動かしている.

Python 3.5.1 :: Anaconda 4.1.0 (x86_64)

1 諸注意
2 今回の目的: お絵描きして遊んでみよう
3 式をもとにグラフを描く
4 次回予告
- 4.1 共有:

諸注意

前回の記事Jupyter のインストール／Python で数学をで Jupyter をインストールしていることを前提にしている.
今回の分は GitHub のこのページを見てもらうと早い.

グラフ描画の前半はそこまでのコードを Jupyter 上で実行して,
ipython にデータを読み込ませていることが前提になっている.
そこだけ実行してもうまく動かないことがあるから注意すること.
該当のコード片は Ctrl+Enter で実行できる.
GitHub 上の ipynb をダウンロードし,
Jupyter を実行してローカルで実行しながら遊ぶのが楽だろう.

またコードの 2 重管理はしんどいので,
いったんここに書いたコードは更新しない.
最新版コードは直接 GitHub を見てほしい.
時間が経つと動かなくなることもあるだろう.
そのときはご指摘頂けると嬉しい.

今回の目的: お絵描きして遊んでみよう

前置きが長くなった.
今回の目的は実際にどんなことができるかを知るために
まずはお絵描きして遊んでみるのが目的だ.

初学者向けには Python のコードの書き方がどうの,
ライブラリの呼び出しがどうのといった細々とした注意は必要だ.
ただそんなのをやっていても面白くない.
とりあえずは目で見てすぐわかるところを見てみようという趣旨だ.

まずはちょっとしたアニメーションサンプルからはじめよう.

%matplotlib nbagg

import numpy as np
import matplotlib.pyplot as plt
import matplotlib.animation as animation

fig = plt.figure()
x = np.arange(0, 10, 0.1)

ims = []
for a in range(50):
    y = np.sin(x - a)
    im = plt.plot(x, y, "r")
    ims.append(im)

ani = animation.ArtistAnimation(fig, ims)
plt.show()

GitHub 上では静止画になってしまっているが,
自分のマシンで Jupyter を立ち上げて実行すればきちんとアニメーションしてくれる.
これなら高校物理のアニメーションも多少はやれそうな希望も湧いてくる.

一応触れておくと numpy の処理が割とすごいことに注目してほしい.
y = np.sin(x - a) と書いただけで y は配列として定義されるのだ.
あとで詳しく見てみよう.

そのあとしばらくいろいろなグラフを並べている.
折れ線を描いたりプロットしたり,
軸を設定したり複数の折れ線を描いたり.

これは単にグラフに関する機能の私の備忘録と思ってほしい.
実際にレポートや論文を書くときには必須の処理だ.
単に遊びたいだけならなくてもいい.

式をもとにグラフを描く

今回の本題はここ.
まず逆 2 乗則をプロットした.

%matplotlib inline
import matplotlib.pyplot as plt

<!-- 距離の変動範囲 -->
r = range(100, 1001, 20)
<!-- 重力定数 -->
G = 6.674*(10**-11)
<!-- 質量 -->
m1 = 0.5
m2 = 1.5
<!-- 重力 -->
F = [G * (m1 * m2) / (dist ** 2) for dist in r]

plt.plot(r, F, marker='o')
<!-- Jupyter 上日本語が文字化けするようなので英語で書いた -->
plt.xlabel('Distance')
plt.ylabel('Newton gravitational force')
plt.title('A graph for the inverse square law')
plt.show()

こんな程度でもはじめて触れると割と感動する.
グラフ描くプログラムを書くのは本当に面倒というイメージがあるし.

重力の値の変化を追うリスト F は numpy を使って書く方がいいのだろう.
書ける人は numpy を使って書き直してみよう.

もう 1 つ, 物理の基本中の基本という感じで放物線に関するグラフも書いている.
ここにコードを載せるのは後者の 3 つの放物線のグラフにしておこう.

v0_list = [20, 40, 60]
theta = 45
tragectories = []

def frange(start, final, interval):
    numbers = []
    # float だと range が使えないのでこの処理
    while start < final:
        numbers.append(start)
        start = start + interval

    return numbers

def get_trajectory(v0, theta):
    # 重力定数の定義
    g = 9.8

    # 入力の角度をラジアンに変換
    theta = math.radians(theta)

    # 着地までの時間の計算
    t_flight = 2 * v0 * math.sin(theta) / g

    # 着地時間までの時間リスト: dt 刻み
    dt = 0.001
    intervals = frange(0, t_flight, dt)

    # x, y 座標の計算
    x = []
    y = []
    for t in intervals:
        x.append(v0*math.cos(theta)*t)
        y.append(v0*math.sin(theta)*t - 0.5*g*t*t)

    return (x,y)

for v0 in v0_list:
    tragectories.append(get_trajectory(v0, theta))

<!-- Add a legend and show the graph -->
plt.plot(tragectories[0][0], tragectories[0][1],
         tragectories[1][0], tragectories[1][1],
         tragectories[2][0], tragectories[2][1])
plt.legend(['20', '40', '60'])
plt.show()

最初のコードは初期速度を入力させてそのグラフを描く処理なので,
コードが面倒くさくなっている.

また上から順番に実行していく前提で Jupyter notebook を作っているので,
本当は上のコードのように get_trajectory や
frange も入れておかなければ動かない.

このコードはまだ numpy をよく知らない状態で書いているので,
frange という関数を作っている.
これは np.arange メソッドでいけるはずだ.
まだよくわかっていないので,
勉強が進み次第書き直そうと思っている.

初回から割とごついコードが出てきてしまった.
とりあえずこんなことができるようになるのか,
という感じで気楽に付き合ってほしい.

次回予告

次回は Jupyter の使い方,
Python や IPython の基礎について触れる予定.
面白くはないがある程度すっきりまとめておいた方が私自身のためでもある.

次回をお楽しみに.