コンテンツにスキップ

Python

ガイド

Pythonの様々なバージョンを導入

Amazon Linuxpyenv

  • Gitの確認:なければsudo yum install git -yで導入
1
git --version
  • pyenvgccをインストール
1
2
git clone https://github.com/pyenv/pyenv.git ~/.pyenv
sudo yum install -y gcc openssl11 openssl11-devel bzip2-devel ncurses-devel libffi-devel readline-devel sqlite-devel.x86_64 xz-devel
  • ~/.bashrcに以下を追記
1
2
export PATH="$HOME/.pyenv/bin:$PATH"
eval "$(pyenv init -)"
  • ~/.bashrcを再読み込み
1
source ~/.bashrc
  • pyenvの確認
1
pyenv --version
  • 導入できるPythonの確認
1
pyenv install -l
  • 適切なバージョンをインストール
1
pyenv install 3.11.7
  • 現状を確認
1
pyenv versions
  • バージョン指定
1
pyenv global 3.11.7
  • 確認
1
python -V

Macasdf

  • URL
  • asdfをインストールする
    • パスを通す
1
brew install asdf
  • 必要なバージョンが入っているか確認
1
asdf list all python
  • Python のプラグインとリポジトリを確認する
1
asdf plugin list all | grep python
  • 上記出力結果を参考にプラグインをインストールする
1
asdf plugin add python https://github.com/danhper/asdf-python.git
  • インストールできたか確認
1
asdf plugin list
  • インストールしたいバージョンに対応しているか確認
1
asdf list all python
  • 適切なバージョンを確認してインストールする
1
asdf install python <version>
  • 使いたいバージョンを設定: 必要に応じてローカルで設定する
1
2
asdf global python <version>
asdf local python <version>
  • アンインストールする場合は次のコマンドを実行する
1
asdf uninstall python <version>

MacHomebrew

  • ChatGPTに聞いた
  • 適当なバージョンをインストールする
    • 下記コマンドではpython3.8にしているが, 適切なバージョンを指定すればよい.
1
brew install python@3.8
  • インストールできたか確認
1
python3.8 --version

venvを作る

  • 適切なPythonのバージョンでvenvを作成
1
python -V
  • 適切なディレクトリに移動
  • .venvを作成
1
python -m venv .venv
  • 仮想環境をアクティベート
1
source .venv/bin/activate
  • venvを解除
1
deactivate

Pythonで数学を

Jupyterのインストール/Pythonで数学を

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

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

まずはインストールから説明したい. この記事は 2016-08 時点でのインストール法を紹介している. 情報が古くなっていたら適宜調べ直してほしい.

Jupyter インストール

2016-08 時点ではQiita のこの記事が参考になる.

詳しくはそちらを見てもらうことにして, ここでは簡単にまとめておく.

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

  • Python 3.5.1 :: Anaconda 4.1.0 (x86_64)

ライブラリのバージョンも書いた方がいいのかもしれないが, それはサボる.

Anaconda でのインストール

Anaconda でのインストールがいいらしい. 私もそれで自機の Mac にインストールしている.

Windows (未検証)

公式サイトからダウンロード, インストール.

Mac (検証済み)

Homebrew で次のコマンド実行でインストールした. brew upgrade でかなり時間がかかった. 注意してほしい.

また次のコマンドリストの # が入っている箇所で # 以下は入力しないようにすること. コメントのつもりで書いている.

1
2
3
4
5
6
7
8
9
$ brew update && brew upgrade
$ brew install pyenv
$ echo 'export PATH="$HOME/.pyenv/bin:$PATH"' >> ~/.bash_profile
$ echo 'eval "$(pyenv init -)"' >> ~/.bash_profile
$ exec $SHELL -l
$ pyenv install -l | grep anaconda # Check the newest version
$ pyenv install anaconda3-4.1.0    # Input the newest
$ pyenv global anaconda3-4.1.0     # Input the newest
$ python --version
Linux (未検証)

Qiita のこの記事を参考にしてほしい. わざわざ Linux を使っている人ならむしろ私よりも詳しいだろうし, 簡単に調べられるだろう.

TODO VPython インストール

2016-08 時点で Python3 対応ができていないようだ. 計算物理の文献で VPython を使っているケースを見かけるので, Python3 に対応次第, こちらも追記したい.

Jupyter 起動
シンプルに

本来はコマンドプロンプトやターミナル (いわゆる「黒いやつ」) で起動しないといけない. しかしこのハードルは高いだろうという気もする. そこでダブルクリックすれば自動で起動できるファイルを作った. 次の URL (GitHub) からファイルをダウンロードしてほしい. Windows の方は bat ファイルを, Mac の方は command ファイルを使うこと.

これを適当なフォルダに置いてダブルクリックで実行すれば, そのフォルダで Jupyter が起動する.

少し詳しく

jupyter/fundamental.ipynb にも書いてある. GitHub 上で確認できる. ただあなたはローカルで確認してみたいと思っているかもしれない. そんなときそもそも起動できなければ確認しようもない.

コマンドを簡単に書いておく.

1
2
$ cd some_directory # 適当なディレクトリに移動
$ jupyter notebook

これでそのディレクトリ上で Jupyter が起動する. ファイルを作ったりすると some_directory にファイルが作られる.

Jupyter Notebook については jupyter/fundamental.ipynb にまとめてあるのでとりあえずはそちら参照. 次回以降で詳しく記事を書いていく. おおよそ次のようなことがまとめてある.

  • Markdown 記法が使えること
  • mathjax が使えるので TeX が使えること.
  • 簡単なキーバインドまとめ.
  • Python の基礎の基礎.
  • IPython の基礎の基礎.

簡単なものでもアニメーションやグラフを描けるのが楽しい. jupyter/fundamental.ipynb を軽く眺めたら, 次に jupyter/math_simple_graph.ipynb を見てほしい. あとは適宜書き進めていく.

あなたにも楽しみにしていてほしいし, 何よりも私が今後を楽しみにしている. 一緒に楽しく勉強していこう.

2016 年現在の数学・物理プログラミング状況

ここからはこの記事シリーズを書き始めた経緯について説明する. 読まなくてもこれから先の記事は読んでいける. しかし Python による数学・物理に興味があるあなたにはぜひ読んでほしい内容をまとめている. 最後, 無料の数学の通信講座の案内もしているから, あなたが大学の数学科の数学にご興味があるならぜひ登録してほしい. 登録ページだけは優先して紹介しておく.

最近 Python が科学技術計算でよく使われているようだ. 少なくとも英語の文献では『Python による計算物理』をテーマにした本もいくつか出ている. さらに Python は統計処理や自然言語処理といった分野でも活躍している. Python での自然言語処理では nltk という定評のあるライブラリがあり, それに関する本も出ている.

私が興味があるのはもちろん数学・物理との関係だ. 特に中高生向けのコンテンツとして, 数学や物理と絡めて遊べる方法をずっと模索していた. 昔から SICP や SICM, さらには Functional Diffential Geometry のように, Scheme を使ったプログラミングと数学・物理というネタはある. 最近そこにさらに Haskell と絡めて物理を勉強しようという意欲的なスタイルが出てきていて, Walck によるプレプリントもある[^mathphysprogramming000001].

[^mathphysprogramming000001]: ちなみに 2016-08 に Walck に問い合わせてみたら, これについての本を執筆中とのこと. 出たら買うしかない.

Haskell に興味があるのでそれはそれでやりたいとは思う. しかし資料の豊富さ, ユーザの多さ, グラフなどの視覚化まで絡めたコーディングの簡単さからすると, Python を使った方がよさそうな気がする. というわけで Python でいろいろ探してみた.

Jupyter という宝具

numpy, scipy, matplotlib は有名だ. メンテや機能追加をはじめとした開発は盛んなようだし, それなりに長い間安定的に使われているようだから, そんなに外れもないだろう.

しかしそうしたことよりもさらにすさまじい良い点として, Jupyter の存在を挙げたい. Jupyter については見てもらうのが早い. まず次のページを見てほしい.

GitHub もすごくて, Jupyter の ipynb をちゃんと表示してくれる. コード実行機能があってグラフをきちんと表示してくれる.

Markdown 形式でテキストも書いていけるし, MathJax との連携もある.

きちんと調べてはいないものの, 拡張子の ipynb の ipy は ipython から, nb は notebook から来ているのだと思う. notebook はもともと Mathematica の機能で, 単なるコード共有だけでなくノート形式での共有を目指した機能のようだ. Jupyter はこれを Python でやっている. 実際には Jupyter プロジェクトは Python だけではなく, Ruby や Haskell をも対象として拡大を続けているらしい.

それはそれとして, 式も入れつつ, コード実行もしつつノート形式で共有できるコンテンツ形式はすごすぎる. Web インタフェース (ブラウザ) で書いていくことになり, エディタや IDE の強力な機能は使いづらい. しかし Jupyter 内でもタブの補完が効くようになっているから, ブラウザ上であってもそれなりに書きやすくなっている. あと ipynb はテキストファイル (json?) になっていて, バージョン管理もやりやすい.

個人的には「数学・物理コンテンツメイキングのためにこんなのがほしかった!」というのが ほぼ完璧に入っていて, さすがに感動せざるを得ない.

何はともあれ, しばらくこれを使って数学・物理プログラミングの勉強をしていくことにした. 本当はもっと体系的にまとめた方がいいのかもしれない. しかしそうやっているといつまで経っても公開できないので, コンテンツというか記事をちょこちょこ出していくことにする.

野望: 中高の数学学習向けコンテンツも作りたい

実は私は現代数学観光ツアーという企画もやっている. これは現代数学の勉強で苦労している人向けに, 物理や工学との絡みを重視しつつ, 特に解析学の学部 4 年までの大きなストーリーを見てもらうことを目的にした無料の通信講座だ.

登録ページで量子力学とかバリバリに謳っているので, かなりその筋というか, けっこうバリバリの人だけが来ると思っていた. しかし実際に蓋を開けてみると中学や高校の数学を復習しているという大人の参加者も多かった.

現代数学観光ツアーはその主旨からして, 彼らのニーズには全く応えられない. ただせっかく集まった人を切り捨てるのも忍びない. というわけでプログラミングと絡めて中高の数学や物理を復習していく方向での コンテンツメイキングを模索している.

そしてこれがある程度まとまったら, 数学や物理に興味がある中高生にも展開していきたい. Kindle などの電子書籍でまとめた形にするのもいいが, やはり通信講座で少しずつ小分けにして配布していくスタイルも良さそう.

何にせよ, 中高の復習がしたい大人, そして知的好奇心に満ちているものの, その捌け口が見つからない中高生に向けてもコンテンツを展開していきたい. そうした目的で記事を書いていく.

GitHub で次のようなまとめページもあった. 今後, 私の勉強の参考にもしていく.

もしあなたが興味を持ってこれらのページの何かで勉強して面白いのがあったら, ぜひコメントで教えてほしい.

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

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

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

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

  • Python 3.5.1 :: Anaconda 4.1.0 (x86_64)

諸注意

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

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

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

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

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

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

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

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
%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 乗則をプロットした.

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
%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 つの放物線のグラフにしておこう.

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
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_trajectoryfrange も入れておかなければ動かない.

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

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

次回予告

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

次回をお楽しみに.

Jupyter操作の基礎, Python, IPython の基礎/Pythonで数学を

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

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

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

  • Python 3.5.1 :: Anaconda 4.1.0 (x86_64)

諸注意

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

GitHub 上の ipynb をダウンロードし, Jupyter を実行してローカルで実行しながら遊ぶのが楽だろう.

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

Jupyter 起動

ターミナルなりコマンドプロンプトなりを立ち上げて次のコマンドを実行すると ブラウザで Jupyter が立ち上がる.

1
2
$ cd some_directory
$ jupyter notebook

some_directory は適切なディレクトリ名を打つこと. ターミナルやコマンドプロンプトが嫌な方は, ここにある runjupyter.bat か runjupyter.command を適切なディレクトリに置いてダブルクリックすれば, そのディレクトリで Jupyter が起動する.

モード

エディタ vi をご存知の方にはお馴染みの機能だ. コマンドモードとエディットモードがある.

Esc でコマンドモードに戻れるので, とりあえず何かあったら Esc でコマンドモードに戻ろう. そのあと h でヘルプを呼び出すといい.

またコマンドモードのときにセルで m を押すと, そのセルは Markdown モードになる. つまり普通のテキストファイルのような扱いだ. Markdown の書式が使えるのが嬉しい.

さらに嬉しいのは MathJax 連携していて LaTeX も書けること.

\begin{align} \int_{\mathbb{R}^d} f(x) dx. \end{align}

数学・物理利用にはこれ以上ないほど助かる.

各モードのキーバインドはコマンドモード時に h を押せば ヘルプが出てくるのでそれを見るのがいい. 基本的には vi 的なコマンド体系だ. 大事なコマンドはまとめておこう.

どちらのモードでも使える
キーバインド 意味
Ctrl + Enter セルを実行
Shift + Enter セルを実行して下のセルを選択
Alt + Enter セルを実行して下に新しいセルを挿入
Ctrl + S notebook を保存
編集モードで使える
コマンド 意味
Esc コマンドモードに移行
Ctrl + Shift + - セルを分割
コマンドモードで使える
コマンド 意味
h ヘルプ
Enter 編集モードに移行
↑ or k 上のセルに移動
↓ or j 次のセルに移動
y / m セルのタイプを変更: y は cell, m は Markdown
a / b 今のセルの上/下にセルを追加
x / c / v 今のセルを cut/copy/paste
dd 今のセルを削除
z 最後の削除を取り消す (undo)
Shift + = 下のセルとマージ

Python のイントロダクション

ここからは GitHub のこのページから ipynb を落としてきて, 実際に Jupyter で実行しながら読むのがいい.

もしかしたら生の ipynb をそのまま落としたいという方もいるかもしれない. そういう場合はここから落とそう. もちろん git clone で落とした方がいいとは思うけれども, あなたが数学サイドから入ってきたなら Git には不慣れなこともあるだろう.

基本的に Jupyter 上で賄う予定なので, このミニ講座内で直接 IPython を実行することは少ないだろう. ただ実践の場で IPython でちょこちょことコード実行することもあるはずだ. Jupyter 利用と大きくは変わらないと思っているが, そういう場合の参考として IPython 用の注意もいくつか書いてある. そちらもぜひここを参照してほしい.

2 系で解説されているもののドットインストールを見て勉強するのもいいだろう.

1
2
 # まずは Hello, World! から
print("Hello, World!") # Jupyter 上で入力するとわかるが括弧やクオートが自動で両方入る

結果は次の通り. ちなみにこの結果は print の部分で出力した内容だ. 以下同じなのでそのように読んでほしい.

1
Hello, World!

数学関係の話をするのだからまずは四則演算を確認したい.

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
 # 計算もできる

 # かけ算
print(2*2)     # 演算子の前後は空けなくてもいい

 # 割り算
print(5 / 2)   # 演算子の前後は空けてもいい

 # 切り捨て割り算
print(5 // 3)

 # あまり
print(5 % 3)

 # a の b 乗
print(5 ** 3)

出力結果は次の通り.

1
2
3
4
5
4
2.5
1
2
125

あまり使うことはないと思うけれども, 文字列のエスケープも簡単に説明しておこう.

1
2
3
4
5
 # 文字列のエスケープ
print("Hello \"world\"")             #  ダブルクオートでくくると文字列の中のダブルクオートをエスケープする必要あり
print("A list:\n* item 1\n* item 2") # \n は改行
print("C:\\path\\on\\windows")       # バックスラッシュ自体のエスケープ
print(r"C:\path\on\windows")         # r"" は生リテラル (raw literal)

出力結果は次の通りだ.

1
2
3
4
5
6
Hello "world"
A list:
* item 1
* item 2
C:\path\on\windows
C:\path\on\windows
リスト

numpy, scipy を使うなら ndarray がメインになるのであまり使わなさそうな印象がある. しかし sympy だとまた違うかもしれない. 何にしろ Python としては基本的で大事なので紹介しておく. 数学での利用法としては数をたくさん叩き込むのに使える. Haskell と違って形式的にでも無限リストを作れない (はず) なので, 無限に続く数列は作れない.

sympy との絡みでどうなっていくか, それは今後 sympy を勉強していく中で見極めていきたいところ. numpy, scipy 利用時はなるべく ndarray で片付けるべき, というところだけ理解している. それ用のコードも順次公開していきたい.

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
items = [3, 4, 5, 7, 6]
print(items)
print(len(items)) # リストの長さ
print(items[0])   # リストの要素を取得: 0 始まり
print(items[-1])  # マイナスにすると後ろから数えて取得
print(items[0:4]) # 複数要素を一気に取得: 第2引数の扱いに注意
items.append(8)   # 要素追加:破壊的に追加
print(items)

print("\nrange の出力") # Python3 ではリストではなくイテレータで返る
print(range(10))
print(range(1, 11))
print(range(0, 30, 5))
print(range(0, 10, 3))
print(range(0, -10, -1))
print(range(0))
print(range(1, 0))

<!-- リストを操作する関数 -->
print(sum(items))
他のコンテナ

タプル, 辞書, 集合がある. 集合は数学利用もあるのでは? と思ったかもしれない. ただ数学用途にはいろいろと足りないことがある. numpy にも集合関係のメソッドはある. しかし私が注目しているのは sympy での集合関係のメソッドだ. こちらはある程度の数学利用に耐える内容という印象がある. もちろんそのうちきちんと紹介する.

まずタプル, 辞書, Python 標準の集合を見てみよう.

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
mytuple = (1, 2, 3)
print(mytuple)
print(mytuple[0])

 # 辞書: キーと値のペア
mydict = {'a': 1, 'b': 2, 'c': 3}
print('a:', mydict['a'])
 # キーのリストを取得
print(mydict.keys())
 # 値のリストを取得
print(mydict.values())

 # 集合: 数学の集合と同じで重複は排除
myset = set([1, 2, 3, 2, 1])
print(myset)

結果は次の通り.

1
2
3
4
5
6
(1, 2, 3)
1
a: 1
dict_keys(['c', 'a', 'b'])
dict_values([3, 1, 2])
{1, 2, 3}
ループ

いわゆる for だ. Python でも for がある. while もあるがここでは省略.

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
a = [5,6,7,8,9]

<!-- 要素だけ取る -->
for x in a:
    print(x)

print("\nenumerate の出力")
<!-- キーと値を両方取得 -->
for k, v in enumerate(a):
    print(k, v)

出力は次の通り.

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
5
6
7
8
9

enumerate の出力
0 5
1 6
2 7
3 8
4 9
リスト内包表記

すぐにわかる必要はない. 数学の集合の記法っぽく書けるのでそれがはまるところでは見やすくなる. for で同じことがやれるので, まずは for で書けるようになることを目指そう.

1
2
3
a = [5,6,7,8,9]
b = [x**2 for x in a if x % 2 == 0] # 偶数だけ取得してさらに 2 乗する
print(b)

結果は次の通り.

1
[36, 64]
条件分岐

これももちろん if 関係. Python では if, elif, else を使う. elif がかなり特殊な書き方なので注意しよう.

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
a = [1,2,3]
for x in a:
    if x % 3 == 1:
        print("3 で割って余り 1")
        print(x)
    elif x % 3 == 2: # elif という気持ち悪い記法
        print("3 で割って余り 2")
        print(x)
    else:
        print("3 で割り切れる")
        print(x)

結果は次の通り.

1
2
3
4
5
6
3 で割って余り 1
1
3 で割って余り 2
2
3 で割り切れる
3
関数

関数を定義できる. 数学の関数よりも広い意味で使われている. 数学の関数をきちんとイメージできるならとりあえずそれでいい. 実際そういう使い方もしていく. 専門家には怒られるだろうが, 上でも少し触れた「メソッド」もとりあえず関数だと思っておいてほしい.

1
2
3
4
5
6
7
8
def is_even(number):
    """Return whether an integer is even or not.""" # 3 つのクオートは複数行コメント
    return number % 2 == 0

print(is_even(1))
print(is_even(2))
print(is_even(3))
print(is_even(4))

結果は次の通り. True, False はいわゆる真偽値. 真偽値という概念があるのだ.

1
2
3
4
False
True
False
True
以下 IPython

IPython を単独では使わない. しかし Jupyter が裏で使っている. Jupyter を活用していく予定なので必然的に使うことになる.

そしてすごいことに Jupyter 上からコマンド実行できる. Unix コマンド前提なので以下のプログラムは素の Windows では実行できないようだ. Windows10 の Ubuntu の bash なり git bash なり msys2 なりを入れておけばおそらく実行可能(未検証).

パスが通っていることが前提だが, ターミナルで ipython と打てば起動する. 「パス」がわからない方は無視しよう. 動くならそれでいい.

1
$ ipython

Jupyter 上では以下の通り直接実行できる。

1
%pwd # 以下の結果は私の環境で実行した場合の結果: 人によって変わる

私の場合は次のような出力結果になった. どうなるかはあなたの実行環境次第だ.

1
'/Users/phasetr/codes/mathcodes/jupyter'

Unix コマンド実行結果.

1
%ls -l

結果は次の通り.

1
2
3
4
total 600
-rw-r--r--  1 phasetr  staff   35589  8 17 22:37 fundamentals.ipynb
-rw-r--r--  1 phasetr  staff    5987  8 14 20:15 math_simple_calculation.ipynb
-rw-r--r--  1 phasetr  staff  260296  8 17 22:38 math_simple_graph.ipynb
コマンド実行結果を Python の変数に入れる

細かく対話型シェルを使って調べているとき, コマンド実行結果を簡単に Python の変数に入れられるのはかなり便利. この辺はある程度プログラミングをやっているとわかる. この Jupyter 講座のメインは数学, 特にまずは中学・高校の数学と sympy を主軸に据えていくので まだあまりご利益を感じるような場面はない.

ただあなたが統計関係で使っていきたいと思うなら, ちょこちょこと IPython でコマンド実行する機会もあるだろうし. そういう場合には便利だろう.

1
2
files = !ls -1 -S | grep test # ! をつけるとコマンド実行してくれる
print(files)

結果は次の通り.

1
['test1', 'test2', 'test3']

よくあるやつだが _ で前の出力が取れる.

1
2
3
 # IPython では _ で前回の出力が取れる
%pwd
print(_*2)

出力結果は次の通り. 「/Users/phasetr/codes/mathcodes/jupyter」が 2 回くり返されていることに注意しよう.

1
/Users/phasetr/codes/mathcodes/jupyter/Users/phasetr/codes/mathcodes/jupyter
他の言語の実行

純粋な IPython 内では次のコードが実行できる: もちろん Haskell (GHCi) インストール済みの場合. IHaskell や IRuby もある.

1
2
#%%script ghci
#putStrLn "Hello, World"

出力結果は次の通り.

1
2
3
GHCi, version 8.0.1: http://www.haskell.org/ghc/  :? for help
Prelude> Hello, World
Prelude> Leaving GHCi.
Jupyter からファイルの読み書きもできる

IPython からやれるのは当然だ. これが Jupyter 上からもやれるのが割とすごい.

1
2
%%writefile myfile.txt
Hello world!

実行結果は次の通り. 既に myfile.txt があったのでこういう出力になっている.

1
Overwriting myfile.txt

ここからいくつか Unix コマンド実行が続く.

1
2
<!-- 前のセルからの続き -->
!more myfile.txt
1
2
Hello world!
[?1l>
1
2
3
<!-- 前のセルからの続き -->
!rm myfile.txt
!ls -l
1
2
3
4
5
6
rm: myfile.txt: No such file or directory
total 32
-rw-r--r--@ 1 phasetr  staff  16035  8 13 21:35 fundamentals.ipynb
-rw-r--r--@ 1 phasetr  staff      0  8 13 21:09 test1
-rw-r--r--@ 1 phasetr  staff      0  8 13 21:09 test2
-rw-r--r--@ 1 phasetr  staff      0  8 13 21:09 test3
1
!head -n5 t # ここで tab を押すと補完候補が出る: Jupyter 上でも出る。

iframe で YouTube の動画の呼出もできてしまう. Jupyter 上で表示してくれる. これが本当に素敵. 本当にいい notebook.

1
2
3
<!-- GUI 系のことだってできてしまう -->
from IPython.display import YouTubeVideo
YouTubeVideo('AyJI5XyMQ9M')

実行結果は次の通り.

1
2
3
4
5
6
7
    <iframe
        width="400"
        height="300"
        src="https://www.youtube.com/embed/AyJI5XyMQ9M"
        frameborder="0"
        allowfullscreen
    ></iframe>
インタラクティブなこともできる

グラフ描画のところでも紹介している。 セルを実行すると Jupyter 上で入力ボックスが出て入力待ちになり, その入力を受け取ってプログラム実行してくれる.

1
2
3
4
from ipywidgets import interact
@interact(x=(0, 10))
def square(x):
    print("The square of %d is %d." % (x, x**2))

実行結果は次の通り. ぜひ Jupyter を起動してやってみてほしい.

1
The square of 5 is 25.
IPython でのファイル読み込み

以下のセルは実行すると本当にファイルができるので注意すること.

1
2
3
4
<!-- 準備のためにディレクトリファイル作成 -->
#!mkdir subdir
!cd subdir; touch test1.txt test2.txt test3.txt test1.org test2.org
!cd subdir; ls -l

実行結果は次の通り.

1
2
3
4
5
6
7
total 0
-rw-r--r--  1 phasetr  staff  0  8 14 12:56 test1.org
-rw-r--r--@ 1 phasetr  staff  0  8 14 12:56 test1.txt
-rw-r--r--  1 phasetr  staff  0  8 14 12:56 test2.org
-rw-r--r--@ 1 phasetr  staff  0  8 14 12:56 test2.txt
-rw-r--r--  1 phasetr  staff  0  8 13 23:16 test3.org
-rw-r--r--@ 1 phasetr  staff  0  8 14 12:56 test3.txt

? でヘルプを出してくれる.

1
2
3
4
import networkx
<!-- 変数名の後に ? をつけるとヘルプ的なものが出る -->
networkx.Graph?
<!-- 他にも %pdef, %pdoc, %psource, %pfile などがある -->

どこまで意味があるかは微妙だが, Jupyter でコードを書いて, そのコードを Python のスクリプトに書き出せる.

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
%%writefile filelist.py
<!-- coding: utf-8 -->
import sys
import os
<!-- コマンドライン呼び出しの第 1 引数として検索するディレクトリを指定 -->
if len(sys.argv) > 1:
    directory = sys.argv[1]

<!-- 該当ディレクトリ内の全てのディレクトリをリストアップ -->
files = os.listdir(directory)
<!-- 拡張子を切り捨ててファイル名だけ取り出す -->
identifiers = [file.split('.')[0] for file in files]
<!-- set で重複を取り除いてソートする -->
names = sorted(set(identifiers))

出力結果は次の通り.

1
Overwriting filelist.py

Jupyter 上から作ったスクリプトを実行できる.

1
2
3
<!-- スクリプト実行 -->
%run filelist.py subdir
print(names)

そして実行結果.

1
['test1', 'test2', 'test3']

次回予定

最初は numpy, scipy を中心にしていこうと思っていた. 微分方程式をゴリゴリ解いてみたいと思っていたのだ.

しかし中高くらいの数学からやり直したい需要が大きいらしいことを知った. その目的にはむしろ Sympy が役に立つ. そこで Sympy と中高数学くらいのところを狙ってやっていこうと思っている. 次回も楽しみにしていてほしい.

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

中高数学のためのSympyのサンプル集/Pythonで数学を

諸注意

ここまでの記事はタグ Pythonでの数学物理 を参考にしてほしい. 特に Jupyter のインストールを前提にしている. 今回の分は GitHub のこのページを見てもらうと早い.

GitHub 上の ipynb をダウンロードし, Jupyter を実行してローカルで実行しながら遊ぶのが楽だろう.

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

はじまる前に

この Web 上の連続記事は現代数学観光ツアーで中高の数学をやってほしいという要望がいくつかあったので, それへの対応としてプログラミングとの連携も絡めた展開としてはじめている. その中で numpy, scipy 以上に sympy がよさそうなことに気付き, それで sympy の勉強をはじめた.

今回は sympy でどんなことができるかのイメージ作りのため, コードサンプルを載せていく回だ. これなら中高の数学との連携が取りやすいことをイメージしてもらえるだろう.

一応 Jupyter 利用を前提にしているので, まずは Jupyter 上での TeX 出力のためのコードを置いておく.

次のコードを読み込ませてほしい. ついでに Sympy で使うシンボルも定義している. シンボルといった言葉も追々説明していく. 要は方程式の未知変数みたいなものだと思ってほしい.

1
2
3
4
5
6
7
from sympy import *
from IPython.display import display
init_printing(use_unicode=True)

 * シンボル定義
x = Symbol('x')
y = Symbol('y')
サンプル
1
2
3
4
5
 # display で出力すると mathjax 連携してくれて素敵
display(x + y + x - 4 * y)

a = Integral(cos(x)*exp(x), x)
display(Eq(a, a.doit()))

結果は次の通り.

$$2 x - 3 y$$

$$\int e^{x} \cos{\left (x \right )}\, dx = \frac{e^{x}}{2} \sin{\left (x \right )} + \frac{e^{x}}{2} \cos{\left (x \right )}$$

参考サイト

何と言ってもまずは公式だろう.

次のサイトは sympy の結果をランダムで色々教えてくれる。 適当に眺めるのにすごくいい。

次のサイトにアクセスすれば sympy をインストールしなくてもオンラインで使える。

単純なサンプルその 2

今回はとにかくコードをずらずらと並べている. まずはこれでどんなことをやれるのか, やれそうか感じ取ってほしい. 細かい解説は今後の記事を楽しみにしていてほしい.

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
import sympy
display(sympy.sqrt(3))

 # 勝手に計算やってくれるので素敵
display(sympy.sqrt(8))

from sympy import symbols
x, y = symbols('x y')
expr = x + 2*y
display(expr)
display(expr + 1)
display(expr - x)
display(x*expr)

from sympy import expand, factor
expanded_expr = expand(x*expr)
display(expanded_expr)
display(factor(expanded_expr))

from sympy import *
x, t, z, nu = symbols('x t z nu')

print('計算: diff(sin(x)*exp(x), x)')
display(diff(sin(x)*exp(x), x))

print('計算: integrate(exp(x)*sin(x) + exp(x)*cos(x), x)')
display(integrate(exp(x)*sin(x) + exp(x)*cos(x), x))

print('計算: integrate(sin(x**2), (x, -oo, oo))')
display(integrate(sin(x**2), (x, -oo, oo)))

print('計算: limit(sin(x)/x, x, 0)')
display(limit(sin(x)/x, x, 0))

print('計算: solve(x**2 - 2, x)')
display(solve(x**2 - 2, x))

print('微分方程式の解: y" - y = e ** t')
f = Function('f')
display(dsolve(Eq(f(t).diff(t, t) - f(t), exp(t)), f(t)))

print('行列の固有値: Matrix([[1, 2], [2, 2]]).eigenvals()')
display(Matrix([[1, 2], [2, 2]]).eigenvals())

print('ベッセル関数を球ベッセル関数で書く')
print('besselj(nu, z).rewrite(jn)')
display(besselj(nu, z).rewrite(jn))

print('LaTeX にする')
print('latex(Integral(cos(x)**2, (x, 0, pi)))')
display(latex(Integral(cos(x)**2, (x, 0, pi))))

結果は次の通りだ.

$$\sqrt{3}$$

$$2 \sqrt{2}$$

$$x + 2 y$$

$$x + 2 y + 1$$

$$2 y$$

$$x \left(x + 2 y\right)$$

$$x^{2} + 2 x y$$

$$x \left(x + 2 y\right)$$

1
計算: diff(sin(x)*exp(x), x)

$$e^{x} \sin{\left (x \right )} + e^{x} \cos{\left (x \right )}$$

1
計算: integrate(exp(x)*sin(x) + exp(x)*cos(x), x)

$$e^{x} \sin{\left (x \right )}$$

1
計算: integrate(sin(x**2), (x, -oo, oo))

$$\frac{\sqrt{2} \sqrt{\pi}}{2}$$

1
計算: limit(sin(x)/x, x, 0)

$$1$$

1
計算: solve(x**2 - 2, x)

$$\left [ - \sqrt{2}, \quad \sqrt{2}\right ]$$

1
微分方程式の解: y" - y = e ** t

$$f{\left (t \right )} = C_{2} e^{- t} + \left(C_{1} + \frac{t}{2}\right) e^{t}$$

1
行列の固有値: Matrix([[1, 2], [2, 2]]).eigenvals()

$$\left { \frac{3}{2} + \frac{\sqrt{17}}{2} : 1, \quad - \frac{\sqrt{17}}{2} + \frac{3}{2} : 1\right }$$

1
2
ベッセル関数を球ベッセル関数で書く
besselj(nu, z).rewrite(jn)

$$\frac{\sqrt{2} \sqrt{z}}{\sqrt{\pi}} j_{\nu - \frac{1}{2}}\left(z\right)$$

1
2
3
4
LaTeX にする
latex(Integral(cos(x)**2, (x, 0, pi)))

'\\int_{0}^{\\pi} \\cos^{2}{\\left (x \\right )}\\, dx'

数値の計算

任意精度での計算ができる. 無限大もシンボリックに表現できるのも便利. あとで紹介するように極限も取れる.

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
display(pi**2)
display(pi.evalf())
display((pi + exp(1)).evalf())

print('sqrt(2) を 100 桁まで評価')
display(sqrt(2).evalf(100))

print('無限大')
display(oo)
display(oo > 99999)
display(oo + 1)

上記のコード片の実行結果は次の通り.

$$\pi^{2}$$

$$3.14159265358979$$

$$5.85987448204884$$

$$1.414213562373095048801688724209698078569671875376948073176679737990732478462107038850387534327641573$$

1
無限大

$$\infty$$

$$\mathrm{True}$$

$$\infty$$

展開

ここでは三角や複素数への展開だけだが, もちろん多項式の展開もできる. テイラー展開や漸近展開もある. これも今後の解説を楽しみにしていてほしい. 三角関数の展開があるというだけでもかなり可能性を感じてもらえると思っている.

1
2
3
4
display(expand(x + y, complex=True))
display(I*im(x) + I*im(y) + re(x) + re(y))
display(expand(cos(x + y), trig=True))
display(cos(x)*cos(y) - sin(x)*sin(y))

上記のコード片の実行結果は次の通り.

$$\Re{x} + \Re{y} + i \Im{x} + i \Im{y}$$

$$\Re{x} + \Re{y} + i \Im{x} + i \Im{y}$$

$$- \sin{\left (x \right )} \sin{\left (y \right )} + \cos{\left (x \right )} \cos{\left (y \right )}$$

$$- \sin{\left (x \right )} \sin{\left (y \right )} + \cos{\left (x \right )} \cos{\left (y \right )}$$

多項式

というわけでお待ちかねの多項式の処理だ. 展開だけではなく因数分解もある. 多項式の割り算や最大公約元などの処理も一通りある.

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
expr = (x + y)**5
display(expr)
print('多項式を展開してくれる: (x + y)**5')
display(expand(expr))

a = (x + x**2)/(x*sin(y)**2 + x*cos(y)**2)
print('式を単純化: (x + x**2)/(x*sin(y)**2 + x*cos(y)**2)')
display(simplify(a))

print('多項式の割り算: div(x**2 - 4 + x, x-2)')
display(div(x**2 - 4 + x, x-2))
print('最大公約数: gcd(2*x**2 + 6*x, 12*x)')
display(gcd(2*x**2 + 6*x, 12*x))
print('最小公倍数: lcm(2*x**2 + 6*x, 12*x)')
display(lcm(2*x**2 + 6*x, 12*x))
print('因数分解: factor(x**4/2 + 5*x**3/12 - x**2/3)')
display(factor(x**4/2 + 5*x**3/12 - x**2/3))
print('多変数の因数分解: factor(x**2 + 4*x*y + 4*y**2)')
display(factor(x**2 + 4*x*y + 4*y**2))
print('根を求める: solve(x**2 + 4*x*y + 4*y**2)')
display(solve(x**2 + 4*x*y + 4*y**2))
print('根を求める: solve(x**2 + 4*x*y + 4*y**2, y)')
display(solve(x**2 + 4*x*y + 4*y**2, y))
print('複雑な根を求める: solve(x**2 + 4*x + 181, x)')
display(solve(x**2 + 4*x + 181, x))
print('無理数解: solve(x**3 + 4*x + 181, x)')
display(solve(x**3 + 4*x + 181, x))
print('多変数多項式を解く: solve_poly_system([y**2 - x**3 + 1, y*x], x, y)')
display(solve_poly_system([y**2 - x**3 + 1, y*x], x, y))

ちょっと長いけれども上記のコード片の実行結果は次の通り.

$$\left(x + y\right)^{5}$$

1
多項式を展開してくれる: (x + y)**5

$$x^{5} + 5 x^{4} y + 10 x^{3} y^{2} + 10 x^{2} y^{3} + 5 x y^{4} + y^{5}$$

1
式を単純化: (x + x**2)/(x*sin(y)**2 + x*cos(y)**2)

$$x + 1$$

1
多項式の割り算: div(x**2 - 4 + x, x-2)

$$\left ( x + 3, \quad 2\right )$$

1
最大公約数: gcd(2*x**2 + 6*x, 12*x)

$$2 x$$

1
最小公倍数: lcm(2*x**2 + 6*x, 12*x)

$$12 x^{2} + 36 x$$

1
因数分解: factor(x**4/2 + 5*x**3/12 - x**2/3)

$$\frac{x^{2}}{12} \left(2 x - 1\right) \left(3 x + 4\right)$$

1
多変数の因数分解: factor(x**2 + 4*x*y + 4*y**2)

$$\left(x + 2 y\right)^{2}$$

1
根を求める: solve(x**2 + 4*x*y + 4*y**2)

$$\left [ \left { x : - 2 y\right }\right ]$$

1
根を求める: solve(x**2 + 4*x*y + 4*y**2, y)

$$\left [ - \frac{x}{2}\right ]$$

1
複雑な根を求める: solve(x**2 + 4*x + 181, x)

$$\left [ -2 - \sqrt{177} i, \quad -2 + \sqrt{177} i\right ]$$

1
無理数解: solve(x**3 + 4*x + 181, x)

$$\left [ \frac{4}{\left(- \frac{1}{2} - \frac{\sqrt{3} i}{2}\right) \sqrt[3]{\frac{4887}{2} + \frac{3 \sqrt{2654409}}{2}}} - \frac{1}{3} \left(- \frac{1}{2} - \frac{\sqrt{3} i}{2}\right) \sqrt[3]{\frac{4887}{2} + \frac{3 \sqrt{2654409}}{2}}, \quad - \frac{1}{3} \left(- \frac{1}{2} + \frac{\sqrt{3} i}{2}\right) \sqrt[3]{\frac{4887}{2} + \frac{3 \sqrt{2654409}}{2}} + \frac{4}{\left(- \frac{1}{2} + \frac{\sqrt{3} i}{2}\right) \sqrt[3]{\frac{4887}{2} + \frac{3 \sqrt{2654409}}{2}}}, \quad - \frac{1}{3} \sqrt[3]{\frac{4887}{2} + \frac{3 \sqrt{2654409}}{2}} + \frac{4}{\sqrt[3]{\frac{4887}{2} + \frac{3 \sqrt{2654409}}{2}}}\right ]$$

1
多変数多項式を解く: solve_poly_system([y**2 - x**3 + 1, y*x], x, y)

$$\left [ \left ( 0, \quad - i\right ), \quad \left ( 0, \quad i\right ), \quad \left ( 1, \quad 0\right ), \quad \left ( - \frac{1}{2} - \frac{\sqrt{3} i}{2}, \quad 0\right ), \quad \left ( - \frac{1}{2} + \frac{\sqrt{3} i}{2}, \quad 0\right )\right ]$$

関数の評価

関数のある点での値を求めることもできる. 代入 (substitution) をしていると思ってほしい.

1
2
3
4
f = x**2 + 3*x + 2
display(f)
f1 = f.subs([(x, 1)])
display(f1)

上記コード片の実行結果は次の通り.

$$x^{2} + 3 x + 2$$

$$6$$

方程式をある変数で解く

ある範囲の方程式をある範囲で解くことができる. いまのところ偏微分方程式の厳密解は 1 階線型までのようだ. 2 階の線型の厳密解表示があると便利ではある. もちろん境界条件の指定などあってかなり難しいとは思っている. Sympy 利用の有限要素法ライブラリもあるようだし, それを紹介していこうとも思っている. そこまで含めれば大学水準の話もかなりできそうで, 私自身勉強を進めていて楽しみでならない.

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
expr = cos(x) * ln(y) + 2/y
display(expr)
display(solve(expr, x))
display(solve(x**4 - 1, x))
display(solve([x + 5*y - 2, -3*x + 6*y - 15], [x, y]))

f = x**4 - 3*x**2 + 1
display(factor(f))
display(factor(f, modulus=5))

print('超越方程式も一部サポート')
display(solve(exp(x) + 1, x))

上記コード片の実行結果は次の通り

$$\log{\left (y \right )} \cos{\left (x \right )} + \frac{2}{y}$$

$$\left [ - \operatorname{acos}{\left (- \frac{2}{y \log{\left (y \right )}} \right )} + 2 \pi, \quad \operatorname{acos}{\left (- \frac{2}{y \log{\left (y \right )}} \right )}\right ]$$

$$\left [ -1, \quad 1, \quad - i, \quad i\right ]$$

$$\left { x : -3, \quad y : 1\right }$$

1
超越方程式も一部サポート

$$\left [ i \pi\right ]$$

$$\left(x^{2} - x - 1\right) \left(x^{2} + x - 1\right)$$

$$\left(x - 2\right)^{2} \left(x + 2\right)^{2}$$

論理式

よく論理的思考を学ぶには数学をやるのがいいと言われる. その有効性に疑問をもってはいるが, そういう勉強にも役立つライブラリではあるだろう.

1
2
display(satisfiable(x & y))
display(satisfiable(x & ~x))

上記コード片の実行結果は次の通り.

1
2
3
{y: True, x: True}

False

連立方程式を解く

連立方程式だって解けてしまう. 大学レベルから言えば行列で何とかしたいところだ. そしてもちろん行列系の演算もある. 楽しみにしていてほしい.

1
2
3
4
5
6
7
8
x, y, z = symbols('x y z')
eq1 = x * y * z + 234
eq2 = x + y + z - 20
eq3 = 5 * x - y + 2 * z - 85
display(eq1)
display(eq2)
display(eq3)
solve([eq1,eq2,eq3], [x,y,z])

上記コード片の実行結果は次の通り.

$$x y z + 234$$

$$x + y + z - 20$$

$$5 x - y + 2 z - 85$$

$$\left [ \left ( 13, \quad -2, \quad 9\right ), \quad \left ( - \frac{9 \sqrt{17}}{4} + \frac{39}{4}, \quad - \frac{9 \sqrt{17}}{4} - \frac{21}{4}, \quad \frac{31}{2} + \frac{9 \sqrt{17}}{2}\right ), \quad \left ( \frac{9 \sqrt{17}}{4} + \frac{39}{4}, \quad - \frac{21}{4} + \frac{9 \sqrt{17}}{4}, \quad - \frac{9 \sqrt{17}}{2} + \frac{31}{2}\right )\right ]$$

線型代数

というわけで線型代数関係.

1
2
3
4
5
6
7
8
9
from sympy import Matrix
display(Matrix([[1,0], [0,1]]))

print('Numpy と違ってシンボルを含む行列が書ける')
x = Symbol('x')
y = Symbol('y')
A = Matrix([[1,x], [y,1]])
display(A)
display(A**2)

上記コード片の実行結果は次の通り.

$$\left[\begin{matrix}1 & 0\0 & 1\end{matrix}\right]$$

1
Numpy と違ってシンボルを含む行列が書ける

$$\left[\begin{matrix}1 & x\y & 1\end{matrix}\right]$$

$$\left[\begin{matrix}x y + 1 & 2 x\2 y & x y + 1\end{matrix}\right]$$

微分積分

極限

高校水準の微分積分をやるのに便利だと思ったのはこれを見たからだ. もはやその筋 (計算機関係の人達) には当たり前のことなのだろう. しかし最初に見るとけっこう感動する. 「こんな便利なものがフリーで使えるのか」と. 中高生はもちろんのこと, 大人もガンガン使ってほしい.

1
2
3
4
5
display(limit(tan(x), x, pi/2))
display(limit(tan(x), x, pi/2, dir="-"))
display(limit(sin(x)/x, x, 0))
 # 無限大は oo で
display(limit(sin(x)/x, x, oo))

上記コード片の実行結果は次の通り.

$$-\infty$$

$$\infty$$

$$1$$

$$0$$

微分

微分の結果も厳密にシンボリックに計算できる. この辺が numpy, scipy よりも中高数学をやるのに優れた点.

1
2
3
4
5
display(diff((sin(x) * x**2) / (1 + tan(cot(x)))))
display(diff(cot(x*y), y))
display(diff(y(x)**2 - 5* sin(x), x))
f = x**2 / y + 2 * x - ln(y)
display(diff(f,x))

$$\frac{x^{2} \cos{\left (x \right )}}{\tan{\left (\cot{\left (x \right )} \right )} + 1} - \frac{x^{2} \sin{\left (x \right )}}{\left(\tan{\left (\cot{\left (x \right )} \right )} + 1\right)^{2}} \left(\tan^{2}{\left (\cot{\left (x \right )} \right )} + 1\right) \left(- \cot^{2}{\left (x \right )} - 1\right) + \frac{2 x \sin{\left (x \right )}}{\tan{\left (\cot{\left (x \right )} \right )} + 1}$$

$$x \left(- \cot^{2}{\left (x y \right )} - 1\right)$$

$$2 y{\left (x \right )} \frac{d}{d x} y{\left (x \right )} - 5 \cos{\left (x \right )}$$

$$\frac{2 x}{y} + 2$$

高階微分

高階微分も一気に計算できる.

1
2
3
4
print('第 3 引数で指定する')
display(diff(sin(2*x), x, 1))
display(diff(sin(2*x), x, 2))
display(diff(sin(2*x), x, 3))

上記コード片の実行結果は次の通り.

1
第 3 引数で指定する

$$2 \cos{\left (2 x \right )}$$

$$- 4 \sin{\left (2 x \right )}$$

$$- 8 \cos{\left (2 x \right )}$$

微分方程式

常微分方程式ならけっこういろいろ解ける. ある程度は非線型もいけるし超幾何関数も扱える. これも楽しみに待っていてほしい.

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
 # 未定義の関数を定義
f, g = symbols('f g', cls=Function)
display(f)
display(f(x))
display(f(x).diff(x, x) + f(x))

diffeq = f(x).diff(x) + f(x)
display(dsolve(diffeq, f(x)))

diffeq = f(x).diff(x, x) + f(x)
display(dsolve(diffeq, f(x)))

diffeq = f(x).diff(x, x, x) + f(x)
display(dsolve(diffeq, f(x)))

 # キーワード引数を使ってヒントを出すこともできる。
 # 方程式が可分離であることを知っていルナらキーワードとして hint='separable'
display(dsolve(sin(x)*cos(f(x)) + cos(x)*sin(f(x))*f(x).diff(x), f(x), hint='separable'))

display(dsolve(x * f(x).diff(x) + f(x) - f(x) ** 2))
display(dsolve(x * f(x).diff(x) + f(x) - f(x) ** 2, hint='Bernoulli'))

上記コード片の実行結果は次の通り. ちょっと長い.

1
f

$$f{\left (x \right )}$$

$$f{\left (x \right )} + \frac{d^{2}}{d x^{2}} f{\left (x \right )}$$

$$f{\left (x \right )} = C_{1} e^{- x}$$

$$f{\left (x \right )} = C_{1} \sin{\left (x \right )} + C_{2} \cos{\left (x \right )}$$

$$f{\left (x \right )} = C_{3} e^{- x} + \left(C_{1} \sin{\left (\frac{\sqrt{3} x}{2} \right )} + C_{2} \cos{\left (\frac{\sqrt{3} x}{2} \right )}\right) \sqrt{e^{x}}$$

$$\left [ f{\left (x \right )} = - \operatorname{asin}{\left (\sqrt{\frac{C_{1}}{\sin^{2}{\left (x \right )} - 1} + 1} \right )} + \pi, \quad f{\left (x \right )} = \operatorname{asin}{\left (\sqrt{\frac{C_{1}}{\sin^{2}{\left (x \right )} - 1} + 1} \right )} + \pi, \quad f{\left (x \right )} = - \operatorname{asin}{\left (\sqrt{\frac{C_{1}}{\sin^{2}{\left (x \right )} - 1} + 1} \right )}, \quad f{\left (x \right )} = \operatorname{asin}{\left (\sqrt{\frac{C_{1}}{\sin^{2}{\left (x \right )} - 1} + 1} \right )}\right ]$$

$$f{\left (x \right )} = - \frac{C_{1}}{- C_{1} + x}$$

$$f{\left (x \right )} = \frac{1}{x \left(C_{1} + \frac{1}{x}\right)}$$

テイラー展開

お待ちかねテイラーて展開.

1
2
display(sin(x).series(x))
display(series(sin(x), x, pi/2))

上記コード片の実行結果は次の通り.

$$x - \frac{x^{3}}{6} + \frac{x^{5}}{120} + \mathcal{O}\left(x^{6}\right)$$

$$1 - \frac{1}{2} \left(x - \frac{\pi}{2}\right)^{2} + \frac{1}{24} \left(x - \frac{\pi}{2}\right)^{4} + \mathcal{O}\left(\left(x - \frac{\pi}{2}\right)^{6}; x\rightarrow\frac{\pi}{2}\right)$$

積分

不定積分も定積分もある. 出せるときには厳密解が出るので, そこは numpy, scipy よりも便利. もちろん実際にはきちんとした使い分けが必要だ. 中高の数学をやるならやはり sympy に軍配が上がるだろう.

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
print('積分')
print('integrate(6*x**5, x)')
display(integrate(6*x**5, x))

display(integrate(sin(x), x))
display(integrate(log(x), x))
display(integrate(2*x + sinh(x), x))

expr = cos(x) * ln(y) + 2/y
display(integrate(expr, x))
display(integrate(tan(x)))
display(integrate(2*x + y, y))
display(integrate(2*x + y, (x, 1, 3)))
display(integrate(2*x + y, (x, 1, 3), (y, 2, 4)))

print('特殊関数もいける')
display(integrate(exp(-x**2)*erf(x), x))

print('定積分')
display(integrate(x**3, (x, -1, 1)))
display(integrate(sin(x), (x, 0, pi/2)))
display(integrate(cos(x), (x, -pi/2, pi/2)))

print('広義積分')
display(integrate(exp(-x), (x, 0, oo)))
display(integrate(exp(-x**2), (x, -oo, oo)))
display(integrate(tan(x), (x, 0, pi/2)))

print('積分の厳密解')
display(integrate(1/(x**2 + 1), (x, 0, oo)))
display(integrate(exp(x) / (1 + exp(2 * x))))
display(integrate((2 * x+3)**7))

上記コード片の実行結果は次の通り.

1
2
積分
integrate(6*x**5, x)

$$x^{6}$$

$$- \cos{\left (x \right )}$$

$$x \log{\left (x \right )} - x$$

$$x^{2} + \cosh{\left (x \right )}$$

$$\frac{2 x}{y} + \log{\left (y \right )} \sin{\left (x \right )}$$

$$- \frac{1}{2} \log{\left (\sin^{2}{\left (x \right )} - 1 \right )}$$

$$2 x y + \frac{y^{2}}{2}$$

$$2 y + 8$$

$$28$$

1
特殊関数もいける

$$\frac{\sqrt{\pi}}{4} \operatorname{erf}^{2}{\left (x \right )}$$

1
定積分

$$0$$

$$1$$

$$2$$

1
広義積分

$$1$$

$$\sqrt{\pi}$$

$$\infty + \frac{i \pi}{2}$$

1
積分の厳密解

$$\frac{\pi}{2}$$

$$\operatorname{RootSum} {\left(4 z^{2} + 1, \left( i \mapsto i \log{\left (2 i + e^{x} \right )} \right)\right)}$$

$$16 x^{8} + 192 x^{7} + 1008 x^{6} + 3024 x^{5} + 5670 x^{4} + 6804 x^{3} + 5103 x^{2} + 2187 x$$

まとめ

まずは Sympy の射程距離や威力を知ってもらおうと思い, サンプル集という形で出した. これだけでも大人の復習というところでは遊び所がわかる人も多いだろう. あなたは十分に可能性を感じてくれたかもしれない.

もちろんまだ「これをどう使っていけばいいの?」と思っている方もいるだろう. あなたもそうかもしれない. まずは私自身の勉強もあるので, 具体的な中高数学の復習などに入る前にライブラリの様子を細かく知る部分で記事を書いていく. ぜひ今後も記事を追いかけてきてほしい.

今回の分は GitHub のこのページに ipython notebook としてまとまっているので, ぜひそちらを Jupyter で実行しながら遊んでほしい.

Python で数学をするシリーズはタグにまとめている. ぜひタグ Pythonでの数学物理を定期的に巡回して更新チェックや面白そうなところのつまみぐいをしてほしい. かなり長いシリーズになりそうなので, まずは面白そうなところをつまみぐいして, そこから徐々に体系的なところに入っていくのがいいだろう. なるべく最初の負荷は楽にするべきだ. 楽しく勉強を続けるために工夫をしてみてほしい.

中高数学のためのSympy: 諸々のコードサンプル

ここまでの記事はタグ Pythonでの数学物理 を参考にしてほしい. 特に Jupyter のインストールを前提にしている. 今回の分は Sympy に関する諸々のコードサンプルを出す.

本当は 1 つ 1 つ出した方がいいのだろうけれども, ちょっと時間が取れないのでコードを一気に出す. GitHub のここの sympy_ 関連のファイルを見てほしい.

何に関するコードかだけでも書いておこう.

0詰め, パディング

  • http://www.lifewithpython.com/2015/10/python-zero-padding.html
1
2
3
4
5
for i in range(32):
    padded = str('{0:02d}'.format(i))
    tmpdate = "2017-08-" + padded
    print(tmpdate)
    print(tmpdate)
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
... ... ... ... 2017-08-00
2017-08-00
2017-08-01
2017-08-01
2017-08-02
2017-08-02
2017-08-03
2017-08-03
2017-08-04
2017-08-04
2017-08-05
2017-08-05
2017-08-06
2017-08-06
2017-08-07
2017-08-07
2017-08-08
2017-08-08
2017-08-09
2017-08-09
2017-08-10
2017-08-10
2017-08-11
2017-08-11
2017-08-12
2017-08-12
2017-08-13
2017-08-13
2017-08-14
2017-08-14
2017-08-15
2017-08-15
2017-08-16
2017-08-16
2017-08-17
2017-08-17
2017-08-18
2017-08-18
2017-08-19
2017-08-19
2017-08-20
2017-08-20
2017-08-21
2017-08-21
2017-08-22
2017-08-22
2017-08-23
2017-08-23
2017-08-24
2017-08-24
2017-08-25
2017-08-25
2017-08-26
2017-08-26
2017-08-27
2017-08-27
2017-08-28
2017-08-28
2017-08-29
2017-08-29
2017-08-30
2017-08-30
2017-08-31
2017-08-31

datetime, 日付

日付の宣言

現在の日付

1
2
3
4
import datetime

print(datetime.date.today()) # 2021-05-15
print(datetime.datetime.today()) # 2021-05-15 23:58:55.230456
  • datetime.date: 年月日だけ
  • datetime.datetime: 時分秒ミリ秒まで持つ

任意の日付

1
2
3
4
import datetime

print(datetime.date(2017, 11, 12)) # 2021-05-15
print(datetime.datetime(2017, 11, 12, 9, 55, 28)) # 2021-05-15 09:55:28
  • datetime.date: 引数に年月日を指定できる
  • datetime.datetime: 引数に年月日, 時分秒ミリ秒, タイムゾーンを指定できる

年月日・時分秒

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
import datetime

now = datetime.datetime(2021, 05, 15, 10, 0, 0)
print(now.year) # 2021
print(now.month) # 5
print(now.day) # 15
print(now.hour) # 10
print(now.minute) # 00
print(now.second) # 00

#### 曜日
`weekday()`を使う.
値は次の通り.

| 数値 | 意味   |
|------|--------|
|    0 | 月曜日 |
|    1 | 火曜日 |
|    2 | 水曜日 |
|    3 | 木曜日 |
|    4 | 金曜日 |
|    5 | 土曜日 |
|    6 | 日曜日 |

```python
import datetime

d = datetime.date(2021, 5, 15)
print(d.weekday()) # 5

d = d + datetime.timedelta(days=1)
print(d.weekday()) # 6

d = d + datetime.timedelta(days=1)
print(d.weekday()) # 0

日付だけ取る

1
2
3
4
import datetime

now = datetime.datetime(2021, 5, 15, 10, 0, 0)
print(now.date()) # 2021-05-15

時間だけ

1
2
3
4
import datetime

now = datetime.datetime(2021, 5, 15, 10, 0, 0)
print(now.time()) # 10:00:00

日付計算

1
2
3
4
5
6
7
8
import datetime

d = datetime.datetime(2021, 5, 15, 10, 0, 0)
print(d + datetime.timedelta(weeks=1))   # 2021-05-22 10:00:00
print(d + datetime.timedelta(days=1))    # 2021-05-16 10:00:00
print(d + datetime.timedelta(hours=1))   # 2021-05-15 11:00:00
print(d + datetime.timedelta(minutes=1)) # 2021-05-15 10:01:00
print(d + datetime.timedelta(seconds=1)) # 2021-05-15 10:00:01
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
import datetime
from dateutil.relativedelta import relativedelta

d = datetime.datetime(2021, 5, 15, 10, 0, 0)
print(d + relativedelta(years=1))   # 2022-05-15 10:00:00
print(d + relativedelta(months=1))  # 2021-06-15 10:00:00
print(d + relativedelta(weeks=1))   # 2021-05-22 10:00:00
print(d + relativedelta(days=1))    # 2021-05-16 10:00:00
print(d + relativedelta(hours=1))   # 2021-05-15 11:00:00
print(d + relativedelta(minutes=1)) # 2021-05-15 10:01:00
print(d + relativedelta(seconds=1)) # 2021-05-15 10:00:01

日付の比較

datetimedate比較できない.

1
2
3
4
5
6
7
8
9
import datetime

d1 = datetime.datetime(2021, 5, 15, 10, 0, 0)
d2 = datetime.datetime(2021, 5, 15, 10, 0, 1)
print((d1 > d2))  # False
print((d1 >= d2)) # False
print((d1 < d2))  # True
print((d1 <= d2)) # True
print((d1 == d2)) # False

日付変換

文字列から日付型

datetime.datetime.strptime を使います. 第一引数は文字列の日付で, 第二引数は日付フォーマットです.

1
2
3
4
5
6
7
8
9
import datetime

d_str = "2021/5/15"
d = datetime.datetime.strptime(d_str, "%Y/%m/%d")
print(d) # 2021-05-15 00:00:00

dt_str = "2021/5/15 10:00:00"
dt = datetime.datetime.strptime(dt_str, "%Y/%m/%d %H:%M:%S")
print(dt) # 2021-05-15 10:00:00

日付型から文字列

1
2
3
4
5
import datetime

df = datetime.datetime(2021, 5, 15, 10, 0, 0)
print(df.strftime("%Y/%m/%d")) # 2021/05/15
print(df.strftime("%Y/%m/%d %H:%M:%S")) # 2021/05/15 10:00:00

フォーマットで使う要素のメモ.

書式 意味
%Y 4 桁西暦 の 10 進表記
%m 0 埋めした 10 進数表記の月
%d 0 埋めした 10 進数表記の日した月中の日にち。
%H 0 埋めした 10 進数表記の 24 時表記時
%I 0 埋めした 10 進数表記の 12 時間表記時
%M 0 埋めした 10 進数表記の分
%S 0 埋めした 10 進数表記の秒
%f 0 埋めした 10進数表記のマイクロ秒
%b ロケール月名の短縮形
%B ロケール月名

dict, 辞書

要素を単純削除: del

すべての要素を削除するclear()メソッド, 個別の要素を削除するpop(), popitem()メソッド, そしてdelがある.

ここではdel文のサンプルを書く. 何となくこれを調べる契機になったyaml読み書きを前提にしたサンプルにしてある. もちろん純粋に辞書だけの問題ならimport yamlはいらない.

1
2
3
4
5
6
7
8
import yaml

input_yaml_file_name = ymlname
output_yaml_file_name = "mkdocs.yml"
orig_data = {}
with open (input_yaml_file_name, "r", encoding="utf-8", newline="") as f:
    orig_data = yaml.load(f, Loader=yaml.FullLoader)
    del orig_data["plugins"][0]["search"]["prebuild_"]

ループのためのメソッド

d.keys(), d.values(), d.items()でキー・値・キーと値のループが作れる.

exception, 例外

1
2
3
4
5
try:
    pass # some process
except Exception as e:
    print(e)
    pass # some process

file/directory

pathlib

ディレクトリの区切り文字をスラッシュにする

1
2
from pathlib import Path
p = Path(fname).as_posix()

ディレクトリ・ファイルを取る

1
2
from pathlib import Path
a = Path("/path/to").glob("**/*.org")

yaml

基本

ここでの記述は PyYAML を前提にしているのでインストールしよう.

1
pip install pyyaml

ファイルの読み込み・書き込み

文字コード指定を忘れるとひどいことになるので encoding='utf-88' をきちんとつけよう.

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
import yaml

input_yaml_file_name = ymlname
output_yaml_file_name = "mkdocs.yml"
orig_data = {}
with open (input_yaml_file_name, "r", encoding="utf-8", newline="") as f:
    orig_data = yaml.load(f, Loader=yaml.FullLoader)
    del orig_data["plugins"][0]["search"]["prebuild_"]
with open(output_yaml_file_name, "w", encoding="utf-8", newline="") as f:
    yaml.dump(orig_data, f, encoding='utf-8', allow_unicode=True, default_flow_style=False)

空のファイルを作る: touch

1
2
3
import pathlib
path = "path/to/1.txt"
pathlib.Path(path).touch()

ディレクトリ生成

os.makedir() より os.makedirs() が便利です. 特に exist_ok=True とすれば対象ディレクトリがあってもエラーになりません.

1
2
import os
os.makedirs(path, exist_ok=True)

ディレクトリの階層番号を取る

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
def get_hnum(p):
    """ディレクトリの階層番号を取る"""
    if os.path.isdir(p):
        # ディレクトリならカウントする
        pstr = str(p.as_posix())
        h_num = len(pstr.split("/"))
    else:
        # ファイルは 1 にする
        h_num = 1
    return h_num - 1

ディレクトリかどうかを判定

1
2
import os
os.path.isdir(fname)

ファイル移動

1
2
import shutil
shutil.move(from, to)

ファイルかどうかを判定

1
2
import os
os.path.isfile(fname)

ファイルコピー

1
2
import shutil
shutil.copyfile(from, to)

ファイル削除

1
2
3
import os
fname = "some file name"
os.remove(fname)

ファイル名取得

1
2
3
import os
filepath = 'dir/subdir/filename.ext'
dirname, basename = os.path.split(filepath)
1
2
from pathlib import Path
Path('path/to').name

ファイル読み書き

"r""w" に変えるとファイルに書き込めます.

1
2
3
4
5
6
7
#読むとき
with open(fname, "r", encoding="utf-8", newline="") as f:
    pass # 何かの処理

#書くとき
with open(fname, "w", encoding="utf-8", newline="") as f:
    pass # 何かの処理

Jupyter jupytext pyipynbに変換する

1
2
3
pip install jupytext
jupytext --to notebook --execute sample.py
jupytext --to notebook --execute arithmetic001.py
  • 変換できないとき: utf-8-unixでBOMなしにしないとエラーになる模様.

Jupyter nbconvert ipynbpyに変換

1
jupyter nbconvert --to script 00-introduction_01_install.ipynb --no-prompt

Jupyter nbconvert ipynbpyに変換時にin []を消す

  • URL
  • --no-promptを追加する
1
jupyter nbconvert --to script 00-introduction_01_install.ipynb --no-prompt

Jupyter nbconvertのテンプレート取得

  • 公式のGitHubshare/templatesから取ってくればよい
  • TODO: 落としてきたtemplateのディレクトリをどこに置くとエラーが出なくなる?

Jupyter Infoのログを出力しない

1
2
3
import logging
logger = logging.getLogger()
logger.setLevel(logging.CRITICAL)

Jupyter 「カーネルが見つからない」と言われたとき

念のため次の二つを確認する.

1
2
import sys
sys.executable
1
jupyter kernelspec list

後者が指定しているPython実行パスが正しいかが問題. WindowsではAppData/Roaming/jupyter/kernelsに対象ディレクトリがある. カーネルごとのディレクトリ内にkernel.jsonがあるから, その中のパスを調べる. 間違っていたら適切な値に修正する. Pythonをバージョンアップするとこの指定がずれていることがあるので, バージョンアップしたあとは特に注意して確認すること.

Jupyter カーネルの確認

1
jupyter kernelspec list

Jupyter カーネルの削除

1
jupyter kernelspec remove kernel_name

Jupyter 関連するツールのバージョン把握

1
jupyter --version

list, リスト

a-zのアルファベットのリストを作りたい

小文字だけ

1
2
3
4
5
letters = [chr(i) for i in range(97, 123)] # 小文字だけ
capitals = [chr(i) for i in range(65, 91)] # 大文字だけ

import string
all = [i for i in string.ascii_letters] # a-Z

リストを逆順にしたい

大文字・小文字を無視してソートしたい: key=str.lower

公式のソート HOW TOに記述があります.

1
2
>>> sorted("This is a test string from Andrew".split(), key=str.lower)
['a', 'Andrew', 'from', 'is', 'string', 'test', 'This']

リストを逆順にしたい

リストを破壊的に変更するreverse()メソッド, 順番を変えるスライスlst[::-1], イテレーターを作るreversed(lst)があります.

reverse()メソッドの返り値はNoneなので, reverse()したあとのオブジェクトを取る必要があります.

1
2
3
4
5
6
a = [1, 2, 3, 4, 5]
print(a[::-1])
print(list(reversed(a)))

a.reverse()
print(a)

リストを文字列にする

join()メソッドを使います.

1
[1, 2, 3].join(",")

リストを連結したい

+を使います.

Matplotlibでquiver

コードサンプル 1

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt

n = 8
X, Y = np.mgrid[0:n, 0:n]
T = np.arctan2(Y - n / 2., X - n/2.)
R = 10 + np.sqrt((Y - n / 2.0) ** 2 + (X - n / 2.0) ** 2)
U, V = R * np.cos(T), R * np.sin(T)

plt.axes([0.025, 0.025, 0.95, 0.95])
plt.quiver(X, Y, U, V, R, alpha=.5)
plt.quiver(X, Y, U, V, edgecolor='k', facecolor='None', linewidth=.5)

plt.xlim(-1, n)
plt.xticks(())
plt.ylim(-1, n)
plt.yticks(())

plt.show()

コードサンプル 2

図 1

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
import matplotlib.pyplot as plt
import numpy as np
from numpy import ma

X, Y = np.meshgrid(np.arange(0, 2 * np.pi, .2), np.arange(0, 2 * np.pi, .2))
U = np.cos(X)
V = np.sin(Y)


plt.figure()
Q = plt.quiver(U, V)
qk = plt.quiverkey(Q, 0.5, 0.92, 2, r'$2 \frac{m}{s}$', labelpos='W',
                   fontproperties={'weight': 'bold'})
l, r, b, t = plt.axis()
dx, dy = r - l, t - b
plt.axis([l - 0.05*dx, r + 0.05*dx, b - 0.05*dy, t + 0.05*dy])

plt.title('picture 1')

図 2

1
2
3
4
5
6
7
8
plt.figure()
Q = plt.quiver(X, Y, U, V, units='width')
qk = plt.quiverkey(Q, 0.9, 0.95, 2, r'$2 \frac{m}{s}$',
                   labelpos='E',
                   coordinates='figure',
                   fontproperties={'weight': 'bold'})
plt.axis([-1, 7, -1, 7])
plt.title('picture 2')

図 3

1
2
3
4
5
6
7
8
plt.figure()
Q = plt.quiver(X[::3, ::3], Y[::3, ::3], U[::3, ::3], V[::3, ::3],
               pivot='mid', color='r', units='inches')
qk = plt.quiverkey(Q, 0.5, 0.03, 1, r'$1 \frac{m}{s}$',
                   fontproperties={'weight': 'bold'})
plt.plot(X[::3, ::3], Y[::3, ::3], 'k.')
plt.axis([-1, 7, -1, 7])
plt.title("picture 3")

図 4

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
plt.figure()
M = np.hypot(U, V)
Q = plt.quiver(X, Y, U, V, M,
               units='x',
               pivot='tip',
               width=0.022,
               scale=1 / 0.15)
qk = plt.quiverkey(Q, 0.9, 1.05, 1, r'$1 \frac{m}{s}$',
                   labelpos='E',
                   fontproperties={'weight': 'bold'})
plt.plot(X, Y, 'k.')
plt.axis([-1, 7, -1, 7])
plt.title("picture 4")

図 5

1
2
3
4
5
6
7
8
plt.figure()
Q = plt.quiver(X[::3, ::3], Y[::3, ::3], U[::3, ::3], V[::3, ::3],
               color='r', units='x',
               linewidths=(2,), edgecolors=('k'), headaxislength=5)
qk = plt.quiverkey(Q, 0.5, 0.03, 1, r'$1 \frac{m}{s}$',
                   fontproperties={'weight': 'bold'})
plt.axis([-1, 7, -1, 7])
plt.title("picture 5")

図 6

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
plt.figure()
M = np.zeros(U.shape, dtype='bool')
M[U.shape[0]//3:2*U.shape[0]//3,
  U.shape[1]//3:2*U.shape[1]//3] = True
U = ma.masked_array(U, mask=M)
V = ma.masked_array(V, mask=M)
Q = plt.quiver(U, V)
qk = plt.quiverkey(Q, 0.5, 0.92, 2, r'$2 \frac{m}{s}$', labelpos='W',
                   fontproperties={'weight': 'bold'})
l, r, b, t = plt.axis()
dx, dy = r - l, t - b
plt.axis([l - 0.05 * dx, r + 0.05 * dx, b - 0.05 * dy, t + 0.05 * dy])
plt.title('picture 6')

plt.show()

Matplotlibでグラフの例

参考情報

具体的にプロットしたいイメージがあるなら公式のギャラリーを見に行くのがいいようだ. IDLE など interactive shell で対話的にやる場合は pylab の方が一気にインポートできて楽. スクリプトに書くなら matplotlib.pyplot (as plt) の方がいい.

org-babel 上だと軸が消えるようだができる図を直接見るときちんと軸が表示されるようだ. Jupyter (または通常のプログラム実行) で見るなら問題はない.

単純なグラフ

1
2
3
4
5
6
7
8
9
%matplotlib inline
import matplotlib.pyplot as plt

x_numbers = [2, 4, 6]
y_numbers = [3, 4, 1]

plt.plot(x_numbers, y_numbers)
plt.savefig("tmp/graph_example001_tmp.png")
plt.show()

プロットしたところを強調するマーカーをつけてみる

1
2
3
4
5
6
7
8
%matplotlib inline
import matplotlib.pyplot as plt

x_numbers = [2, 4, 5, 6]
y_numbers = [3, 4, 2, 1]

plt.plot(x_numbers, y_numbers, marker='o') # marker は 'o', '*', 'x', +' から選べる
plt.show()

marker を省くとプロットだけになる

点だけなのでものすごく見づらいがきちんと描画されている.

1
2
3
4
5
6
7
8
%matplotlib inline
import matplotlib.pyplot as plt

x_numbers = [2, 4, 5, 6]
y_numbers = [3, 4, 2, 1]

plt.plot(x_numbers, y_numbers, 'o') # marker は 'o', '*', 'x', +' から選べる
plt.show()

直線を引いてみる

1
2
3
4
5
6
7
8
%matplotlib inline
import matplotlib.pyplot as plt

x_numbers = [1, 2, 3]
y_numbers = [2, 4, 6]

plt.plot(x_numbers, y_numbers) # marker は 'o', '*', 'x', +' から選べる
plt.show()

適当な折れ線を引く

1
2
3
4
5
%matplotlib inline
import matplotlib.pyplot as plt

lst = [53.9, 56.3, 56.4, 53.4, 54.5, 55.8, 56.8, 55.0, 55.3, 54.0, 56.7, 56.4, 57.3]
plt.plot(lst, marker='o')

横軸を適当に設定

org-babel だと横軸が見えないので微妙. ここを見てみよう.

1
2
3
4
5
6
%matplotlib inline
import matplotlib.pyplot as plt

num = 18
h_axis = range(num, num + len(lst))
plt.plot(h_axis, lst, marker='o')

複数のグラフを同時に載せる

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
%matplotlib inline
import matplotlib.pyplot as plt

h_axis_multiple = range(1, 13)
lst1 = [31.3, 37.3, 47.2, 51.0, 63.5, 71.3, 72.3, 72.7, 66.0, 57.0, 45.3, 31.1]
lst2 = [40.9, 35.7, 43.1, 55.7, 63.1, 71.0, 77.9, 75.8, 66.6, 56.2, 51.9, 43.6]
lst3 = [37.3, 40.9, 50.9, 54.8, 65.1, 71.0, 78.8, 76.7, 68.8, 58.0, 43.9, 41.5]

plt.plot(h_axis_multiple, lst1, h_axis_multiple, lst2, h_axis_multiple, lst3, marker='o')


plt.legend(['graph1', 'graph2', 'graph3'])

グラフに色々追加

$x$ 軸, $y$ 軸に軸の名前がつくのだが org-babel で作った図だと見えない模様. ここを見るといい.

1
2
3
4
5
6
7
8
9
%matplotlib inline
import matplotlib.pyplot as plt

plt.plot(h_axis_multiple, lst1, h_axis_multiple, lst2, h_axis_multiple, lst3)

plt.title('graph title')
plt.xlabel('x axis')
plt.ylabel('y axis')
plt.legend(['graph1', 'graph2', 'graph3'])

軸の表示を変える 1

1
2
3
4
5
6
7
8
9
%matplotlib inline
import matplotlib.pyplot as plt

lst = [53.9, 56.3, 56.4, 53.4, 54.5, 55.8, 56.8, 55.0, 55.3, 54.0, 56.7, 56.4, 57.3]
plt.plot(lst, marker='o')
plt.axis(xmin=-1)
plt.axis(xmax=13)
plt.axis(ymin=50)
plt.axis(ymax=60)

軸の表示を変える 2

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
%matplotlib inline
import matplotlib.pyplot as plt


time_of_day = ['4 AM', '7 AM', '10 AM', '1 PM', '4 PM', '7PM', '10 PM']
forecast_temp = [10, 6, 15, 17, 20, 19, 16]
time_interval = range(1, len(time_of_day) + 1)

plt.plot(time_interval, forecast_temp, 'o-')

plt.xticks(time_interval, time_of_day)
plt.show()

縦横比を合わせる

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt

def draw_quadrant():
    """四分円の円周を描画"""
    circle_x = np.arange(0,1,0.001)
    circle_y = np.sqrt(1- circle_x * circle_x)
    plt.plot(circle_x, circle_y)
    plt.gca().set_aspect('equal', adjustable='box')
    plt.show()

draw_quadrant()

ファイルとして保存

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
%matplotlib inline
import matplotlib.pyplot as plt
x_numbers = [1, 2, 3]
y_numbers = [2, 4, 6]
plt.plot(x_numbers, y_numbers)




plt.savefig('tmp/graph_example012_tmp.png')
plt.savefig('tmp/graph_example012_tmp.svg')
plt.savefig('tmp/graph_example012_tmp.pdf')

円を描く

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
%matplotlib inline
import matplotlib.pyplot as plt

def create_circle():
    circle = plt.Circle((0, 0), radius = 0.5)
    return circle

def show_shape(patch):
    ax = plt.gca()
    ax.add_patch(patch)
    plt.axis('scaled')
    plt.show()

c = create_circle()
show_shape(c)

式をもとにグラフを描く: ニュートンの逆 2 乗則

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
%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')

plt.xlabel('Distance')
plt.ylabel('Newton gravitational force')
plt.title('A graph for the inverse square law')
plt.show()

式をもとにグラフを描く: ファンデルワールスの状態方程式

\begin{align} (P + \frac{a}{V^2})(V-b) = RT. \end{align}

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
%matplotlib inline
import matplotlib
import numpy as np
import matplotlib.pyplot as plt

matplotlib.rcParams['xtick.direction'] = 'out'
matplotlib.rcParams['ytick.direction'] = 'out'


a = 3.592
b = 0.0427
R = 0.08206


pointsnum = 1000



v = np.linspace(0.064, 0.18, pointsnum)
p = np.linspace(0, 100.0, pointsnum)
X, Y = np.meshgrid(v, p)


Z = (Y + (a / X ** 2)) * (X - b) / R


plt.figure()
CS = plt.contour(X, Y, Z, 25, linewidth=1.0, colors='k')
zc = CS.collections[0:8]
plt.setp(zc, linewidth=0.2, linestyle='dashed')

plt.clabel(CS, inline=1, fontsize=10, fmt='%4.0f')
plt.title('T (Kelvin)')
plt.xlabel('V (litres)')
plt.ylabel('P (atm)')
plt.show()

物体の投射 (放物線)

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
%matplotlib inline
import matplotlib.pyplot as plt
import math

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)

v0 = 20    # m/s
theta = 30 # 度数法で入力
tragectory = get_trajectory(v0, theta)

plt.plot(tragectory[0], tragectory[1])
plt.xlabel('x')
plt.ylabel('y')
plt.title('Projectile motion of a ball')
plt.show()

複数の放物線を重ねる

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
v0_list = [20, 40, 60]
theta = 45
tragectories = []

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

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()

フィボナッチ数列の隣項比が黄金比に近づくことを数値的に検証

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
%matplotlib inline
import matplotlib.pyplot as plt


def make_fibonacci(n):
    # n = 1, 2 のとき
    if n == 1:
        return [1]
    if n == 2:
        return [1, 1]

    # n > 2 の場合
    a = 1
    b = 1

    # 数列の最初の 2 要素
    series = [a, b]
    for i in range(n):
        c = a + b
        series.append(c)
        # 再計算をせずに済むように次の計算に備えて前 2 項を待機: 時間節約
        a = b
        b = c

    return series

def plot_ratio(series):
    ratios = []
    for i in range(len(series)-1):
        ratios.append(series[i + 1] / series[i])
    plt.plot(ratios)
    plt.title('Ratio between Fibonacci numbers & Golden ratio')
    plt.ylabel('Ratio')
    plt.xlabel('No.')
    plt.show()


num = 100

series = make_fibonacci(num)

plot_ratio(series)

Numpy も使って三角関数

1
2
3
4
5
6
7
8
%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt

x = np.arange(-3, 3, 0.1)
y = np.sin(x)
plt.plot(x, y)
plt.show()

Numpy も使って散布図

1
2
3
4
5
6
7
8
%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
x = np.random.randn(30)
y = np.sin(x) + np.random.randn(30)

plt.plot(x, y, "o")  # "o"は小さい円(circle marker)
plt.show()

Numpy も使って色を変更

1
2
3
4
5
6
7
8
%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt

x = np.random.randn(30)
y = np.random.randn(30)
plt.plot(x, y, "ro")  # "r"はredの省略
plt.show()

Numpy も使って正規分布のヒストグラム

1
2
3
4
5
6
%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt

plt.hist(np.random.randn(1000))
plt.show()

Numpy も使って 100x100 行列 1

1
2
3
4
5
6
%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt

plt.imshow(np.random.randn(100, 100))
plt.show()

Numpy も使って 100x100 行列 2

1
2
3
4
5
6
%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt

plt.imshow(np.random.randn(100, 100), interpolation='nearest')
plt.show()

Numpy も使って指数関数を描いてみる

1
2
3
4
5
6
7
8
9
%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt

x = np.arange(0, 10, 0.1)
y = np.exp(x)
plt.plot(x, y)
plt.title("exponential function: $y = e^x$") # タイトルを TeX で書いてくれる
plt.ylim(0, 5000)  # yを0-5000の範囲に限定

破線を入れる

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt

xmin, xmax = -np.pi, np.pi
x = np.arange(xmin, xmax, 0.1)
y_sin = np.sin(x)
y_cos = np.cos(x)
plt.plot(x, y_sin)
plt.plot(x, y_cos)
plt.hlines([-1, 1], xmin, xmax, linestyles="dashed")  # y=-1, 1に破線を描画
plt.title(r"$\sin(x)$ and $\cos(x)$")
plt.xlim(xmin, xmax)
plt.ylim(-1.3, 1.3)
plt.show()

別々に描画して 1 つの図にする

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt

xmin, xmax = -np.pi, np.pi
x = np.arange(xmin, xmax, 0.1)
y_sin = np.sin(x)
y_cos = np.cos(x)


plt.subplot(2, 1, 1)
plt.plot(x, y_sin)
plt.title(r"$\sin x$")
plt.xlim(xmin, xmax)
plt.ylim(-1.3, 1.3)


plt.subplot(2, 1, 2)
plt.plot(x, y_cos)
plt.title(r"$\cos x$")
plt.xlim(xmin, xmax)
plt.ylim(-1.3, 1.3)

plt.tight_layout()  # タイトルの被りを防ぐ

プロットの順番

*plt.subplot() を使う.

1
plt.subplot(行数, 列数, 何番目のプロットか)
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt

def plot_function_name(name, x=0, y=0):
    plt.text(x, y, name, alpha=0.3, size=25, ha="center", va="center")

x = np.arange(-3, 3.01, 0.01) # 3.01

plt.subplot(231)
plot_function_name("1: sin")
plt.plot(x, np.sin(x))

plt.subplot(232)
plot_function_name("2: cos")
plt.plot(x, np.cos(x))

plt.subplot(233)
plot_function_name("3: tan")
plt.plot(x, np.tan(x))
plt.ylim(-1, 1)  # y軸が無限まで行ってしまうので制限

plt.subplot(234)
plot_function_name("4: sinh")
plt.plot(x, np.sinh(x))

plt.subplot(235)
plot_function_name("5: cosh", x=0, y=6)
plt.plot(x, np.cosh(x))

plt.subplot(236)
plot_function_name("6: tanh")
plt.plot(x, np.tanh(x))

plt.show()

オブジェクト指向でプロット

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
x = np.arange(-3, 3, 0.01)
y_sin = np.sin(x)
y_cos = np.cos(x)


fig = plt.figure()
print("type(fig): {}".format(type(fig)))


ax1 = fig.add_subplot(211)
ax2 = fig.add_subplot(212)
print("type(ax1): {}".format(type(ax1)))


ax1.plot(x, y_sin)
ax2.plot(x, y_cos)


ax1.set_ylim(-1.3, 1.3)
ax2.set_ylim(-1.3, 1.3)


ax1.set_title("$\sin x$")
ax2.set_title("$\cos x$")


ax1.set_xlabel("x")
ax2.set_xlabel("x")

fig.tight_layout()  # タイトルとラベルが被るのを解消

複数のグラフを描く時のポイント

普通に描くと縦軸ずれる

org-babel 上では縦軸が見えないのでわかりづらいが, 実際には縦軸の刻みが違う.

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt

n = 300
x = np.random.randn(n)
y1 = np.exp(x) + np.random.randn(n)
y2 = np.exp(x) * 3 + np.random.randn(n)

fig = plt.figure()
ax1 = fig.add_subplot(121)
ax2 = fig.add_subplot(122)

ax1.plot(x, y1, ".")
ax2.plot(x, y2, ".")

sharex, sharey

縦軸を合わせるには sharey を指定する.

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt

fig = plt.figure()
ax1 = fig.add_subplot(121)
ax2 = fig.add_subplot(122, sharey=ax1)

ax1.plot(x, y1, ".")
ax2.plot(x, y2, ".")

1 つの図の中に分割してデータを入れる.

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
n = 300
x = np.random.randn(n)
y = np.sin(x) + np.random.randn(n) * 0.3

fig = plt.figure()


ax1 = fig.add_axes((0, 0.2, 1, 0.8))
ax2 = fig.add_axes((0, 0, 1, 0.2), sharex=ax1)


ax1.tick_params(labelbottom="off")
ax2.tick_params(labelleft="off")

ax1.plot(x, y, "x")
ax2.hist(x, bins=20)

凡例

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt

x = np.linspace(0, 10)
y_sin = np.sin(x)
y_cos = np.cos(x)

plt.plot(x, y_sin, label="sin")
plt.plot(x, y_cos, label="cos")
plt.legend()  # 凡例をグラフにプロット
plt.show()

legend で凡例とラベルを明示的に対応させる

1
2
3
4
5
6
7
8
%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt

p1, = plt.plot(x, y_sin)
p2, = plt.plot(x, y_cos)
plt.legend([p1, p2], ["sin", "cos"])
plt.show()

凡例の位置調節

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt

x = np.random.randn(100)
y1 = x + np.random.randn(100)
y2 = x * 3 + np.random.randn(100)
plt.subplot(121)
plt.plot(x, y1, "r.", label="label 1")
plt.plot(x, y2, "k.", label="label 2")
plt.legend(loc="upper left")

plt.subplot(122)
plt.plot(x, y1, "r.", label="label 1")
plt.plot(x, y2, "k.", label="label 2")
plt.legend(loc="lower right")

plt.show()

Matplotlibでグラフの例: 3Dグラフ

参考情報

その 1: ファンデルワールスの状態方程式

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
%matplotlib inline
from mpl_toolkits.mplot3d import Axes3D
import matplotlib.pyplot as plt
import numpy as np

fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.w_xaxis.set_pane_color((1.0, 1.0, 1.0, 1.0))
ax.w_yaxis.set_pane_color((1.0, 1.0, 1.0, 1.0))
ax.w_zaxis.set_pane_color((1.0, 1.0, 1.0, 1.0))


a = 3.592
b = 0.0427
R = 0.08206


pointsnumber = 500

v = np.linspace(0.064, 0.15, pointsnumber)

t = np.linspace(260, 320, pointsnumber)
X, Y = np.meshgrid(v, t)

Z = (R * Y / (X - b)) - (a / X ** 2)

ax.plot_wireframe(X, Y, Z, rstride=25, cstride=25, lw=0.5, color='b')

plt.xlabel('V')

plt.ylabel('T')
plt.show()

その 2: ワイヤフレーム

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
%matplotlib inline
from mpl_toolkits.mplot3d import Axes3D
import matplotlib.pyplot as plt
import numpy as np

def three_d_data_set():
    # メッシュに分割
    x = np.arange(-3, 3, 0.25)
    y = np.arange(-3, 3, 0.25)

    # meshgrid に入れるのがポイント
    X, Y = np.meshgrid(x, y)

    # 値の計算は Z に入れる: ベクトル計算を使っている
    Z = np.sin(X)+ np.cos(Y)

    return (X, Y, Z)

data3 = three_d_data_set()


fig = plt.figure()

ax = Axes3D(fig)

ax.plot_wireframe(data3[0], data3[1], data3[2])

plt.show()

サーフェス表示

1
2
3
4
5
data3 = three_d_data_set()
fig = plt.figure()
ax = Axes3D(fig)
ax.plot_surface(data3[0], data3[1], data3[2], rstride=1, cstride=1)
plt.show()

3 次元プロット

1
2
3
4
5
data3 = three_d_data_set()
fig = plt.figure()
ax = Axes3D(fig)
ax.plot3D(np.ravel(data3[0]),np.ravel(data3[1]),np.ravel(data3[2]))
plt.show()

等高線 1

1
2
3
4
5
data3 = three_d_data_set()
fig = plt.figure()
ax = Axes3D(fig)
ax.contour3D(data3[0], data3[1], data3[2])
plt.show()

等高線 2

1
2
3
4
5
data3 = three_d_data_set()
fig = plt.figure()
ax = Axes3D(fig)
ax.contourf3D(data3[0], data3[1], data3[2])
plt.show()

Numpyの基本

SymPy だと LaTeX 形式で出力できるが NumPy だとできないようだ. 参考までに SymPy の出力を見てみる.

1
2
3
4
5
6
7
8
9
from sympy import *
from IPython.display import display
init_printing(use_latex=True)

#シンボル定義
a, b, x, y, z, t = symbols ('a b x y z t')
f, g, h = symbols ('f g h', cls=Function)

display(latex(Integral(sqrt(1/x), x)))

基本

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
import numpy as np


print(np.pi)

print(str(np.pi))


print('{:10.4f}'.format(np.pi))
print('{:10.8f}'.format(np.pi))
1
2
3
4
: 3.141592653589793
: 3.141592653589793
:     3.1416
: 3.14159265

多項式

順番に見ていこう.

多項式の作り方

1
2
3
4
5
6
7
8
9
import numpy as np

p1 = np.poly1d([1, 2, 3])
p2 = np.poly1d([2, 3, 4, 5, 6])
p3 = np.poly1d([0, 1], True)
print('配列の 0 番目が最高階の係数となる形式で順次指定されていく')
print(p1)
print('第 2 引数を True にすると指定配列を根とする多項式になる')
print(p3)
1
2
3
4
5
6
: 配列の 0 番目が最高階の係数となる形式で順次指定されていく
:    2
: 1 x + 2 x + 3
: 第 2 引数を True にすると指定配列を根とする多項式になる
:    2
: 1 x - 1 x

多項式の係数

1
2
3
4
5
print(p1)
print('多項式の係数の配列')
print(p1.c)
print('多項式の k 番目の係数')
print(p1[1])
1
2
3
4
5
6
:    2
: 1 x + 2 x + 3
: 多項式の係数の配列
: [1 2 3]
: 多項式の k 番目の係数
: 2

多項式の次数

1
2
print(p1)
print(p1.order)
1
2
3
4
:    2
: 1 x + 2 x + 3
: 多項式の次数
: 2

多項式の評価

1
2
3
print(p1)
print(np.polyval(p1, 5))
print(p1(5))
1
2
3
4
:    2
: 1 x + 2 x + 3
: 38
: 38

多項式の演算

和, 差, 積.

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
print(p1)
print(p2)
print("")

print('多項式の和')
print(np.polyadd(p1, p2))

print('多項式の差: 第1引数 - 第2引数')
print(np.polysub(p1, p2))

print('多項式の積')
print(np.polymul(p1, p2))
print('普通の積記号でも大丈夫:他も同じ')
print(p1 * p2)

べきと割り算関係.

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
print('多項式のべき')
print(p1 ** 2)

q = np.polydiv(p2, p1)
print('多項式の割り算')
print(q)
print('多項式の商')
print(q[0])
print('多項式の割り算:あまり')
print(q[1])
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
多項式のべき
   4     3      2
1 x + 4 x + 10 x + 12 x + 9
多項式の割り算
(poly1d([ 2., -1.,  0.]), poly1d([ 8.,  6.]))
多項式の商
   2
2 x - 1 x
多項式の割り算:あまり

8 x + 6

多項式の微積分

1
2
3
4
5
6
7
print('多項式の導多項式')
for i in range(0, 6):
    print('第 {} 導関数'.format(i))
    print(np.polyder(p2, m = i))

print('多項式の積分')
print(np.polyint(p1))

多項式の根

1
2
3
4
print('多項式の根')
print(np.sort(np.roots(p2)))
print(np.sort(p2.r))
print(np.sort(p3.r))
1
2
3
4
5
6
: 多項式の根
: [-1.08691531-0.76526803j -1.08691531+0.76526803j  0.33691531-1.25867457j
:   0.33691531+1.25867457j]
: [-1.08691531-0.76526803j -1.08691531+0.76526803j  0.33691531-1.25867457j
:   0.33691531+1.25867457j]
: [ 0.  1.]

配列

配列のすべての要素に対して同じ処理ができることがある. 高速な上に凄まじく便利なので使い慣れる一択しかない.

1
2
3
4
a = np.array([1.7, 233.3234235, 21, 234.23421, 2.621321])
print(np.trunc(a))  # 配列の全ての数に対して切り捨て
print(np.ceil(a))   # 配列の全ての数に対して切り上げ
print(np.around(a)) # 配列の全ての数に対して四捨五入
1
2
3
: [   1.  233.   21.  234.    2.]
: [   2.  234.   21.  235.    3.]
: [   2.  233.   21.  234.    3.]

関数

階乗

1
2
3
4
5
6
import numpy as np
from scipy.special import factorial

args = np.arange(10)
print(np.math.factorial(6))
print(factorial(args)) # scipy.special の factorial は引数に配列を当てられる: numpy だと駄目らしい
1
2
3
4
: 720
: [  1.00000000e+00   1.00000000e+00   2.00000000e+00   6.00000000e+00
:    2.40000000e+01   1.20000000e+02   7.20000000e+02   5.04000000e+03
:    4.03200000e+04   3.62880000e+05]

指数

1
2
3
exponents = np.arange(0, 10, 0.1)
powers = np.power(np.e, exponents)
print(powers)
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
[  1.00000000e+00   1.10517092e+00   1.22140276e+00   1.34985881e+00
   1.49182470e+00   1.64872127e+00   1.82211880e+00   2.01375271e+00
   2.22554093e+00   2.45960311e+00   2.71828183e+00   3.00416602e+00
   3.32011692e+00   3.66929667e+00   4.05519997e+00   4.48168907e+00
   4.95303242e+00   5.47394739e+00   6.04964746e+00   6.68589444e+00
   7.38905610e+00   8.16616991e+00   9.02501350e+00   9.97418245e+00
   1.10231764e+01   1.21824940e+01   1.34637380e+01   1.48797317e+01
   1.64446468e+01   1.81741454e+01   2.00855369e+01   2.21979513e+01
   2.45325302e+01   2.71126389e+01   2.99641000e+01   3.31154520e+01
   3.65982344e+01   4.04473044e+01   4.47011845e+01   4.94024491e+01
   5.45981500e+01   6.03402876e+01   6.66863310e+01   7.36997937e+01
   8.14508687e+01   9.00171313e+01   9.94843156e+01   1.09947172e+02
   1.21510418e+02   1.34289780e+02   1.48413159e+02   1.64021907e+02
   1.81272242e+02   2.00336810e+02   2.21406416e+02   2.44691932e+02
   2.70426407e+02   2.98867401e+02   3.30299560e+02   3.65037468e+02
   4.03428793e+02   4.45857770e+02   4.92749041e+02   5.44571910e+02
   6.01845038e+02   6.65141633e+02   7.35095189e+02   8.12405825e+02
   8.97847292e+02   9.92274716e+02   1.09663316e+03   1.21196707e+03
   1.33943076e+03   1.48029993e+03   1.63598443e+03   1.80804241e+03
   1.99819590e+03   2.20834799e+03   2.44060198e+03   2.69728233e+03
   2.98095799e+03   3.29446808e+03   3.64095031e+03   4.02387239e+03
   4.44706675e+03   4.91476884e+03   5.43165959e+03   6.00291222e+03
   6.63424401e+03   7.33197354e+03   8.10308393e+03   8.95529270e+03
   9.89712906e+03   1.09380192e+04   1.20883807e+04   1.33597268e+04
   1.47647816e+04   1.63176072e+04   1.80337449e+04   1.99303704e+04]

対数

自然対数、常用対数、2 進対数の関数がある.

1
2
3
4
args = np.arange(0.1, 10, 0.1)
print(np.log(np.e))
print(np.log10(10))
print(np.log2(2))
1
2
3
: 1.0
: 1.0
: 1.0

$\log$ も配列を食わせて配列を返せる.

1
2
3
print(np.log(args))
print(np.log10(args))
print(np.log2(args))
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
[-2.30258509 -1.60943791 -1.2039728  -0.91629073 -0.69314718 -0.51082562
 -0.35667494 -0.22314355 -0.10536052  0.          0.09531018  0.18232156
  0.26236426  0.33647224  0.40546511  0.47000363  0.53062825  0.58778666
  0.64185389  0.69314718  0.74193734  0.78845736  0.83290912  0.87546874
  0.91629073  0.95551145  0.99325177  1.02961942  1.06471074  1.09861229
  1.13140211  1.16315081  1.19392247  1.22377543  1.25276297  1.28093385
  1.30833282  1.33500107  1.36097655  1.38629436  1.41098697  1.43508453
  1.45861502  1.48160454  1.5040774   1.5260563   1.54756251  1.56861592
  1.58923521  1.60943791  1.62924054  1.64865863  1.66770682  1.68639895
  1.70474809  1.7227666   1.74046617  1.75785792  1.77495235  1.79175947
  1.80828877  1.82454929  1.84054963  1.85629799  1.87180218  1.88706965
  1.90210753  1.91692261  1.93152141  1.94591015  1.96009478  1.97408103
  1.98787435  2.00148     2.01490302  2.02814825  2.04122033  2.05412373
  2.06686276  2.07944154  2.09186406  2.10413415  2.11625551  2.12823171
  2.14006616  2.1517622   2.16332303  2.17475172  2.18605128  2.19722458
  2.20827441  2.21920348  2.2300144   2.24070969  2.2512918   2.2617631
  2.27212589  2.28238239  2.29253476]
[-1.         -0.69897    -0.52287875 -0.39794001 -0.30103    -0.22184875
 -0.15490196 -0.09691001 -0.04575749  0.          0.04139269  0.07918125
  0.11394335  0.14612804  0.17609126  0.20411998  0.23044892  0.25527251
  0.2787536   0.30103     0.32221929  0.34242268  0.36172784  0.38021124
  0.39794001  0.41497335  0.43136376  0.44715803  0.462398    0.47712125
  0.49136169  0.50514998  0.51851394  0.53147892  0.54406804  0.5563025
  0.56820172  0.5797836   0.59106461  0.60205999  0.61278386  0.62324929
  0.63346846  0.64345268  0.65321251  0.66275783  0.67209786  0.68124124
  0.69019608  0.69897     0.70757018  0.71600334  0.72427587  0.73239376
  0.74036269  0.74818803  0.75587486  0.76342799  0.77085201  0.77815125
  0.78532984  0.79239169  0.79934055  0.80617997  0.81291336  0.81954394
  0.8260748   0.83250891  0.83884909  0.84509804  0.85125835  0.8573325
  0.86332286  0.86923172  0.87506126  0.88081359  0.88649073  0.8920946
  0.89762709  0.90308999  0.90848502  0.91381385  0.91907809  0.92427929
  0.92941893  0.93449845  0.93951925  0.94448267  0.94939001  0.95424251
  0.95904139  0.96378783  0.96848295  0.97312785  0.97772361  0.98227123
  0.98677173  0.99122608  0.99563519]
[-3.32192809 -2.32192809 -1.73696559 -1.32192809 -1.         -0.73696559
 -0.51457317 -0.32192809 -0.15200309  0.          0.13750352  0.26303441
  0.37851162  0.48542683  0.5849625   0.67807191  0.76553475  0.84799691
  0.92599942  1.          1.07038933  1.13750352  1.20163386  1.26303441
  1.32192809  1.37851162  1.43295941  1.48542683  1.5360529   1.5849625
  1.63226822  1.67807191  1.72246602  1.76553475  1.80735492  1.84799691
  1.88752527  1.92599942  1.96347412  2.          2.03562391  2.07038933
  2.10433666  2.13750352  2.169925    2.20163386  2.23266076  2.26303441
  2.29278175  2.32192809  2.35049725  2.37851162  2.40599236  2.43295941
  2.45943162  2.48542683  2.51096192  2.5360529   2.56071495  2.5849625
  2.60880924  2.63226822  2.65535183  2.67807191  2.70043972  2.72246602
  2.7441611   2.76553475  2.78659636  2.80735492  2.82781902  2.84799691
  2.86789646  2.88752527  2.9068906   2.92599942  2.94485845  2.96347412
  2.98185265  3.          3.01792191  3.03562391  3.05311134  3.07038933
  3.08746284  3.10433666  3.1210154   3.13750352  3.15380534  3.169925
  3.18586655  3.20163386  3.21723072  3.23266076  3.24792751  3.26303441
  3.27798475  3.29278175  3.30742853]

その他

sin, cos, tan, arcsin, arccos, arctan などもある.

度数とラジアンの変換

arange* で一気に配列を作り degrees** を表示する.

1
2
degrees = np.arange(0, 181, 0.1)
print(degrees)
1
2
: [  0.00000000e+00   1.00000000e-01   2.00000000e-01 ...,   1.80700000e+02
:    1.80800000e+02   1.80900000e+02]

度数をラジアンに変換.

1
2
radians = np.deg2rad(degrees)
print(radians)
1
2
: [  0.00000000e+00   1.74532925e-03   3.49065850e-03 ...,   3.15380996e+00
:    3.15555529e+00   3.15730062e+00]

ラジアンを度数に変換.

1
print(np.rad2deg(radians))
1
2
: [  0.00000000e+00   1.00000000e-01   2.00000000e-01 ...,   1.80700000e+02
:    1.80800000e+02   1.80900000e+02]

NumPyでアニメーション

NumPyでアニメーション: 数学編

参考 1

リーマン積分, 導関数と接線, フーリエ級数に関するアニメーションがある.

NumPyで常微分方程式

1 階の線型常微分方程式

一番単純でしかも実際に使われる微分方程式としてまずは 1 階の線型常微分方程式を考えよう. ちょっと不吉な例であるが放射性物質の崩壊の方程式を紹介する. 導出をしたければちゃんと物理を勉強してもらう必要がある. ここでは物理は省略して数学に集中する.

\begin{align} \frac{dx}{dt} = - c x. \end{align}

厳密解は $x = C_0 e^{-ct}$ だ. 初期値を設定すれば $C_0$ はそこから決まる.

微分を単純に離散化すると次のようになる.

\begin{align} \frac{x_{n+1} - x_{n}}{\Delta t} = -c x_{n}. \end{align}

$\Delta t$ は $h$ と書くこともある. 整理すると次の通り.

\begin{align} x_{n+1} = x_{n} - c (\Delta t) x_{n}. \end{align}

これに沿って計算したのがいわゆるオイラー法. これをコードに落としてみよう.

オイラー法

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
from matplotlib import animation

def radioactive_euler(nt, init = 10):
    dt = 2 / (nt - 1)
    # 初期条件設定
    x = np.zeros(nt)
    x[0] = init

    for i in range(1, nt):
        x[i] = x[i-1] - c * dt * x[i-1]

    # ベクトル計算で書き直したい

    return x


c = 1
nt = 101
init = 5
x1 = radioactive_euler(nt, init)
plt.plot(np.linspace(0, 2, nt), x1)


t = np.linspace(0, 2, nt)
x2 = init * np.exp(- c * t)
plt.plot(t, x2)


plt.legend(['approximation', 'exact'])

file_name = "tmp/ode1.tmp.png"
plt.savefig(file_name)

plt.show()

ルンゲ-クッタ法

少なくとも上の図示の限りではほとんど違いが見られないくらい精度がいい. ただいろいろやっているとオイラー法の限界もいろいろ見える. ちょっと先走ってルンゲ-クッタ法で実装してみよう. 結果は次の通り: 面倒なので導出は省略する.

とりあえず一般的な次の次の初期値問題を考えよう. \begin{align} \frac{dx}{dt} = f(t, x), \quad x(t_0) = x_{0}. \end{align}

とりあえず右辺は $f$ で一般化しておいた. 上の例でいうと $f(t, x) = x$ だ.

ではルンゲ-クッタ法による式 (数列から見れば漸化式) を書いてみる.

\begin{align} x_{n+1} &= x_{n} + \frac{h}{6} (k_{1} + 2 k_{2} + 2 k_{3} + k_{4}), \ t_{n+1} &= t_{n} + h. \end{align}

$h$ はさっき $\Delta t$ と書いた量だ. そして $k_{i}$ は次のように定義している.

\begin{align} k_{1} &= f \left(t_{n}, x_{n} \right), \ k_{2} &= f \left(t_{n} + \frac{h}{2}, x_{n} + \frac{h}{2} k_{1} \right), \ k_{3} &= f \left(t_{n} + \frac{h}{2}, x_{n} + \frac{h}{2} k_{2} \right), \ k_{4} &= f \left(t_{n} + h, x_{n}+ h k_{3} \right). \end{align}

放射性物質の崩壊の方程式で書き直す

最初の式は変わらない. $k_{i}$ を書き直せばいいだけだ.

\begin{align} x_{n+1} &= x_{n} + \frac{h}{6} (k_{1} + 2 k_{2} + 2 k_{3} + k_{4}), \ t_{n+1} &= t_{n} + h, \ k_{1} &= x_{n}, \ k_{2} &= x_{n} + \frac{h}{2} k_{1}, \ k_{3} &= x_{n} + \frac{h}{2} k_{2}, \ k_{4} &= x_{n}+ h k_{3}. \end{align}

ではコードに落としてみよう.

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
from matplotlib import animation

def radioactive_rk(nt, init = 10):
    dt = 2 / (nt - 1)
    # 初期条件設定
    x = np.zeros(nt)
    x[0] = init

    # ベクトル計算で書き直したい
    for i in range(1, nt):
        k1 = x[i-1]
        k2 = x[i-1] + k1 * dt / 2
        k3 = x[i-1] + k2 * dt / 2
        k4 = x[i-1] + k3 * dt
        x[i] = x[i-1] - c * dt / 6 * (k1 + 2 * k2 + 2 * k3 + k4)

    return x


c = 1
nt = 101
init = 5
x1 = radioactive_rk(nt, init)
t = np.linspace(0, 2, nt)
plt.plot(np.linspace(0, 2, nt), x1)


x2 = init * np.exp(- c * t)
plt.plot(t, x2)


plt.legend(['approximation', 'exact'])

plt.savefig("tmp/ode2.tmp.png")
plt.show()

TODO 連立の 1 階常微分方程式

ローレンツ方程式を見てみよう.

カオスで有名な次の方程式. 連立でしかも非線型なので厳密解は書けない (知られていない).

\begin{align} \frac{dx}{dt} &= -px + py, \ \frac{dy}{dt} &= -xz + rx - y, \ \frac{dz}{dt} &= xy - bz. \end{align}

いったん scipyodeint を使った実装を書いておく.

*TODO オイラー, ルンゲ-クッタで書き直す.

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
import numpy as np
from scipy.integrate import odeint
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D


def func(v, t, p, r, b):
    return [-p*v[0]+p*v[1], -v[0]*v[2]+r*v[0]-v[1], v[0]*v[1]-b*v[2]]

p = 10
r = 28
b = 8/3
v0 = [0.1, 0.1, 0.1]
t = np.arange(0, 100, 0.01)

v = odeint(func, v0, t, args=(p, r, b))

fig = plt.figure()
ax = fig.gca(projection='3d')
ax.plot(v[:, 0], v[:, 1], v[:, 2])

2 階の常微分方程式

次に 2 階の常微分方程式を紹介しよう. 高校の物理で出てくるばねの振動(単振動)がまさにこの例だ. 項を増やすと減衰振動になったり、外力をつけたりといろいろなケースがある. まずは一番単純な式を考えよう.

\begin{align} \frac{d^2 x}{dt^2} = - \omega^2 x. \end{align}

さっきのオイラー法なりルンゲ-クッタ法なりは 1 階の方程式に対する計算法なので直接は使えない. この場合は中間処理として $v = dx/dt$ を置いて計算すればいい. これは単なる数値計算の便法ではない. 速度の意味もあるから, という表面的な理由ではなくもっと深く解析力学の文脈で物理としても大事な視点だ. もっといえばシンプレクティック計算法などもっといい計算法にも発展する.

とりあえずオイラー法で計算したい. まずは微分方程式自体を書き直す.

\begin{align} \frac{dx}{dt} = v, \quad \frac{dv}{dt} = - \omega^2 x. \end{align}

これをオイラー法で近似しよう.

\begin{align} x_{n+1} = x_{n} + h v_{n}, \quad v_{n+1} = v_{n} - h \omega^{2} x_{n} \end{align}

オイラー法をコードに落とす.

オイラー法のコード

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
from matplotlib import animation

def harmonic_euler(nt, init = (5, 0)):
    dt = t_range / (nt - 1)
    # 初期条件設定
    x = np.zeros(nt)
    v = np.zeros(nt)
    x[0] = init[0]
    v[0] = init[1]

    for i in range(1, nt):
        x[i] = x[i-1] + dt * v[i-1]
        v[i] = v[i-1] - dt * (omega ** 2) * x[i-1]

    # ベクトル計算で書き直したい

    return (x, v)


omega = 2 * np.pi
nt = 101
t_range = 2
init = (5, 0)
harm = harmonic_euler(nt, init)
t = np.linspace(0, 2, nt)

x_exact = init[0] * np.cos(- omega * t)
v_exact = - omega * init[0] * np.sin(omega * t)


plt.subplot(3, 1, 1)
plt.title('x-t graph')
plt.plot(np.linspace(0, 2, nt), harm[0])
plt.plot(t, x_exact)
plt.legend(['x approximation', 'x exact'])

plt.subplot(3, 1, 2)
plt.title('v-t graph')
plt.plot(np.linspace(0, 2, nt), harm[1])
plt.plot(t, v_exact)
plt.legend(['v approximation', 'v exact'])

plt.subplot(3, 1, 3)
plt.title('phase space')
plt.plot(harm[0], harm[1])


plt.tight_layout()
plt.savefig("tmp/ode5.tmp.png")

見ての通り時間が進むごとに誤差が大きくなる. nt を大きくすると少しはましになる. 実際に上のコードで nt を大きくして再計算してみてほしい.

あとまずいのは phase space の図だ. この系はエネルギーが保存する系だから相空間内の軌道が閉じてほしいのにそうなっていない. シンプレクティックにやれば解消できるようだが, とにかくここではよろしくない.

この方程式でオイラー法はよろしくないことがわかった. とりあえずルンゲ-クッタでやってみよう.

ルンゲ-クッタ法

とりあえず近似式を書く.

\begin{align} x_{n+1} &= x_{n} + \frac{h}{6} (k_{1} + 2 k_{2} + 2 k_{3} + k_{4}), \ t_{n+1} &= t_{n} + h, \ k_{1} &= v_{n}, \ k_{2} &= v_{n} + \frac{h}{2} k_{1}, \ k_{3} &= v_{n} + \frac{h}{2} k_{2}, \ k_{4} &= v_{n} + h k_{3}. \end{align}

次が $v$ の式.

\begin{align} v_{n+1} &= v_{n} + \frac{h}{6} (k_{1} + 2 k_{2} + 2 k_{3} + k_{4}), \ t_{n+1} &= t_{n} + h, \ k_{1} &= - \omega^2 x_{n}, \ k_{2} &= x_{n} - \frac{h}{2} \omega^2 k_{1}, \ k_{3} &= x_{n} - \frac{h}{2} \omega^2 k_{2}, \ k_{4} &= x_{n} - h \omega^2 k_{3}. \end{align}

ルンゲ-クッタ法のコード

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
import numpy as np
import matplotlib.pyplot as plt
from matplotlib import animation

def harmonic_rk(nt, init = 10):
    dt = t_range / (nt - 1)
    # 初期条件設定
    x = np.zeros(nt)
    v = np.zeros(nt)
    x[0] = init[0]
    v[0] = init[1]

    def fx(t, x, v):
        return v

    def fv(t, x, v):
        return - (omega ** 2) * x

    # ベクトル計算で書き直したい
    for i in range(1, nt):
        xk1 = fx(dt * (i - 1), x[i-1], v[i-1])
        vk1 = fv(dt * (i - 1), x[i-1], v[i-1])

        xk2 = fx(dt * (i - 1/2), x[i-1] + xk1 * dt / 2, v[i-1] + vk1 * dt / 2)
        vk2 = fv(dt * (i - 1/2), x[i-1] + xk1 * dt / 2, v[i-1] + vk1 * dt / 2)

        xk3 = fx(dt * (i - 1/2), x[i-1] + xk2 * dt / 2, v[i-1] + vk2 * dt / 2)
        vk3 = fv(dt * (i - 1/2), x[i-1] + xk2 * dt / 2, v[i-1] + vk2 * dt / 2)

        xk4 = fx(dt * (i - 1),   x[i-1] + xk3 * dt,     v[i-1] + vk3 * dt)
        vk4 = fv(dt * (i - 1),   x[i-1] + xk3 * dt,     v[i-1] + vk3 * dt)

        x[i] = x[i-1] + dt / 6 * (xk1 + 2 * xk2 + 2 * xk3 + xk4)
        v[i] = v[i-1] + dt / 6 * (vk1 + 2 * vk2 + 2 * vk3 + vk4)

    return (x, v)


omega = 2 * np.pi
nt = 101
t_range = 2
init = (5, 0)
harm = harmonic_rk(nt, init)
t = np.linspace(0, 2, nt)


x_exact = init[0] * np.cos(- omega * t)
v_exact = - omega * init[0] * np.sin(omega * t)


plt.subplot(3, 1, 1)
plt.title('x-t graph')
plt.plot(np.linspace(0, 2, nt), harm[0])
plt.plot(t, x_exact)
plt.legend(['x approximation', 'x exact'])

plt.subplot(3, 1, 2)
plt.title('v-t graph')
plt.plot(np.linspace(0, 2, nt), harm[1])
plt.plot(t, v_exact)
plt.legend(['v approximation', 'v exact'])

plt.subplot(3, 1, 3)
plt.title('phase space')
plt.plot(harm[0], harm[1])


plt.tight_layout()

NumPyで確率論

整数の乱数を発生させるモジュール randint

1
2
3
4
%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
from numpy.random import randint

1-6 の乱数を 10 個発生させる

1
2
from numpy.random import randint
print(randint(1, 7, 10))
1
: [1 2 2 3 2 3 4 4 1 4]

1-6 の乱数を 1000 回発生させてヒストグラムに表示

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
from numpy.random import randint


dice = randint(1, 7, 1000)

fig = plt.figure(figsize=(3, 3))
subplot = fig.add_subplot(1, 1, 1)
subplot.set_xlim(0.5, 6.5)
subplot.hist(dice, bins=np.linspace(0.5, 6.5, 7))
plt.show()

2 個サイコロを振った合計を 1000 回分計算

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
from numpy.random import randint

dice1 = randint(1, 7, 1000)
dice2 = randint(1, 7, 1000)

dice_total = dice1 + dice2

fig = plt.figure(figsize=(5, 3))
subplot = fig.add_subplot(1, 1, 1)
subplot.set_xlim(1.5, 12.5)
subplot.hist(dice_total, bins=np.linspace(1.5, 12.5, 12))
plt.show()

他のモジュール: random, normal

この 2 つも試してみよう.

random で 0-1 の範囲の乱数を 5 個作る

1
2
3
4
5
6
7
%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
from numpy.random import random
from numpy.random import normal

print(random(5))

100-150 の範囲の乱数

1
2
3
4
5
6
7
%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
from numpy.random import random
from numpy.random import normal

print(100 + random(5) * 50)
1
: [ 125.5122101   146.17318573  100.94123889  114.62349853  123.51461347]

正規分布の乱数 normal

標準偏差 1 の場合と標準偏差 2 の場合をグラフで比較してみる.

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
from numpy.random import random
from numpy.random import normal

fig = plt.figure(figsize=(5, 6))


subplot = fig.add_subplot(2, 1, 1)
data = normal(5, 1, 1000)
subplot.hist(data, bins=np.linspace(0, 10, 20))


subplot = fig.add_subplot(2, 1, 2)
data = normal(5, 2, 1000)
subplot.hist(data, bins=np.linspace(0, 10, 20))

plt.show()

大数の法則

1000 回分のサイコロの目を用意して各回の平均を順に計算してグラフに描く.

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
from numpy.random import randint


dice = randint(1,7,1000)
result = []
for i in range(1, 1001):
    sample = dice[:i]
    result.append(np.mean(sample))

fig = plt.figure(figsize=(5, 3))
subplot = fig.add_subplot(1, 1, 1)
subplot.set_ylim(2, 5)
subplot.plot(result)
plt.show()

コイン投げ

1 回分のゲームをシミュレーションする関数を用意する. 手持ちのポイントと賭けるポイントを渡すと, ゲームの結果に応じて増減した後の手持ちのポイントが返るようにする. 毎回 10 ポイント賭けるという戦略で 100 回分ゲームを繰り返す. その結果で手持ちのポイントがどう変わるかをグラフに描く.

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
from numpy.random import randint

def coin_game(money, bet):
    coin = randint(2)
    if coin == 0:
        money += bet
    else:
        money -= bet
    return money

money = 1000
result = []
for i in range(100):
    # 10 ポイントごと賭ける
    money = coin_game(money, 10)
    result.append(money)

fig = plt.figure(figsize=(5,3))
subplot = fig.add_subplot(1,1,1)
subplot.plot(range(len(result)), result)
plt.show()

毎回手持ちのポイントの半分を賭けるという戦略を取って 100 回ゲームする. その結果手持ちのポイントがどう変わるかをグラフに描く.

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
money = 1000
result1 = []
for i in range(100):
    # 半分ごと賭ける
    money = coin_game(money, money/2)
    result1.append(money)

fig = plt.figure()
ax1 = fig.add_subplot(121)
ax2 = fig.add_subplot(122)
ax1.plot(range(len(result1)), result1)


money = 1000
result2 = []
for i in range(100):
    result2.append(money)
    # 半分ごと賭ける
    money = coin_game(money, money/2)

ax2.plot(range(len(result2)), result2)

株の値動きのシミュレーションのようなもの

ちょっと遊んでみよう.

1 回分のゲームをシミュレーション

手持ちのポイントは賭けたポイントの分だけ増えるか, 賭けたポイントの半分だけ減る.

毎回 10 ポイント賭ける戦略で 100 回ゲームを繰り返す. 手持ちのポイントがどう変わるかグラフで描く.

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
from numpy.random import randint

def coin_game(money, bet):
    coin = randint(2)
    if coin == 0:
        money += bet
    else:
        money -= bet / 2
    return money

money = 1000
result = []
for i in range(100):
    money = coin_game(money, 10)
    result.append(money)

fig = plt.figure(figsize=(5,3))
subplot = fig.add_subplot(1,1,1)
subplot.plot(range(len(result)), result)
plt.show()

毎回手持ちのポイントを半分賭ける

100 回ゲームを繰り返す. 手持ちのポイントがどう変わるかをグラフで描く.

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
money = 1000
result = []
for i in range(100):
    result.append(money)
    money = coin_game(money, money / 2)

fig = plt.figure(figsize=(5, 3))
subplot = fig.add_subplot(1, 1, 1)
subplot.plot(range(len(result)), result)
plt.show()

モンテカルロ法

参考

一様乱数

numpy.random.rand() で 0-1 の一様乱数を生成する. 引数を指定すれば複数の乱数を生成できる. 乱数の範囲を変えたい場合は後からベクトル演算する.

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
from numpy.random import *

print('0-1の乱数を1個生成')
print(rand())
print('')

print('0-1の乱数を100個生成')
print(rand(100))
print('')

print('0-1の乱数で 10x10 の行列を生成')
print(rand(10,10))
print('')

print('30-70の乱数を100個生成')
print(rand(100) * 40 + 30)
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
0-1の乱数を1個生成
0.2991300250393705

0-1の乱数を100個生成
[ 0.96008499  0.4876809   0.2032633   0.3515062   0.09395832  0.82046741
  0.09894689  0.88181204  0.91717783  0.67966579  0.70159484  0.77883058
  0.36860837  0.1836394   0.59647458  0.94684601  0.77424186  0.19880258
  0.63949375  0.8190393   0.94407509  0.10071247  0.26399755  0.56802998
  0.45412078  0.57331073  0.49222259  0.55373428  0.58883185  0.04028661
  0.84054433  0.91416415  0.62588206  0.04036401  0.95063797  0.91012057
  0.8069295   0.02557208  0.71591708  0.02961248  0.93197455  0.21329549
  0.31917112  0.46084521  0.47293933  0.17530845  0.31995878  0.40964597
  0.61651931  0.78196877  0.04742495  0.38418352  0.75048599  0.38486768
  0.89891665  0.61948902  0.44826018  0.40118431  0.30269675  0.7331052
  0.60302038  0.26506878  0.98582242  0.05995304  0.00127997  0.87175875
  0.00854238  0.59077752  0.33181066  0.68507039  0.84599931  0.05243945
  0.16143193  0.96565688  0.85844458  0.64651156  0.70840907  0.79623183
  0.61208582  0.19972924  0.28507776  0.62411044  0.82616062  0.20850016
  0.55577483  0.87397947  0.67298962  0.47055815  0.23906173  0.74671261
  0.54771724  0.28208002  0.36276463  0.59117694  0.33044094  0.57519493
  0.53522023  0.86243342  0.22243464  0.04603091]

0-1の乱数で 10x10 の行列を生成
[[ 0.52351035  0.77780342  0.06046234  0.05099905  0.53615257  0.25082531
   0.93969339  0.39506631  0.4688679   0.0323669 ]
 [ 0.40057593  0.12485139  0.22740708  0.94114815  0.79743915  0.9309769
   0.38248251  0.28966811  0.06485752  0.43783663]
 [ 0.67432477  0.14535369  0.63298396  0.33801197  0.51658446  0.37371677
   0.61763152  0.76407791  0.22566605  0.96020584]
 [ 0.39898091  0.16531665  0.37862637  0.41239287  0.08324925  0.06900096
   0.48383583  0.40779061  0.40727017  0.89574793]
 [ 0.35067508  0.41216795  0.11103469  0.47683551  0.26370249  0.84601269
   0.70482798  0.95529879  0.52109438  0.598296  ]
 [ 0.76916897  0.06869007  0.49815425  0.31708733  0.74498014  0.75691811
   0.52486371  0.82246753  0.738168    0.2484319 ]
 [ 0.02245602  0.43706805  0.75552175  0.37756288  0.35889601  0.85616382
   0.64400511  0.69913794  0.84546117  0.39521931]
 [ 0.08119396  0.88051937  0.96311977  0.48251766  0.61977593  0.33504418
   0.33484545  0.06191795  0.81446698  0.75770015]
 [ 0.34372958  0.33347022  0.84614375  0.76819626  0.90859494  0.94151852
   0.86645554  0.19259208  0.12609405  0.13717528]
 [ 0.77963946  0.50745395  0.80544665  0.44550312  0.1705092   0.11841045
   0.73157783  0.96556873  0.47356432  0.9708871 ]]

30-70の乱数を100個生成
[ 44.66413134  36.19492985  42.41025844  69.47761869  53.69453721
  69.51140937  52.97038786  41.5271023   40.85007111  67.84417929
  68.75254461  53.92350523  34.84851189  67.44000515  58.47783332
  35.06566833  44.47175201  32.88240286  40.22975356  56.27435134
  62.96654307  60.68728179  30.33974845  64.44747346  32.76055753
  56.52732791  36.30078828  53.56613501  67.25497887  56.18455251
  51.72082103  30.89339084  66.34685313  50.50820848  55.90118935
  66.32162376  52.25309297  48.36118134  51.40277549  67.43573023
  42.52509628  63.78551886  31.35449976  52.46083215  38.03023251
  60.15890544  64.9335873   37.6473036   36.64335116  55.93439082
  34.81244331  69.8234886   53.35277006  35.7300399   59.61860678
  61.44432019  58.04037214  46.14783698  67.93400212  62.1192039
  36.04040751  64.84002325  41.54577584  55.50281339  31.43213752
  52.62696038  49.06379487  60.49771853  36.49773853  53.54037806
  40.15673115  58.47193921  53.0984042   31.51360816  42.63096431
  46.35753943  64.04877554  30.25394998  65.08223575  52.53882757
  34.48577627  48.80941172  66.26757892  54.53815429  35.60388037
  30.70455326  59.5685379   55.11094793  38.25641108  64.82042875
  53.23847495  64.48403816  68.72632037  44.77839931  40.06353377
  44.57931327  69.86805955  55.66530533  41.18298176  56.24740705]

モンテカルロ法

半径 1 の四分円の面積を求めることで円周率を計算する. まずは関数定義.

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
from datetime import datetime

def draw_quadrant():
    """四分円の円周を描画"""
    circle_x = np.arange(0,1,0.001)
    circle_y = np.sqrt(1- circle_x * circle_x)
    plt.plot(circle_x, circle_y)

def calculate_pi(N):
    """モンテカルロ法で円周率を近似計算する"""

    # 的に当たった階数
    X = 0

    # 試行時間 測定開始
    start_time = datetime.now()

    # 試行開始
    for i in range(1, N):
        score_x = np.random.rand()
        score_y = np.random.rand()
        if score_x * score_x + score_y * score_y < 1:
            # 的に入ったら赤
            plt.plot(score_x, score_y,"ro")
            X = X + 1
        else:
            # 的から外れたら青
            plt.plot(score_x, score_y,"bo")

        # pi の近似値計算
        pi = 4 * float(X) / float(N)

    # 試行時間 測定終了
    end_time = datetime.now()
    time = end_time - start_time

    print("time = {}".format(time))
    print("$\pi$ = {}".format(pi))

    # 描画
    draw_quadrant()
    plt.grid(True)
    plt.xlabel('X')
    plt.ylabel('Y')
    plt.gca().set_aspect('equal', adjustable='box')
    plt.show()

試行回数 100

1
calculate_pi(100)

試行回数 1000

1
calculate_pi(1000)

試行回数 10000

少し時間がかかるのでそのつもりで実行しよう.

1
calculate_pi(10000)

NumPyで統計学

確率論に入れてもいいのかもしれない. まあいいだろう.

平均 (mean)

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
meanlst = [1,2,3]


lstsum = sum(meanlst)
print(lstsum)


lengthlst = len(meanlst)
print(lengthlst)


mean = lstsum / lengthlst
print(mean)
1
2
3
: 6
: 3
: 2.0

中央値

前提としてリストがソートされているとする. このリストの中央値は次のように定義する: 計算すべきリストの項数が奇数ならちょうど半分のところの値. 項数が偶数なら中央 2 項の平均として定義する.

ソート

1
2
3
4
5
medianlst = [4, 1, 3]
medianlst.sort()


print(medianlst)
1
: [1, 3, 4]

中央値を計算する

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
def calculate_median(numbers):
    u"""中央値を計算する関数.
    引数としては数のリストを取る.
    """

    N = len(numbers)
    # ソートしておく
    numbers.sort()

    if N % 2 == 0:
        # N が偶数

        # 中央の 2 つの値を取得するための項番号を取得
        m1 = N/2
        m2 = (N/2) + 1
        # 整数に変換
        m1 = int(m1) - 1
        m2 = int(m2) - 1

        # 平均を取る
        median = (numbers[m1] + numbers[m2])/2
    else:
        # 中央の項番
        m = (N+1)/2
        # 整数に変換
        m = int(m) - 1
        # 中央値
        median = numbers[m]

    return median

donations = [100, 65, 73, 934, 110, 218, 510, 518, 532, 643, 1128, 1215]
median = calculate_median(donations)
N = len(donations)


print('要素数 {0} で中央値は {1}.'.format(N, median))
1
: 要素数 12 で中央値は 514.0.

最頻値 (mode) と度数分布表 (frequency table)

文章でコードの説明を書くよりも結果を見た方がわかりやすい.

1
2
3
4
5
6
7
8
9
import collections
lst = [4, 2, 1, 3, 4, 2, 2, 2, 1, 5, 6]
c = collections.Counter(lst)
print(c.most_common())
print(c.most_common(1))
print(c.most_common(2))

mode = c.most_common(1)
print(mode)
1
2
3
4
: [(2, 4), (1, 2), (4, 2), (3, 1), (5, 1), (6, 1)]
: [(2, 4)]
: [(2, 4), (1, 2)]
: [(2, 4)]

別のリストでも最頻値を計算してみよう.

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
import collections
def calculate_mode(numbers):
    c = collections.Counter(numbers)
    mode = c.most_common(1)
    return mode[0][0]

scores = [7, 8, 9, 2, 10, 9, 9, 9, 9, 4, 5, 6, 1, 5, 6, 7, 8, 6, 1, 10]
mode = calculate_mode(scores)


print('リストの最頻値は {0}'.format(mode))
1
: リストの最頻値は 9

最頻値が複数ある場合を見てみよう.

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
import collections
def calculate_mode(numbers):
    c = collections.Counter(numbers)
    numbers_freq = c.most_common()
    max_count = numbers_freq[0][1]

    # 最頻値が複数ある場合もある
    modes = [x for x in numbers_freq if x[1] == max_count]
    return modes

scores = [5, 5, 5, 4, 4, 4, 9, 1, 3]
modes = calculate_mode(scores)
print('リストの最頻値は {0}'.format(modes))
1
: リストの最頻値は [(4, 3), (5, 3)]

度数分布表

1
2
3
4
5
6
7
8
9
import collections
def frequency_table(numbers):
    table = collections.Counter(numbers)
    tablelst = ["Number,Frequency"] + ["{0},{1}".format(x[0], x[1]) for x in table.most_common()]
    return "\n".join(tablelst)

scores = [7, 8, 9, 2, 10, 9, 9, 9, 9, 4, 5, 6, 1, 5, 6, 7, 8, 6, 1, 10]
table_string = frequency_table(scores)
print(table_string)
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
Number,Frequency
9,5
6,3
1,2
5,2
7,2
8,2
10,2
2,1
4,1

何となく別コードを書いてみる.

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
import collections
import pprint
pprint.PrettyPrinter(indent=2)
pp = pprint.pprint

scores = [7, 8, 9, 2, 10, 9, 9, 9, 9, 4, 5, 6, 1, 5, 6, 7, 8, 6, 1, 10]
table = collections.Counter(scores)
tablelst = [["Number", "Frequency"]] + [[x[0], x[1]] for x in table.most_common()]

pp(tablelst)
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
[['Number', 'Frequency'],
 [9, 5],
 [6, 3],
 [1, 2],
 [5, 2],
 [7, 2],
 [8, 2],
 [10, 2],
 [2, 1],
 [4, 1]]

上の度数分布表をソートしよう.

1
2
3
4
5
6
7
8
9
import collections
def sorted_frequency_table(numbers):
    table = collections.Counter(numbers)
    tablelst = ["Number,Frequency"] + ["{0},{1}".format(x[0], x[1]) for x in sorted(table.most_common())]
    return "\n".join(tablelst)

scores = [7, 8, 9, 2, 10, 9, 9, 9, 9, 4, 5, 6, 1, 5, 6, 7, 8, 6, 1, 10]
table_string = sorted_frequency_table(scores)
print(table_string)
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
Number,Frequency
1,2
2,1
4,1
5,2
6,3
7,2
8,2
9,5
10,2

分散

範囲 (range)

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
def find_range(numbers):
    lowest = min(numbers)
    highest = max(numbers)
    # 範囲の計算
    r = highest-lowest
    return lowest, highest, r

donations = [100, 60, 70, 900, 100, 200, 500, 500, 503, 600, 1000, 1200]
lowest, highest, r = find_range(donations)


print('最小: {0} 最大: {1} 範囲: {2}'.format(lowest, highest, r))
1
: 最小: 60 最大: 1200 範囲: 1140

偏差 (variance) と標準偏差 (standard deviation)

偏差 $V$ は次式で定義する. \begin{align} V = \frac{\sum_{k=1}^{n} (x_k - m)}{n}. \end{align}

ここで $m$ は平均. 標準偏差は $\sqrt{V}$ で定義する.

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
def calculate_mean(numbers):
    """渡ってきた数のリストの平均を計算する"""
    s = sum(numbers)
    N = len(numbers)
    mean = s/N
    return mean

def find_differences(numbers):
    """平均との差を取った数列を返す"""
    # 平均
    mean = calculate_mean(numbers)
    # 平均との差を取ったリストを返却
    diff = [x - mean for x in numbers]
    return diff

def calculate_variance(numbers):
    """分散を計算する"""
    # 平均との差を取ったリスト
    diff = find_differences(numbers)
    # 平均との差の 2 乗を取る
    squared_diff = [x**2 for x in diff]

    # 分散を計算する
    sum_squared_diff = sum(squared_diff)
    variance = sum_squared_diff / len(numbers)
    return variance

donations = [100, 60, 70, 900, 100, 200, 500, 500, 503, 600, 1000, 1200]
variance = calculate_variance(donations)
print('リストの分散は {0}'.format(variance))


std = variance**0.5
print('リストの標準偏差は {0}'.format(std))
1
2
: リストの分散は 141047.35416666666
: リストの標準偏差は 375.5627166887931

相関係数

定義は次の通り.

\begin{align} \frac{n \sum_{k}^{n} x_{k} y_{k} - \sum_{k=1}^{n} x_{k} \sum_{k=1}^{n} y_{k}}{\sqrt{(n \sum_{k=1}^{n} x_{k}^2 - (\sum_{k=1}^{n} x_{k})^2) (n \sum_{k=1}^{n} y_{k}^2 - (\sum_{k=1}^{n} x_{k})^2)}}. \end{align}

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
def find_corr_x_y(x,y):
    """x と y の相関係数を調べる."""
    n = len(x)

    # 2 数の積
    prod = [xi * yi for (xi, yi) in zip(x,y)]

    sum_prod_x_y = sum(prod)
    sum_x = sum(x)
    sum_y = sum(y)
    squared_sum_x = sum_x ** 2
    squared_sum_y = sum_y ** 2
    # 2 乗する
    x_square = [xi**2 for xi in x]
    # 和を計算する
    x_square_sum = sum(x_square)

    y_square=[yi**2 for yi in x]
    # 和を計算する
    y_square_sum = sum(y_square)
    # 相関係数の公式に代入する変数の準備

    numerator = n*sum_prod_x_y - sum_x*sum_y
    denominator_term1 = n*x_square_sum - squared_sum_x
    denominator_term2 = n*y_square_sum - squared_sum_y
    denominator = (denominator_term1*denominator_term2)**0.5
    correlation = numerator/denominator
    return correlation

NumPyで物理

アニメーションもあるのでとりあえず Jupyter 管理だけにしておく.

Pandas

  • CSV 読み込み
1
2
3
4
5
6
7
message_csv = "メッセージ定義.csv"


df = pd.read_csv(message_csv, encoding="utf-8").fillna("")


last_row_number = len(df.index)
  • 行数を取得
1
2
3
4
5
6
7
message_csv = "メッセージ定義.csv"


df = pd.read_csv(message_csv, encoding="utf-8").fillna("")


last_row_number = len(df.index)
  • CSV のループ
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
message_csv = "メッセージ定義.csv"


df = pd.read_csv(message_csv, encoding="utf-8").fillna("")
last_row_number = len(df.index)
java_values = []
js_values = []

for key, row in df.iterrows():
    # Java 側の内容
    java_separator = ";" if (key + 1) == last_row_number else ","
    # read_csv で header を指定しないと 1 行目の内容で項目が取れる
    condition = "" if row["条件"] == "" else row["条件"] + " "
  • Excel 読み込み
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
import pandas as pd
import pprint
pp = pprint.PrettyPrinter(indent=2)


book = "1.3.1_テーブル定義書.xlsx"
newfilename = "現基幹システムテーブル・カラム一覧.csv"
csvdata = [["種別", "論理テーブル名", "物理テーブル名", "カラム名"]]

i = 0

file = pd.ExcelFile(book)
for i, sheet_name in enumerate(file.sheet_names):
    sheet = file.parse(sheet_name, header=None)
    pass_flag = False

    if i > 0:
        data1 = ["論理名"]
        data2 = ["物理名"]
        for key, row in sheet.iterrows():
            if key == 4 or key == 5:
                data1.append(row[2])
                data2.append(row[2])
            elif key > 12:
                if row[0] == "インデックス情報" or row[0] == "登録ユーザID":
                    pass_flag = True
                if pass_flag:
                    break
                data1.append(row[1])
                data2.append(row[2])

        csvdata.append(data1)
        csvdata.append(data2)

pd.DataFrame(csvdata).to_csv(newfilename)

pip

pip自体のアップデート

1
pip install -U pip

requirements.txtによるライブラリのインストール

1
pip install -r requirements.txt

アップデートが必要なパッケージのリスト

1
pip list -o

インストールされているライブラリの確認

1
pip list

インストールされているライブラリをfreeze形式で表示

1
pip freeze

パッケージのバージョン指定インストール

1
pip install <package-name>==<version>

パッケージをパッケージ名==バージョンで表示

requirements.txtを作るときに便利.

1
pip freeze > requirements.txt

パッケージのアップデート

1
pip install -U <package-name>

依存関係の確認

1
pip check

print

即時表示

flush=True をつけます.

1
print("出力の即時表示", flush=True)

re, 正規表現

compile

1
2
3
4
5
6
import re
pattern = re.compile(r"\n\*\* ")

fname = "a.org"
with open(fname, "r", encoding="utf-8") as f:
    cont = sorted(pattern.split(f.read()), key=str.lower)

Sphinxでのdocstringドキュメント生成

ドキュメント生成

注意

ドキュメント生成時にファイルを実行する模様. 実行部分はif __name__ == "__main__":をつけて隔離する必要がある.

初回だけ

1
2
3
sphinx-quickstart docs

> Separate source and build directories (y/n) [n]:
  • あとは適当に設定.
  • docs/conf.pyを修正.
  • 以下をアンコメント.
1

  • 以下のように修正.
1
sys.path.insert(0, os.path.abspath('../'))
  • conf.pyextensionsを以下のように修正
1
2
3
4
5
extensions = [
    'sphinx.ext.autodoc',
    'sphinx.ext.viewcode',
    'sphinx.ext.todo',
]
  • docs/index.pymain(実際のソースファイル)を追加.

ドキュメント生成コマンド

  • 上記「初回だけ」のディレクトリ構成を前提にする.
  • プロジェクトルートから次のコマンドを実行.
1
sphinx-apidoc -f -o ./docs . & sphinx-build -b singlehtml ./docs ./docs/_build
  • -b singlehtmlは適当なオプションを指定してよい.
  • 参考: 実行時にWARNING: autodoc: failed to import module 'plot' from module 'xxxx'; the following exception was raised: No module named 'matplotlib'のようなエラーが出たら, 次のコードをconf.pyに追加.
1
autodoc_mock_imports = ["some_module_name"]

string, 文字列

いろいろな都合で文字列の項目であってもリストなど他の項目に入っていることもあります.

ファイル名からファイル名と拡張子を取得

os.path.split() でファイル名とフォルダ名のペアが取れます. 拡張子はドットを含むこと, つまり a.txt に対して .txt であることに注意しましょう.

1
2
3
4
5
import os
fname = "string.org"
froot, fext = os.path.split()
#froot: string
#fext:  .org

日時の文字列を取る

datetimeから作ります. 適切なバージョンの datetime オブジェクトのドキュメントを見るといいでしょう.

1
2
3
4
import datetime
dtstr = datetime.datetime.now().strftime('%Y%m%d-%H%M%S')

dt_now = datetime.datetime.now()

文字詰め・パディング

右寄せゼロ埋め zfill

1
s = "1".zfill(3) # 001

右寄せ・左寄せ・中央寄せ rjust, ljust, center

1
2
3
4
5
6
7
8
9
s = '1234'.rjust(8, '0')  # 00001234
s = '1234'.ljust(8, '0')  # 12340000
s = '1234'.center(8, '0') # 00123400

import math
for n in [10**k for k in range(0, 10)]:
    print(f"n      = {str(n).rjust(10, ' ')}")
    print(f"sqrt n = {str(int(math.sqrt(n))).rjust(10, ' ')}")
    print()

文字の連番を取得する

1
map(chr, range(ord('B'), ord('R')+1))

全体コードは次の通り.

1
2
3
for char in map(chr, range(ord('A'), ord('C')+1)):
    s = "\\newcommand{\scr%s}{\mathscr{%s}}" % (char, char)
    print(s)
1
2
3
\newcommand{\scrA}{\mathscr{A}}
\newcommand{\scrB}{\mathscr{B}}
\newcommand{\scrC}{\mathscr{C}}

見ての通りTeXのコマンドを一気に定義したかった.

文字列をリストにする

文字列のsplit()メソッドを使います.

1
ss = "a,b,c".sprit(",")

SymPyことはじめ

参考サイト

プロットのサンプル (関数のグラフを描く) は別項目でまとめる.

サンプル

まずは式絡みの計算を

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
from sympy import *
from IPython.display import display
init_printing()


x = Symbol('x')
y = Symbol('y')


print(x + y + x - 4 * y)
print(latex(x + y + x - 4 * y))

a = Integral(cos(x)*exp(x), x)
print(Eq(a, a.doit()))
print(latex(Eq(a, a.doit())))
1
2
3
4
: 2*x - 3*y
: 2 x - 3 y
: Eq(Integral(exp(x)*cos(x), x), exp(x)*sin(x)/2 + exp(x)*cos(x)/2)
: \int e^{x} \cos{\left (x \right )}\, dx = \frac{e^{x}}{2} \sin{\left (x \right )} + \frac{e^{x}}{2} \cos{\left (x \right )}

\begin{gather} 2 x - 3 y \ \int e^{x} \cos{\left (x \right )}\, dx = \frac{e^{x}}{2} \sin{\left (x \right )} + \frac{e^{x}}{2} \cos{\left (x \right )} \end{gather}

ルート

Jupyter なら display を使った方が適切に LaTeX (MathJax) 表示できる.

1
2
3
4
5
import sympy
from IPython.display import display
init_printing(use_unicode=True)

display(sympy.sqrt(3))
1
2
import sympy
print(latex(sympy.sqrt(3)))
1
: \sqrt{3}

\begin{gather} \sqrt{3} \end{gather}

勝手に計算もやってくれる. 状況によってはよくない.

1
2
3
4
from sympy import symbols
x, y = symbols('x y')

print(latex(sympy.sqrt(8)))
1
: 2 \sqrt{2}

\begin{align} 2 \sqrt{2} \end{align}

多項式のかけ算

1
2
3
4
5
6
7
from sympy import symbols
x, y = symbols('x y')
expr = x + 2*y
print(latex(expr))
print(latex(expr + 1))
print(latex(expr - x))
print(latex(x * expr))
1
2
3
4
: x + 2 y
: x + 2 y + 1
: 2 y
: x \left(x + 2 y\right)

\begin{gather} x + 2 y \ x + 2 y + 1 \ 2 y \ x \left(x + 2 y\right) \end{gather}

微分

1
2
3
4
5
from sympy import *
x, t, z, nu = symbols('x t z nu')

print(latex(sin(x)*exp(x)))
print(latex(diff(sin(x)*exp(x), x)))
1
2
: e^{x} \sin{\left (x \right )}
: e^{x} \sin{\left (x \right )} + e^{x} \cos{\left (x \right )}

\begin{gather} e^{x} \sin{\left (x \right )} \ e^{x} \sin{\left (x \right )} + e^{x} \cos{\left (x \right )} \end{gather}

積分その 1

1
2
print(latex(exp(x)*sin(x) + exp(x)*cos(x)))
print(latex(integrate(exp(x)*sin(x) + exp(x)*cos(x), x)))
1
2
: e^{x} \sin{\left (x \right )} + e^{x} \cos{\left (x \right )}
: e^{x} \sin{\left (x \right )}

\begin{gather} e^{x} \sin{\left (x \right )} + e^{x} \cos{\left (x \right )} \ e^{x} \sin{\left (x \right )} \end{gather}

積分その 2

1
2
print(latex(exp(x)*sin(x) + exp(x)*cos(x)))
print(latex(integrate(exp(x)*sin(x) + exp(x)*cos(x), x)))
1
2
: e^{x} \sin{\left (x \right )} + e^{x} \cos{\left (x \right )}
: e^{x} \sin{\left (x \right )}

\begin{gather} e^{x} \sin{\left (x \right )} + e^{x} \cos{\left (x \right )} \ e^{x} \sin{\left (x \right )} \end{gather}

積分その 3

$\sin x^{2}$ を $\mathbb{R}$ 上で積分する.

1
print(latex(integrate(sin(x**2), (x, -oo, oo))))
1
: \frac{\sqrt{2} \sqrt{\pi}}{2}

\begin{gather} \frac{\sqrt{2} \sqrt{\pi}}{2} \end{gather}

極限

\begin{align} \lim_{x \to 0} \frac{\sin x}{x} = 1. \end{align}

1
print(latex(limit(sin(x)/x, x, 0)))
1
: 1

多項式の根

\begin{align} x^{2} - 2 = 0. \end{align}

1
print(latex(solve(x**2 - 2, x)))
1
: \left [ - \sqrt{2}, \quad \sqrt{2}\right ]

\begin{gather} \left [ - \sqrt{2}, \quad \sqrt{2}\right ] \end{gather}

微分方程式の解

\begin{align} y'' - y = e^{t}. \end{align}

1
2
f = Function('f')
print(latex(dsolve(Eq(f(t).diff(t, t) - f(t), exp(t)), f(t))))
1
: f{\left (t \right )} = C_{2} e^{- t} + \left(C_{1} + \frac{t}{2}\right) e^{t}

\begin{gather} f{\left (t \right )} = C_{2} e^{- t} + \left(C_{1} + \frac{t}{2}\right) e^{t} \end{gather}

行列の固有値

\begin{align} \begin{pmatrix} 1 & 2 \ 2 & 2 \end{pmatrix}. \end{align}

1
print(latex(Matrix([[1, 2], [2, 2]]).eigenvals()))
1
: \left \{ \frac{3}{2} + \frac{\sqrt{17}}{2} : 1, \quad - \frac{\sqrt{17}}{2} + \frac{3}{2} : 1\right \}

\begin{gather} \left { \frac{3}{2} + \frac{\sqrt{17}}{2} : 1, \quad - \frac{\sqrt{17}}{2} + \frac{3}{2} : 1\right } \end{gather}

ベッセル関数を球ベッセル関数で書く

1
print(latex(besselj(nu, z).rewrite(jn)))
1
: \frac{\sqrt{2} \sqrt{z}}{\sqrt{\pi}} j_{\nu - \frac{1}{2}}\left(z\right)

\begin{align} \frac{\sqrt{2} \sqrt{z}}{\sqrt{\pi}} j_{\nu - \frac{1}{2}}\left(z\right) \end{align}

LaTeX コードを出力する

1
print(latex(Integral(cos(x)**2, (x, 0, pi))))
1
: \int_{0}^{\pi} \cos^{2}{\left (x \right )}\, dx

\begin{gather} \int_{0}^{\pi} \cos^{2}{\left (x \right )}\, dx \end{gather}

数値の計算

円周率

1
2
3
display(pi**2)
display(pi.evalf())
display((pi + exp(1)).evalf())

$\sqrt{2}$ を 100 桁まで評価

1
print(latex(sqrt(2).evalf(100)))
1
: 1.414213562373095048801688724209698078569671875376948073176679737990732478462107038850387534327641573

無限大

1
2
3
print(latex(oo))
print(latex(oo > 99999))
print(latex(oo + 1))
1
2
3
: \infty
: \mathrm{True}
: \infty

展開

1
2
3
4
print(latex(expand(x + y, complex=True)))
print(latex(I*im(x) + I*im(y) + re(x) + re(y)))
print(latex(expand(cos(x + y), trig=True)))
print(latex(cos(x)*cos(y) - sin(x)*sin(y)))
1
2
3
4
: \Re{x} + \Re{y} + i \Im{x} + i \Im{y}
: \Re{x} + \Re{y} + i \Im{x} + i \Im{y}
: - \sin{\left (x \right )} \sin{\left (y \right )} + \cos{\left (x \right )} \cos{\left (y \right )}
: - \sin{\left (x \right )} \sin{\left (y \right )} + \cos{\left (x \right )} \cos{\left (y \right )}

\begin{gather} \Re{x} + \Re{y} + i \Im{x} + i \Im{y} \ \Re{x} + \Re{y} + i \Im{x} + i \Im{y} \ - \sin{\left (x \right )} \sin{\left (y \right )} + \cos{\left (x \right )} \cos{\left (y \right )} \ - \sin{\left (x \right )} \sin{\left (y \right )} + \cos{\left (x \right )} \cos{\left (y \right )} \end{gather}

TODO 多項式

多項式本編に統合する?

SymPy 基礎

SymPy との連携

次のコードはエラーになる.

1
2
3
4
from sympy import Symbol
import math
theta = Symbol('theta')
math.sin(theta) + math.sin(theta)

次のようにすると正しく処理できる.

1
2
3
4
5
6
7
from sympy import Symbol
import sympy
import math
theta = Symbol('theta')

print(latex(sympy.sin(math.pi/2)))
print(latex(sympy.sin(theta) + sympy.sin(theta)))
1
2
: 1.0
: 2 \sin{\left (\theta \right )}

\begin{gather} 1.0 \ 2 \sin{\left (\theta \right )} \end{gather}

変数とシンボルのアルファベットは同じでなくてもいい

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
from sympy import Symbol

a = Symbol('x')
print(latex(a + a + 1))

x = Symbol('x')

print(latex(a.name))
print(latex(x.name))


b = Symbol('test')
print(latex(b.name))
print(latex(b**2 + 2 * b + 1))
1
2
3
4
5
: 2 x + 1
: x
: x
: test
: test^{2} + 2 test + 1

\begin{gather} 2 x + 1 \ x \ x \ test \ test^{2} + 2 test + 1 \end{gather}

複数のシンボルを一気に定義できる

1
2
3
4
from sympy import symbols
x,y,z = symbols('x,y,z')
s = x * y + x * y
print(latex(s))
1
: 2 x y

\begin{gather} 2 x y \end{gather}

簡単な計算なら自動処理してくれる

1
2
p = x * (x + x)
print(latex(p))
1
: 2 x^{2}

\begin{gather} 2 x^{2} \end{gather}

少し複雑になると展開などの処理はしない

1
2
p = (x + 2)*(x + 3)
print(latex(p))
1
: \left(x + 2\right) \left(x + 3\right)

\begin{gather} \left(x + 2\right) \left(x + 3\right) \end{gather}

代入

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
from sympy import Symbol
x = Symbol('x')
y = Symbol('y')

expr = x*x + x*y + x*y + y*y
print(latex(expr))
res = expr.subs({x:1, y:2})
print(latex(res))


print(latex(expr.subs({x:1-y})))
1
2
3
: y^{2} + 2 x y + x^{2}
: 9
: \left(1 - y\right)^{2} + 2 y \left(1 - y\right) + y^{2}

\begin{gather} y^{2} + 2 x y + x^{2} \ 9 \ \left(1 - y\right)^{2} + 2 y \left(1 - y\right) + y^{2} \end{gather}

文字列を式に変換

1
2
3
4
5
6
7
8
from sympy import sympify
init_printing(order='lex')
expr = "x**2 + 3*x + x**3 + 2*x"
expr = sympify(expr)
print(latex(expr))


print(latex(2 * expr))
1
2
: x^{3} + x^{2} + 5 x
: 2 x^{3} + 2 x^{2} + 10 x

\begin{align} x^{3} + x^{2} + 5 x \ 2 x^{3} + 2 x^{2} + 10 x \end{align}

簡単な方程式を解く

1 次方程式

1
2
3
4
5
from sympy import Symbol, solve, latex
x = Symbol('x')
expr = x - 5 - 7
sol = solve(expr)
print(latex(sol))
1
: \left [ 12\right ]

\begin{align} \left [ 12\right ] \end{align}

連立 1 次方程式

1
2
3
4
5
x = Symbol('x')
y = Symbol('y')
expr1 = 2*x + 3*y - 6
expr2 = 3*x + 2*y - 12
print(latex(solve((expr1, expr2), dict=True)))
1
: \left [ \left \{ x : \frac{24}{5}, \quad y : - \frac{6}{5}\right \}\right ]

\begin{gather} \left [ \left { x : \frac{24}{5}, \quad y : - \frac{6}{5}\right }\right ] \end{gather}

2 次方程式

1
2
3
4
5
6
7
8
from sympy import solve, latex
x = Symbol('x')
expr = x**2 + 5*x + 4
print(latex(solve(expr, dict=True)))

x = Symbol('x')
expr = x**2 + x + 1
print(latex(solve(expr, dict=True)))
1
2
: \left [ \left \{ x : -4\right \}, \quad \left \{ x : -1\right \}\right ]
: \left [ \left \{ x : - \frac{1}{2} - \frac{\sqrt{3} i}{2}\right \}, \quad \left \{ x : - \frac{1}{2} + \frac{\sqrt{3} i}{2}\right \}\right ]

\begin{gather} \left [ \left { x : -4\right }, \quad \left { x : -1\right }\right ] \ \left [ \left { x : - \frac{1}{2} - \frac{\sqrt{3} i}{2}\right }, \quad \left { x : - \frac{1}{2} + \frac{\sqrt{3} i}{2}\right }\right ] \end{gather}

不等式を解く

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
from sympy import Symbol, sympify, SympifyError, latex
from sympy import solve_poly_inequality, solve_rational_inequalities
from sympy import solve_univariate_inequality, Poly
from sympy.core.relational import Relational, Equality
def isolve(ineq_obj):
    """不等式を解く"""
    x = Symbol('x')
    # 不等式の左辺を取る
    expr = ineq_obj.lhs
    # 関係演算子の取得
    rel = ineq_obj.rel_op

    if expr.is_polynomial():
        # 式が多項式の場合
        p = Poly(expr, x)
        return solve_poly_inequality(p, rel)
    elif expr.is_rational_function():
        # 式が有理式の場合
        p1, p2 = expr.as_numer_denom()
        num = Poly(p1)
        denom = Poly(p2)
        return solve_rational_inequalities([[((num, denom), rel)]])
    else:
        return solve_univariate_inequality(ineq_obj , x, relational=False)

ineq = "x**2 - 4 > 0"
ineq_obj = sympify(ineq)
print(latex(isolve(ineq_obj)))
1
: \left [ \left(-\infty, -2\right), \quad \left(2, \infty\right)\right ]

\begin{gather} \left [ \left(-\infty, -2\right), \quad \left(2, \infty\right)\right ] \end{gather}

物体の投射の式も解ける

1
2
3
4
5
6
from sympy import sin, solve, Symbol
u = Symbol('u')
t = Symbol('t')
g = Symbol('g')
theta = Symbol('theta')
print(latex(solve(u*sin(theta)-g*t, t)))
1
: \left [ \frac{u}{g} \sin{\left (\theta \right )}\right ]

\begin{gather} \left [ \frac{u}{g} \sin{\left (\theta \right )}\right ] \end{gather}

条件指定

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
x = Symbol('x', positive=True)

if x > 0:
    print('Do Something')
else:
    print('Do Something else')

if (x + 5) > 0:
    print('Do Something')
else:
    print('Do Something else')
1
2
: Do Something
: Do Something

極限

まずはこれ。

\begin{align} \lim_{x \to \infty} \frac{1}{x} = 0. \end{align}

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
from sympy import Limit, Symbol, S
x = Symbol('x')


l = Limit(1/x, x, S.Infinity)
print(latex(l))



print(latex(l.doit()))

print(latex(Limit(1/x, x, 0).doit()))
print(latex(Limit(1/x, x, 0, dir='+').doit()))
print(latex(Limit(1/x, x, 0, dir='-').doit()))
1
2
3
4
5
: \lim_{x \to \infty} \frac{1}{x}
: 0
: \infty
: \infty
: -\infty

\begin{gather} \lim_{x \to \infty} \frac{1}{x} \ 0 \ \infty \ \infty \ -\infty \ \end{gather}

不定形

1
2
from sympy import Symbol, sin
print(latex(Limit(sin(x) / x, x, 0).doit()))
1
: 1

\begin{gather} 1 \end{gather}

自然対数の底 $e$

1
2
3
from sympy import Limit, Symbol, S
n = Symbol('n')
print(latex(Limit((1 + 1 / n) ** n, n, S.Infinity).doit()))
1
: e

\begin{gather} e \end{gather}

次の式を評価してみよう.

\begin{align} A = P (1 + \frac{r}{n})^{nt}. \end{align}

1
2
3
4
5
from sympy import Symbol, Limit, S
p = Symbol('p', positive  =True)
r = Symbol('r', positive = True)
t = Symbol('t', positive = True)
print(latex(Limit(p * (1 + r / n) ** (n * t), n, S.Infinity).doit()))
1
: p e^{r t}

\begin{gather} p e^{r t} \end{gather}

微分: 極限を取る

次の関数を微分する.

\begin{align} S(t) = 5 t^2 + 2 t + 8. \end{align}

導関数は次の通り. \begin{align} S'(t) = 10 t + 2. \end{align}

1
2
3
4
5
6
7
8
9
from sympy import Symbol, Limit
t = Symbol('t')
St = 5 * t ** 2 + 2 * t + 8
t1 = Symbol('t1')
delta_t = Symbol('delta_t')
St1 = St.subs({t: t1})
St1_delta = St.subs({t: t1 + delta_t})

print(latex(Limit((St1_delta - St1) / delta_t, delta_t, 0).doit()))
1
: 10 t_{1} + 2

\begin{gather} 10 t_{1} + 2 \end{gather}

微分: 微分するメソッド

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
from sympy import Symbol, Derivative
t = Symbol('t')
St = 5*t**2 + 2*t + 8
print(latex(Derivative(St, t)))

d = Derivative(St, t)
print(latex(d.doit()))


print(latex(d.doit().subs({t:t1})))
print(latex(d.doit().subs({t:1})))
1
2
3
4
: \frac{d}{d t}\left(5 t^{2} + 2 t + 8\right)
: 10 t + 2
: 10 t_{1} + 2
: 12

\begin{gather} \frac{d}{d t}\left(5 t^{2} + 2 t + 8\right) \ 10 t + 2 \ 10 t_{1} + 2 \ 12 \end{gather}

導関数を求めるプログラム

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
from sympy import Symbol, Derivative, sympify, pprint
from sympy.core.sympify import SympifyError

def derivative(f, var):
    """
    @var: 微分に使う変数
    @f: 微分する関数. べきは ^ ではなく ** で入力すること.
    """
    var = Symbol(var)
    d = Derivative(f, var).doit()
    print(latex(d))

f = "x ** 2 + 3 * x + 10"
var = "x"
derivative(f, var)

f = "e ** x + log(x)"
var = "x"
derivative(f, var)
1
2
: 2 x + 3
: e^{x} \log{\left (e \right )} + \frac{1}{x}

\begin{gather} 2 x + 3 \ e^{x} \log{\left (e \right )} + \frac{1}{x} \end{gather}

高階の導関数や極大極小

元の関数と導関数.

1
2
3
4
5
6
from sympy import Symbol, solve, Derivative, pprint
x = Symbol('x')
f = x**5 - 30*x**3 + 50*x
print(latex(f))
d1 = Derivative(f, x).doit()
print(latex(d1))
1
2
: x^{5} - 30 x^{3} + 50 x
: 5 x^{4} - 90 x^{2} + 50

\begin{align} x^{5} - 30 x^{3} + 50 x \ 5 x^{4} - 90 x^{2} + 50 \end{align}

臨界点.

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
critical_points = solve(d1)
print(latex(critical_points))

A = critical_points[2]
B = critical_points[0]
C = critical_points[1]
D = critical_points[3]

d2 = Derivative(f, x, 2).doit()
print(latex(d2))
1
2
: \left [ - \sqrt{- \sqrt{71} + 9}, \quad \sqrt{- \sqrt{71} + 9}, \quad - \sqrt{\sqrt{71} + 9}, \quad \sqrt{\sqrt{71} + 9}\right ]
: 20 x \left(x^{2} - 9\right)

\begin{align} \left [ - \sqrt{- \sqrt{71} + 9}, \quad \sqrt{- \sqrt{71} + 9}, \quad - \sqrt{\sqrt{71} + 9}, \quad \sqrt{\sqrt{71} + 9}\right ] \ 20 x \left(x^{2} - 9\right) \end{align}

極値の評価.

1
2
3
4
print(d2.subs({x:B}).evalf())
print(d2.subs({x:C}).evalf())
print(d2.subs({x:A}).evalf())
print(d2.subs({x:D}).evalf())
1
2
3
4
: 127.661060789073
: -127.661060789073
: -703.493179468151
: 703.493179468151

端点との比較.

1
2
3
4
5
6
x_min = -5
x_max = 5
print(f.subs({x:A}).evalf())
print(f.subs({x:C}).evalf())
print(f.subs({x:x_min}).evalf())
print(f.subs({x:x_max}).evalf())
1
2
3
4
: 705.959460380365
: 25.0846626340294
: 375.000000000000
: -375.000000000000

最急降下法

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
import math
from sympy import Derivative, Symbol, sin
def grad_ascent(x0, f1x, x):
    epsilon = 1e-6
    step_size = 1e-4
    x_old = x0
    x_new = x_old + step_size * f1x.subs({x: x_old}).evalf()
    while abs(x_old - x_new) > epsilon:
        x_old = x_new
        x_new = x_old + step_size * f1x.subs({x: x_old}).evalf()
    return x_new

def find_max_theta(R, theta):
    # Calculate the first derivative
    R1theta = Derivative(R, theta).doit()
    theta0 = 1e-3
    theta_max = grad_ascent(theta0, R1theta, theta)
    return theta_max


g = 9.8


u = 25


theta = Symbol('theta')
R = u ** 2 * sin(2 * theta) / g

theta_max = find_max_theta(R, theta)
print(latex('Theta: {0}'.format(math.degrees(theta_max))))
print(latex('Maximum Range: {0}'.format(R.subs({theta:theta_max}))))
1
2
: Theta: 44.997815081691805
: Maximum Range: 63.7755100185965

一般的な最急降下法

*それなりに実行時間がかかるので注意しよう.

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
from sympy import Derivative, Symbol, sympify

def grad_ascent(x0, f1x, x):
    epsilon = 1e-6
    step_size = 1e-4
    x_old = x0
    x_new = x_old + step_size*f1x.subs({x:x_old}).evalf()
    while abs(x_old - x_new) > epsilon:
        x_old = x_new
        x_new = x_old + step_size*f1x.subs({x:x_old}).evalf()
    return x_new




f = "cos(y)"
var = "y"
var0 = 1

try:
    f = sympify(f)
except SympifyError:
    print('Invalid function entered')
else:
    var = Symbol(var)
    d = Derivative(f, var).doit()
    var_max = grad_ascent(var0, d, var)
    print('{0}: {1}'.format(var.name, var_max))
    print('Maximum value: {0}'.format(f.subs({var:var_max})))
1
2
: y: 0.00999904956272288
: Maximum value: 0.999950009920428

積分

不定積分.

1
2
3
4
5
6
from sympy import Integral, Symbol
x = Symbol('x')
k = Symbol('k')

i = Integral(k*x, x)
print(latex(i.doit()))
1