2018-01-10

調和振動子と平面波の比較

量子力学

零点振動している調和振動子と同じエネルギーを持つ平面波とで波動関数を比較する。
$\displaystyle l_0 = \sqrt{ \frac{ \hbar }{ m \omega } } \\ \displaystyle \xi = \frac{ x }{ l_0 } \\ \displaystyle \psi_n( x ) = A_n H_n\left( \frac{ x }{ l_0 } \right) e^{ - \frac{1}{2} \left( \frac{ x }{ l_0 } \right)^2 } = A_n H_n( \xi ) e^{ - \frac{ \xi^2 }{ 2 } } \\ \displaystyle A_n = \sqrt{ \frac{ 1 }{ 2^n \; n! } } \; \sqrt{ \frac{ l_0 }{ \sqrt{ \pi } } } \\ \displaystyle H_0 = 1 \\ \displaystyle E_0 = \frac{ \hbar \omega }{ 2 }$
最後ら辺、サボりました。

ここで、 $E_0 = E_k$ として自由平面波のエネルギーを代入して、波数と有効距離 $l_0$ を結びつけると、
$\displaystyle l_0 = \sqrt{ \frac{ \hbar }{ m } \frac{ \hbar }{ 2 E_k } } = \sqrt{ \frac{ \hbar^2 }{ m } \frac{ 2 m }{ 2 \hbar^2 k^2_0 } } = \frac{ 1 }{ k_0 }$
綺麗に有効距離の逆数になる。

これで求めた $k_0$ を持つ平面波（ $cos( k_0 x )$ ）を $n=0$ の調和振動子波動関数と比較する（規格化定数は無視）と、全体で一周期分をちょっと超えたかなくらいの波が調和振動子の中に入っていることがわかる。
f:id:koideforest:20180110005333p:plain

もう少し厳密に見てみると、 $k_0$ を逆に固定して、 $k_0 x_0 = \pi$ となる $x_0$ における調和振動子(係数は無視)の値を求めると、
$\displaystyle e^{ - \frac{1}{2} \left( \frac{ x }{ l_0 } \right)^2 } = e^{ - \frac{1}{2} ( k_0 x )^2 } \rightarrow e^{ - \frac{1}{2} ( k_0 x_0 )^2 } = e^{ - \frac{\pi^2}{2} } = 0.00719... \sim 0.01 \ll 1$
したがって、十分小さいと見なせるところまでの範囲で、平面波一周期分が入っていることがわかる。

思っていたよりも調和振動子と平面波の対応関係がしっかりしていてビックリした。
原点近傍でポテンシャルがゼロに近くなるから、その近傍で自由解に近づくから当然といえば当然であるが。
自由解だけ波動関数があるような領域は、運動エネルギー（波動関数）ではなくポテンシャル項によってエネルギーが補填されているんだなぁとしみじみ感じた。

2017-12-14

PythonのLaguerre陪関数について

プログラミング数学量子力学

Wikipediaには日頃からよくお世話になっているが、自分で確認するのは大事だなと感じた事件。

水素原子の波動関数の可視化とかカッコイイなぁと思い、ネットサーフィン中に以下のサイト様を訪問。
scipyの特殊関数 - 篠突く雨の日記

自分も早速コピーしてやってみて、ふむふむと散布図の便利さに舌鼓を打っていた。
ちょっとアレンジして、波動関数の正負で色を変えるようにしたら気付いてしまった。
「ん？1s波動関数に節があんぞ？？」
波動関数の部分はscipy.special.assoc_laguerreで実装されていて、原因はこいつであることがわかった。
何か引数がおかしいのかと思い、Wikipediaで水素原子の解析解におけるLaguerre陪関数を確認。
水素原子におけるシュレーディンガー方程式の解 - Wikipedia
$\displaystyle \rho \frac{ d^{ 2 } u( \rho ) }{ d\rho^{ 2 } } + ( 2l + 2 - \rho ) \frac{ d u( \rho ) }{ dx } + ( n - l - 1 ) u( \rho ) = 0 \\ \displaystyle \rho \frac{ d^{ 2 } u( \rho ) }{ d\rho^{ 2 } } + ( ( 2l + 1 ) + 1 - \rho ) \frac{ d u( \rho ) }{ dx } + ( ( n + l ) - ( 2l + 1 ) ) u( \rho ) = 0 \\ \displaystyle u( \rho ) = c L^{ 2 l + 1 }_{ n + l }( \rho ) \\ \displaystyle L^m_k( \rho ) = \frac{ d^m }{ d\rho^m } e^{ \rho } \frac{ d^k }{ d\rho^k } ( e^{ -\rho } \rho^k ) \\$
Legendreではlとmは直接打ち込めば良かったけど、Laguerreはそういえばややこしかったんだよなぁ、と遠い記憶が蘇ってきた。
しかし、可視化スクリプトの引数はこれと対応したものが入っている。
実は上付きと下付きは逆の表記になっているとかあるか？？と思い、他も参照。

ラゲールの陪多項式 - Wikipedia
$\displaystyle \left( x \frac{ d^{ k + 2 } }{ dx^{ k + 2 } } + ( k + 1 - x ) \frac{ d^{ k + 1 } }{ dx^{ k + 1 } } + ( n - k ) \frac{ d^k }{ dx^k } \right) L^k_n( x ) = 0 \\ \displaystyle L^k_n( x ) = \frac{ d^{ k } }{ dx^{ k } } L_n( x ) = \frac{ d^{ k } }{ dx^{ k } } \left( e^x \frac{ d^{ n } }{ dx^{ n } } ( x^n e^{ -x } ) \right)$
上がkだったり下がkだったりしてめちゃくちゃややこしいが、とりあえずLaguerre多項式（下付き添字のみ）を上付き添字の分だけ微分しまくるのが陪関数であることがわかった。
その意味で、水素原子におけるシュレーディンガー方程式の解 - Wikipediaの上付き下付きの意味はラゲールの陪多項式 - Wikipediaと同じ。

scipy.special.assoc_laguerreのサイトに行ってみた。
scipy.special.assoc_laguerre — SciPy v1.1.0 Reference Guide
見る限り、スクリプトの引数の順番も間違っていなさそうだった。

しかし、

scipy.special.assoc_laguerre( x, n + l, 2 * l + 1 )

はn=1, l=0で $-x + 1$ のプロットを示し、意味不明。

そして英語版のwikipediaに辿り着く。
Laguerre polynomials - Wikipedia
$\displaystyle xy'' + ( \alpha + 1 - x ) y' + \beta y = 0 \\ \displaystyle L^{ ( \alpha ) }_0( x ) = 1 \\ \displaystyle L^{ ( \alpha ) }_1( x ) = 1 + \alpha - x$
これや、、、これやったんや。。。

つまり、scipy.special.assoc_laguerreを水素様原子の波動関数に使うには、
$\displaystyle L^{ \alpha }_{ \beta } = L^{ 2 l + 1 }_{ n - l - 1 }$
ということになる。
確かにこれで1s (n=1, l=0)の時 $L^1_0 = 1$ となり、めでたく波動関数の節が無くなる。

お手軽に特殊関数が使える様になったとはいえ、油断は禁物ですな。

2017-12-14

モンテカルロ積分で水素原子のエネルギー期待値を確認

プログラミング量子力学

「スレーター型軌道の数値積分は大変」と呪文のように聞かされて続けて心にハードルが出来てしまったので、ここらで打破することを試みた。
前回、多次元のモンテカルロ積分をまとめた。
koideforest.hatenadiary.com
今回はスクリプトを組んで、どれくらいちゃんと動くのか見てみることにした。
目標は、非相対論水素原子の解析解を使って、ハミルトニアンの期待値を出す。
以下、スクリプトのポイントをまとめる。

運動エネルギーを出すところの微分は、中心差分法を使った。簡単便利。

積分範囲は、水素原子波動関数の値が0.01を下回る距離をrmaxと定義し、その範囲で水素原子の8分の1が入る箱を考えた。対称性から、その箱の中の積分値を8倍すれば全体の積分値が求まる。

そんなことをしたのは、乱数の範囲が $[ 0, 1 ]$ の0始まりなので、全体を平行移動したりして水素原子全体を箱の中に無理矢理入れるのは余り得策ではなく感じたから。今回は素直に水素原子の原点を0にした。

箱の中の全体で積分値を求めるから、乱数点に関する制限は必要ない。下手にif分を入れると遅くなるらしいので、そこまで一般化してスクリプトを作るのは余り良くないかもしれない。今回はそこの一般化は行っていない。

以下ソース。

import numpy as np
import time

def MonteCalro3d_ws( func, scales ): # ws: whole space
    M = 10**6
    vec_rand = np.random.rand( M * 3 ).reshape( M, 3 )
    # Caution: the number of samplings is still M, not M * 3
    vec_rand = np.array( scales ) * vec_rand
    tot = np.sum( [ func( v[0], v[1], v[2] ) for v in vec_rand ] )
    return np.prod( scales ) * tot / float( M )

def h1s( x, y, z ):
    vec = np.array( [ x, y, z] )
    r   = np.linalg.norm( vec )
    return 2. * np.exp( - r ) / np.sqrt( 4. * np.pi )

def d2func_dx2( func, x ):
    h = 1.0e-7
    return ( func( x + 2.*h ) - 2.*func(x) + func( x - 2.*h ) ) / (2.*h)**2

def integrand( x, y, z ):
    v = np.array( [ x, y, z] )
    r = np.linalg.norm( v )
    deriv2x = d2func_dx2( lambda x: h1s( x,     v[1], v[2] ), v[0] )
    deriv2y = d2func_dx2( lambda y: h1s( v[0],     y, v[2] ), v[1] )
    deriv2z = d2func_dx2( lambda z: h1s( v[0],  v[1],    z ), v[2] )
    ke = - ( deriv2x + deriv2y + deriv2z ) / 2.
    pot_n = - h1s( v2[0], v2[1], v2[2] ) / r
    return h1s( v[0], v[1], v[2] ) * ( ke + pot_n )

rmax = max( [ x_ for x_ in np.linspace( 0, 5, 1000 ) \
              if h1s( x_, 0., 0. ) / h1s( 0., 0., 0. ) > 0.01 ] )
scales = [ rmax, rmax, rmax ]

start = time.time()
result = 8. * MonteCalro3d_ws( integrand, scales )
print( "Hydrogen atom energy = {}".format( result )  )
end = time.time()
print( end - start )

結果：
-0.500215806919
148.051234961

原子単位系は原子単位系でも、Hartree unitで計算しているので、理論値は-0.5でちゃんと計算できている。

ちなみにHartreeかRydbergかは運動エネルギーが $1/2$ されているかどうかで判断できる。

$1/2$ の無いRydbergの場合には、核ポテンシャルの項も二倍して理論値は-1になる。

細かい値を忘れたが、ざっくり 1 Hartree ~ 27 eV, 1 Rydberg ~ 13 eV, 1 Hartree = 2 Rydberg である。最後の等式は厳密。

値の安定性はサンプリング数 $M$ に依存するが、 $M = 10^6$ であれば有効数字一桁は安心できると思う。

cythonとか使えばもっと速くなる。

2017-12-13

N次元のモンテカルロ積分

数学

一般的なアイデアとして、積分を体積 x 密度（割合）に焼き直すところがスタート。
$\displaystyle I = \int_{ \Omega } d \textbf{ r } \; f( \textbf{ r } ) \ = \left( \int_{ \Omega } d \textbf{ r } \right) \ \left( \frac{ \int_{ \Omega } d \textbf{ r } \; f( \textbf{ r } ) }{ \int_{ \Omega } d \textbf{ r } } \right) \ = V_{\Omega} < f >_{\Omega}$

密度（割合）は単位量当たりにどれだけ対象となる事物があるかを表し、いわゆる期待値として焼き直せる。連続量の期待値なので、（確率）密度 or 重みをかけた積分に対応。それを離散量で近似しましょうというノリ。このあたりは本当は測度論とかが顔を出すんだろうなぁと妄想していますが、不勉強なので回避。
$\displaystyle < f >_{\Omega} = \int_{ \Omega } d \textbf{ r } \; \rho_{ \Omega }( \textbf{ r } ) f( \textbf{ r } ) \ \sim \sum^{ m }_{ \mu_1 } \cdots \sum^{ m }_{ \mu_N } p_{\Omega}( r^{ (1) }_{ \mu_1 }, \cdots, r^{ (N) }_{ \mu_N } ) f( r^{ (1) }_{ \mu_1 }, \cdots, r^{ (N) }_{ \mu_N } )$

$m$ はサンプリング数またはメッシュの細かさ。トータルのメッシュ数は $M \equiv m^N$ で定義する。
一次元と三次元に関して先ほどの和を解して書いてみると、
$\displaystyle \sum^{ m }_{ \mu } p_{\Omega}( x_{ \mu } ) f( x_{ \mu } ) = p_{\Omega}( x_1 ) f( x_1 ) + p_{\Omega}( x_2 ) f( x_2 ) + \cdots + p_{\Omega}( x_m ) f( x_m ) \\ \displaystyle \sum^{ m }_{ \mu_x } \sum^{ m }_{ \mu_y } \sum^{ m }_{ \mu_z } p_{\Omega}( x_{ \mu_x }, y_{ \mu_y }, z_{ \mu_z } ) f( x_{ \mu_x }, y_{ \mu_y }, z_{ \mu_z } ) \\ \displaystyle \qquad = p_{\Omega}( x_1, y_1, z_1 ) f( x_1, y_1, z_1 ) + p_{\Omega}( x_2, y_1, z_1 ) f( x_2, y_1, z_1 ) \\ \displaystyle \qquad + \cdots + p_{\Omega}( x_m, y_m, z_{ m-1 } ) f( x_m, y_m, z_{ m-1 } ) + p_{\Omega} ( x_m, y_m, z_m ) f( x_m, y_m, z_m )$

一次元において、 $0 \leq x \leq a$ における積分を、等間隔のメッシュを切った長方形で近似する場合で書いてみれば、もう少し期待値の描像が明確になる。
$\displaystyle \int^{a}_{0} dx \; f(x) \sim \sum^{ m }_{ i = 1 } \frac{ a }{ m } f( ( i - 1 ) a / m ) \\ \displaystyle \int^{a}_{0} dx \; f(x) = a < f >_a \\ \displaystyle \therefore \; < f >_a \sim \sum^{ m }_{ i = 1 } \frac{ 1 }{ m } f( ( i - 1 ) a / m ) \ = \sum^{ m }_{ i = 1 } p f( a ( i - 1 ) / m ) \\ \displaystyle \qquad \sim \sum^{ m }_{ i = 1 } p f( a R_i )$

全空間において一様に離散化・サンプリングすると、それは平均値を取りましょうということになる。
そして、これの近似として、 $0 < R_i < 1$ で規格化された昇順の乱数列 $\{ R \}$ を定義して置き換えたものが、積分領域に制約のない場合のモンテカルロ積分となる。
近似関係が分かりやすいようにあえて昇順と言ったが、結局は和なので順番は関係ない。
とりあえず多重積分したいという場合にはこっちの考え方になるだろう。いちいちメッシュの切り方を自分で考えなくて良いので、手軽で良い。
$p$ が単純にサンプリング数の逆数なので、モンテカルロ積分の例題で有名な $\pi$ を求める問題（円の中に入った点の数を数える）は、この考えでは出来ない。

これまでの話は領域 $\Omega$ の中で閉じていた。
なので、今までの $p_{\Omega}$ はあまり確率の意味はなく、単純に重みの役割が大きかった。
この領域 $\Omega$ を扱い易い任意の箱 $\Box$ に拡張するが、その際に余分な領域での値がゼロになるような関数 $\Theta$ を定義する。
$\displaystyle \int_{ \Omega } d \textbf{ r } \; f( \textbf{ r } ) \ = \int_{ \Box } d \textbf{ r } \; \Theta( \textbf{ r } ) f( \textbf{ r } ) = V_{ \Box } < \Theta f >_{ \Box } \\ \Theta( \textbf{ r } ) = \left\{ \ \begin{array}{ ll } 1 & ( \textbf{ r } \in \Omega ) \\ 0 & (otherwise) \end{array} \right.$

これを先ほどと同様に密度関数で表したのちに離散化すると、
$\displaystyle < \Theta f >_{\Box} = \int_{ \Box } d \textbf{ r } \; \rho_{ \Box }( \textbf{ r } ) \Theta( \textbf{ r } ) f( \textbf{ r } ) \ = \int_{\Box} d \textbf{ r } \; \rho( \textbf{ r } ) f( \textbf{ r } ) \\ \displaystyle \sim \sum^{ m }_{ \mu_1 } \cdots \sum^{ m }_{ \mu_N } p( r^{ (1) }_{ \mu_1 }, \cdots, r^{ (1) }_{ \mu_N } ) f( r^{ (1) }_{ \mu_1 }, \cdots, r^{ (N) }_{ \mu_N } )$

上と同様に一次元の例で考えれば、
$\displaystyle \int^{a}_{0} dx \; f(x) = \int^{b}_{0} dx \; \theta( a - x ) f(x) \sim b \sum^{ m }_{ i = 1 } \frac{ \sigma( ( i - 1 ) b / m ) }{ m } f( ( i - 1 ) b / m ) \\ \displaystyle \int^{a}_{0} dx \; f(x) = b < \theta_a f >_b \\ \displaystyle \therefore \; < \theta_a f >_b \sim \sum^{ m }_{ i = 1 } \frac{ \sigma( ( i - 1 ) b / m ) }{ m } f( ( i - 1 ) b / m ) \ = \sum^{ m }_{ i = 1 } p( b ( i - 1 ) / m ) f( b ( i - 1 ) / m ) \\ \displaystyle \qquad \sim \sum^{ m }_{ i = 1 } p( b R_i ) \; f( b R_i ) \\$

ここで、
$\displaystyle \sigma( x ) = \left\{ \begin{array}{ll} 1 & ( x \in [ 0, \; a ] ) \\ 0 & (otherwise) \end{array} \right.$
と定義した。

$\sigma$ 及び $p$ で、積分領域の調整を乱数点が領域に入った入らないかで行うことが出来る。積分領域が複雑な時には、この方法が威力を発揮することになる。
乱数自体は $[ 0, 1 ]$ の間でしか値が出ないでの、計算したい系に合わせてスケール因子を掛ける必要がある。
一つの次元を $c$ 倍する度に結果を $c$ 倍すれば正しい積分結果が得られる。

2017-11-16

ラグランジュの未定乗数法と等高線プロット

プログラミング数学

ラグランジュの未定乗数法を直感的に非常に分かりやすく描いてあるサイト様を発見。
ラグランジュの未定乗数法の解説と直感的な証明

せっかくなので、実際に等高線プロット使って図示してみた。
f:id:koideforest:20171116003805p:plain

図で考えているのは、

$f(x,y) = ( x - 3 )^2 + y^2$
$g(x,y) = x^2 + y^2 - 4$

のような、何かポテンシャル的な場が $f(x,y)$ で与えられていて、軌道 $g(x,y)=0$ の上で $f(x,y)$ を最大もしくは最小にする $(x,y)$ は何でしょう？という問題。
数学の問題っぽく言えば、束縛条件 $g(x,y) = 0$ の下で $f(x,y)$ を最大（最小）にする $(x,y)$ を求める問題。
文章を読むと、最大最小化する $f(x,y)$ が主役な感じがするが、図を見るとわかる様に、（当たり前だが） $g(x,y)=0$ の軌道の中の点で良さそうなのを探す方がイメージに近い。
ちなみに、二次元の今の問題において得られる最大（最小）の座標群がなぜ零次元（つまり点）なのかというと、 $g(x,y)=0$ の時点で $x$ と $y$ は独立ではなく、既に一次元の問題になっているところが味噌である。そこに $f(x,y)$ を介して条件式が加わると、零次元になるというわけである。三次元でやると得られる座標群は線（つまり軌道・一次元）になる。
単純に変数 $x, y$ それぞれに対して微分するから、方程式が未定乗数を含めた変数の数だけ出てきて、最終的に値がカッチリ決まるというもの。未定乗数の微分では $g(x,y) = 0$ がそのまま出てくる。
この意味で、主役はなんとなく $g(x,y)=0$ のような気がしてならない。

図を見れば明らかに
最大： $( -2, 0 )$
最小： $( 2, 0 )$
であるが、ラグランジュの未定乗数法を使うと、

$L = ( x - 3 )^2 + y^2 - \lambda [ x^2 + y^2 - 4$ ]
$\frac{ \partial L }{ \partial x } = 2 ( 1 - \lambda ) x - 6 = 0$
$\frac{ \partial L }{ \partial y } = 2 ( 1 - \lambda ) y = 0$
$\frac{ \partial L }{ \partial \lambda } = - ( x^2 + y^2 - 4 )= 0$

ここから、

（2018年10月26日　訂正）
$x = 3 / ( 1 - \lambda )$
$( y, \lambda ) = ( 0, a ) \; or \; ( a, 1 ) \; for \; \forall a \in \mathbb{R}$
$x^2 + y^2 = 4$

三つ目のは元々の束縛条件まんまなので、当たり前は当たり前。
$y = 0$ をその三つ目に入れれば、 $x = \pm 2$ が求まって、図から類推した答えと合っている。
$\lambda = 1$ は最初の式を満たさない（発散する）ので、棄却出来る。

以下、ソースの中身。

import numpy as np
import matplotlib.pyplot as pet

N = 100
x = y = np.linspace( -5, 5, N )
X, Y = np.meshgrid( x, y )
F = ( X - 3 )**2 + Y**2
G = X**2 + Y**2 - 4

plt.text( 0, 0, 'g(x,y)=0' )
plt.text( -2.5, 3.5, 'f(x,y)', color='#0068b7')

plt.contour( X, Y, G, [0] )
plt.gray()
plt.pcolor( X, Y, F )
plt.cool()
plt.gca().set_aspect( 'equal', adjustable='box' )
plt.show()

以下のサイト様にお世話になりました。

等高線プロット
matplotlibカラープロット - 週末はいつも晴れ

陰関数
http://tokyoneet.com/how-to-use-matplotlib-2/

2017-11-15

matplotlibの等高線プロットでCDのジャケットを作る

プログラミング

等高線プロットを練習してたら何かカッコイイのが出来たので、CDジャケット風にしてみた。
f:id:koideforest:20171115232319p:plain

import numpy as np
import matplotlib.pyplot as plt

N = 5

x = y = np.arange( N )
X, Y = np.meshgrid( x, y )

z = np.random.rand( N * N )
Z = z.reshape( N, N )

plt.figure( figsize = ( 4.72, 4.72 ) )
plt.tick_params( left='off', right='off', bottom='off', top='off' \
                 labelleft='off', labelright='off', labelbottom='off', labeltop='off' )
plt.contour( X, Y, Z )
plt.viridis()

plt.show()

CDジャケットのサイズは約4.72平方インチらしい。
敢えて小さめにして枠を残して印刷してもカッコイイかも。

2017-11-15

python: matplotlib.pyplot: インスタンスを使わない一通りのまとめ

プログラミング

いつも忘れるので、個人的なメモ

参考にしたサイトの皆様
グラフの体裁を整える — matplotlib 1.0 documentation
matplotlibでの凡例(ラベル)の表示場所・形式を変更する - ゆるふわめも
 matplotlibで忘れがちな文法まとめ - Qiita
【Python】matplotlibによるグラフ描画時のColormapのカスタマイズ - Qiita
matplotlibの凡例を半透明に - Qiita
matplotlib入門 - りんごがでている
 [Python] matplotlib: 論文用に図の体裁を整える - Qiita
【Python】matplotlibのフォントを変えた際にクリアする必要のあるキャッシュ - 俺言語。
matplotlib と Seaborn の軸の日本語設定 - Qiita
api example code: font_family_rc.py — Matplotlib 2.0.2 documentation
matplotlibでグラフの文字サイズを大きくする - Qiita
自分用メモ: matplotlibを使ったグラフの調整 - Qiita
Matplotlib 超入門(3)文字の大きさ，グリッド幅
 matplotlibで忘れがちな文法まとめ - Qiita
matplotlibで色付きのベクトルを作成 - ゴルゾンの日記
 matplotlib/plot - memoring
matplotlibで、系列の凡例を枠外に表示する - Symfoware
matplotlibでrectangleを描く | mwSoft
matplotlib.axes.Axes.tick_params — Matplotlib 2.1.0.post851+g2fd479a documentation
Customizing matplotlib — Matplotlib 2.0.2 documentation

import matplotlib as mpl

# change font
mpl.get_cachedir()
mpl.rcParams[ "font.family" ] = [ "Arial" ]

# change font size
mpl.rcParams[ "font.size" ] = 18

# change the line width of the frame/
mpl.rcParams[ "axes.linewidth" ] = 2.0

# adjust the positon of ticks
mpl.rcParams[ "xtick.major.pad" ] = 3.5
mpl.rcParams[ "ytick.major.pad" ] = 3.5


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

# set figure size
plt.figure( figsize = ( 8, 6 ) )

# plot data
for id, data in enumerate( data_list ):
  x, y = data[0], data[1]
  #cy = "#00a0e9"
  #ma = "#e4007f"
  #gr = "#009944"
  #or = "#f39800"
  #bl = "#0068b7"
  #clolar_list = [ cy, ma, gr, or, bl ]
  #color = cm.tab10( id  )
  #color = cm.viridis( id / len( data_list )  )
  color = cm.cool( id / len( data_list )  )
  plt.plot( x, y, lw = 2, c = color, label = label_list[ id ] )

# add comments
plt.text( xt, yt, "comment", fontsize = 10 )

# add arrows
plt.quiver( x_position, y_position, \
            x_component, y_component, \
            angles = "xy", scale_units = "xy" )

# add rectangular
plt.gca().add_patch( plt.rectangle( [ x_left, y_bottom ], width, height, fill=False ) )

# add horizontal and/or vertical lines
# when you use hlines or vlines, added lines are behind the plot.
#plt.hlines( y_list, xmin, xmax, linestyle = "dashed" )
#plt.vlines( x_list, ymin, ymax, linestyle = "dashed" )
plt.plot( [ xmin, xmax ], [ y, y ], linestyle = "dashed"  )
plt.plot( [ x, x ], [ ymin, ymax ], linestyle = "dashed"  )

# change representation of labels of ticks
plt.xticks( [ x1, x2, x3 ], [ 'X1', 'X2', 'X3' ] )
plt.yticks( np.linspace( ymin, ymax, dividing_number ) )

# change parameters of ticks
plt.tick_params( direction = 'out', width = 1, length = 8, labelsize = 12)

# set each minimum and maximum
plt.xlim( xmin, xmax )
plt.ylim( ymin, ymax )

# set labels for each axis
plt.xlabel( 'horizontal', fontsize = 14 )
plt.ylabel( 'vertical', fontsize = 14 )

# allow to appear legends
#plt.legend( loc = 'upper left', bbox_to_anchor = ( 1.05, 1.00 ),  frameon = False, fontsize = 10 )
#plt.subplots_adjust( right = 0.7 )
plt.legend( loc = 'best', frameon = False, fontsize = 10 )

# avoid lacking of characters of the labels
plt.tight_layout()

# save figure
plt.savefig( 'filename.eps', format='eps', transparent=True, dpi=600 )

ちなみに調べたところ、tickという単語そのものには目盛りの意味は無く、「瞬間」という意味らしい。
tick markで目盛りという意味になるらしいので、そこから来ていると思われる。

カラーマップの実際の色は下記のサイトで見れる。
Choosing Colormaps — Matplotlib 2.0.2 documentation

色の選択には下記のサイトを参照
統計グラフの色

f:id:koideforest:20171115080506p:plain
左が通常、右がphotoshop仕様

でも matplotlib が version 2 になってからdefaultの色の設定でも十分見易くなったと思う。color map の viridis は default で設定されているが、全然悪くない。
今度論文書くときは安直だがそれで図を作ってみることにする。

nano_exit

基礎的なことこそ、簡単な例が必要だと思うのです。