Python

ガイド

Pythonで数学を

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

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

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

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

Jupyter インストール

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

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

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

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

Anaconda でのインストール

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

Windows (未検証)

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

Mac (検証済み)

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

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

$ 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 上で確認できる. ただあなたはローカルで確認してみたいと思っているかもしれない. そんなときそもそも起動できなければ確認しようもない.

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

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

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

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

簡単なものでもアニメーションやグラフを描けるのが楽しい. 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 を書いた時点では次のバージョンで動かしている.

諸注意

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

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

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

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

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

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

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

%matplotlib nbagg

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

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

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

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

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

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

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

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

式をもとにグラフを描く

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

%matplotlib inline
import matplotlib.pyplot as plt

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

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

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

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

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

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

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

    return numbers

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

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

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

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

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

    return (x,y)

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

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

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

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

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

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

次回予告

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

次回をお楽しみに.

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

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

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

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

諸注意

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

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

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

Jupyter 起動

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

$ 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 系で解説されているもののドットインストールを見て勉強するのもいいだろう.

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

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

Hello, World!

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

 # 計算もできる

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

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

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

 # あまり
print(5 % 3)

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

出力結果は次の通り.

4
2.5
1
2
125

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

 # 文字列のエスケープ
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)

出力結果は次の通りだ.

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 で片付けるべき, というところだけ理解している. それ用のコードも順次公開していきたい.

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 標準の集合を見てみよう.

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)
1
a: 1
dict_keys(['c', 'a', 'b'])
dict_values([3, 1, 2])
{1, 2, 3}
ループ

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

a = [5,6,7,8,9]

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

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

出力は次の通り.

5
6
7
8
9

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

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

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

結果は次の通り.

[36, 64]
条件分岐

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

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)

結果は次の通り.

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

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

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 はいわゆる真偽値. 真偽値という概念があるのだ.

False
True
False
True
以下 IPython

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

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

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

$ ipython

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

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

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

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

Unix コマンド実行結果.

%ls -l

結果は次の通り.

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 でコマンド実行する機会もあるだろうし. そういう場合には便利だろう.

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

結果は次の通り.

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

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

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

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

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

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

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

出力結果は次の通り.

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

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

%%writefile myfile.txt
Hello world!

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

Overwriting myfile.txt

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

<!-- 前のセルからの続き -->
!more myfile.txt
Hello world!
[?1l>
<!-- 前のセルからの続き -->
!rm myfile.txt
!ls -l
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
!head -n5 t # ここで tab を押すと補完候補が出る: Jupyter 上でも出る。

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

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

実行結果は次の通り.

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

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

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

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

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

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

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

実行結果は次の通り.

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

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

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

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

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

出力結果は次の通り.

Overwriting filelist.py

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

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

そして実行結果.

['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 で使うシンボルも定義している. シンボルといった言葉も追々説明していく. 要は方程式の未知変数みたいなものだと思ってほしい.

from sympy import *
from IPython.display import display
init_printing(use_unicode=True)

 * シンボル定義
x = Symbol('x')
y = Symbol('y')
サンプル
 # 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

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

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

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

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

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

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

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

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

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

$$1$$

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

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

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

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

行列の固有値: 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 }$$

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

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

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

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

数値の計算

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

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$$

無限大

$$\infty$$

$$\mathrm{True}$$

$$\infty$$

展開

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

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 )}$$

多項式

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

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}$$

多項式を展開してくれる: (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}$$

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

$$x + 1$$

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

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

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

$$2 x$$

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

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

因数分解: 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)$$

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

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

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

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

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

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

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

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

無理数解: 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 ]$$

多変数多項式を解く: 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) をしていると思ってほしい.

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

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

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

$$6$$

方程式をある変数で解く

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

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 }$$

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

$$\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}$$

論理式

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

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

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

{y: True, x: True}

False

連立方程式を解く

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

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 ]$$

線型代数

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

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]$$

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]$$

微分積分

極限

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

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 よりも中高数学をやるのに優れた点.

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$$

高階微分

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

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

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

第 3 引数で指定する

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

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

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

微分方程式

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

 # 未定義の関数を定義
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'))

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

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)}$$

テイラー展開

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

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 に軍配が上がるだろう.

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

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

積分
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$$

特殊関数もいける

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

定積分

$$0$$

$$1$$

$$2$$

広義積分

$$1$$

$$\sqrt{\pi}$$

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

積分の厳密解

$$\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詰め, パディング

+BEGIN_SRC python :session :results output

for i in range(32): padded = str('{0:02d}'.format(i)) tmpdate = "2017-08-" + padded print(tmpdate) print(tmpdate)

+END_SRC

+RESULTS:

+begin_example

... ... ... ... 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

+end_example

datetime, 日付

日付の宣言

現在の日付

import datetime

print(datetime.date.today()) # 2021-05-15
print(datetime.datetime.today()) # 2021-05-15 23:58:55.230456

任意の日付

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

年月日・時分秒

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

日付だけ取る

import datetime

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

時間だけ

import datetime

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

日付計算

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
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比較できない.

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 を使います. 第一引数は文字列の日付で, 第二引数は日付フォーマットです.

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

日付型から文字列

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はいらない.

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, 例外

try:
    pass # some process
except Exception as e:
    print(e)
    pass # some process

file/directory

pathlib

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

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

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

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

yaml

基本

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

pip install pyyaml

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

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

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

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

ディレクトリ生成

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

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

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

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

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

import os
os.path.isdir(fname)

ファイル移動

import shutil
shutil.move(from, to)

ファイルかどうかを判定

import os
os.path.isfile(fname)

ファイルコピー

import shutil
shutil.copyfile(from, to)

ファイル削除

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

ファイル名取得

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

ファイル読み書き

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

#読むとき
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に変換する

pip install jupytext
jupytext --to notebook --execute sample.py
jupytext --to notebook --execute arithmetic001.py

Jupyter nbconvert ipynbpyに変換

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

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

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

Jupyter nbconvertのテンプレート取得

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

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

import sys
sys.executable
jupyter kernelspec list

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

Jupyter カーネルの確認

jupyter kernelspec list

Jupyter カーネルの削除

jupyter kernelspec remove kernel_name

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

jupyter --version

list, リスト

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

小文字だけ

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に記述があります.

>>> 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()したあとのオブジェクトを取る必要があります.

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

a.reverse()
print(a)

リストを文字列にする

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

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

リストを連結したい

+を使います.

Matplotlibでquiver

コードサンプル 1

+BEGIN_SRC ipython :session :file tmp/quiver1.tmp.png

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

+END_SRC

+RESULTS:

[[file:tmp/quiver1.tmp.png]]

コードサンプル 2

図 1

+BEGIN_SRC ipython :session :file tmp/quiver2.tmp.png

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.05dx, r + 0.05dx, b - 0.05dy, t + 0.05dy])

plt.title('picture 1')

+END_SRC

+RESULTS:

[[file:tmp/quiver2.tmp.png]]

図 2

+BEGIN_SRC ipython :session :file tmp/quiver3.tmp.png

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

+END_SRC

+RESULTS:

[[file:tmp/quiver3.tmp.png]]

図 3

+BEGIN_SRC ipython :session :file tmp/quiver4.tmp.png

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

+END_SRC

+RESULTS:

[[file:tmp/quiver4.tmp.png]]

図 4

+BEGIN_SRC ipython :session :file tmp/quiver5.tmp.png

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

+END_SRC

+RESULTS:

[[file:tmp/quiver5.tmp.png]]

図 5

+BEGIN_SRC ipython :session :file tmp/quiver6.tmp.png

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

+END_SRC

+RESULTS:

[[file:tmp/quiver6.tmp.png]]

図 6

+BEGIN_SRC ipython :session :file tmp/quiver7.tmp.png

plt.figure() M = np.zeros(U.shape, dtype='bool') M[U.shape[0]//3:2U.shape[0]//3, U.shape[1]//3:2U.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()

+END_SRC

+RESULTS:

[[file:tmp/quiver7.tmp.png]]

Matplotlibでグラフの例

参考情報

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

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

単純なグラフ

+BEGIN_SRC ipython :session :file tmp/graph_example001_tmp.png

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

+END_SRC

+RESULTS:

[[file:tmp/graph_example001_tmp.png]]

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

+BEGIN_SRC ipython :session :file tmp/graph_example002_tmp.png

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

+END_SRC

+RESULTS:

[[file:tmp/graph_example002_tmp.png]]

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

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

+BEGIN_SRC ipython :session :file tmp/graph_example003_tmp.png

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

+END_SRC

+RESULTS:

[[file:tmp/graph_example003_tmp.png]]

直線を引いてみる

+BEGIN_SRC ipython :session :file tmp/graph_example004_tmp.png

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

+END_SRC

+RESULTS:

[[file:tmp/graph_example004_tmp.png]]

適当な折れ線を引く

+BEGIN_SRC ipython :session :file tmp/graph_example005_tmp.png

%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')

+END_SRC

+RESULTS:

[[file:tmp/graph_example005_tmp.png]]

横軸を適当に設定

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

+BEGIN_SRC ipython :session :file tmp/graph_example006_tmp.png

%matplotlib inline import matplotlib.pyplot as plt

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

+END_SRC

+RESULTS:

[[file:tmp/graph_example006_tmp.png]]

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

+BEGIN_SRC ipython :session :file tmp/graph_example007_tmp.png

%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'])

+END_SRC

+RESULTS:

[[file:tmp/graph_example007_tmp.png]]

グラフに色々追加

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

+BEGIN_SRC ipython :session :file tmp/graph_example008_tmp.png

%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'])

+END_SRC

+RESULTS:

[[file:tmp/graph_example008_tmp.png]]

軸の表示を変える 1

+BEGIN_SRC ipython :session :file tmp/graph_example009_tmp.png

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

+END_SRC

+RESULTS:

[[file:tmp/graph_example009_tmp.png]]

軸の表示を変える 2

+BEGIN_SRC ipython :session :file tmp/graph_example010_tmp.png

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

+END_SRC

+RESULTS:

[[file:tmp/graph_example010_tmp.png]]

縦横比を合わせる

+BEGIN_SRC ipython :session :file tmp/graph_example011_tmp.png

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

+END_SRC

+RESULTS:

[[file:tmp/graph_example011_tmp.png]]

ファイルとして保存

+BEGIN_SRC ipython :session :results silent

%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')

+END_SRC

円を描く

+BEGIN_SRC ipython :session :file tmp/graph_example013_tmp.png

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

+END_SRC

+RESULTS:

[[file:tmp/graph_example013_tmp.png]]

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

+BEGIN_SRC ipython :session :file tmp/graph_example014_tmp.png

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

+END_SRC

+RESULTS:

[[file:tmp/graph_example014_tmp.png]]

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

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

+BEGIN_SRC ipython :session :file tmp/graph_example035_tmp.png

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

+END_SRC

+RESULTS:

[[file:tmp/graph_example035_tmp.png]]

物体の投射 (放物線)

+BEGIN_SRC ipython :session :file tmp/graph_example015_tmp.png

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

+END_SRC

+RESULTS:

[[file:tmp/graph_example015_tmp.png]]

複数の放物線を重ねる

+BEGIN_SRC ipython :session :file tmp/graph_example016_tmp.png

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

+END_SRC

+RESULTS:

[[file:tmp/graph_example016_tmp.png]]

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

+BEGIN_SRC ipython :session :file tmp/graph_example017_tmp.png

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

+END_SRC

+RESULTS:

[[file:tmp/graph_example017_tmp.png]]

Numpy も使って三角関数

+BEGIN_SRC ipython :session :file tmp/graph_example018_tmp.png

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

+END_SRC

+RESULTS:

[[file:tmp/graph_example018_tmp.png]]

Numpy も使って散布図

+BEGIN_SRC ipython :session :file tmp/graph_example019_tmp.png

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

+END_SRC

+RESULTS:

[[file:tmp/graph_example019_tmp.png]]

Numpy も使って色を変更

+BEGIN_SRC ipython :session :file tmp/graph_example020_tmp.png

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

+END_SRC

+RESULTS:

[[file:tmp/graph_example020_tmp.png]]

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

+BEGIN_SRC ipython :session :file tmp/graph_example021_tmp.png

%matplotlib inline import numpy as np import matplotlib.pyplot as plt

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

+END_SRC

+RESULTS:

[[file:tmp/graph_example021_tmp.png]]

Numpy も使って 100x100 行列 1

+BEGIN_SRC ipython :session :file tmp/graph_example022_tmp.png

%matplotlib inline import numpy as np import matplotlib.pyplot as plt

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

+END_SRC

+RESULTS:

[[file:tmp/graph_example022_tmp.png]]

Numpy も使って 100x100 行列 2

+BEGIN_SRC ipython :session :file tmp/graph_example023_tmp.png

%matplotlib inline import numpy as np import matplotlib.pyplot as plt

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

+END_SRC

+RESULTS:

[[file:tmp/graph_example023_tmp.png]]

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

+BEGIN_SRC ipython :session :file tmp/graph_example024_tmp.png

%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の範囲に限定

+END_SRC

+RESULTS:

[[file:tmp/graph_example024_tmp.png]]

破線を入れる

+BEGIN_SRC ipython :session :file tmp/graph_example025_tmp.png

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

+END_SRC

+RESULTS:

[[file:tmp/graph_example025_tmp.png]]

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

+BEGIN_SRC ipython :session :file tmp/graph_example026_tmp.png

%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() # タイトルの被りを防ぐ

+END_SRC

+RESULTS:

[[file:tmp/graph_example026_tmp.png]]

プロットの順番

*plt.subplot() を使う.

plt.subplot(行数, 列数, 何番目のプロットか)

+BEGIN_SRC ipython :session :file tmp/graph_example027_tmp.png

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

+END_SRC

+RESULTS:

[[file:tmp/graph_example027_tmp.png]]

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

+BEGIN_SRC ipython :session :file tmp/graph_example028_tmp.png

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() # タイトルとラベルが被るのを解消

+END_SRC

+RESULTS:

[[file:tmp/graph_example028_tmp.png]]

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

普通に描くと縦軸ずれる

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

+BEGIN_SRC ipython :session :file tmp/graph_example029_tmp.png

%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, ".")

+END_SRC

+RESULTS:

[[file:tmp/graph_example029_tmp.png]]

sharex, sharey

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

+BEGIN_SRC ipython :session :file tmp/graph_example030_tmp.png

%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, ".")

+END_SRC

+RESULTS:

[[file:tmp/graph_example030_tmp.png]]

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

+BEGIN_SRC ipython :session :file tmp/graph_example031_tmp.png

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)

+END_SRC

+RESULTS:

[[file:tmp/graph_example031_tmp.png]]

凡例

+BEGIN_SRC ipython :session :file tmp/graph_example032_tmp.png

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

+END_SRC

+RESULTS:

[[file:tmp/graph_example032_tmp.png]]

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

+BEGIN_SRC ipython :session :file tmp/graph_example033_tmp.png

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

+END_SRC

+RESULTS:

[[file:tmp/graph_example033_tmp.png]]

凡例の位置調節

+BEGIN_SRC ipython :session :file tmp/graph_example034_tmp.png

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

+END_SRC

+RESULTS:

[[file:tmp/graph_example034_tmp.png]]

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

参考情報

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

+BEGIN_SRC ipython :session :file tmp/graph_example_3d_001_tmp.png

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

+END_SRC

+RESULTS:

[[file:tmp/graph_example_3d_001_tmp.png]]

その 2: ワイヤフレーム

+BEGIN_SRC ipython :session :file tmp/graph_example_3d_002_tmp.png

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

+END_SRC

+RESULTS:

[[file:tmp/graph_example_3d_002_tmp.png]]

サーフェス表示

+BEGIN_SRC ipython :session :file tmp/graph_example_3d_003_tmp.png

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

+END_SRC

+RESULTS:

[[file:tmp/graph_example_3d_003_tmp.png]]

3 次元プロット

+BEGIN_SRC ipython :session :file tmp/graph_example_3d_004_tmp.png

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

+END_SRC

+RESULTS:

[[file:tmp/graph_example_3d_004_tmp.png]]

等高線 1

+BEGIN_SRC ipython :session :file tmp/graph_example_3d_005_tmp.png

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

+END_SRC

+RESULTS:

[[file:tmp/graph_example_3d_005_tmp.png]]

等高線 2

+BEGIN_SRC ipython :session :file tmp/graph_example_3d_006_tmp.png

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

+END_SRC

+RESULTS:

[[file:tmp/graph_example_3d_006_tmp.png]]

Numpyの基本

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

+BEGIN_SRC ipython :session

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

+END_SRC

+RESULTS:

\int \sqrt{\frac{1}{x}}\, dx

基本

+BEGIN_SRC ipython :session :results output

import numpy as np

print(np.pi)

print(str(np.pi))

print('{:10.4f}'.format(np.pi)) print('{:10.8f}'.format(np.pi))

+END_SRC

+RESULTS:

3.141592653589793
3.141592653589793
3.1416
3.14159265

多項式

順番に見ていこう.

多項式の作り方

+BEGIN_SRC ipython :session :results output

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)

+END_SRC

+RESULTS:

配列の 0 番目が最高階の係数となる形式で順次指定されていく
2
1 x + 2 x + 3
第 2 引数を True にすると指定配列を根とする多項式になる
2
1 x - 1 x

多項式の係数

+BEGIN_SRC ipython :session :results output

print(p1) print('多項式の係数の配列') print(p1.c) print('多項式の k 番目の係数') print(p1[1])

+END_SRC

+RESULTS:

2
1 x + 2 x + 3
多項式の係数の配列
[1 2 3]
多項式の k 番目の係数
2

多項式の次数

+BEGIN_SRC ipython :session :results output

print(p1) print(p1.order)

+END_SRC

+RESULTS:

2
1 x + 2 x + 3
多項式の次数
2

多項式の評価

+BEGIN_SRC ipython :session :results output

print(p1) print(np.polyval(p1, 5)) print(p1(5))

+END_SRC

+RESULTS:

2
1 x + 2 x + 3
38
38

多項式の演算

和, 差, 積.

+BEGIN_SRC ipython :session :results output

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)

+END_SRC

べきと割り算関係.

+BEGIN_SRC ipython :session :results output

print('多項式のべき') print(p1 ** 2)

q = np.polydiv(p2, p1) print('多項式の割り算') print(q) print('多項式の商') print(q[0]) print('多項式の割り算:あまり') print(q[1])

+END_SRC

+RESULTS:

+begin_example

多項式のべき 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

+end_example

多項式の微積分

+BEGIN_SRC ipython :session :results output

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

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

+END_SRC

多項式の根

+BEGIN_SRC ipython :session :results output

print('多項式の根') print(np.sort(np.roots(p2))) print(np.sort(p2.r)) print(np.sort(p3.r))

+END_SRC

+RESULTS:

多項式の根
[-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.]

配列

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

+BEGIN_SRC ipython :session :results output

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

+END_SRC

+RESULTS:

[ 1. 233. 21. 234. 2.]
[ 2. 234. 21. 235. 3.]
[ 2. 233. 21. 234. 3.]

関数

階乗

+BEGIN_SRC ipython :session :results output

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 だと駄目らしい

+END_SRC

+RESULTS:

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]

指数

+BEGIN_SRC ipython :session :results output

exponents = np.arange(0, 10, 0.1) powers = np.power(np.e, exponents) print(powers)

+END_SRC

+RESULTS:

+begin_example

[ 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]

+end_example

対数

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

+BEGIN_SRC ipython :session :results output

args = np.arange(0.1, 10, 0.1) print(np.log(np.e)) print(np.log10(10)) print(np.log2(2))

+END_SRC

+RESULTS:

1.0
1.0
1.0

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

+BEGIN_SRC ipython :session :results output

print(np.log(args)) print(np.log10(args)) print(np.log2(args))

+END_SRC

+RESULTS:

+begin_example

[-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]

+end_example

その他

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

度数とラジアンの変換

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

+BEGIN_SRC ipython :session :results output

degrees = np.arange(0, 181, 0.1) print(degrees)

+END_SRC

+RESULTS:

[ 0.00000000e+00 1.00000000e-01 2.00000000e-01 ..., 1.80700000e+02
1.80800000e+02 1.80900000e+02]

度数をラジアンに変換.

+BEGIN_SRC ipython :session :results output

radians = np.deg2rad(degrees) print(radians)

+END_SRC

+RESULTS:

[ 0.00000000e+00 1.74532925e-03 3.49065850e-03 ..., 3.15380996e+00
3.15555529e+00 3.15730062e+00]

ラジアンを度数に変換.

+BEGIN_SRC ipython :session :results output

print(np.rad2deg(radians))

+END_SRC

+RESULTS:

[ 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}

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

オイラー法

+BEGIN_SRC ipython :session :file tmp/ode1.tmp.png

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

+END_SRC

+RESULTS:

[[file:tmp/ode1.tmp.png]]

ルンゲ-クッタ法

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

とりあえず一般的な次の次の初期値問題を考えよう. \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}

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

+BEGIN_SRC ipython :session :file tmp/ode2.tmp.png

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

+END_SRC

+RESULTS:

[[file:tmp/ode2.tmp.png]]

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

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

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

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

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

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

+BEGIN_SRC ipython :session :file tmp/ode4.tmp/png

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 [-pv[0]+pv[1], -v[0]v[2]+rv[0]-v[1], v[0]v[1]-bv[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])

+END_SRC

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}

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

オイラー法のコード

+BEGIN_SRC ipython :session :file tmp/ode5.tmp.png

%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")

+END_SRC

+RESULTS:

[[file: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}

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

+BEGIN_SRC ipython :session :file tmp/ode6.tmp.png

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

+END_SRC

+RESULTS:

[[file:tmp/ode6.tmp.png]]

NumPyで確率論

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

+BEGIN_SRC ipython :session :results none

%matplotlib inline import numpy as np import matplotlib.pyplot as plt from numpy.random import randint

+END_SRC

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

+BEGIN_SRC ipython :session :results output

from numpy.random import randint print(randint(1, 7, 10))

+END_SRC

+RESULTS:

[1 2 2 3 2 3 4 4 1 4]

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

+BEGIN_SRC ipython :session :file tmp/probability.tmp.png

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

+END_SRC

+RESULTS:

[[file:tmp/probability.tmp.png]]

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

+BEGIN_SRC ipython :session :file tmp/probability2.tmp.png

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

+END_SRC

+RESULTS:

[[file:tmp/probability2.tmp.png]]

他のモジュール: random, normal

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

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

+BEGIN_SRC ipython :session :results output

%matplotlib inline import numpy as np import matplotlib.pyplot as plt from numpy.random import random from numpy.random import normal

print(random(5))

+END_SRC

+RESULTS:

[ 0.66619195 0.29721425 0.07389146 0.58977243 0.97906712]

100-150 の範囲の乱数

+BEGIN_SRC ipython :session :results output

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

+END_SRC

+RESULTS:

[ 125.5122101 146.17318573 100.94123889 114.62349853 123.51461347]

正規分布の乱数 normal

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

+BEGIN_SRC ipython :session :file tmp/probability3.tmp.png

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

+END_SRC

+RESULTS:

[[file:tmp/probability3.tmp.png]]

大数の法則

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

+BEGIN_SRC ipython :session :file tmp/probability4.tmp.png

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

+END_SRC

+RESULTS:

[[file:tmp/probability4.tmp.png]]

コイン投げ

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

+BEGIN_SRC ipython :session :file tmp/probability5.tmp.png

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

+END_SRC

+RESULTS:

[[file:tmp/probability5.tmp.png]]

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

+BEGIN_SRC ipython :session :file tmp/probability6.tmp.png

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)

+END_SRC

+RESULTS:

[[file:tmp/probability6.tmp.png]]

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

ちょっと遊んでみよう.

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

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

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

+BEGIN_SRC ipython :session :file tmp/probability7.tmp.png

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

+END_SRC

+RESULTS:

[[file:tmp/probability7.tmp.png]]

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

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

+BEGIN_SRC ipython :session :file tmp/probability8.tmp.png

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

+END_SRC

+RESULTS:

[[file:tmp/probability8.tmp.png]]

モンテカルロ法

参考

一様乱数

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

+BEGIN_SRC ipython :session :results output

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)

+END_SRC

+RESULTS:

+begin_example

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]

+end_example

モンテカルロ法

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

+BEGIN_SRC ipython :session :results silent

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

+END_SRC

試行回数 100

+BEGIN_SRC ipython :session :file tmp/probability9.tmp.png

calculate_pi(100)

+END_SRC

+RESULTS:

[[file:tmp/probability9.tmp.png]]

試行回数 1000

+BEGIN_SRC ipython :session :file tmp/probability10.tmp.png

calculate_pi(1000)

+END_SRC

+RESULTS:

[[file:tmp/probability10.tmp.png]]

試行回数 10000

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

+BEGIN_SRC ipython :session :file tmp/probability11.tmp.png

calculate_pi(10000)

+END_SRC

+RESULTS:

[[file:tmp/probability11.tmp.png]]

NumPyで統計学

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

平均 (mean)

+BEGIN_SRC ipython :session :results output

meanlst = [1,2,3]

lstsum = sum(meanlst) print(lstsum)

lengthlst = len(meanlst) print(lengthlst)

mean = lstsum / lengthlst print(mean)

+END_SRC

+RESULTS:

6
3
2.0

中央値

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

ソート

+BEGIN_SRC ipython :session :results output

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

print(medianlst)

+END_SRC

+RESULTS:

[1, 3, 4]

中央値を計算する

+BEGIN_SRC ipython :session :results output

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

+END_SRC

+RESULTS:

要素数 12 で中央値は 514.0.

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

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

+BEGIN_SRC ipython :session :results output

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)

+END_SRC

+RESULTS:

[(2, 4), (1, 2), (4, 2), (3, 1), (5, 1), (6, 1)]
[(2, 4)]
[(2, 4), (1, 2)]
[(2, 4)]

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

+BEGIN_SRC ipython :session :results output

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

+END_SRC

+RESULTS:

リストの最頻値は 9

最頻値が複数ある場合を見てみよう.

+BEGIN_SRC ipython :session :results output

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

+END_SRC

+RESULTS:

リストの最頻値は [(4, 3), (5, 3)]

度数分布表

+BEGIN_SRC ipython :session :results output

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)

+END_SRC

+RESULTS:

+begin_example

Number,Frequency 9,5 6,3 1,2 5,2 7,2 8,2 10,2 2,1 4,1

+end_example

何となく別コードを書いてみる.

+BEGIN_SRC ipython :session :results output

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)

+END_SRC

+RESULTS:

+begin_example

[['Number', 'Frequency'], [9, 5], [6, 3], [1, 2], [5, 2], [7, 2], [8, 2], [10, 2], [2, 1], [4, 1]]

+end_example

上の度数分布表をソートしよう.

+BEGIN_SRC ipython :session :results output

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)

+END_SRC

+RESULTS:

+begin_example

Number,Frequency 1,2 2,1 4,1 5,2 6,3 7,2 8,2 9,5 10,2

+end_example

分散

範囲 (range)

+BEGIN_SRC ipython :session :results output

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

+END_SRC

+RESULTS:

最小: 60 最大: 1200 範囲: 1140

偏差 (variance) と標準偏差 (standard deviation)

偏差 $V$ は次式で定義する. \begin{align} V = \frac{\sum_{k=1}^{n} (x_k - m)}{n}. \end{align}

ここで $m$ は平均. 標準偏差は $\sqrt{V}$ で定義する.

+BEGIN_SRC ipython :session :results output

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

+END_SRC

+RESULTS:

リストの分散は 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}

+BEGIN_SRC ipython :session :results output

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

+END_SRC

NumPyで物理

アニメーションもあるのでとりあえず Jupyter 管理だけにしておく.

Pandas

+BEGIN_SRC python :session :results output

message_csv = "メッセージ定義.csv"

df = pd.read_csv(message_csv, encoding="utf-8").fillna("")

last_row_number = len(df.index)

+END_SRC

+BEGIN_SRC python :session :results output

message_csv = "メッセージ定義.csv"

df = pd.read_csv(message_csv, encoding="utf-8").fillna("")

last_row_number = len(df.index)

+END_SRC

+BEGIN_SRC python :session :results output

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["条件"] + " "

+END_SRC

+BEGIN_SRC python :session :results output

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)

+END_SRC

pip

pip 自体のアップデート

pip install -U pip

アップデートが必要なパッケージのリスト

pip list -o

パッケージのバージョン指定インストール

pip install <package-name>==<version>

パッケージをパッケージ名==バージョンで表示

requirements.txtを作るときに便利.

pip freeze > requirements.txt

パッケージのアップデート

pip install -U <package-name>

依存関係の確認

pip check

print

即時表示

flush=True をつけます.

print("出力の即時表示", flush=True)

re, 正規表現

compile

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__":をつけて隔離する必要がある.

初回だけ

sphinx-quickstart docs

> Separate source and build directories (y/n) [n]:

sys.path.insert(0, os.path.abspath('../'))
extensions = [
    'sphinx.ext.autodoc',
    'sphinx.ext.viewcode',
    'sphinx.ext.todo',
]

ドキュメント生成コマンド

sphinx-apidoc -f -o ./docs . & sphinx-build -b singlehtml ./docs ./docs/_build
autodoc_mock_imports = ["some_module_name"]

string, 文字列

いろいろな都合で文字列の項目であってもリストなど他の項目に入っていることもあります.

ファイル名からファイル名と拡張子を取得

os.path.split() でファイル名とフォルダ名のペアが取れます. 拡張子はドットを含むこと, つまり a.txt に対して .txt であることに注意しましょう.

import os
fname = "string.org"
froot, fext = os.path.split()
#froot: string
#fext:  .org

日時の文字列を取る

datetimeから作ります. 適切なバージョンの datetime オブジェクトのドキュメントを見るといいでしょう.

import datetime
dtstr = datetime.datetime.now().strftime('%Y%m%d-%H%M%S')

dt_now = datetime.datetime.now()

文字詰め・パディング

右寄せゼロ埋め zfill

s = "1".zfill(3) # 001

右寄せ・左寄せ・中央寄せ rjust, ljust, center

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

文字の連番を取得する

map(chr, range(ord('B'), ord('R')+1))

全体コードは次の通り.

for char in map(chr, range(ord('A'), ord('C')+1)):
    s = "\\newcommand{\scr%s}{\mathscr{%s}}" % (char, char)
    print(s)
\newcommand{\scrA}{\mathscr{A}}
\newcommand{\scrB}{\mathscr{B}}
\newcommand{\scrC}{\mathscr{C}}

見ての通りTeXのコマンドを一気に定義したかった.

文字列をリストにする

文字列のsplit()メソッドを使います.

ss = "a,b,c".sprit(",")

SymPyことはじめ

参考サイト

プロットのサンプル (関数のグラフを描く) は別項目でまとめる.

サンプル

まずは式絡みの計算を

+BEGIN_SRC ipython :session :results output

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

+END_SRC

+RESULTS:

2x - 3y
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) 表示できる.

+BEGIN_SRC ipython :session :results output

import sympy from IPython.display import display init_printing(use_unicode=True)

display(sympy.sqrt(3))

+END_SRC

+BEGIN_SRC ipython :session :results output

import sympy print(latex(sympy.sqrt(3)))

+END_SRC

+RESULTS:

\sqrt{3}

\begin{gather} \sqrt{3} \end{gather}

勝手に計算もやってくれる. 状況によってはよくない.

+BEGIN_SRC ipython :session :results output

from sympy import symbols x, y = symbols('x y')

print(latex(sympy.sqrt(8)))

+END_SRC

+RESULTS:

2 \sqrt{2}

\begin{align} 2 \sqrt{2} \end{align}

多項式のかけ算

+BEGIN_SRC ipython :session :results output

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

+END_SRC

+RESULTS:

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}

微分

+BEGIN_SRC ipython :session :results output

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

+END_SRC

+RESULTS:

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

+BEGIN_SRC ipython :session :results output

print(latex(exp(x)sin(x) + exp(x)cos(x))) print(latex(integrate(exp(x)sin(x) + exp(x)cos(x), x)))

+END_SRC

+RESULTS:

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

+BEGIN_SRC ipython :session :results output

print(latex(exp(x)sin(x) + exp(x)cos(x))) print(latex(integrate(exp(x)sin(x) + exp(x)cos(x), x)))

+END_SRC

+RESULTS:

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}$ 上で積分する.

+BEGIN_SRC ipython :session :results output

print(latex(integrate(sin(x**2), (x, -oo, oo))))

+END_SRC

+RESULTS:

\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}

+BEGIN_SRC ipython :session :results output

print(latex(limit(sin(x)/x, x, 0)))

+END_SRC

+RESULTS:

1

多項式の根

\begin{align} x^{2} - 2 = 0. \end{align}

+BEGIN_SRC ipython :session :results output

print(latex(solve(x**2 - 2, x)))

+END_SRC

+RESULTS:

\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}

+BEGIN_SRC ipython :session :results output

f = Function('f') print(latex(dsolve(Eq(f(t).diff(t, t) - f(t), exp(t)), f(t))))

+END_SRC

+RESULTS:

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}

+BEGIN_SRC ipython :session :results output

print(latex(Matrix([[1, 2], [2, 2]]).eigenvals()))

+END_SRC

+RESULTS:

\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}

ベッセル関数を球ベッセル関数で書く

+BEGIN_SRC ipython :session :results output

print(latex(besselj(nu, z).rewrite(jn)))

+END_SRC

+RESULTS:

\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 コードを出力する

+BEGIN_SRC ipython :session :results output

print(latex(Integral(cos(x)**2, (x, 0, pi))))

+END_SRC

+RESULTS:

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

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

数値の計算

円周率

+BEGIN_SRC ipython :session :results output

display(pi**2) display(pi.evalf()) display((pi + exp(1)).evalf())

+END_SRC

$\sqrt{2}$ を 100 桁まで評価

+BEGIN_SRC ipython :session :results output

print(latex(sqrt(2).evalf(100)))

+END_SRC

+RESULTS:

1.414213562373095048801688724209698078569671875376948073176679737990732478462107038850387534327641573

無限大

+BEGIN_SRC ipython :session :results output

print(latex(oo)) print(latex(oo > 99999)) print(latex(oo + 1))

+END_SRC

+RESULTS:

\infty
\mathrm{True}
\infty

展開

+BEGIN_SRC ipython :session :results output

print(latex(expand(x + y, complex=True))) print(latex(Iim(x) + Iim(y) + re(x) + re(y))) print(latex(expand(cos(x + y), trig=True))) print(latex(cos(x)cos(y) - sin(x)sin(y)))

+END_SRC

+RESULTS:

\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 との連携

次のコードはエラーになる.

+BEGIN_SRC ipython :session :results

from sympy import Symbol import math theta = Symbol('theta') math.sin(theta) + math.sin(theta)

+END_SRC

次のようにすると正しく処理できる.

+BEGIN_SRC ipython :session :results output

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

+END_SRC

+RESULTS:

1.0
2 \sin{\left (\theta \right )}

\begin{gather} 1.0 \ 2 \sin{\left (\theta \right )} \end{gather}

変数とシンボルのアルファベットは同じでなくてもいい

+BEGIN_SRC ipython :session :results output

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

+END_SRC

+RESULTS:

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}

複数のシンボルを一気に定義できる

+BEGIN_SRC ipython :session :results output

from sympy import symbols x,y,z = symbols('x,y,z') s = x * y + x * y print(latex(s))

+END_SRC

+RESULTS:

2 x y

\begin{gather} 2 x y \end{gather}

簡単な計算なら自動処理してくれる

+BEGIN_SRC ipython :session :results output

p = x * (x + x) print(latex(p))

+END_SRC

+RESULTS:

2 x^{2}

\begin{gather} 2 x^{2} \end{gather}

少し複雑になると展開などの処理はしない

+BEGIN_SRC ipython :session :results output

p = (x + 2)*(x + 3) print(latex(p))

+END_SRC

+RESULTS:

\left(x + 2\right) \left(x + 3\right)

\begin{gather} \left(x + 2\right) \left(x + 3\right) \end{gather}

代入

+BEGIN_SRC ipython :session :results output

from sympy import Symbol x = Symbol('x') y = Symbol('y')

expr = xx + xy + xy + yy print(latex(expr)) res = expr.subs({x:1, y:2}) print(latex(res))

print(latex(expr.subs({x:1-y})))

+END_SRC

+RESULTS:

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}

文字列を式に変換

+BEGIN_SRC ipython :session :results output

from sympy import sympify init_printing(order='lex') expr = "x2 + 3*x + x3 + 2*x" expr = sympify(expr) print(latex(expr))

print(latex(2 * expr))

+END_SRC

+RESULTS:

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 次方程式

+BEGIN_SRC ipython :session :results output

from sympy import Symbol, solve, latex x = Symbol('x') expr = x - 5 - 7 sol = solve(expr) print(latex(sol))

+END_SRC

+RESULTS:

\left [ 12\right ]

\begin{align} \left [ 12\right ] \end{align}

連立 1 次方程式

+BEGIN_SRC ipython :session :results output

x = Symbol('x') y = Symbol('y') expr1 = 2x + 3y - 6 expr2 = 3x + 2y - 12 print(latex(solve((expr1, expr2), dict=True)))

+END_SRC

+RESULTS:

\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 次方程式

+BEGIN_SRC ipython :session :results output

from sympy import solve, latex x = Symbol('x') expr = x*2 + 5x + 4 print(latex(solve(expr, dict=True)))

x = Symbol('x') expr = x**2 + x + 1 print(latex(solve(expr, dict=True)))

+END_SRC

+RESULTS:

\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}

不等式を解く

+BEGIN_SRC ipython :session :results output

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

+END_SRC

+RESULTS:

\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}

物体の投射の式も解ける

+BEGIN_SRC ipython :session :results output

from sympy import sin, solve, Symbol u = Symbol('u') t = Symbol('t') g = Symbol('g') theta = Symbol('theta') print(latex(solve(usin(theta)-gt, t)))

+END_SRC

+RESULTS:

\left [ \frac{u}{g} \sin{\left (\theta \right )}\right ]

\begin{gather} \left [ \frac{u}{g} \sin{\left (\theta \right )}\right ] \end{gather}

条件指定

+BEGIN_SRC ipython :session :results output

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

+END_SRC

+RESULTS:

Do Something
Do Something

極限

まずはこれ。

\begin{align} \lim_{x \to \infty} \frac{1}{x} = 0. \end{align}

+BEGIN_SRC ipython :session :results output

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

+END_SRC

+RESULTS:

\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}

不定形

+BEGIN_SRC ipython :session :results output

from sympy import Symbol, sin print(latex(Limit(sin(x) / x, x, 0).doit()))

+END_SRC

+RESULTS:

1

\begin{gather} 1 \end{gather}

自然対数の底 $e$

+BEGIN_SRC ipython :session :results output

from sympy import Limit, Symbol, S n = Symbol('n') print(latex(Limit((1 + 1 / n) ** n, n, S.Infinity).doit()))

+END_SRC

+RESULTS:

e

\begin{gather} e \end{gather}

次の式を評価してみよう.

\begin{align} A = P (1 + \frac{r}{n})^{nt}. \end{align}

+BEGIN_SRC ipython :session :results output

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

+END_SRC

+RESULTS:

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}

+BEGIN_SRC ipython :session :results output

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

+END_SRC

+RESULTS:

10 t_{1} + 2

\begin{gather} 10 t_{1} + 2 \end{gather}

微分: 微分するメソッド

+BEGIN_SRC ipython :session :results output

from sympy import Symbol, Derivative t = Symbol('t') St = 5t2 + 2t + 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})))

+END_SRC

+RESULTS:

\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}

導関数を求めるプログラム

+BEGIN_SRC ipython :session :results output

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)

+END_SRC

+RESULTS:

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}

高階の導関数や極大極小

元の関数と導関数.

+BEGIN_SRC ipython :session :results output

from sympy import Symbol, solve, Derivative, pprint x = Symbol('x') f = x5 - 30*x3 + 50*x print(latex(f)) d1 = Derivative(f, x).doit() print(latex(d1))

+END_SRC

+RESULTS:

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}

臨界点.

+BEGIN_SRC ipython :session :results output

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

+END_SRC

+RESULTS:

\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}

極値の評価.

+BEGIN_SRC ipython :session :results output

print(d2.subs({x:B}).evalf()) print(d2.subs({x:C}).evalf()) print(d2.subs({x:A}).evalf()) print(d2.subs({x:D}).evalf())

+END_SRC

+RESULTS:

127.661060789073
-127.661060789073
-703.493179468151
703.493179468151

端点との比較.

+BEGIN_SRC ipython :session :results output

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

+END_SRC

+RESULTS:

705.959460380365
25.0846626340294
375.000000000000
-375.000000000000

最急降下法

+BEGIN_SRC ipython :session :results output

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

+END_SRC

+RESULTS:

Theta: 44.997815081691805
Maximum Range: 63.7755100185965

一般的な最急降下法

*それなりに実行時間がかかるので注意しよう.

+BEGIN_SRC ipython :session :results output

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_sizef1x.subs({x:x_old}).evalf() while abs(x_old - x_new) > epsilon: x_old = x_new x_new = x_old + step_sizef1x.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})))

+END_SRC

+RESULTS:

y: 0.00999904956272288
Maximum value: 0.999950009920428

積分

不定積分.

+BEGIN_SRC ipython :session :results output

from sympy import Integral, Symbol x = Symbol('x') k = Symbol('k')

i = Integral(k*x, x) print(latex(i.doit()))

+END_SRC

+RESULTS:

\frac{k x^{2}}{2}

\begin{gather} \frac{k x^{2}}{2} \end{gather}

定積分.

+BEGIN_SRC ipython :session :results output

d = Integral(k*x, (x, 0, 2)) print(latex(d.doit())) print(latex(Integral(x, (x, 2, 4)).doit()))

+END_SRC

+RESULTS:

2 k
6

\begin{gather} 2 k, \quad 6 \end{gather}

誤差関数を得る積分計算: 実行に少し時間がかかる.

+BEGIN_SRC ipython :session :results output

from sympy import Symbol, exp, sqrt, pi, Integral x = Symbol('x') p = exp(-(x - 10)*2/2)/sqrt(2pi)

print(latex(Integral(p, (x, 11, 12)).doit()))

+END_SRC

+RESULTS:

  • \frac{1}{2} \operatorname{erf}{\left (\frac{\sqrt{2}}{2} \right )} + \frac{1}{2} \operatorname{erf}{\left (\sqrt{2} \right )}

\begin{gather} - \frac{1}{2} \operatorname{erf}{\left (\frac{\sqrt{2}}{2} \right )} + \frac{1}{2} \operatorname{erf}{\left (\sqrt{2} \right )} \end{gather}

関数の値を具体的に評価: 実行に少し時間がかかる.

+BEGIN_SRC ipython :session :results output

print(latex(Integral(p, (x, 11, 12)).doit().evalf()))

+END_SRC

+RESULTS:

0.135905121983278

$\bbR$ 全体での積分.

+BEGIN_SRC ipython :session :results output

from sympy import Symbol, exp, sqrt, pi, Integral, S x = Symbol('x') p = exp(-(x - 10) ** 2 / 2) / sqrt(2 * pi) print(latex(Integral(p, (x, S.NegativeInfinity, S.Infinity)).doit().evalf()))

+END_SRC

+RESULTS:

1.0

SymPy 多項式

基本的な計算や書き方

+BEGIN_SRC ipython :session :results output

from sympy import Symbol, latex x = Symbol('x') print(latex(x + x +1))

+END_SRC

+RESULTS:

2 x + 1

\begin{gather} 2 x + 1 \end{gather}

+BEGIN_SRC ipython :session :results output

from sympy import Symbol x = Symbol('x') print(x + x + 1) print(x * x + 2 * x + 1) print(x3 + 3 * x2 + 3x + 1) print((x+1)3)

+END_SRC

+RESULTS:

2*x + 1
x*2 + 2x + 1
3x + x3 + 3x*2 + 1
(x + 1)**3

\begin{gather} 2x + 1 \ x2 + 2x + 1 \ 3x + x3 + 3x2 + 1 \ (x + 1)*3 \end{gather}

factor, expand

+BEGIN_SRC ipython :session :results output

from sympy import Symbol from sympy import factor, expand x = Symbol('x') y = Symbol('y')

expr = x2 - y2 f_tmp = factor(expr) print(latex(f_tmp))

e_tmp = expand(f_tmp) print(latex(e_tmp))

expr = x3 + 3*x2y + 3xy2 + y*3 factors = factor(expr) print(latex(factors)) print(latex(expand(factors)))

+END_SRC

+RESULTS:

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

\begin{gather} \left(x - y\right) \left(x + y\right) \ x^{2} - y^{2} \ \left(x + y\right)^{3} \ x^{3} + 3 x^{2} y + 3 x y^{2} + y^{3} \end{gather}

prettyprint

http://docs.sympy.org/dev/tutorial/printing.html を見る限り, tex が入っていると tex のタイプセットで出力してくれる模様.

org-babel だと prettyprint が文字化けした. org-babel だと print(latex()) の方がよさそう.

+BEGIN_SRC ipython :session :results silent

from sympy import pprint from sympy import Symbol from sympy import factor, expand x = Symbol('x') y = Symbol('y')

expr = xx + 2xy + yy print(expr) pprint(expr)

expr = 1 + 2x + 2x**2 pprint(expr)

+END_SRC

多項式の順序反転

ここを見るとオプションがいろいろ書いてある.

+BEGIN_SRC ipython :session :results output

from sympy import init_printing init_printing(order='rev-lex') print(latex(expr))

+END_SRC

+RESULTS:

1 + 2 x + 2 x^{2}

\begin{gather} 1 + 2 x + 2 x^{2} \end{gather}

適当な多項式の表示

+BEGIN_SRC ipython :session :results output

from sympy import Symbol, pprint, init_printing def get_some_series(n): """n 次の適当な多項式"" x = Symbol('x') series = x for i in range(2, n + 1): series = series + (x ** i) / i

return series

n = 10 series = get_some_series(int(n))

init_printing(order='rev-lex') print(latex(series))

+END_SRC

+RESULTS:

x + \frac{x^{2}}{2} + \frac{x^{3}}{3} + \frac{x^{4}}{4} + \frac{x^{5}}{5} + \frac{x^{6}}{6} + \frac{x^{7}}{7} + \frac{x^{8}}{8} + \frac{x^{9}}{9} + \frac{x^{10}}{10}

\begin{gather} x + \frac{x^{2}}{2} + \frac{x^{3}}{3} + \frac{x^{4}}{4} + \frac{x^{5}}{5} + \frac{x^{6}}{6} + \frac{x^{7}}{7} + \frac{x^{8}}{8} + \frac{x^{9}}{9} + \frac{x^{10}}{10} \end{gather}

簡素に表示する

+BEGIN_SRC ipython :session :results output

from sympy import Symbol, simplify x = Symbol('x') y = Symbol('y')

expr = xx + xy + xy + yy

expr_subs = expr.subs({x:1-y}) print(latex(expr_subs)) print(latex(simplify(expr_subs)))

+END_SRC

+RESULTS:

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

\begin{gather} \left(1 - y\right)^{2} + 2 y \left(1 - y\right) + y^{2} \ 1 \end{gather}

多項式の評価

+BEGIN_SRC ipython :session :results output

from sympy import Symbol, pprint, init_printing def get_some_series(n, x_value): """ある有限級数の表示""" x = Symbol('x') series = x for i in range(2, n + 1): series = series + (x ** i) / i

return series

init_printing(order='rev-lex')

n = 5

x_value = 3

series = get_some_series(n, x_value) print(latex(series)) series_value = series.subs({x:x_value}) print('x = {0} での級数の値: {1}'.format(x_value, series_value))

+END_SRC

+RESULTS:

x + \frac{x^{2}}{2} + \frac{x^{3}}{3} + \frac{x^{4}}{4} + \frac{x^{5}}{5}
x = 3 での級数の値: 1707/20

\begin{gather} x + \frac{x^{2}}{2} + \frac{x^{3}}{3} + \frac{x^{4}}{4} + \frac{x^{5}}{5} \end{gather}

SymPy とプロット: 関数のグラフを描く

裏では matplotlib を使っているらしい.

グラフ 1

まずやってみるのは次の関数.

\begin{align} f(x) = x^5 - 30 x^3 + 50 x. \end{align}

定義域は適切に制限する.

+BEGIN_SRC ipython :session :file tmp/sympy_fundamental1.tmp.png

%matplotlib inline from sympy.plotting import plot from sympy import Symbol x = Symbol('x') plot(x 5 - 30 * x 3 + 50 * x, (x, -5, 5), title='A curve', xlabel='x', ylabel='y')

+END_SRC

+RESULTS:

[[file:tmp/sympy_fundamental1.tmp.png]]

グラフ 2

+BEGIN_SRC ipython :session :file tmp/sympy_fundamental2.tmp.png

%matplotlib inline from sympy.plotting import plot from sympy import Symbol x = Symbol('x') plot(2 * x + 3)

+END_SRC

+RESULTS:

[[file:tmp/sympy_fundamental2.tmp.png]]

定義域を制限してみる

+BEGIN_SRC ipython :session :file tmp/sympy_fundamental3.tmp.png

%matplotlib inline from sympy.plotting import plot from sympy import Symbol x = Symbol('x') plot(2 * x + 3, (x, -5, 5), title='A Line', xlabel='x', ylabel='2x+3')

+END_SRC

+RESULTS:

[[file:tmp/sympy_fundamental3.tmp.png]]

陰関数の表示 1

2 変数の場合に限るようだ. 計算のためにそれなりに時間がかかる.

+BEGIN_SRC ipython :session :file tmp/sympy_fundamental4.tmp.png

%matplotlib inline from sympy import Symbol, sympify, solve from sympy.plotting import plot_implicit

expr = "x2 + y3 + 2" expr = sympify(expr) plot_implicit(expr)

+END_SRC

+RESULTS:

[[file:tmp/sympy_fundamental4.tmp.png]]

陰関数の表示 2

計算のためにそれなりに時間がかかる. 太く表示されてしまうのは特異点まわりっぽいのでそこでの計算が難しいからだろうか. と思ったが特異性うんぬん関係なさそうなところでも線の太い細いがある. 何だろう?

+BEGIN_SRC ipython :session :file tmp/sympy_fundamental5.tmp.png

from sympy import symbols from sympy.plotting import plot_implicit

x, y = symbols("x y") f = x 6 + 3 * x 4 * y 2 + 6 * x 4 * y - 2 * x 4 + 3 * x 2 * y 4 \ - 2 * x 2 * y 3 - 6 * x 2 * y 2 - 6 * x 2 * y + 3 * x 2 + y 6 \ - 3 * y 4 + 3 * y 2 - 1 plot_implicit(f, (x, -2, 2), (y, -2, 2))

+END_SRC

+RESULTS:

[[file:tmp/sympy_fundamental5.tmp.png]]

複数の関数のプロット 1

+BEGIN_SRC ipython :session :file tmp/sympy_fundamental6.tmp.png

%matplotlib inline from sympy.plotting import plot from sympy import Symbol x = Symbol('x') plot(2x+3, 3x+1)

+END_SRC

+RESULTS:

[[file:tmp/sympy_fundamental6.tmp.png]]

複数の関数のプロット 2

色をつけたりもできる.

+BEGIN_SRC ipython :session :file tmp/sympy_fundamental7.tmp.png

from sympy.plotting import plot from sympy import Symbol x = Symbol('x') p = plot(2x+3, 3x+1, legend=True, show=False) p[0].line_color = 'b' p[1].line_color = 'r' p.show()

+END_SRC

+RESULTS:

[[file:tmp/sympy_fundamental7.tmp.png]]

Sympy 01 symbol eval printing

この notebook を読み進める前に

TeX 出力のためには次のコードを読み込ませること.

#+BEGIN_SRC ipython :session :exports none 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) #+END_SRC

#+RESULTS:

シンボルを定義する必要がある

定義すればそれをシンボルとして使える.

#+BEGIN_SRC ipython :session a = symbols ('a') display(latex(a + 1)) #+END_SRC

#+RESULTS:
a + 1

\begin{align} a + 1 \end{align}

SymPy での等号

SymPy の等号は exact structural equality を見ている. 例を見た方が早いので見てみよう.

#+BEGIN_SRC ipython :session display((x + 1)2 == x2 + 2*x + 1) #+END_SRC

#+RESULTS:
False

展開すればイコール (True) になる.

#+BEGIN_SRC ipython :session display(expand ((x+1)2) == x2 + 2*x + 1) #+END_SRC

#+RESULTS:
True

真偽判定をするには $a - b = 0$ を示すことを考えればいい. 合わせて simplify も使う.

#+BEGIN_SRC ipython :session a = (x + 1)2 b = x2 + 2*x + 1 display(simplify (a - b)) #+END_SRC

#+RESULTS:
0

equals はランダムに点を選んでその値を評価することで等号成立確認する.

#+BEGIN_SRC ipython :session a = cos(x)2 - sin(x)2 b = cos(2*x) display(a.equals (b)) #+END_SRC

#+RESULTS:
True

^/ に関する注意

#+BEGIN_SRC ipython :session display(latex(x^y)) #+END_SRC

#+RESULTS:
x \veebar y

\begin{align} x \veebar y \end{align}

#+BEGIN_SRC ipython :session display(type (Integer(1) / Integer(3))) #+END_SRC

#+RESULTS:
sympy.core.numbers.Rational

#+BEGIN_SRC ipython :session display(latex(Integer(1) / Integer(3))) #+END_SRC

#+RESULTS:
\frac{1}{3}

\begin{align} \frac{1}{3} \end{align}

SymPy での関数の評価

#+BEGIN_SRC ipython :session :exports none #まずは変数定義 from sympy import * x, y, z = symbols ("x y z") #+END_SRC

#+RESULTS:

代入

#+BEGIN_SRC ipython :session expr = cos(x) + 1 #x に y を代入 expr.subs(x, y) display(latex(cos(y) + 1)) #+END_SRC

#+RESULTS:
\cos{\left (y \right )} + 1

\begin{align} \cos{\left (y \right )} + 1 \end{align}

代入には 2 つの用法がある. 1 つはある点での値の評価.

#+BEGIN_SRC ipython :session expr.subs(x, 0) display(latex(expr.subs(x, 0))) #+END_SRC

#+RESULTS:
2

\begin{align} 2 \end{align}

ある式に他の式を入れて置き換えることも値の評価に含めよう. 例えば $x^{x}^{x}^{x}$ のように適当な対称性のある式を作りたいとき. これは次のようにして作れる.

#+BEGIN_SRC ipython :session expr = x**y display(latex(expr)) #+END_SRC

#+RESULTS:
x^{y}

\begin{align} x^{y} \end{align}

#+BEGIN_SRC ipython :session expr = xy expr = expr.subs(y, xy) display(latex(expr)) #+END_SRC

#+RESULTS:
x^{x^{y}}

\begin{align} x^{x^{y}} \end{align}

#+BEGIN_SRC ipython :session expr = xy expr = expr.subs(y, xx) display(latex(expr)) #+END_SRC

#+RESULTS:
x^{x^{x}}

\begin{align} x^{x^{x}} \end{align}

もう 1 つは複雑な式を単純にしたいとき. 例えば次のようなケース.

#+BEGIN_SRC ipython :session expr = sin(2x) + cos(2x) display(latex(expr)) #+END_SRC

#+RESULTS:
\sin{\left (2 x \right )} + \cos{\left (2 x \right )}

\begin{align} \sin{\left (2 x \right )} + \cos{\left (2 x \right )} \end{align}

$\sin 2x$ だけ展開する.

#+BEGIN_SRC ipython :session expr = sin(2x) + cos(2x) display(latex(expr.subs(sin(2x), 2sin(x)*cos(x)))) #+END_SRC

#+RESULTS:
2 \sin{\left (x \right )} \cos{\left (x \right )} + \cos{\left (2 x \right )}

\begin{align} 2 \sin{\left (x \right )} \cos{\left (x \right )} + \cos{\left (2 x \right )} \end{align}

特殊なケースについては sympy にメソッドがある. 後で紹介する予定. 例えば三角関数の 2 倍角の公式に関する展開とか.

#+BEGIN_SRC ipython :session expr = sin(2x) + cos(2x) display(latex(expand_trig(expr))) #+END_SRC

#+RESULTS:
2 \sin{\left (x \right )} \cos{\left (x \right )} + 2 \cos^{2}{\left (x \right )} - 1

\begin{align} 2 \sin{\left (x \right )} \cos{\left (x \right )} + 2 \cos^{2}{\left (x \right )} - 1 \end{align}

subs メソッドの注意

subs () は新しい式を返すことに注意しよう: SymPy のオブジェクトは不変だから. あなたが不変というのがよくわからないなら次の評価を見て意味をつかんでほしい.

#+BEGIN_SRC ipython :session expr = cos(x) display(latex(expr.subs(x, 0))) #+END_SRC

#+RESULTS:
1

subs() しても expr は cos(x) のまま.

#+BEGIN_SRC ipython :session expr = cos(x) expr.subs(x, 0) display(latex(expr)) #+END_SRC

#+RESULTS:
\cos{\left (x \right )}

\begin{align} \cos{\left (x \right )} \end{align}

$x$ に 0 を代入していても $x$ はシンボル $x$ のまま.

#+BEGIN_SRC ipython :session expr = cos(x) expr.subs(x, 0) display(latex(x)) #+END_SRC

#+RESULTS:
x

\begin{align} x \end{align}

一度にたくさん代入したいときはタプルのリストを渡すといい.

#+BEGIN_SRC ipython :session expr = x3 + 4xy - z display(latex(expr.subs([(x, 2), (y, 4), (z, 0)]))) #+END_SRC

#+RESULTS:
40

\begin{align} 40 \end{align}

リスト内包表記で一気に作ってもいい.

まずはオリジナルの式を TeX で書く.

#+BEGIN_SRC ipython :session expr = x4 - 4*x3 + 4x2 - 2x + 3 display(latex(expr)) #+END_SRC

#+RESULTS:
x^{4} - 4 x^{3} + 4 x^{2} - 2 x + 3

\begin{align} x^{4} - 4 x^{3} + 4 x^{2} - 2 x + 3 \end{align}

$x$ の偶数べきを $y$ の偶数べきに変換する.

#+BEGIN_SRC ipython :session expr = x4 - 4*x3 + 4x2 - 2x + 3 replacements = [(xi, yi) for i in range (5) if i % 2 == 0] display(latex(expr.subs (replacements))) #+END_SRC

#+RESULTS:
  • 4 x^{3} - 2 x + y^{4} + 4 y^{2} + 3

\begin{align} - 4 x^{3} - 2 x + y^{4} + 4 y^{2} + 3 \end{align}

sympify() を使うと Python の string を sympy のシンボルを使った式に変換できる. 一応書いておくが simplify() とは違う.

ただし注意! sympify では eval を使っているので, 得体の知れない文字列を通すと危険.

#+BEGIN_SRC ipython :session str_expr = "x*2 + 3x - 1/2" expr = sympify(str_expr) display(latex(expr)) #+END_SRC

#+RESULTS:
x^{2} + 3 x - \frac{1}{2}

\begin{align} x^{2} + 3 x - \frac{1}{2} \end{align}

subs() してみる.

#+BEGIN_SRC ipython :session str_expr = "x*2 + 3x - 1/2" expr = sympify(str_expr) display(latex(expr.subs(x, 2))) #+END_SRC

#+RESULTS:
\frac{19}{2}

\begin{align} \frac{19}{2} \end{align}

evalf()

数を float にするときに使う.

#+BEGIN_SRC ipython :session expr = sqrt (8) display(latex(expr.evalf())) #+END_SRC

#+RESULTS:
2.82842712474619

\begin{align} 2.82842712474619 \end{align}

evalf() は何桁までも指定できる.

#+BEGIN_SRC ipython :session expr = sqrt (2) display(latex(expr.evalf(100))) #+END_SRC

#+RESULTS:
1.414213562373095048801688724209698078569671875376948073176679737990732478462107038850387534327641573

\begin{align} 1.414213562373095048801688724209698078569671875376948073176679737990732478462107038850387534327641573 \end{align}

シンボルからなる式をある点で評価する時は次のように評価するのが効率的だし安定的な評価になる.

#+BEGIN_SRC ipython :session expr = cos(2*x) display(latex(expr.evalf(subs={x: 2.4}))) #+END_SRC

#+RESULTS:
0.0874989834394464

\begin{align} 0.0874989834394464 \end{align}

評価の時に丸め誤差が一定より小さい時, オプションで切り落とすこともできる.

#+BEGIN_SRC ipython :session one = cos(1)2 + sin(1)2 display((one - 1).evalf()) #+END_SRC

#+RESULTS:
-0.e-124

#+BEGIN_SRC ipython :session one = cos(1)2 + sin(1)2 display((one - 1).evalf(chop=True)) #+END_SRC

#+RESULTS:
0

lambdify()

複数の点での評価がほしい時の効率的な評価法を考えたい. 例えば 1000 個の点で評価したい時 SymPy ではものすごい遅い時がある. 特に任意精度の評価がほしいわけでもない場合は NumPy や SciPy を使う方がいい.

SymPy の式を数値評価可能な式に変換するのに楽には lamdify() を使うことだ.

#+BEGIN_SRC ipython :session import numpy a = numpy.arange (10) expr = sin(x) f = lambdify(x, expr, "numpy") display(f(a)) #+END_SRC

#+RESULTS:
array([ 0. , 0.84147098, 0.90929743, 0.14112001, -0.7568025 ,

-0.95892427, -0.2794155 , 0.6569866 , 0.98935825, 0.41211849])

NumPy の魔力で f (a) で一気に関数値の配列が作れるのがとても素敵. NumPy 以外に Python 標準の math モジュールも使える. 他にも自分で作ったモジュールを使うこともできる:これについては公式のチュートリアル参照.

#+BEGIN_SRC ipython :session import numpy a = numpy.arange(10) expr = sin(x) f = lambdify(x, expr, "math") display(latex(f(0.1))) #+END_SRC

#+RESULTS:
0.09983341664682815

出力

いわゆる pretty print に関わる話. SymPy オブジェクトを C, Fortran, Javascript, Theano, Python コードに変換することもできるけれども, ここでは触れない.

出力設定

その時々で最適な出力設定をしたければ init_printing() を使おう. 環境ごとに適切な設定を選んでくれる.

#+BEGIN_SRC ipython :session from sympy import init_printing init_printing() #+END_SRC

#+RESULTS:

init_session()

init_session () もかなり使える. 必要なインポート, 便利なシンボル定義, プロットに関するセットアップまでやってくれる. 物理で使う想定の場合, $k$, $m$ はばね定数, 質量という実数を割り当てたいシンボルに整数指定が入ってしまうのでちょっと微妙かもしれない.

#+BEGIN_SRC ipython :session from sympy import init_session init_session() #+END_SRC

#+RESULTS:

#+BEGIN_SRC ipython IPython console for SymPy 1.0 (Python 3.5.1-64-bit) (ground types: python)

from future import division from sympy import * x, y, z, t = symbols('x y z t') k, m, n = symbols('k m n', integer=True) f, g, h = symbols('f g h', cls=Function) init_printing()

Documentation can be found at http://docs.sympy.org/1.0/ #+END_SRC

init_printing() でどうなるか

基本的には LaTeX コンパイルされているのが一番読みやすい. LaTeX をインストールしていなくても Jupyter なら MathJax の恩恵が受けられて便利.

どうしても LaTeX を使いたくないなら use_latex=False にしよう. use_unicode = False オプションもある.

出力関数

いろいろ紹介.

str()

式の string 形式が欲しいなら str (expr) を使えばいい. これは print(expr) でも同じ結果が出せる. String 形式は読みやすさだけでなく Python の正しい構文としてコピー・ペーストできるように配慮された形式だ.

#+BEGIN_SRC ipython :session from sympy import * x, y, z = symbols ('x y z') display(str(Integral (sqrt(1/x), x))) #+END_SRC

#+RESULTS:
Integral(sqrt(1/x), x)

#+BEGIN_SRC ipython :session from sympy import * x, y, z = symbols ('x y z') print(Integral(sqrt(1/x), x)) #+END_SRC

#+RESULTS: Integral(sqrt(1/x), x)

srepr()

式の srepr 形式は式の厳密な形を表している:理解のためには出力を実際に見てみた方が早い. 詳しくは ここ 参照. 内部的な式の作り方を見たいなら参考になる形式.

#+BEGIN_SRC ipython :session srepr(Integral(sqrt(1/x), x)) #+END_SRC

#+RESULTS:
Integral(Pow(Pow(Symbol('x'), Integer(-1)), Rational(1, 2)), Tuple(Symbol('x')))
pprint (): ASCII Pretty Printer

ASCII pretty printer は pprint (). よほどの何かのこだわりがない限り, 見にくいだけだから使う必要はない.

#+BEGIN_SRC ipython :session pprint (Integral (sqrt (1/x), x), use_unicode=False) #+END_SRC

#+RESULTS:

文字列形式で欲しいなら pretty () を使おう. ここで積極的に使いたい場面は特にない.

#+BEGIN_SRC ipython :session pretty(Integral(sqrt(1/x), x), use_unicode=False) #+END_SRC

#+RESULTS:
/ \n | \n | ___ \n | / 1 \n | / - dx\n | \/ x \n | \n/

#+BEGIN_SRC ipython :session print(pretty(Integral(sqrt(1/x), x), use_unicode=False)) #+END_SRC

#+RESULTS:

/

|

| ___

| / 1

| / - dx

: \/ x
: /
pprint(), pretty (): Unicode Pretty Printer

Unicode pretty printer にも pprint()pretty () を使う.

#+BEGIN_SRC ipython :session pprint (Integral (sqrt (1/x), x), use_unicode=True) #+END_SRC

#+RESULTS:

⎮ ___

⎮ ╱ 1

⎮ ╱ ─ dx

⎮ ╲╱ x

latex(): LaTeX

latex() を使えばいい. SymPy は使えるが LaTeX が使えないというケースでは表記の参考になることもあるだろう. オプションがたくさんあるので必要ならドキュメントを見よう.

#+BEGIN_SRC ipython :session display(latex(Integral(sqrt(1/x), x))) #+END_SRC

#+RESULTS:
\int \sqrt{\frac{1}{x}}\, dx

\begin{align} \int \sqrt{\frac{1}{x}}\, dx \end{align}

print_mathml() を使う. sympy.printing.mathml からインポートすることに注意する. 対応する文字列が欲しいなら mathml() を使おう.

#+BEGIN_SRC ipython :session from sympy.printing.mathml import print_mathml print_mathml(Integral(sqrt(1/x), x)) #+END_SRC

#+RESULTS:

  <ci>x</ci>

  <root/>
  <apply>
      <power/>
      <ci>x</ci>
      <cn>-1</cn>
  </apply>

dotprint(): dot

dotprint()sympy.printing.dot prints をインポートすることで使えるようになる. dot 自体はグラフ構造を記述する言語で Graphviz で画像にレンダリングする. 詳しくはここを見よう. Graphviz はここ.

#+BEGIN_SRC ipython :session from sympy.printing.dot import dotprint from sympy.abc import x display(dotprint(x+2)) #+END_SRC

#+RESULTS:
digraph{\n\n# Graph style\n"ordering"="out"\n"rankdir"="TD"\n\n#########\n# Nodes #\n#########\n\n"Add(Integer(2), Symbol(x))()" ["color"="black", "label"="Add", "shape"="ellipse"];\n"Integer(2)(0,)" ["color"="black", "label"="2", "shape"="ellipse"];\n"Symbol(x)(1,)" ["color"="black", "label"="x", "shape"="ellipse"];\n\n#########\n# Edges #\n#########\n\n"Add(Integer(2), Symbol(x))()" -> "Integer(2)(0,)";\n"Add(Integer(2), Symbol(x))()" -> "Symbol(x)_(1,)";\n}

#+BEGIN_SRC digraph{

#Graph style "ordering"="out" "rankdir"="TD"

######### #Nodes # #########

"Add (Integer (2), Symbol (x))()" ["color"="black", "label"="Add", "shape"="ellipse"]; "Integer (2)(0,)" ["color"="black", "label"="2", "shape"="ellipse"]; "Symbol (x)_(1,)" ["color"="black", "label"="x", "shape"="ellipse"];

######### #Edges # #########

"Add (Integer (2), Symbol (x))()" -> "Integer (2)(0,)"; "Add (Integer (2), Symbol (x))()" -> "Symbol (x)(1,)"; } #+END_SRC

Sympy quaternion, 四元数

from sympy.algebras.quaternion import Quaternion
q = Quaternion(2,3,1,4)
a = Symbol('a')
p = Quaternion(a**2, a**3, a)
p*q

SymPy 行列の要素にsimplifyをあてる

import sympy
from sympy import pprint
from sympy.abc import *
import sys
sys.displayhook = pprint
from sympy.matrices import *
A = Matrix([[a, b],[c, d]])
B = A.inv()
print((B*A).applyfunc(simplify))

Windowsで実行時にutf-8文字列をprintしたときにエラーにしない

import io, sys
sys.stdout = io.TextIOWrapper(sys.stdout.buffer, encoding='utf-8')

Webサーバー立ち上げ

以下のコマンドで実行ディレクトリをルートに web サーバーが立ち上がる.

python -m http.server 80

YAML書き込み

フロースタイルではなくブロックスタイルで書き込むにはdump()default_flow_style=False を設定する.

+BEGIN_SRC python :session s1 :results output

import yaml yaml_file_name = "1.tmp.yaml" envs = {"test1": "test1", "test2": ["test2-1", "test2-2"]}

with open(yaml_file_name, "w") as f: f.write(yaml.dump(envs, default_flow_style=False))

print("END!")

+END_SRC

+RESULTS:

END!

環境変数

自宅と会社での環境設定の違いを見たかったので, 環境変数を取得するプログラムをPythonで書いた.

+BEGIN_SRC python :session s1 :results output

import subprocess

envlist = subprocess.run(["set", ">", "nul"], shell=True, stdout=subprocess.PIPE).stdout.decode("utf-8").split("\r\n") envlist = filter(lambda x: x != "", envlist) envs = {}

for elem in envlist: name, cont = elem.split("=", 1) cont_split = cont.split(";") cont = cont if len(cont_split) == 1 else cont_split envs[name] = cont

""" with open(yaml_file_name, "w") as f: f.write(yaml.dump(envs, default_flow_style=False)) """

+END_SRC

実際にはこれをyam化している. #でコメントアウトするとmarkdownのセクション記法になってしまうので, ここでは上のようなコメントアウト手法を使った: 細かいことは気にしない.

+RESULTS:

(出力していないから出力なし)

Windows コマンドを実行したいなら shell=True を指定しなければいけない模様. あと > nul を指定するといわゆる > /dev/null ができる. Windows の場合 null ではなく nul なので注意.

Windows10 で少し緩和されたとか何とか聞いてはいるが, とりあえず Windows7 だと環境変数の管理がめんどいので, 環境変数は都度指定するようにした方がいいのかもしれない.

ただ, それだと Emacs でコマンド実行させたりするときにどうなのか, 大丈夫なのかという気はする. あと会社の Windows と自宅の Mac (と Windows) で同じように動かすのがまためんどい.

三項演算子

+BEGIN_SRC python :session :results output

a = "test" b = "TRUE" if a == "test" else "FALSE" c = "TRUE" if a == "tes" else "FALSE"

print(b) print(c)

+END_SRC

+RESULTS:

Python 3.5.1 |Anaconda 4.1.0 (64-bit)| (default, Jun 15 2016, 15:29:36) [MSC v.1900 64 bit (AMD64)] on win32
Type "help", "copyright", "credits" or "license" for more information.

TRUE

FALSE

'org_babel_python_eoe'

辞書

指定回数ループ range

+BEGIN_SRC python :session s1 :results output

for i in range(655): tmp = "実績" if i % 2 == 0 else "予定" print(tmp)

+END_SRC

数学モジュール

四則演算

+BEGIN_SRC ipython :session :results output

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

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

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

print ("あまり: 5 % 3") print (5 % 3)

print ("a の b 乗: 5 3") print (5 3)

+END_SRC

+RESULTS:

+begin_example

かけ算: 2 * 2 4 割り算: 5 / 2 2.5 切り捨て割り算: 5 // 3 1 あまり: 5 % 3 2 a の b 乗: 5 ** 3 125

+end_example

特殊な演算記号

元の値に 2 を足す.

+BEGIN_SRC ipython :session :results output

x = 1 y = 100 * x print(y) # 100 y += 2 # y = y + 2 print(y)

+END_SRC

+RESULTS:

100
102

元の値から引く.

+BEGIN_SRC ipython :session :results output

y -= 22 print(y) # 80

+END_SRC

+RESULTS:

80

元の値にかける.

+BEGIN_SRC ipython :session :results output

y *= 3 print(y) # 240

+END_SRC

+RESULTS:

240

元の値を割る.

+BEGIN_SRC ipython :session :results output

y /= 10 print(y) # 24.0 実数になる

+END_SRC

+RESULTS:

24.0

元の値のあまりを取る.

+BEGIN_SRC ipython :session :results output

y %= 5 print(y) # 4.0

+END_SRC

+RESULTS:

4.0

+BEGIN_SRC ipython :session :results output

print(type(3)) # 整数 print(type(3.5)) # 浮動小数

print(int(3.5)) # 整数にする print(float(3)) # 浮動小数にする

+END_SRC

+RESULTS:

3
3.0

分数

+BEGIN_SRC ipython :session :results output

from fractions import Fraction # インポートが必要 f = Fraction(3, 4) print(f)

print(Fraction(3, 4) + 1 + 1.5) # 浮動小数になってしまう print(Fraction(3, 4) + 1 + Fraction(3, 2)) # 分数で計算してくれる

+END_SRC

+RESULTS:

3/4
3.25
13/4

複素数

定義の仕方は次の通り.

+BEGIN_SRC ipython :session :results output

a = 2 + 3j # 虚数単位には i ではなく j を使う print(type(a))

a = complex(2, 3) # Fraction のように書いてもいい print(a)

+END_SRC

+RESULTS:

(2+3j)

複素数の計算.

+BEGIN_SRC ipython :session :results output

b = 3 + 3j

print(a + b) print(a * b) print(a / b)

+END_SRC

実部・虚部を取る.

+BEGIN_SRC ipython :session :results output

print(a.real) print(a.imag)

+END_SRC

特殊な演算: 共役, 絶対値.

+BEGIN_SRC ipython :session :results output

print(a.conjugate())

print(abs(a))

+END_SRC

+RESULTS:

(2-3j)
3.605551275463989

フォーマット

+BEGIN_SRC ipython :session :results output

print('{0}'.format(1.25456)) # そのまま出力 print('{0:.2f}'.format(1.25456)) # 小数点以下 2 桁まで print('{0:.2f}'.format(1.25356)) # 3 桁目は四捨五入 print('{0:.2f}'.format(1)) # 不足分は 0 詰め

+END_SRC

+RESULTS:

1.25456
1.25
1.25
1.00

ディレクトリまたはファイルの存在確認

+BEGIN_SRC python :session s1 :results output

import os.path home = os.path.expanduser('~').replace("\", "/") print(os.path.exists(f"{home}/.skk-jisyo")) print(os.path.exists(f"{home}/nothing.txt")) print(os.path.exists(f"{home}/junk")) print(os.path.exists(f"{home}/nothing"))

+END_SRC

+RESULTS:

:

True

False

True

False

ディレクトリ内のファイルを走査してリネーム

+BEGIN_SRC python :session :results output

import os

def fild_all_files(directory): for root, dirs, files in os.walk(directory): yield root for file in files: yield os.path.join(root, file)

for file in fild_all_files('./'): print(file) name = file.replace("/", ".").replace(" ", "_").replace("..", "") print(name) # 実際のコードを実行するときはコメントアウトを解除すること # if os.path.isfile(file): # os.rename(file, name)

+END_SRC

入門

Jupyter 上の記録として fundamentals.ipynb がある.

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

Hello, World!

まずは Hello, World! から.

+BEGIN_SRC ipython :session :results output

print ("Hello, World!")

+END_SRC

+RESULTS:

Hello, World!

四則演算

別途, 数学パートで説明する.

文字列のエスケープ

+BEGIN_SRC ipython :session :results output

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)

+END_SRC

+RESULTS:

Hello "world"
A list:
  • item 1
  • item 2

C:\path\on\windows

C:\path\on\windows

リスト

シンプルなリスト.

+BEGIN_SRC ipython :session :results output

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)

+END_SRC

+RESULTS:

[3, 4, 5, 7, 6]
5
3
6
[3, 4, 5, 7]
[3, 4, 5, 7, 6, 8]

Python3 で range はリストではなくイテレータで返る.

+BEGIN_SRC ipython :session :results output

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

+END_SRC

+RESULTS:

range (0, 10)
range (1, 11)
range (0, 30, 5)
range (0, 10, 3)
range (0, -10, -1)
range (0, 0)
range (1, 0)
33

他のコンテナ: タプル, 辞書, 集合

タプル.

+BEGIN_SRC ipython :session :results output

mytuple = (1, 2, 3) print(mytuple) print(mytuple[0])

+END_SRC

+RESULTS:

(1, 2, 3)
1

辞書: キーと値のペア.

+BEGIN_SRC ipython :session :results output

mydict = {'a': 1, 'b': 2, 'c': 3} print('a:', mydict['a'])

print(mydict.keys ())

print(mydict.values ())

+END_SRC

+RESULTS:

a: 1
dict_keys (['a', 'c', 'b'])
dict_values ([1, 3, 2])

集合: 数学の集合と同じで重複は排除.

+BEGIN_SRC ipython :session :results output

myset = set([1, 2, 3, 2, 1]) print(myset)

+END_SRC

+RESULTS:

{1, 2, 3}

ループ

for が基本.

+BEGIN_SRC ipython :session :results output

a = [5,6,7,8,9]

for x in a: print (x)

+END_SRC

+RESULTS:

5
6
7
8
9

enumerate でキーと値を両方取得.

+BEGIN_SRC ipython :session :results output

for k, v in enumerate (a): print (k, v)

+END_SRC

+RESULTS:

0 5
1 6
2 7
3 8
4 9

辞書のループは次のとおり.

+BEGIN_SRC ipython :session :results output

mydict = {'a': 1, 'b': 2, 'c': 3} for k, v in mydict.items(): print (k, v)

+END_SRC

+RESULTS:

a 1
b 2
c 3

ファイル書き込み

+BEGIN_SRC python :session :results output

with open(js_staff_source_path, "w", encoding="utf-8") as f: js_code = """$(function() {{ _.extend(clmsg, {{ {} }}); }}); """.format("\n ".join(js_values)) f.write(js_code)

+END_SRC

ファイル読み込み

例外処理まで含めて with を使うこと. ファイルについては f.read() または f.readlines() で, f.readlines() はファイルの内容を配列にしてくれる.

+BEGIN_SRC python :session :results output

filepath = "create.sql" with open(filepath, "r", encoding="utf-8_sig") as f: filecont = f.readlines() for data in filecont: data = data.rstrip("\n") if data.find("CREATE") >= 0: print(data) elif data.find("is_") >= 0: print(data) elif data.find("_type") >= 0: print(data)

+END_SRC

フォーマット文字列

+BEGIN_SRC python :session :results output

js_staff_source_path = "test.tmp.js" with open(js_staff_source_path, "w", encoding="utf-8") as f: js_code = """abc {} def;""".format("\n ".join(js_values)) f.write(js_code)

+END_SRC

文字列の検索

find または re.search.

+BEGIN_SRC python :session :results output

filepath = "create.sql" with open(filepath, "r", encoding="utf-8_sig") as f: filecont = f.readlines() for data in filecont: data = data.rstrip("\n") if data.find("CREATE") >= 0: print(data) elif data.find("is_") >= 0: print(data) elif data.find("_type") >= 0: print(data)

+END_SRC

文字列の分割

文字列オブジェクトの split() メソッドを使う. 第 2 引数にはいくつめまでを分割対象にするかを指定する. rsplit() メソッドは後ろから分割する.

+BEGIN_SRC python :session s1 :results output

test = "orange,apple,banana,strawberry" print(test.split(",", 2))

+END_SRC

+RESULTS:

['orange', 'apple', 'banana,strawberry']

+BEGIN_SRC python :session s1 :results output

test = "orange,apple,banana,strawberry" print(test.rsplit(",", 2))

+END_SRC

+RESULTS:

['orange,apple', 'banana', 'strawberry']

リストから空の要素を削除 2018-01-11

Python3 だと filter の結果はイテレータになるので, 表示させるためには list() をかませる必要がある.

+BEGIN_SRC python :session s1 :results output

lst = ["a", "b", "", "c", "", "d"] newlst = filter(lambda s:s != "", lst) print(list(newlst))

+END_SRC

+RESULTS:

['a', 'b', 'c', 'd']

同じ処理を Haskell でも書いてみた.

+BEGIN_SRC haskell :session s1

main = do a = ["a", "b", "", "c", "", "d"] print $ filter (\x -> x /= "") a

+END_SRC

+RESULTS:

Prelude> ["a","b","c","d"]

リスト内包表記

+BEGIN_SRC ipython :session :results output

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

+END_SRC

+RESULTS:

[36, 64]

条件分岐

if, elif, else. elif という特殊な書き方に注意する.

+BEGIN_SRC ipython :session :results output

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

print ("")

+END_SRC

+RESULTS:

+begin_example

1 3 で割って余り 1 1

2 3 で割って余り 2 2

3 3 で割り切れる 3

+end_example

関数

+BEGIN_SRC ipython :session :results output

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

+END_SRC

+RESULTS:

False
True
False
True

文字列の連番を取得する方法

次のように教えて頂いた.

コードも転記しておこう.

+BEGIN_SRC ipython :session :results output

for char in map(chr, range(ord('A'), ord('C')+1)): s = "\newcommand{\scr%s}{\mathscr{%s}}" % (char, char) print(s)

+END_SRC

+RESULTS:

\newcommand{\scrA}{\mathscr{A}}
\newcommand{\scrB}{\mathscr{B}}
\newcommand{\scrC}{\mathscr{C}}

見ての通り, TeX のコマンドを一気に定義したかったのだ. Haskell だとどう書けばいいのだろう.