sikino のすべての投稿

時間依存しないシュレーディンガー方程式と変分原理 1/2 (計算手法の説明)

時間依存しないシュレーディンガー方程式を数値的に解くための、変分原理について述べていきます。
変分原理はエネルギー固有値の上限を与えます(厳密な解よりも変分原理で求めた解は必ず大きくなる)。
基底関数を適切に選べば高精度の結果が与えられます。ですが、基底関数をどのように選べばよいか?の指標は無く、経験則に頼らざるを得ない状況が出てきます。

ここでは、変分原理によってなぜエネルギー固有値が求まるのかを述べて行きます。

もしも変分原理より正しい解を得たければ、shooting methodと呼ばれる方法も推奨します。
恐らくこちらの方が色々なポテンシャルに対して有効ですが、
・計算速度が遅い
・発散の可能性がある
・変分原理よりも経験値が必要
です。

まとめ


区間[a,b]で定義された時間依存しないシュレーディンガー方程式
\(
\displaystyle \left[-\frac{1}{2}\frac{d^2}{dx^2}+V(x)\right]\psi(x)=E\psi(x)
\)

の解である固有関数を
\(
\displaystyle \psi(x)=\sum_n c_n \varphi_n(x) \ \ …(6)
\)

の形で固定端条件(\(\psi(a)=\psi(b)=0\))の下で探します。
ここで、\(\varphi_n(x)\)は区間\([a,b]\)で正規直交化されたsin関数系
\(
\displaystyle \varphi_n(x)=\sqrt{\frac{2}{b-a}}\sin\left(\frac{n\pi}{b-a}(x-a)\right),\;\;\;(n=1,2,\cdots,N)
\)

です。この時、行列要素\(H_{ij}=\langle\varphi_i|-\frac{1}{2}\frac{d^2}{dx^2}+V(x)|\varphi_j\rangle\)は
\(
\begin{align}
H_{ij}&=K_{ij}+V_{ij} \\
&\hspace{1em}K_{ij}=\frac{1}{2}\left(\frac{j\pi}{b-a}\right)^2\delta_{ij} \\
&\hspace{1em}V_{ij}=\frac{1}{j}\sum_{m=2}^{M-1}\omega_m^{(GL)} \cos\left(\frac{\pi}{2}x_m^{(GL)}\right)
\sum_{n=0}^{j-1}(-1)^n \sin\left(\frac{i\pi}{2j}(x_m^{GL}+2n+1)\right)V(X_n^{(m)}) \\
&\hspace{2em}X_{n}^{(m)}=\frac{h}{2}(x_m^{(GL)}+2n+1)+a,\;\;h=\frac{b-a}{j}
\end{align}
\)
と求められます(\(j\ge i\)を推奨)。ここで、\(\omega_m^{(GL)},\ x_m^{(GL)}\)はM点のGauss-Lobatto求積法の区間[-1,1]における重みと求積点で、
\(
\displaystyle \omega_m=\frac{2}{M(M-1)[P_{M-1}(x_m)]^2},\ \ x_m\ne \pm 1
\)

であり、求積点は関数
\(
\displaystyle \frac{P_{M-1}(x)+P_{M}(x)}{x-1}
\)

のゼロ点です。\(P_{n}(x)\)は\(n\)次Legendre多項式を表します。

  1. 解説
  2. 時間依存から時間依存しないシュレーディンガー方程式へ
  3. 間接法
    1. 行列要素の計算
    2. K_{nm}の計算
    3. V_{nm}の計算
    4. Gauss-Lobatto求積法による積分の計算
    5. ゼロ点に挟まれる区間の積分
  4. 数値計算プログラム
  5. 正規直交系を用いない場合
  6. 研究分野での数値計算法について
  7. 参考文献

計算の具体例は
時間依存しないシュレーディンガー方程式と変分原理 2/2 (計算例)をご覧ください。

[adsense1]

解説


対象とする問題を導出します。何はともあれ、時間依存シュレーディンガー方程式から出発します。

時間依存から時間依存しないシュレーディンガー方程式へ


時間依存シュレーディンガー方程式は
\(
\displaystyle i\frac{\partial \Psi(x,t)}{\partial t}=\hat{H}(x)\Psi(x,t)
\)

と記述されます。ここで、簡単ために原子単位系(\(\hbar=1, m=1\))を考えます。
変数分離の形で解を仮定します。
\(
\displaystyle \Psi(x,t)=\psi(x)U(t)
\)

代入して移行すると以下の式を得ます。

\(
\begin{align}
i\psi(x)\frac{\partial U(t)}{\partial t}=U(t)\hat{H}\psi(x) \\
i\frac{1}{U(t)}\frac{\partial U(t)}{\partial t}=\frac{1}{\psi(x)}\hat{H}\psi(x) = E
\end{align}
\)

2番目の式の右辺\(E\)は何らかの定数です。
左辺は時間のみを含む式、真ん中の式は位置のみを含む式なので、いつどんな時でもこの関係式が成り立っていないといけないという条件からこれらはなんらかの定数でなければなりません。その定数を\(E\)と置きました。

すなわち、2つの式
\(
\begin{align}
i\frac{\partial U(t)}{\partial t}=E U(t)\;\;\; \cdots(1)\\
\hat{H}\psi(x) = E\psi(x)\;\;\; \cdots(2)
\end{align}
\)
が導けます。

\(t\)に関する式は\(\hat{H}(x)\)に依らず解けてしまい、
\(U(t)=e^{-i E t}\)(1)
です。位置に関する式は
\(
\hat{H}\psi(x) = E\psi(x)
\)
(2)
を満たさなければなりません。この(2)が時間依存しないシュレーディンガー方程式と呼ばれているものです。
これを数値計算で解いていくことを考えます。

左から\(\psi^*\)を掛けて全空間で積分しますと、
\(
\begin{align}
\int \psi^* \hat{H} \psi dx – E \int \psi^* \psi dx = 0\;\;\; \cdots(3)
\end{align}
\)
という式を得ます。この問題は
条件\(\int \psi^* \psi dx=1\)の下で、\(\int \psi^* \hat{H} \psi dx\)の極値を求めよ、
という問いになります。
上記問題を考えた時、(3)と同一の式になります。

すなわち、原理的にはどんな\(\psi(x)\)を与えても\(E(=E[\psi])\)の値を計算できるわけです。
その無数に考えられる\(\psi\)の中で、特に(3)の解が停留する\(\psi\)、すなわち
\(
\displaystyle \delta \int \psi^* (\hat{H}-E) \psi dx = 0
\)
となる\(\psi\)を探せ、という問題になります。

↑がもしも停留する\(\psi\)が見付かれば、それはシュレーディンガー方程式と同一になるということは以下の通りに示せます。

\(
\displaystyle \delta \int \psi^* (\hat{H}-E) \psi dx =0 \;\;\;\cdots (4)
\)
厳密な\(\psi\)の変分に対して停留、すなわち\(\psi\to \psi+\delta \psi\)、もしくは複素平面で\(\frac{\pi}{2}\)ずれた方向\(\psi\to \psi+i\delta \psi\)に対して0でなければなりません。
(4)に対して2つの方向に対して変分を取ると、
\(
\begin{align}
\int (\delta \psi)^* (\hat{H}-E)\psi dx &+ \int \psi^* (\hat{H}-E)(\delta\psi) dx =0 \\
-i\int (\delta \psi)^* (\hat{H}-E)\psi dx &+ i\int \psi^* (\hat{H}-E)(\delta\psi) dx =0
\end{align}
\)
と2式が導けます。2式は変分をとった方向にかかわらず同時に満たされていなければなりません。同時に成立するためには
\(
(\hat{H}-E)\psi=0,\;\; \text{or}\;\; \psi^*(\hat{H}^{\dagger}-E^*)=0
\)
を満たしている必要があります。これは時間依存しないシュレーディンガー方程式と一致します。

以上をまとめますと、

厳密な\(\psi\)の変分はシュレーディンガー方程式に一致する

近似的に\(\psi\)を変分するとシュレーディンガー方程式の近似的な解に帰着されるだろう

近似的に\(\psi\)を構成し、変分問題を解く

代表的な方法として、
手法1)直接法…\(\psi\)の具体的な形を変分パラメータと共に与える
手法2)間接法…変分パラメータとしてある正規直交基底関数系を採用する
の2つがありますが、数値計算で使われる間接法を説明します。

間接法


間接法による変分原理は、

  • 固有値の上限を与える(変分原理によるエネルギー固有値は必ず厳密なエネルギー固有値よりも大きくなる)
  • 基底関数の数を増やすほど精度は上がる(完全系の場合)

という特性があります。

さて、今解きたい問題は
条件
\(
\displaystyle \int \psi^* \psi dx =1
\)

の下で
\(
\displaystyle \int \psi^* \hat{H} \psi dx
\)

の極値を与える\(\psi\)を探すこと、です。

すなわち
\(
\displaystyle \int \psi^* \hat{H} \psi dx – E \int \psi^* \psi dx =0 \ \ …(5)
\)

を考えて変分します。求めたい波動関数は
\(
\displaystyle \psi(x)=\sum_n c_n \varphi_n(x) \ \ …(6)
\)

と展開できる、と考えます。ここで、\(c_n\)は変分パラメータ、\(\varphi_n(x)\)は区間[a,b]で定義されている正規直交関数系です。正規直交関数系\(\{\varphi_n(x)\}\)は以下の性質を持ちます。
\(
\displaystyle \int_a^b \varphi_n(x) \varphi_m(x)dx =\delta_{nm}
\)

この制約が無くても変分は可能ですが、以降の式が簡略になるので、正規直交関数だけを考えます。

(6)を(5)に代入すると、
\(
\begin{align}
\sum_{n}\sum_{m} c_n^* c_m \left[H_{nm}-E \delta_{nm}\right]=0 \\
\end{align}\ \ …(7)
\)

なります。ここで、
\(
\begin{align}
H_{nm}&=\int_a^b \varphi_n(x) \hat{H} \varphi_m(x)dx
\end{align}
\)

とおきました。(7)式を満たす係数\(\{c_n\}\)の内、自明な解(\(c_n=0\ (for\ all\ n)\))はすぐ見つかります。しかし、この解は\(\psi(x)=0\)であり、条件\(
\displaystyle \int \psi^* \psi dx =1\)を満たしません。なので解としては不適切です。

では非自明な解を見つけましょう。(7)式の解を探すため、極値を探します。今、係数の組を変分パラメータとみなしているので、(変分パラメータで微分したもの)=0を計算します。すなわち
\(
\displaystyle \frac{\partial}{\partial c_n^*} (for\ all\ n)
\)

を(7)式に作用させます。すると、
\(
\begin{align}
\sum_{n}\sum_{m} \left[H_{nm}-E\delta_{nm}\right]c_m=0
\end{align}
\)

となります。\(n,m\)を行列の行、列ととらえれば、

\(
\begin{align}
\left(
\begin{array}{cccc}
H_{11} & H_{12} & \cdots & H_{1N}\\
H_{21} & H_{22} & \cdots & H_{2N}\\
\ddots & \ddots & \ldots & \ddots \\
H_{N1} & H_{N2} & \cdots & H_{NN}
\end{array}
\right)
\left(
\begin{array}{c}
c_{1} \\
c_{2} \\
\ddots \\
c_{N}
\end{array}
\right)
=
E
\left(
\begin{array}{c}
c_{1} \\
c_{2} \\
\ddots \\
c_{N}
\end{array}
\right)
\end{align}
\)

となるわけです。ということは、左辺の正方行列を対角化すれば右辺のE,すなわち固有値が求まる、という事になります。

行列要素の計算


行列を対角化するためには行列要素\(H_{nm}\)が分かっていないとなりません。
これからこの行列要素を計算していきます。

ハミルトニアンが
\(
\displaystyle \hat{H}=-\frac{1}{2}\frac{d^2 }{d x^2} + V(x)
\)

として表現されている場合、行列要素は
\(
\begin{align}
H_{nm}&=K_{nm}+V_{nm} \\
&\hspace{3em} K_{nm}=-\frac{1}{2}\int_a^b \varphi_n(x)\frac{d^2}{dx^2} \varphi_m(x)dx \\
&\hspace{3em} V_{nm}=\int_a^b \varphi_n(x)V(x) \varphi_m(x)dx
\end{align}
\)
と書き表され、これらを数値計算で求めます。

\(K_{nm}\)の計算


まず、微分を含んだ
\(
\displaystyle K_{nm}=-\frac{1}{2}\int_a^b \varphi_n(x)\frac{d^2}{dx^2} \varphi_m(x)dx
\)

の項を考えます。微分は数値計算で行いません。
というのも正規直交関数系はこちらで用意する知っている関数系なので、求めることが出来てしまうからです。
多くの正規直交関数系には漸化式と微分方程式が存在します。一番簡単な基底としてsin関数系を考えてみましょう。
区間\([a,b]\)で正規直交化されたsin関数系は
\(
\displaystyle \varphi_n(x)=\sqrt{\frac{2}{b-a}}\sin\left(\frac{n\pi}{b-a}(x-a)\right),\;\;\;(n=1,2,\cdots,N)
\)

と書き表せます。sin関数なので、2階微分はもちろん
\(
\displaystyle \frac{d^2}{dx^2}\varphi_n(x)=-\left(\frac{n\pi}{b-a}\right)^2\varphi_n(x)
\)

を満たします。これを積分の中に代入して、
\(
\begin{align}
K_{nm}&=-\frac{1}{2}\int_a^b \varphi_n(x)\frac{d^2}{dx^2} \varphi_m(x)dx \\
&=\frac{1}{2}\left(\frac{n\pi}{b-a}\right)^2\int_a^b \varphi_n(x) \varphi_m(x)dx \\
&=\frac{1}{2}\left(\frac{n\pi}{b-a}\right)^2\delta_nm
\end{align}
\)

となります。積分を計算することなく求める事が出来ました。

より一般には一回微分の存在しない形、もしくは漸化式によって手計算で微分を含む積分を片付けてしまったり、部分積分を用いて何とか計算します。
正規直交系であれば、うまく整理することによって誤差が全く入らず、かつ高速に計算できます。

\(V_{nm}\)の計算


さて、残る問題は
\(
\displaystyle V_{nm}=\int_a^b \varphi_n(x)V(x) \varphi_m(x)dx
\)

の計算です。基底関数の次数が低い場合、この計算は何も問題がありません。
調和ポテンシャル\(V(x)=x^2\)の時、\(n=1,m=2\)の時の被積分関数は

と非常になだらかですが、\(n=5,m=99\)の時の被積分関数は

となります。なので細かい分点をとることになり、非常に積分が難しいです。

離散フーリエ変換によって求められそうな気もしますが、周期境界条件が満たされていないことは明らかなのであきらめましょう。
有限区間で早く振動する関数の積分は難しいです。二重指数関数型数値積分も振動する関数は苦手です。

こういった場合の計算方法として、
ゼロ点の位置が分かっている場合、そこで区間を分割して積分する
という直観的な方法があります[2]。
ゼロ点の位置が分かっていれば、その間の被積分関数の振る舞いはシンプルな形になるので滑らかな関数として考えられます。

また、端点での値がゼロだと分かっている場合、Gauss-Lobatto求積法が非常に有用です。

Gauss-Lobatto求積法による積分の計算


Gauss-Lobatto求積法はGauss-Legendre求積法と同じく被積分関数を多項式で近似します[1]。
Gauss-Legendre求積法と異なる点は端点で分点を持つということです。端点を明示的に書けば、積分は
\(
\displaystyle \int_{-1}^1 f(x)dx\sim \frac{2}{N(N-1)}(f(-1)+f(1))+\sum_{i=2}^{N-1}\omega_i f(x_i)
\)

と書けます。ここで、\(\omega_i,x_i\)はGauss-Lobatto求積法における重み、求積点です。
Gauss-Lobatto求積法の重みは
\(
\displaystyle \omega_i=\frac{2}{N(N-1)[P_{N-1}(x_i)]^2},\ \ x_i\ne \pm 1
\)

であり、求積点は
\(
\displaystyle \frac{P_{N-1}(x)+P_{N}(x)}{x-1}
\)

のゼロ点です。\(P_{n}(x)\)は\(n\)次Legendre多項式を表します。

Gauss-Lobatto求積法は
\(f(\pm 1)\)が既知の値(特に0)を持つ場合
に便利なことが分かります。
関数\(f(\pm 1)\)がゼロであれば、その点を計算する必要性がなくなりますし、見かけ上の特異点であってもこちらで値を指定すれば参照されることはなく、計算できないという状況に陥りません。

ゼロ点に挟まれる区間の積分


ゼロ点に挟まれる区間の積分は、結局以下の式になります[2]。
\(
\displaystyle \varphi_j(x)=\sqrt{\frac{2}{b-a}}\sin\left(\frac{j\pi}{b-a}\right)
\)

とすると、積分は
\(
\begin{align}
&\int_a^b \varphi_i(x) V(x) \varphi_j(x) dx \\
&=\sum_{n=0}^{j-1} \int_{x_n}^{x_{n+1}} \varphi_i(x) V(x) \varphi_j(x)dx \\
&\approx \sum_{n=0}^{j-1}\sum_{m=2}^{M-1}\frac{h}{2} \omega_m^{(GL)}\varphi_i(X_n^{(m)}) V(X_n^{(m)}) \varphi_j(X_n^{(m)}) \\
&=\frac{h}{2}\sum_{m=2}^{M-1}\omega_m^{(GL)}\varphi_j(X_0^{(m)})S_j^{(m)} \\
\end{align}
\)
ここで、
\(
\begin{align}
&x_n=nh+a,\ h=\frac{b-a}{j},\ (n=0,1,\cdots,j) \\
&X_n^{(m)}=\frac{h}{2}x_m^{(GL)}+\frac{x_{n+1}+x_n}{2} \\
&X_{n+1}^{(m)}=X_n^{(m)}+h \\
&S_j^{(m)}=\sum_{n=0}^{j-1}(-1)^n \varphi_i(X_n^{(m)})V(X_n^{(m)})
\end{align}
\)
と置きました。
\(x_m^{(GL)},\ \omega_m^{(GL)}\)は区間\([-1,1]\)でのM点Gauss-Lobatto求積法の求積点と重みです。
計算を可能な限り減らすため、\(\varphi_j(x)\)の周期性を用いています。

数値計算プログラム


行列要素を組み立てた後は適当なルーチンを用いて対角化します。
今回の場合、基底関数をsin関数に選んだため、強制的に固定端条件が課されます。
なので、この基底関数系を用いるといかなる問題も端点\(a,b\)に無限に高いポテンシャルの壁がある、と読み替えることが出来ます。
よって行列は対称行列、エルミートです。故に固有値は必ず実数です。

この問題を解く場合、実対称行列を対角化するルーチン、lapackのdsyevdを用いるのが良いでしょう。
対角化のプログラムに関する詳細は対角化へどうぞ。

以下のプログラムは行列要素を上記手法で作り、対角化し、得られた波動関数を規格化し、出力するものです。

[adsense2]

コンパイルと実行は

$ ifort -mkl main.f90
$ ./a.out

などとしてください。lapackを用いますので、それと一緒にコンパイルしてください。

ラグランジュの未定乗数について


ラグランジュの未定乗数は特に物理的な意味を持たない魔法の変数、として置きますが、今回の場合は未定乗数がエネルギー固有値というちゃんとした意味を持つのです。

正規直交系を用いない場合


\(
\begin{align}
\sum_{n}\sum_{m} c_n^* c_m \left[H_{nm}-E I_{nm}\right]=0 \\
\sum_{n}\sum_{m} c_n^* c_m \left[J_{nm} H_{nm}-E \right]=0
\end{align}\ \ …(7)
\)

となります。ここで、
\(
\begin{align}
H_{nm}&=\int_a^b \varphi_n(x) \hat{H} \varphi_m(x)dx\\
J_{nm}&=I_{nm}^{-1}\\
I_{nm}&=\int_a^b \varphi_n(x) \varphi_m(x)dx
\end{align}
\)
と置きました。行列の形では、
\(
\begin{align}
\left(
\begin{array}{cccc}
J_{11}H_{11} & J_{12}H_{12} & \cdots & J_{1N}H_{1N}\\
J_{21}H_{21} & J_{22}H_{22} & \cdots & J_{2N}H_{2N}\\
\ddots & \ddots & \ldots & \ddots \\
J_{N1}H_{N1} & J_{N2}H_{N2} & \cdots & J_{NN}H_{NN}
\end{array}
\right)
\left(
\begin{array}{c}
c_{1} \\
c_{2} \\
\ddots \\
c_{N}
\end{array}
\right)
=
E
\left(
\begin{array}{cccc}
1 & 0 & \cdots & 0\\
0 & 1 & \cdots & 0\\
\ddots & \ddots & \ldots & \ddots \\
0 & 0 & \cdots & 1
\end{array}
\right)
\left(
\begin{array}{c}
c_{1} \\
c_{2} \\
\ddots \\
c_{N}
\end{array}
\right)
\end{align}
\)
と書けます。

\(
\begin{align}
\sum_{n}\sum_{m} c_n^* c_m \left[H_{nm}-E I_{nm}\right]=0 \\
\sum_{n}\sum_{m} c_n^* c_m \left[I_{nm}^{-1} H_{nm}-E \right]=0
\end{align}
\)

研究分野での数値計算法について


実際の研究の場では対角化はしますが、上記方法によって行列要素を作ることはしません。遅いからです。
特別な基底関数を用いると行列要素\(H_{nm}\)が1~2項、しかも非常に簡単な式に置き換えることが出来ます。精度も速度も上記方法以上です。
この計算方法に関してはこれ以上は述べることはしませんが、こういった方法も存在する、とだけ述べておきます。

参考文献


[1]P. J. Davis, P. Rabinowitz著, 森正武訳, 「計算機による数値積分法」, Gauss-Lobatto求積法 p.88
[2]P. J. Davis, P. Rabinowitz著, 森正武訳, 「計算機による数値積分法」, ゼロ点に挟まれる区間の積分 p.131
[3]ランダウ=リフシッツ,「量子力学1」第3刷, p81-84

ガンマ関数の数値計算

複素平面で有効なガンマ関数を数値計算します。
目指すは大浦様が公開してくださっている複素ガンマ関数cdgamma,cqgammaの計算速度を上回ることです。

結果として、本稿の複素ガンマ関数の計算コードは大浦様のコードの0.85~0.9倍ほどの時間で計算できるコードになりました。

大浦様のコードは、Lanczos近似という方法に基づいています。
Lanczos近似はガンマ関数を

\(
\begin{align}
\Gamma(z)=\frac{(z+g-0.5)^{z-0.5}}{e^{z+g-0.5}}L_g(z) \\
L_g(z)=C_0+\sum_{k=1}^{N-1}\frac{C_N}{z+k-1}
\end{align}
\)
の形で展開します[6]。係数\(C_k\)は定数\(g\)の値に依存します。
定数\(g\)は魔法の数のようなもので、値によって収束具合が異なります。このあたりの理論に関しては[5-6]が詳しいのでそちらをご参考ください。
おおよそ10-20項程度で収束し、非常にシンプルなコードになっているのが特徴です。

同じ考えでは超えられるはずがないので、ガンマ関数の特性をフルに使って計算していきます。

実数のみのガンマ関数が欲しい場合、gfortranでは組み込み関数dgamma, gammaを最適化オプション  -Oまたは-O3 とともにコンパイルすると500分の1くらいの時間で計算できるため、圧倒的に早いです。

[adsense1]

方針


方針としては

  • 漸化式
  • 相反公式
  • 逆ガンマ関数のTaylor展開
  • 対数ガンマ関数の漸近展開

を用いて計算していきます。

詳細は以下に折りたたんでおきます。

領域ごとに使っている関係式をまとめると、このような感じです。

※範囲変更(cgamma,円の大きさ8→9)

考え方は以上です。
実際にFortran90でプログラムしてみます。

fortran90によるプログラム


用意してあるのは、

複素倍精度用のガンマ関数cgamma
(本当のガンマ関数値との相対誤差 < 2e-14 [-10,10]の範囲, 8e-14 [-40,40]の範囲)


4倍精度複素ガンマ関数zgamma
(本当のガンマ関数値との相対誤差 < 3e-32 [-10,10]の範囲, 8e-32 [-40,40]の範囲)
です。

2017/01/10 逆ガンマ関数の項数の最適化
2017/10/27 正の整数の時、1ずれることの修正

function cgamma(z)
  ! sikinote (http://slpr.sakura.ne.jp/qp)
  ! author : sikino
  ! date   : 2017/01/09
  !          2017/01/10 optimize the number of a(M)
  !          2017/10/27 fix for positive integer
 
  implicit none
  complex(kind(0d0)),intent(in)::z
  complex(kind(0d0))::cgamma
 
  integer::i,n,M
  double precision,parameter::pi=3.14159265358979324d0
  double precision,parameter::e1=0.3678794411714423216d0
  double precision,parameter::ln2pi2=0.91893853320467274d0
  double precision,parameter::sq2pi=2.5066282746310005d0
  double precision::a(1:30),AB(1:14),dz,iz
  complex(kind(0d0))::q,s,w,r,q1,q2

  dz=dble(z)
  iz=imag(z)

  if(dz.eq.nint(dz).and.iz.eq.0d0)then
     if(dz.gt.0d0)then
        s=dcmplx(1d0,0d0)
        n=nint(dz)
        do i=2,n-1
           s=s*i
        enddo
     else
        s=1d+300
     endif
     cgamma=s
  else
     q=z
     if(dz.lt.0d0)q=1d0-z

     if(abs(iz).lt.1.45d0)then
        ! For x=arb, |y|\lt 1.45
        n=int(dble(q))
        s=1d0
        do i=1,n
           s=s*(q-i)
        enddo
        w=q-dble(n)

        a(1) = 1d0
        a(2) = 0.57721566490153286d0
        a(3) =-0.65587807152025388d0
        a(4) =-0.42002635034095236d-1
        a(5) = 0.16653861138229149d0
        a(6) =-0.42197734555544337d-1
        a(7) =-0.96219715278769736d-2
        a(8) = 0.72189432466630995d-2
        a(9) =-0.11651675918590651d-2
        a(10)=-0.21524167411495097d-3
        a(11)= 0.12805028238811619d-3
        a(12)=-0.20134854780788239d-4
        a(13)=-0.12504934821426707d-5
        a(14)= 0.11330272319816959d-5
        a(15)=-0.20563384169776071d-6
        a(16)= 0.61160951044814158d-8
        a(17)= 0.50020076444692229d-8
        a(18)=-0.11812745704870201d-8
        a(19)= 0.10434267116911005d-9
        a(20)= 0.77822634399050713d-11
        a(21)=-0.36968056186422057d-11
        a(22)= 0.51003702874544760d-12
        a(23)=-0.20583260535665068d-13
        a(24)=-0.53481225394230180d-14
        a(25)= 0.12267786282382608d-14
        a(26)=-0.11812593016974588d-15
        a(27)= 0.11866922547516003d-17
        a(28)= 0.14123806553180318d-17
        a(29)=-0.22987456844353702d-18
        a(30)= 0.17144063219273374d-19

        M=int(11.3*abs(w)+13)
        if(M.gt.30)M=30
        r=a(M)
        do i=M-1,1,-1
           r=r*w+a(i)
        enddo
       
        cgamma=s/(r*w)
     else
        ! For x=arb, |y|\gt 1.45
        s=1d0
        if(abs(q).lt.9d0)then
           do i=0,7
              s=s*(q+i)
           enddo
           s=1d0/s
           q=q+8d0
        endif

        AB(1) = 0.83333333333333333d-1
        AB(2) =-0.27777777777777778d-2
        AB(3) = 0.79365079365079365d-3
        AB(4) =-0.59523809523809524d-3
        AB(5) = 0.84175084175084175d-3
        AB(6) =-0.19175269175269175d-2
        AB(7) = 0.64102564102564103d-2
        AB(8) =-0.29550653594771242d-1

        q1=1d0/q
        q2=q1*q1
       
        r=AB(8)
        do i=7,1,-1
           r=r*q2+AB(i)
        enddo
        cgamma=s*exp((q-0.5d0)*log(q)-q+ln2pi2+r*q1)
     endif
     if(dz.lt.0d0)cgamma=pi/(cgamma*sin(pi*z))
  endif

  return
end function cgamma

function zgamma(z)
  ! sikinote (http://slpr.sakura.ne.jp/qp)
  ! author : sikino
  ! date   : 2017/01/09
  !          2017/01/10
  !          2017/10/27 fix for positive integer
  implicit none
  complex*32,intent(in)::z
  complex*32::zgamma
 
  integer::i,n,M
  real*16,parameter::    pi=3.14159265358979323846264338327950288q0
  real*16,parameter::    e1=0.36787944117144232159552377016146087q0
  real*16,parameter::ln2pi2=0.91893853320467274178032973640561764q0
  real*16,parameter:: sq2pi=2.50662827463100050241576528481104525q0
  real*16::a(1:50),AB(1:30),dz,iz,aw
  complex*32::q,s,w,r,q1,q2

  dz=real(z); iz=imag(z)
  if(dz.eq.nint(dz).and.iz.eq.0q0)then
     if(dz.gt.0q0)then
        s=dcmplx(1q0,0q0)
        n=nint(dz)
        do i=2,n-1
           s=s*i
        enddo        
     else
        s=1q+3000
     endif
     zgamma=s
  else
     q=z
     if(dz.lt.0q0)q=1q0-z  

     if(abs(imag(q)).lt.1.35q0)then
        ! For x=arb, |y|\lt 1.35
        n=int(real(q))
        s=1q0
        do i=1,n
           s=s*(q-i)
        enddo
        w=q-real(n)
       
        a(1) = 1q0
        a(2) = 0.577215664901532860606512090082402431q0
        a(3) =-0.655878071520253881077019515145390481q0
        a(4) =-0.420026350340952355290039348754298187q-1
        a(5) = 0.166538611382291489501700795102105236q0
        a(6) =-0.421977345555443367482083012891873913q-1
        a(7) =-0.962197152787697356211492167234819898q-2
        a(8) = 0.721894324666309954239501034044657271q-2
        a(9) =-0.116516759185906511211397108401838867q-2
        a(10)=-0.215241674114950972815729963053647806q-3
        a(11)= 0.128050282388116186153198626328164323q-3
        a(12)=-0.201348547807882386556893914210218184q-4
        a(13)=-0.125049348214267065734535947383309224q-5
        a(14)= 0.113302723198169588237412962033074494q-5
        a(15)=-2.056338416977607103450154130020572837q-7
        a(16)= 6.116095104481415817862498682855342867q-9
        a(17)= 5.002007644469222930055665048059991303q-9
        a(18)=-1.181274570487020144588126565436505578q-9
        a(19)= 1.043426711691100510491540332312250191q-10
        a(20)= 7.782263439905071254049937311360777226q-12
        a(21)=-3.696805618642205708187815878085766237q-12
        a(22)= 5.100370287454475979015481322863231803q-13
        a(23)=-2.058326053566506783222429544855237420q-14
        a(24)=-5.348122539423017982370017318727939949q-15
        a(25)= 1.226778628238260790158893846622422428q-15
        a(26)=-1.181259301697458769513764586842297831q-16
        a(27)= 1.186692254751600332579777242928674071q-18
        a(28)= 1.412380655318031781555803947566709037q-18
        a(29)=-2.298745684435370206592478580633699260q-19
        a(30)= 1.714406321927337433383963370267257067q-20
        a(31)= 1.337351730493693114864781395122268023q-22
        a(32)=-2.054233551766672789325025351355733797q-22
        a(33)= 2.736030048607999844831509904330982015q-23
        a(34)=-1.732356445910516639057428451564779799q-24
        a(35)=-2.360619024499287287343450735427531031q-26
        a(36)= 1.864982941717294430718413161878666894q-26
        a(37)=-2.218095624207197204399716913626860384q-27
        a(38)= 1.297781974947993668824414486330595180q-28
        a(39)= 1.180697474966528406222745415509988346q-30
        a(40)=-1.124584349277088090293654674261653627q-30
        a(41)= 1.277085175140866203990206677750563204q-31
        a(42)=-7.391451169615140823461289330001078362q-33
        a(43)= 1.134750257554215760954165252401478595q-35
        a(44)= 4.639134641058722029944804894433875948q-35
        a(45)=-5.347336818439198875077418068570325956q-36
        a(46)= 3.207995923613352622861238213439296964q-37
        a(47)=-4.445829736550756882101576639628564322q-39
        a(48)=-1.311174518881988712901096702293607192q-39
        a(49)= 1.647033352543813886818056435068597372q-40
        a(50)=-1.056233178503581218582902184722885523q-41

        aw=abs(w)
        if(aw.le.0.42q0)then
           M=int(35q0*aw+20q0)+1
        else
           M=int(15.5q0*aw+28.2q0)+1
        endif
        if(M.gt.50)M=50
        r=a(M)
        do i=M-1,1,-1
           r=r*w+a(i)
        enddo

        zgamma=s/(r*w)
     else
        ! For x=arb, |y|\ge 1.35
        s=1q0
        if(abs(q).lt.18q0)then
           do i=0,16
              s=s*(q+i)
           enddo
           s=1q0/s
           q=q+17q0
        endif

        AB(1) = 0.83333333333333333333333333333333333q-1
        AB(2) =-0.27777777777777777777777777777777778q-2
        AB(3) = 0.79365079365079365079365079365079365q-3
        AB(4) =-0.59523809523809523809523809523809524q-3
        AB(5) = 0.84175084175084175084175084175084175q-3
        AB(6) =-0.19175269175269175269175269175269175q-2
        AB(7) = 0.64102564102564102564102564102564103q-2
        AB(8) =-0.29550653594771241830065359477124183q-1
        AB(9) = 0.17964437236883057316493849001588940q0
        AB(10)=-1.3924322169059011164274322169059011q0
        AB(11)= 13.402864044168391994478951000690131q0
        AB(12)=-156.84828462600201730636513245208897q0
        AB(13)= 2193.1033333333333333333333333333333q0
        AB(14)=-36108.771253724989357173265219242231q0
        AB(15)= 691472.268851313067108395250775673468q0
        AB(16)=-1.52382215394074161922833649588867805q7

        q1=1q0/q; q2=q1*q1
       
        r=AB(16)
        do i=15,1,-1
           r=r*q2+AB(i)
        enddo

        zgamma=s*exp((q-0.5q0)*log(q)-q+ln2pi2+r*q1)  
     endif
     if(dz.lt.0q0)zgamma=pi/(zgamma*sin(pi*z))
  end if

  return
end function zgamma

[adsense2]

計算速度の比較


大浦様のコードと計算速度を比較します。
比較方法は、
乱数を[-20,20]の範囲で40000点振り、それを10回繰り返した時の平均cpu時間を計測します。
上の計算を20回繰り返し、その平均を取ることによって計算時間を比較します。
その結果は以下のようになりました。
大浦様のコードによる結果は、私のコードによる結果はです。

大浦様のコードと比べて、私のプログラムの方が計算時間が0.85~0.9倍になっていることが分かります。
嬉しい結果となりました。満足です。

コードの利用に関しては連絡先とサイトについてをご覧下さい。

ガンマ関数の概形


最後に、ガンマ関数の概形です。

参考文献

[1] Gamma Function -mathworld
[2]Reciprocal gamma function -wikipedia
[3]Gamma function – OeisWiki
[4] Wrench, J.W. (1968). Concerning two series for the gamma function. Mathematics of Computation, 22, 617–626
[5]Wrench, J.W. (1973). Erratum: Concerning two series for the gamma function. Mathematics of Computation, 27, 681–682
※[4]の係数表には間違いがあり、同著者によって[4]で訂正されています。
[6]The Lanczos Approximation -Boost
[7]Lanczos_approximation -wikipedia
[8] Abramowitz & Stegun (1965), ISBN 978-0486612720 (Numerical Method, 6.1.40, Page 257)
[8] Abramowitz & Stegun (1965), ISBN 978-0486612720 (Numerical Method, 6.1.34, Page 256)

Javaによる任意精度演算

javaによる任意精度演算のプログラムです。
fortranでやろうとしましたがサポートされてないようです。

私はjavaに精通していません。とりあえず何のプログラムでもいいので任意精度演算を使おうとした際、javaが一番手軽だっただけです。

ここでは、ベルヌーイ数を計算することを考えます。
擬コードとして、wikipediaより、

for m in 0, 1, …, n do
     A[m] ← 1 / (m + 1)
     for j in m, m − 1, …, 1 do
         A[j − 1] ← j × (A[j − 1] - A[j])
     end for
end for
return A[0]    # which is Bn

でベルヌーイ数が得られます。
倍精度演算(10進数で16桁)では桁落ちのため\(B_{14}\)で6桁しか一致しません。

計算過程で100桁まで用いて(割り算で100桁目で四捨五入)、\(B_{50}\)を計算してみましょう。

import java.math.BigDecimal;
 
public class bernoulli {

  final static int divideN=100;
 
  public static void main(String[] args) {
    int m,j;
    BigDecimal b  = new BigDecimal("0");
    BigDecimal bm = new BigDecimal("0");
    BigDecimal a[]= new BigDecimal[51];

    a[0]=new BigDecimal("0");

    for(m=0;m<=50;m++){
       bm = BigDecimal.valueOf(m);
         a[m]=(new BigDecimal("1")).divide(bm.add(new BigDecimal("1")), divideN, BigDecimal.ROUND_HALF_UP);
       for(j=m;j>=1;j--){
             a[j-1]=(BigDecimal.valueOf(j)).multiply(a[j-1].subtract(a[j]));
         }
       b=a[0];
       String text=String.format("%.36e",b);
       System.out.println(text);
      }
    }
}

コンパイルと実行は

javac bernoulli.java
java bernoulli

です。

出力は36桁としています。上記プログラムによる結果は、
7.500866746076964366855720075757575758e+24
となり、wolframによる値Table[bernoulli b(i), {i,0,50}] with 36 digits
7.50086674607696436685572007575757576e+24
です。全桁一致しました。
[adsense1]


以下、javaの任意精度演算に関するメモです。

メモ

変数の定義方法
  • 整数型
    int i,j;
  • 倍精度実数型
    double x;
  • 任意精度型
    BigDecimal a,b,c;
配列の定義方法
BigDecimal r[]= new BigDecimal[51];

↑で配列r[0]~r[50]が確保される。

配列の初期化

for文でまわしましょう。

for(i=0;i<=N;i++){
    r[i]=new BigDecimal("0");
}
変数へ値の代入

aに任意精度型で\(1.2000000000\ldots\)を代入する。

a=new BigDecimal("1.2");
足し算
\(
a=b+1.2
\)
a=b.add(new BigDecimal("1.2"))

(整数型の型変換)
\(
a=b+i
\)
\(i\)は整数型

a=b.add(BigDecimal.valueOf(i));

引き算
\(
a=b-1.2
\)
a=b.subtract(new BigDecimal("1.2"));

掛け算
\(
a=b\times c
\)
a=b.multiply(c);

割り算

\(
a=b/c
\)
丸めの桁を100桁にする時

a=b.divide(c, 100, BigDecimal.ROUND_HALF_UP);

足し算と割り算

\(
a=b/(c+9.92)
\)
丸めの桁を100桁にする時

a=b.divide(c.add(new BigDecimal("9.92")), 100, BigDecimal.ROUND_HALF_UP);

累乗
\(
a=b^c
\)
a=b.pow(c);

ループ

\(
\displaystyle s=\sum_{i=0}^{N-1}a_i
\)

b=new BigDecimal("0");
for(i=0;i<=N-1;i++){
    b=b.add(r[i]);
}

端末への指数表記,有効桁数36桁で出力
String text=String.format("%.36e",a);
System.out.println(text);

[adsense2]

参考文献


【Java】任意精度整数(多倍長整数)で演算する -ヽ|∵|ゝ(Fantom) の 開発blog?
Java BigDecimalの四則演算とフォーマット処理 -ITのおもちゃ箱
intとBigDecimalの型変換 -Javaプログラミング講座

逆行列の数値計算

Lapackを用いて数値計算で逆行列を求めます。
基本的に数値的にあらわに逆行列を求める意味はそんなにありません。
連立一次方程式を解きたければ、行列のLU分解を通じて解けば良いです。この場合、逆行列を求めることはしません。

そうはいっても使ってみたい時があります。そんな状況を考えましょう。

Lapack, もしくはインテルのマス・カーネル・ライブラリ(MKL)を用います。

用いるルーチンは

?getrf, ?getri

の2つです。

ここに用意したのは配列の大きさNと配列Aを入れると逆行列を出力するプログラムです。
倍精度、複素倍精度のルーチンを書きました。

ifort -mkl main.f90

でコンパイルしてください。

program main
  implicit none
  integer::i,j,N
  complex(kind(0d0)),allocatable::A(:,:),B(:,:)
 
  N=3
  allocate(A(1:N,1:N),B(1:N,1:N))
  A(1,1:3)=(/(2d0,0d0),(3d0,0d0),(0d0,0d0)/)
  A(2,1:3)=(/(2d0,0d0),(9d0,0d0),(1d0,3d0)/)
  A(3,1:3)=(/(3d0,4d0),(5d0,0d0),(3d0,1d0)/)

  B=A
 
  do i=1,N
     do j=1,N
        write(6,'(2f15.10,$)')A(i,j)
     enddo
     write(6,*)
  enddo
  write(6,*)

  call zinvmat(size(A,1),A)
 
  do i=1,N
     do j=1,N
        write(6,'(2f15.10,$)')A(i,j)
     enddo
     write(6,*)
  enddo
  write(6,*)

  A=matmul(B,A)
 
  do i=1,N
     do j=1,N
        write(6,'(2f15.10,$)')A(i,j)
     enddo
     write(6,*)
  enddo
  write(6,*)

  stop
end program main
   
subroutine invmat(N,A)
  implicit none
  integer,intent(in)::N
  double precision,intent(inout)::A(1:N)

  integer::info,lwork,ipiv(1:N)
  double precision::tw(1:1)
  double precision,allocatable::work(:)
 
  info=0
  call dgetrf(N,N,A,N,ipiv,info)
  if(info.ne.0)then
     write(6,*)"Error at inverse matrix ?getrf, info",info
     stop
  endif
 
  !size query
  call dgetri(N,A,N,ipiv,tw,-1,info)
  lwork=nint(tw(1))
  allocate(work(1:lwork+1))
  work=0d0
 
  call dgetri(N,A,N,ipiv,work,lwork,info)
  if(info.ne.0)then
     write(6,*)"Error at inverse matrix ?dgetri, info",info
     stop
  endif
 
  return
end subroutine invmat


subroutine zinvmat(N,A)
  implicit none
  integer,intent(in)::N
  complex(kind(0d0)),intent(inout)::A(1:N)

  integer::info,lwork,ipiv(1:N)
  complex(kind(0d0))::tw(1:1)
  complex(kind(0d0)),allocatable::work(:)
 
  info=0
  call zgetrf(N,N,A,N,ipiv,info)
  if(info.ne.0)then
     write(6,*)"Error at inverse matrix ?getrf, info",info
     stop
  endif
 
  !size query
  call zgetri(N,A,N,ipiv,tw,-1,info)
  lwork=nint(dble(tw(1)))
  allocate(work(1:lwork+1))
  work=dcmplx(0d0,0d0)
 
  call zgetri(N,A,N,ipiv,work,lwork,info)
  if(info.ne.0)then
     write(6,*)"Error at inverse matrix ?dgetri, info",info
     stop
  endif
 
  return
end subroutine zinvmat

ニュートン写像に現れる綺麗な画像

ニュートン写像はニュートン=ラフソン法によって、ある初期値から適当な値へ収束するまでの回数でその関数を特徴づけるものです。
フラクタルと似ている考えであることを明記しておきます。

ニュートン写像


\(
\displaystyle N_f(z)=z-\frac{f(z)}{f'(z)}
\)

の繰り返しによって適当な初期値からスタートし、収束させていきます。このページでは画像が主です。説明の方はリーマンのゼータ関数の複素力学系をご覧ください。

基本はこれだけです。簡単な複素関数でも相当複雑になります。
早速適当な関数を用意して、ニュートン写像を見ていきましょう。

[adsense1]

まずは基本的な多項式を見ます。
\(
f(z)=z^2-1
\)

\(
f(z)=z^3-1
\)

\(
f(z)=z^4-1
\)

\(
\displaystyle f(z)=z^5-1, z^6-1, z^7-1, z^10-1
\)

少し複雑な関数の場合
\(
\displaystyle f(z)=\ln z, e^z-1
\)

\(
\displaystyle f(z)=e^{1/z}-0.01
\)

\(
\displaystyle f(z)=ze^{z}\ln z-0.01
\)

\(
\displaystyle f(z)=z^5e^z-0.1
\)

\(
\displaystyle f(z)=z^4e^z\ln z-0.1
\)

\(
\displaystyle f(z)=sin(z)/z
\)

\(
\displaystyle f(z)=z^2 e^{-z^2}-1
\)

\(
\displaystyle f(z)=ln(z) e^{-z^2}-1
\)

最後にリーマンのゼータ関数

複素関数\(f(z)\)が多くの根を持つとき、非常に複雑で綺麗な画像が得られるようです。
恐らく遠く離れたところからたまたま根のあたりに到達することが多くなるからでしょう。

[adsense2]

fortran90によるプログラム


fortran90によるプログラムを置いておきます。入力するべき場所は!=====の間と、関数fと導関数dfです。

program main
  implicit none
  integer,parameter::kmax=49
  double precision,parameter::eps=1d-8
  integer::i,j,k,Nx,Ny
  double precision::xa,xb,ya,yb,x,y,hx,hy
  complex(kind(0d0))::f,df,z,z1
  external::f,df

  !==================
  xa=-5d0; xb=5d0
  Nx=201

  ya=-5d0; yb=5d0
  Ny=201
  !==================
   
  hx=(xb-xa)/dble(Nx-1)
  hy=(yb-ya)/dble(Ny-1)
  do i=0,Nx-1
     x=xa+i*hx
     do j=0,Ny-1
        y=ya+j*hy
        z=dcmplx(x,y)
        do k=0,kmax
           z1=z-f(z)/df(z)
           if(abs((z-z1)/z).le.eps)exit
           z=z1
        enddo
        write(10,'(2e18.8e3,i5)')x,y,k
     enddo
     write(10,*)
     if(mod(i,100).eq.0)write(6,'(A,i6,A,i6)')"  ",i," / ",Nx
  enddo
 
  stop
end program main

function f(z)
  implicit none
  complex(kind(0d0))::f,z
 
  f=z**5*exp(z)-0.1d0
 
  return
end function f

function df(z)
  implicit none
  complex(kind(0d0))::df,z
 
  df=5*z**4*exp(z)+z**5*exp(z)
 
  return
end function df

リーマンのゼータ関数の複素力学系

本稿の目的はリーマンのゼータ関数を、ニュートン法による写像を用いて別の視点から見ることです。
綺麗な絵が得られる[1]ので自分のコンピュータ上で見てみたいと思いました。

複素力学系とは?


複素力学系とは、簡単に言えば写像のことです。
実関数による写像の繰り返しで得られる系を力学系と呼び、
複素関数による写像の繰り返しで得られる系を複素力学系と呼びます。
”力学”の言葉の由来はニュートン力学から来ています。

”ニュートン力学”は目に見えるくらいの物体を扱う学問ですが、その正体は時間に関する2階の微分方程式です。
なので、時刻\(t\)の物体の振る舞い\(x(t)\)からちょっとの時間\(\Delta t\)秒後の物体の振る舞い\(x(t+\Delta t)\)は、適当な時間発展をさせる演算子\(f\)によって
\(
x(t)\overset{f}{\longrightarrow} x(t+\Delta t)
\)
と移り変わります。これを何度も繰り返し作用させて任意の時刻\(t’\)まで時間発展させる事が出来ます。

\(
x(t)\overset{f}{\longrightarrow} x(t+\Delta t)\overset{f}{\longrightarrow} x(t+2\Delta t)\overset{f}{\longrightarrow}x(t+3\Delta t)\overset{f}{\longrightarrow} \cdots\overset{f}{\longrightarrow} x(t’)
\)
となると、ニュートン力学を性質を知ろうとするときに重要な情報は、演算子\(f\)初期値に集約されます。

演算子\(f\)の空間へどんどん写像していくわけですので、ある写像の繰り返しも力学のようだ、だから”力学系”という名前が付けられたわけです。
もう少し詳しく知りたければ、まずは[3]を読むと良いでしょう。

本稿では複素力学系の一つであるニュートン=ラフソン法によって得られる写像(ここでは便宜的にニュートン写像と呼ぶことにします)を考えます。

複素関数\(f(z)\)の具体的なニュートン写像の数式は以下で与えられます。
\(
\displaystyle N_f(z)=z-\frac{f(z)}{f'(z)}
\)

ある複素平面の初期値からスタートして、この写像を何回繰り返したら収束したか?をその系を特徴づける値とします。
要は、ニュートン=ラフソン法の初期値を複素平面上の全ての点を取り、収束するまでの回数をカウントするのです。

[adsense1]

ニュートン写像の具体例


本稿の目的はゼータ関数に対して行うことですが、まずは簡単な複素関数のニュートン写像を見てみましょう。
まず、複素関数
\(
\displaystyle f(z)=z^2-1
\)

を考えます。根の周り(\(z=\pm 1\))では収束回数は少なくなることが予想できます。実際にやってみますとこんな感じに。

非常に単調です。実部の値が正か負かでのみ特徴づけられるようです。

続いて、複素関数
\(
\displaystyle f(z)=z^3-1
\)

を考えます。根は\(z=1, e^{\pm i2\pi/3}\)。

3つは確かに分けられていますが、その境界辺りでは非常に奇怪な模様になります。

さらに、複素関数
\(
\displaystyle f(z)=z^4-1
\)

を考えます。根は\(z=\pm 1, \pm i\)。

ゼータ関数のニュートン写像


では、ゼータ関数に適応してみましょう。
ゼータ関数の導関数は導出方法とプログラムは[2],[9]/リーマンのゼータ関数の数値計算コード(複素平面)に追記したのでそちらをご覧ください。
するとこんなニュートン写像になります。

ゼロ点付近では収束が早くなっています。ゼロ点のある場所がマークしたところです。

カラーを対数にしますとこんな図に。

(2017/01/13 追記)
もっとほかの範囲のニュートン写像を見ていきます。

全体像

拡大その1(全体像の左側の四角)

拡大その2(全体像の右側の四角)

拡大その2, 高解像度版 highres_zetanewton.png (15MB)

元の画像はこちらからダウンロードできます。

[adsense2]

ゼータ関数のニュートン写像を得るfortran90のプログラムはこちらです。

ifort -O3 -qopenmp main.f90

計算範囲はxa,xb,ya,ybで指定されています。

参考文献


[1]Tomoki Kawahira “Riemann’s zeta function and Newton’s method: Numerical experiments from a complex-dynamical viewpoint”

[2]Tom M. Apóstol “Formulas for Higher Derivatives of the Riemann Zeta Function”, Math. of. Comp vol 44, 169, 1985, pp. 223-232.

[3]Chapter 14 力学系(補講1)

[4]ジュリア集合の色付けを工夫して芸術的なフラクタル図形を描く

[5]複素力学系, ニュートン法, ジュリア集合

[6]西沢清子, 藤村雅代 “多項式のニュートン写像における鉢の幅について(カオスをめぐる力学系の諸問題)”

[7]Riemann Zeta Function -wolfram mathworld

[8]ディガンマ関数の数値計算 -シキノート
[9]リーマンのゼータ関数の数値計算コード(複素平面) -シキノート

二重指数関数型数値積分

二重指数関数型数値積分公式(DE公式またはtanh-sinh公式)とは、変数変換によって被積分関数を別の関数に変換し積分する方法のことです。
計算アイデアは

  1. 変数変換によって台形則が良く働く被積分関数に変換する
  2. 台形則によって数値積分する

という考えに立脚しています。端点に特異点が存在する場合に最適、という特徴を持ちます。
1974年に、高橋 秀俊と森 正武によって提案されました[1]。

二重指数関数型数値積分とは、関数\(f(x)\)の\([-1,1]\)に渡る積分を

\(
\begin{align}
\int_{-1}^{1} f(x) dx &\approx \sum_{i=-N_-}^{N^+} w_i f(x_i),\\
x_i&=\tanh\left(\frac{\pi}{2}\sinh(ih)\right), \\
w_i&=\frac{h\frac{\pi}{2}\cosh(ih)}{\cosh^2\left(\frac{\pi}{2}\sinh(ih)\right)}
\end{align}
\)

として近似します。\(N_-, N_+\)は離散化誤差と打ち切り誤差が等しくなるように決められます。
\(N_-, N_+\)はプログラム上では、\(x_i\)がコンピュータの扱える桁数を超えず、\(w_i\)がアンダーフローを起こさない範囲で決められます。


[adsense1]

計算アイデア


DE公式は変数変換型の数値積分公式と呼ばれます。
DE公式の計算は、

  1. 変数変換によって台形則が良く働く被積分関数に変換する
  2. 台形則によって数値積分する

という考えが元になっています。

大きな特徴として、

・端点特異性に強い
・応用範囲が広い

ことが挙げられます。
例えば、端点で発散してしまう積分
\(
\displaystyle \int_0^1 \frac{1}{\sqrt{x}} dx
\)

や、厄介な積分
\(
\displaystyle \int_0^1 \sin\left(\frac{1}{\sqrt{x}}\right)/\sqrt{x} dx
\)

の計算も少ない分点数で実行できます(後者の厄介な積分はDE公式を用いて4桁程度一致します。しかし、ほかの補間型の積分公式では1桁合うかどうかでしょう)。

なぜ変数変換でうまく計算できるのか、これを理解するにはまず、台形則が良く働く場合とは何かを知らなければなりません。

台形則が良く働く関数とは、関数の端点での微分の値が近い事です。
台形則の刻み幅を\(h\), 区間\([a,b]\)を台形則によって関数\(f(x)\)を数値積分する時、本来の積分値と台形則で求めた値の誤差を表す第1項は、
\(
\displaystyle \frac{h^2}{12}\{f'(b)-f'(a)\}
\)

と表されます。刻み幅を小さくすれば誤差が小さくなるのは直観的に分かりますが、誤差を小さくできる要因は\(\{f'(b)-f'(a)\}\)にもあります。これは、
・両端で微分の値が等しい(周期関数の1周期に渡る積分、\(a,b\)に近づくつれて関数が定数に近づいていく積分)
を考えるとき、高精度の結果を与えると言っています。

上記条件が満たされるとき、他の高次の誤差項\(f^{(n)}(b)-f^{(n)}(a)\)もほとんどゼロになることが期待されるため、数値計算誤差が限りなく小さくなっていきます。

では変数変換によって端点で減衰する被積分関数に変換することを考えましょう。
高橋-森によって以下の変数変換が提案され、その特性が調べられました[1]。
\(
\displaystyle x=\varphi(t)=\tanh\left(\frac{\pi}{2}\sinh(t)\right)
\)

この変数変換をすることによって端点で”急速に“ゼロに近づく被積分関数に変換されるのです。

この”急速に“の速度はすさまじく、tが大きくなる時、被積分関数は
\(
\displaystyle f(\varphi(t))\varphi'(t)\sim A\exp\{-c\exp(-|t|)\}
\)

の形を持ちます。指数の指数の形で減衰するため、二重指数関数型という名前が付けられているのです。

まぁ、なぜ簡単に計算できるのかと言えば、\(1/\sqrt{x}|_{x\to 0}\)で発散していく時、これよりも早く減衰する関数で割れば収束しますよね、という事です。

※三重、四重指数関数型があるのでは、と考えるのは自然です。しかし、このような関数は(多分ですが)正則であるという条件の下では存在しないことが示されています[3,4]。一重指数関数型数値積分公式は存在します(例えばIMT公式)。

ニュートン・コーツ型数値積分、ガウス求積法は補間型の数値積分公式と呼ばれます。
これらの公式は
与えられたデータ点とデータ点を(重み関数)×(多項式)として補間し、その補間された関数を積分する
というアイデアの元考えられています。
しかし、重み関数は方程式から人間が知っていないとならず、上の形以外の関数を無理やり積分しようとしても積分精度は悪いです。
そのため、問題ごとに合わせたアルゴリズムの変更が必要不可欠になります。

数値計算で気を付けるべき点


端点特異点がある場合の収束

端点特異点がある場合、変数変換後の空間で計算区間が足らなくなり、大きな打切り誤差が発生します。
端点で発散している時に顕著です。
例えば、積分
\(
\displaystyle \int_0^1 \sqrt{x}dx
\)

は高精度(16桁近く一致)の結果を与えますが、
積分
\(
\displaystyle \int_0^1 \frac{1}{\sqrt{x}}dx
\)

は8桁程度の一致しかしません。この原因は打切り誤差です。桁落ちが起きないように処理が必要な問題となります。

積分区間の変更


積分
\(
\displaystyle \int_a^b f(x) dx
\)

を行いたい場合、変数変換
\(
\displaystyle x=\frac{b-a}{2}t+\frac{b+a}{2}
\)

を行うことで区間\([-1,1]\)の積分に置き換えることが出来ます。

広義積分に対する二重指数関数型数値積分


半無限区間\([0,\infty]\)の場合、

\(
\begin{align}
\int_{0}^{\infty} f(x) dx &\approx \sum_{i=-N_-}^{N^+} w_i f(x_i),\\
x_i&=\exp\left(\frac{\pi}{2}\sinh(t)\right), \\
w_i&=h\frac{\pi}{2}\cosh(ih)\exp\left(\frac{\pi}{2}\sinh(ih)\right)
\end{align}
\)

無限区間\([-\infty,\infty]\)の場合、

\(
\begin{align}
\int_{-\infty}^{\infty} f(x) dx &\approx \sum_{i=-N_-}^{N^+} w_i f(x_i),\\
x_i&=\sinh\left(\frac{\pi}{2}\sinh(t)\right), \\
w_i&=h\frac{\pi}{2}\cosh(ih)\cosh\left(\frac{\pi}{2}\sinh(ih)\right)
\end{align}
\)

のように変数変換を行うことで二重指数で減衰する関数に変わります。

また、半無限区間の積分の多くの場合は\(exp(-x)\)で減衰します。
これを考慮すると
\(
x=\exp(t-\exp(-t))
\)

という変数変換が有効だということが分かります[1]。本稿では数値計算のアルゴリズムの都合上載せません。

計算の工夫


被積分関数の評価回数を可能な限り減らすため、アルゴリズムを工夫します。

変数変換後の空間で刻み幅\(h\)で求めた数値積分値\(I_h\)は
\(
\displaystyle I_h=h\sum_{i=-N}^N f(\varphi(ih))\varphi'(ih)
\)

と表されます。刻み幅を半分\(h/2\)にすると、分点数は\(4N+1\)になり、
\(
\displaystyle I_{h/2}={h/2}\sum_{i=-2N}^{2N} f(\varphi(ih/2))\varphi'(ih/2)
\)

と表されます。
\(I_h\)を利用して、\(I_{h/2}\)を表現すると、
\(
\displaystyle I_{h/2}=\frac{1}{2}\left\{I_h + h\sum_{i=-N}^{N-1} f(\varphi(ih+h/2))\varphi'(ih+h/2)\right\}
\)

と表されます。

この表式の利点は以下の通りです。
本来、\(I_{h/2}\)を計算するためには\(4N+1\)回右辺を計算しなければなりません。
しかし、\(I_{h}\)と\(I_{h/2}\)の分点の位置が重なっているときがちょうど2N+1点有ります。よって重なっていない残りの\(2N\)点だけを計算してやれば\(I_{h/2}\)を計算できるのです。評価回数が\(4N+1\)回から\(2N\)回に減ったことで単純に計算時間が半分になります。

[adsense2]

fortranプログラム


関数fの数値積分を実行します。
対応しているのは、

・1,2,3次元
・実関数/複素関数(実/複素引数)
・有限区間/半無限区間/無限区間



・極座標、球面座標での全空間積分

となっています。
数値計算精度は有効桁数の問題から、8桁程度と思うのが良いと思います。

複素関数の積分は、複素数点と複素数点を直線で結んだ線積分です。
また、半無限区間、無限区間の複素関数の積分では始点\(ca,wa,va\)と角度\(thx,thy,thz\)によって経路を指定することが出来ます。

厳密に全てのルーチンが正しく動作するかは確認していません。
利用する際は連絡先とサイトについてをご覧下さい。

モジュール↓(2000行近くあります)

具体的なプログラム例


積分
\(
\displaystyle \int_0^{5}\sin(\sqrt(x))
\)

を計算。

厳密値
4.3340264879445362
数値計算結果
4.33402648794427

プログラム

全具体例


実装してある計算を全ルーチンを利用したプログラムを書きます。
詳細はコメントを確認してください。

基本的にaを含む変数は積分の始点、bは積分の終点です。

program main
  use DE_integration
  implicit none
  integer::info
  double precision::xa,xb,ya,yb,za,zb,eps,s,thx,thy,thz
  complex(kind(0d0))::ca,cb,wa,wb,va,vb
  double precision,parameter::pi=dacos(-1d0)
  complex(kind(0d0))::cs
 
  double precision::f1,f2,f3,f4,p2,p3,fh1,fi1,fi2,fp2,fi3,fp3
  complex(kind(0d0))::g1,g2,g3,g4,g5,g6,h1,h2,h3,gh1,hh1,&
       gi1,hi1,gi2,hi2,gp2,gi3,hi3,gp3
  external::f1,f2,f3,f4,g1,g2,g3,g4,g5,g6,p2,p3,h1,h2,h3,&
       fh1,gh1,hh1,fi1,gi1,hi1,fi2,gi2,hi2,fp2,gp2,gi3,hi3,gp3,fi3,fp3
 
 
  ! +------- 1D Integration -------+

  ! real f, real x, 1D, finite-range
  xa=0d0;  xb=1d0;  
  eps=1d-12; s=0d0  
  call dde1d(f1,xa,xb,eps,s,info)
  write(6,*)info,s

  ! complex f, real x, 1D, finite-range
  xa=0d0;  xb=1d0;  
  eps=1d-12; s=0d0  
  call cde1d(g1,xa,xb,eps,cs,info)
  write(6,*)info,cs

  ! complex f, complex x, 1D, finite-range
  ca=dcmplx(0d0,1d0);  cb=dcmplx(1d0,3d0)  
  eps=1d-12; s=0d0  
  call zde1d(h1,ca,cb,eps,cs,info)
  write(6,*)info,cs

  ! real f, real x, 1D, semi-infinite
  xa=2d0;  
  eps=1d-12; s=0d0  
  call dde1d_hinf(fh1,xa,eps,s,info)
  write(6,*)info,s

  ! complex f, real x, 1D, semi-infinite
  xa=2d0;
  eps=1d-12; s=0d0
  call cde1d_hinf(gh1,xa,eps,cs,info)
  write(6,*)info,cs

  ! complex f, complex x, 1D, semi-infinite
  ca=dcmplx(0d0,0d0);  thx=pi/4d0
  eps=1d-12; s=0d0
  call zde1d_hinf(hh1,ca,thx,eps,cs,info)
  write(6,*)info,cs

  ! real f, real x, 1D, infinite
  xa=2d0;  
  eps=1d-12; s=0d0  
  call dde1d_inf(fi1,eps,s,info)
  write(6,*)info,s

  ! complex f, real x, 1D, infinite
  xa=2d0;
  eps=1d-12; s=0d0
  call cde1d_inf(gi1,eps,cs,info)
  write(6,*)info,cs

  ! complex f, complex x, 1D, infinite
  ca=dcmplx(0d0,0d0);  thx=pi/4d0
  eps=1d-12; s=0d0
  call zde1d_inf(hi1,ca,thx,eps,cs,info)
  write(6,*)info,cs
 
 
  ! +------- 2D Integration -------+
 
  ! real f, real x,y, 2D, finite-range
  xa=0d0; xb=1d0; ya=0d0; yb=2d0
  eps=1d-12; s=0d0  
  call dde2d(f2,xa,xb,ya,yb,eps,s,info)
  write(6,*)info,s

  ! complex f, real x,y, 2D, finite-range
  xa=0d0; xb=1d0; ya=0d0; yb=2d0
  eps=1d-12; s=0d0  
  call cde2d(g2,xa,xb,ya,yb,eps,cs,info)
  write(6,*)info,cs

  ! complex f, complex x,y, 2D, finite-range
  ca=dcmplx(0d0,1d0);  cb=dcmplx(1d0,3d0)  
  wa=dcmplx(1d0,-1d0);  wb=dcmplx(2d0,0d0)  
  eps=1d-12; s=0d0  
  call zde2d(h2,ca,cb,wa,wb,eps,cs,info)
  write(6,*)info,cs

  ! real f, real x,y, 2D, infinite-range
  eps=1d-8; s=0d0
  call dde2d_inf(fi2,eps,s,info)
  write(6,*)info,s

  ! complex f, real x,y, 2D, infinite-range
  eps=1d-12; s=0d0  
  call cde2d_inf(gi2,eps,cs,info)
  write(6,*)info,cs

  ! complex f, complex x,y, 2D, infinite-range
  ca=dcmplx(0d0,0d0);  thx=0d0
  wa=dcmplx(0d0,0d0);  thy=0d0
  eps=1d-12; s=0d0
  call zde2d_inf(hi2,ca,thy,wa,thx,eps,cs,info)
  write(6,*)info,cs

  ! real f, real x,y, 2D polar, infinite-range
  eps=1d-12; s=0d0
  call dde_polar(fp2,eps,s,info)
  write(6,*)info,s

  ! complex f, real x,y, 2D polar, infinite-range
  eps=1d-12; s=0d0  
  call cde_polar(gp2,eps,cs,info)
  write(6,*)info,cs
 
 
  ! +------- 3D Integration -------+

  ! real f, real x,y,z, 3D, finite-range
  xa=0d0; xb=1d0; ya=0d0; yb=2d0; za=2d0; zb=3d0
  eps=1d-12; s=0d0
  call dde3d(f3,xa,xb,ya,yb,za,zb,eps,s,info)
  write(6,*)info,s

  ! complex f, real x,y,z, 3D, finite-range
  xa=0d0; xb=1d0; ya=0d0; yb=2d0; za=2d0; zb=3d0
  eps=1d-12; s=0d0  
  call cde3d(g3,xa,xb,ya,yb,za,zb,eps,cs,info)
  write(6,*)info,cs

  ! complex f, complex x,y,z, 3D, finite-range  
  ca=dcmplx(0d0,1d0);  cb=dcmplx(1d0,3d0)  
  wa=dcmplx(1d0,-1d0);  wb=dcmplx(2d0,0d0)  
  va=dcmplx(1d0,1d0);  vb=dcmplx(0d0,1d0)
  eps=1d-12; s=0d0  
  call zde3d(h3,ca,cb,wa,wb,va,vb,eps,cs,info)
  write(6,*)info,cs

  ! real f, real x,y,z, 3D, finite-range
  eps=1d-12; s=0d0
  call dde3d_inf(fi3,eps,s,info)
  write(6,*)info,s

  ! complex f, real x,y,z, 3D, finite-range
  eps=1d-12; s=0d0  
  call cde3d_inf(gi3,eps,cs,info)
  write(6,*)info,cs

  ! complex f, complex x,y,z, 3D, finite-range
  ca=dcmplx(0d0,1d0);  thx=0d0
  wa=dcmplx(1d0,-1d0); thy=pi/2d0
  va=dcmplx(1d0,1d0);  thz=0d0
  eps=1d-12; s=0d0  
  call zde3d_inf(hi3,ca,thx,wa,thy,va,thz,eps,cs,info)
  write(6,*)info,cs

  ! real f, real x,y, 3D spherical, infinite-range
  eps=1d-12; s=0d0
  call dde_polar(fp3,eps,s,info)
  write(6,*)info,s

  ! complex f, real x,y, 3D spherical, infinite-range
  eps=1d-12; s=0d0  
  call cde_polar(gp3,eps,cs,info)
  write(6,*)info,cs
 
  stop
end program main

function f1(x)
  implicit none
  double precision,intent(in)::x
  double precision::f1
 
  f1=sin(sqrt(x))
 
  return
end function f1

function g1(x)
  implicit none
  double precision,intent(in)::x
  complex(kind(0d0))::g1
 
  g1=dcmplx(sin(sqrt(x)),exp(-x))
 
  return
end function g1

function h1(z)
  implicit none
  complex(kind(0d0)),intent(in)::z
  complex(kind(0d0))::h1
 
  h1=sin(z)
 
  return
end function h1

function fh1(x)
  implicit none
  double precision,intent(in)::x
  double precision::fh1
 
  fh1=exp(-x)
 
  return
end function fh1

function gh1(x)
  implicit none
  double precision,intent(in)::x
  complex(kind(0d0))::gh1
 
  gh1=dcmplx(sqrt(x)*exp(-x),exp(-x))
 
  return
end function gh1

function hh1(z)
  implicit none
  complex(kind(0d0)),intent(in)::z
  complex(kind(0d0))::hh1
  double precision,parameter::pi=dacos(-1d0)
 
  hh1=exp(dcmplx(0d0,1d0)*0.5d0*pi*z**2)
 
  return
end function hh1

function fi1(x)
  implicit none
  double precision,intent(in)::x
  double precision::fi1
 
  fi1=exp(-x*x)
 
  return
end function fi1

function gi1(x)
  implicit none
  double precision,intent(in)::x
  complex(kind(0d0))::gi1
 
  gi1=dcmplx(1d0/(1d0+x*x),exp(-x*x))
 
  return
end function gi1

function hi1(z)
  implicit none
  complex(kind(0d0)),intent(in)::z
  complex(kind(0d0))::hi1
  double precision,parameter::pi=dacos(-1d0)
 
  hi1=1d0/(1d0+z*z)
 
  return
end function hi1


function f2(x,y)
  implicit none
  double precision,intent(in)::x,y
  double precision::f2
 
  f2=sin(sqrt(x))*exp(-y)
 
  return
end function f2

function g2(x,y)
  implicit none
  double precision,intent(in)::x,y
  complex(kind(0d0))::g2
 
  g2=dcmplx(sin(sqrt(x))/sqrt(y),exp(-x)*y)
 
  return
end function g2

function h2(z,w)
  implicit none
  complex(kind(0d0)),intent(in)::z,w
  complex(kind(0d0))::h2
 
  h2=sin(z)*w
 
  return
end function h2

function fi2(x,y)
  implicit none
  double precision,intent(in)::x,y
  double precision::fi2
 
  fi2=1d0/(1d0+x**4+y**4)
 
  return
end function fi2

function gi2(x,y)
  implicit none
  double precision,intent(in)::x,y
  complex(kind(0d0))::gi2
 
  gi2=dcmplx(exp(-y*y-x*x),exp(-y*y-x*x))
 
  return
end function gi2

function hi2(z,w)
  implicit none
  complex(kind(0d0)),intent(in)::z,w
  complex(kind(0d0))::hi2
 
  hi2=exp(-z*z-w*w)
 
  return
end function hi2

function fp2(r,theta)
  implicit none
  double precision,intent(in)::r,theta
  double precision::fp2
 
  fp2=r*exp(-r)*sin(theta)**2
  fp2=r*fp2 ! Jacobian
 
  return
end function fp2

function gp2(r,theta)
  implicit none
  double precision,intent(in)::r,theta
  complex(kind(0d0))::gp2
 
  gp2=dcmplx(0d0,1d0)*r*exp(-r)*sin(theta)**2
  gp2=r*gp2 ! Jacobian
 
  return
end function gp2

function f3(x,y,z)
  implicit none
  double precision,intent(in)::x,y,z
  double precision::f3
 
  f3=sin(x*z)*exp(-y)
 
  return
end function f3

function g3(x,y,z)
  implicit none
  double precision,intent(in)::x,y,z
  complex(kind(0d0))::g3
 
  g3=dcmplx(sin(sqrt(x))/sqrt(y),exp(-x)*y*sqrt(z))
 
  return
end function g3

function h3(z,w,v)
  implicit none
  complex(kind(0d0)),intent(in)::z,w,v
  complex(kind(0d0))::h3
 
  h3=sin(z)*w*v**2
 
  return
end function h3

function fi3(x,y,z)
  implicit none
  double precision,intent(in)::x,y,z
  double precision::fi3
 
  fi3=exp(-x*x-y*y-z*z)
 
  return
end function fi3

function gi3(x,y,z)
  implicit none
  double precision,intent(in)::x,y,z
  complex(kind(0d0))::gi3
 
  gi3=dcmplx(exp(-y*y-x*x-z*z),2d0*exp(-y*y-x*x-z*z))
 
  return
end function gi3

function hi3(z,w,v)
  implicit none
  complex(kind(0d0)),intent(in)::z,w,v
  complex(kind(0d0))::hi3
 
  hi3=exp(-abs(z)**2-abs(w)**2-abs(v)**2)
 
  return
end function hi3

function fp3(r,theta,phi)
  implicit none
  double precision,intent(in)::r,theta,phi
  double precision::fp3
 
  fp3=r*exp(-r)*sin(theta)**2
  fp3=r*r*sin(theta)*fp3 ! Jacobian
 
  return
end function fp3

function gp3(r,theta,phi)
  implicit none
  double precision,intent(in)::r,theta,phi
  complex(kind(0d0))::gp3
 
  gp3=dcmplx(0d0,1d0)*r*exp(-r)*sin(theta)**2
  gp3=r*r*sin(theta)*gp3 ! Jacobian
 
  return
end function gp3

参考文献


[1] H. Takahasi and M. Masatake, “Double Exponential Formulas for Numerical Integration” (1974)
[2]渡辺二太, 二重指数関数型数値積分公式について(1990)

[3] 森 正武、「数値解析における二重指数関数型変換の最適性」
[4] M. Sugihara, Optimality of the double exponential formula-functional analysis approach,
Numer. Math. 75 (1997) 379-395.

数値積分法の入門とDE公式の説明として、
[5] 戸田 英雄, 小野 令美 数値解析における一つの話題

[6] 森 正武「二重指数関数型変換のすすめ」

[7]森 正武, “FORTRAN 77 数値計算プログラミング”、p.174-176岩波書店、(1986)第一刷

ロシアに行く前に

ロシア、特にモスクワ近辺しか行っていないので余り言えませんが、準備しておくと良い点を挙げます。

・水はけが悪い→トレッキングシューズ推奨
・トイレットペーパーを水に流してはいけない→数日程度なら日本から持参するべき
・室外と室内で温度差がすごい→多少暖かい服にとても厚いコート類10月初旬当り(1℃-10℃前後)

水はけが悪いため、道が絶えずぬかるんでいたり水たまりができていたりします。このため、トレッキングシューズなどを履いていくのが良いでしょう。

トイレットペーパーを水と一緒に流してはいけません。ペーパーが水に溶けないため詰まります
拭いた後のペーパーは横にあるバケツに捨てます。日本では考えられないかもしれませんが、首都モスクワであってもこれは普通のことです。
ホテルではバケツはないかもしれません。その時はきっとペーパーが流せるタイプなのでしょう。
ペーパーを日本から持っていくのが良いでしょう。
一応ロシアのスーパーマーケット等にも水で流せるタイプのペーパーは売っています。

室外は冬服でも寒く感じます。室外では帽子やフードをかぶるべきです。
半面、室内は半そででいられるほど暑いです。どうやら薄い服の上に厚手のコート類を着るのがこちらのスタイルのようです。

また、トイレは至る所にあるわけではなく、有料のトイレがかなりあります。
行きたいとなってから探しても見つけるのは難しいです。その時は、トイレを使わせてもらうためにお金を払うのを覚悟で近くのcafeや、レストラン等に行くのが良いでしょう。モスクワ周辺ならマクドナルドが分かりやすくて良いと思います。
他、大き目のショッピングセンターであれば無料で使えるトイレがあります。

畑政義による写像

畑政義による写像によって得られる点の集合は、簡単な式で書かれるにもかかわらず、その綺麗さに惹きつけられます
このページでは、畑政義による写像の定義、fortranコード、図を掲載します。

  1. 畑政義による写像の数式
  2. 代表的なパラメータ
  3. パラメータを予想する(2つの複素定数がゼロの場合)
  4. 考察
  5. 参考文献


ある日、こんなツイートを見ました。


とても綺麗で感動しました。また、


というツイートに触発されました。自分でも作ってみよう、と。

畑政義による写像の数式


畑政義によって考えられた写像は、簡単な数式で綺麗な図が得られます。

このページでは畑政義の写像で得られる画像の”綺麗さ”を主とします。

畑政義による写像の数式は論文[1]より、
\(
\begin{align}
F_1(z)&=\alpha z+\beta \bar{z} \\
F_2(z)&=\gamma (z-1)+\delta (\bar{z}-1)+1
\end{align}
\)

です(\(\alpha,\beta,\gamma,\delta\)は複素定数,\(\bar{z}\)は\(z\)の複素共役を意味します)。

ある1点に対して写像を行うと、\(F_1,F_2\)によっての2つ点に写像されます。
ある1点からスタートして、\(n\)回作用させていきます。\(n\)回作用させた結果、最後に得られる点の数は\(2^n\)点になります。

具体的にどういう風に計算をすればいいかといいますと
1.初期値\(z_0\)を用意
2.値\(F_1(z_0),F_2(z_0)\)を計算
3.値\(F_1(F_1(z_0)),F_1(F_2(z_0)),F_2(F_1(z_0)),F_2(F_2(z_0))\)を計算
4…
という具合に計算していくのです。そうして得られた最後の複素数の組を複素平面上に打っていくのです。

fortranコードはこちら。重複して計算するのが嫌なので少し工夫しています。

program main
  implicit none
  integer,parameter::N=12
  integer::i
  complex(kind(0d0))::a,b,c,d,z0,h(1:2**N)
 
  a=dcmplx(0.7d0,0.2d0)
  b=dcmplx(0.0d0,0.0d0)
  c=dcmplx(0d0,0d0)
  d=dcmplx(2d0/3d0,0d0)
  z0=dcmplx(1d0,0d0)

  h=dcmplx(0d0,0d0)

  call hatamap(N,a,b,c,d,z0,h)
  do i=1,2**N
     write(10,'(2f10.6)')h(i)
  enddo
 
  stop
end program main

subroutine hatamap(N,a,b,c,d,z0,h)
  implicit none
  integer,intent(in)::N
  complex(kind(0d0)),intent(in)::a,b,c,d,z0
  complex(kind(0d0)),intent(out)::h(1:2**N)

  integer::i,j,k,l,m
  complex(kind(0d0))::z,F,h0(1:2**N)
  external::F
 
  h=dcmplx(0d0,0d0); h0=dcmplx(0d0,0d0)  
 
  h(1)=z0;  h0(1:2**N)=h(1:2**N)
  do j=1,N
     k=1-2**(N-j)
     l=1-2**(N-j+1)
     do i=1,2**j
        if(mod(i,2).eq.1)l=l+2**(N-j+1)
        k=k+2**(N-j)
        m=mod(i+1,2)+1
        h(k)=F(m,a,b,c,d,h0(l))
     enddo
     h0(1:2**N)=h(1:2**N)
  enddo  
 
  return
end subroutine hatamap

function F(n,a,b,c,d,z)
  implicit none
  integer,intent(in)::n
  complex(kind(0d0)),intent(in)::a,b,c,d,z
  complex(kind(0d0))::F
 
  if(n.eq.1)then
     F=a*z+b*conjg(z)
  elseif(n.eq.2)then
     F=c*(z-1d0)+d*(conjg(z)-1d0)+1d0
  else
     write(6,*)"***error"; stop
  endif
 
  return
end function F

さて、プログラムが正しく動いていることを確かめるため、論文のパラメータで写像を再現してみます。
論文[1]では\((\alpha,\beta,\gamma,\delta)\)のセットとして描いています。
論文[1]ではこんな感じに紹介されています。
hata_paper_c

私のプログラムでは、
(左上) \((0.4614+i0.4614, 0, 0.622-i0.196, 0)\)
(右上) \( (0, 0.3+i0.3, 0, 41/50)\)
(左下) \((0, 0.5+i0.5, 0, -0.5+i0.5)\)
(右下) \((0.4614+i0.4614, 0, 0, 0.2896-i0.585)\)
hatamap_reconst_c
となり、ちゃんと正しく計算されていることが分かります。綺麗ですね。

[adsense1]

代表的なパラメータ


さて綺麗な画像を得るために探さなければならないパラメータは\(\alpha,\beta,\gamma,\delta,z_0\)の全部で5つ。各々複素定数なので、実部と虚部で合計10個のパラメータをいじって綺麗な画像を探さなければなりません。これは多いです。
まずは[2],[3],[4]に掲載されているパラメータで計算を行って見ましょう。
ずらっと並べます。
hatamap2_c

hatamap3_c

hatamap4_c

hatamap5_c

hatamap6_c

hatamap1_c

[adsense2]

パラメータを予想する


さて、綺麗な画像として紹介されているのは大体\(\alpha,\beta,\gamma,\delta\)のうち2つがゼロであり、\(0,\pm0.25,\pm 0.5, \pm0.75, \pm1\)のどれかです。この組み合わせだけでも甚大で、
\( _4C_2\times 9^8=258 280 326\)通りの組み合わせがあります。
2億個の画像を見てこれは綺麗、これは違うを判断することはできません。

注記しておきますが、写像を作る際に\(\alpha,\beta,\gamma,\delta\)の絶対値が1以下である必要はありません。
この条件は恐らく指数関数の発散を最低限防ぐためという予想だと思います。

おおよその傾向として以下の推測が出来ました。
\(\kappa, \lambda\)を\(\alpha,\beta,\gamma,\delta\)のゼロではないどれかだとすると、

  • \({\rm Re}{\kappa}={\rm Re}{\lambda}=0\)の時、つまらない画像になる
  • \({\rm Im}{\kappa}={\rm Im}{\lambda}=0\)の時、つまらない画像になる
  • \(\kappa=0\)または\(\lambda=0\)はつまらない画像になる
  • \(\alpha=\beta=0\)または\(\gamma=\delta=0\)はつまらない画像になる

だと感覚的にわかったので、これらを除外します。
さらに、\(0,\pm 0.5, \pm1\)のみを考えて僕自身が綺麗だと思う画像のみを集めてみました。
都合上、\(h=0.5\)と置いています。
また、範囲を固定するために点の撮りうる最大の値を固定してあります。

hatamap001_c3

hatamap002_c3

hatamap003_c3

hatamap004_c3

hatamap005_c3

hatamap006_c3

hatamap007_c3

考察


簡単に分かったことがあります。図の傾向は0.5の数や1の数に密接に関係しているようです。
ここで考えている4つのパラメータは、\(\kappa, \lambda\)のそれぞれ実部、虚部です。
ここで、この4つのパラメータを\((i,j,k,l)\)の組と考えましょう。括弧内の順番は特に関係なく、4つのパラメータの内どれか、を意味するとします。

例えば、初めの方にある、長方形の頂点に点が置かれたような図が得られるとき、4つのパラメータの組は\((0,0,\pm 1,\pm h)\)の組み合わせが多いことに気が付きます。
h100_c

続いて、これぞフラクタルだ、と思わせる形は\((\pm h, \pm h, \pm h, \pm h)\)、
hhhh_c

最後の方の紋様のような模様では\((\pm 1, \pm 1, \pm 1, \pm 1)\)で埋められているようです。
1111_c

恐らく、フラクタル的な雰囲気を持つ画像は、\(\kappa, \lambda\)の各々の値が有限で、絶対値が1とか、そんな時に出てくるのではないでしょうか。

今回調べた範囲のは、\(\alpha,\beta,\gamma,\delta\)のうち2つがゼロ、しかも0.5刻みしか許さないという非常に強い条件を入れた範囲です。
これだけでも膨大な数になることが、上の画像の量を見ただけで分かるでしょう。

本来は\(\alpha,\beta,\gamma,\delta\)がゼロではなくて良く、0.5刻みである必要もなく、実数ですから、今回調べたのはごくごく一部です。

綺麗な画像の定義があいまいなのもいけませんね。どうにかしてみたいものです。

また、同じく@snowy_treeさんが、


という綺麗なパラメータを集めたgifを公開してくれました。凄く綺麗なパラメータがたくさんあります!
とんでもなく綺麗な写像を見てみたいものです。

参考文献


[1]Dimensión, Análisis Multifractal y Aplicaciones
内の「3. Fractals in Mathematics (Hata)」
※中を見ると掲載されている元の論文はPatterns and Waves pp.259-278(1986)という論文のようです。
畑さんは、写像を拡張し統一的に扱おうとしたようです。なので、畑さんがこの写像を1から全て作った、というのは語弊があるかと思います。

[2]カオス #3,#6,#10,#11,#12,#13,#14,#15,#27畑政義写像 -主にコーディング

Web上で畑政義による写像が書けるサイトがありますので紹介いたします。
[2]畑政義写像で遊ぼう -救済に紹介されています。
(Play with Hata-map / 畑政義写像で遊ぼう)

[3]畑政義写像(1) -閃光的網站・弛緩複合体 Blog
※1. ただし、若干異なっています。
n回繰り返した時のみの写像の集まりのはずですが、どうやらn回に至るまでの点もプロットしているようです。n+1回繰り返しても似た画像しか出てこないのでまぁ、問題ないと思いますが…。

※2. また、どちらも上下が逆さまになっています。時々この上の二つで作成した値を私のプログラムに入れても再現できない場合があります。これに起因する問題なのか、別の要因なのかははっきりしません。

ロシアの食べ物と飲み物

ロシアのモスクワに行く機会がありまして、そこで食べた美味しい食べ物を紹介します。

カフェ、レストランの二つに分けられるでしょう。
カフェはそこそこの味、価格でご飯が食べられます。
レストランは美味しくて高いご飯が食べられます。

基本的なロシアの食事スタイルは学食に近いです(MyMy(ムームー)というカフェが典型的)。
列に沿って進んでいき、食べたいものを注文、最後に会計というシステムが多いです。

典型的な食事スタイルは
スープ、サラダ、メイン、黒パン、飲み物
の5つです。

メインは肉/魚類と、ライス/蕎麦の実/マカロニ/マッシュポテト から選んで構成されています。
スープはボルシチや、サリャンカ等のことです。
サラダは種類がいろいろありますが、ビーツや野菜、豆から作られているものが主流です。

上のスタイルを選んでいくと、モスクワ郊外では300rub前後、モスクワ市内では400-600rub程度かかります。

例えばこんな感じの献立になります。
meal1_c


borsh_c
ボルシチборщ
※←実際にお皿に入っている画像を撮ってなかったので家で食べる用のパッケージです。
スープです。”ボルシチ”と発音しても通じません。”ボ↑ル↓シ↓”で通じます。
好みによってサワークリーム(スメタナсметана)を入れる。
この組み合わせは美味しい。どうやら牛の骨髄をダシに使うようで日本では作るのが難しいと思います。
また、スープにはディル(イノンド)と呼ばれるハーブが大体乗っかっています。
生でディルをむしゃむしゃ食べましたが、ローズマリーに近い雰囲気がします。
※なんと、生のディルが西友で売っていました。驚きです。

サリャンカСолянка
スープです。魚、肉、キノコ等種類は豊富にあるようですが、私は魚のサリャンカを食べました。
とても美味しいです。味噌の代わりに別の何かが入ったあら汁に近い感じですかね。

syarma_c
シャルマШаурма
ファストフードです。ケバブのチキン版とでもいえばいいのでしょうか。ケバブとは別物です。
アラブ、イスラエル地方のファストフードのようです。
非常に美味です。値段も安く量もたっぷり。ビールと合います。
ケバブのように焼かれている鶏肉の塊に指さしながらこれください、といえばきっと通じます。
110-120rub
Шаурма-google画像検索


黒パン
酸っぱいパン。
正直に言うとなんでみんなこれを食べてるのかわからない、そんな美味しいものではないと思います。
子供の時からずっと食べてきているものなのでしょう。
安く、どこに行ってもあるので食べることが出来ます。面白い味なので試してみてはどうでしょう。
10rub

sprat_c
スプラット(шпроты)
スプラットという魚のオイル漬け缶詰です。
燻製されているためとってもおいしい。鰹節の風味がします。
レモンを添えてウォッカと一緒に食べるそうです。
100rub



飲み物
gohan2_c
コンポートкомпот
果物を煮て作ったジュースです。
図の中央上にある飲み物がそうです。
ジャムのような一般的なコンポートとは違います。
ご飯を食べるときに一緒に飲みます。売られているのは基本的に甘いです。
とある家庭にお邪魔して頂いたのですが、砂糖の甘味は余りなくスッキリ爽やかな味です。各家では8月ごろに50ビン位作っておくそうです。
10-20rub


kbac
クバスквас
コーラ、ドクターペッパーとは違った方向を向いた微炭酸飲料です。
本格的なシーズンは夏。外にクバス売りが来るそうです。
ペットボトルのクバスは美味しくありません…紹興酒に近い匂いがしました。
手作りのクバスはとても美味しいです。味の表現が難しいですね。
一応缶でしたら、←のラベルの物がまぁまぁの味でした。


baikal
バイカルбайкал
コカ・コーラをやさしくしてドクターペッパーのチェリーをなくした感じの炭酸飲料です。
個人的にコーラよりも好き。
2リットル60rub



tarkhuna
タルフンтархун
エストラゴンというハーブから作られた微炭酸飲料です。
すごい緑色しているけど着色はしていない模様。
2種類飲みましたが、このメーカー(左の画像)のタルフン(ASCANIA)が植物感があって良いです。
実際にエストラゴンを生で食べましたが、まさに飲み物と同じ味でした。
※「タラゴン」の名前で乾燥させたエストラゴンが西友で売られていました。たしかS&Bだった気がします。
また、Amazonでも普通に売られているんですね。

もう一種類飲みましたが、ただのメロンソーダで面白くありませんでした。
60rub



kefir
ケフィアкефир
あのケフィア。爽やかに酸っぱい。牛乳とヨーグルトの中間位でトロトロしてる飲み物。
とってもとってもおいしい。ちなみに牛乳はモロコмолокоと呼ばれています。
日本では種菌しか買えませんが、ロシアではペットボトルで売られています。
80rub前後



エナジードリンク
gorilla0_c
Gorilla URL:http://gorillaenergy.com
緑(original):普通のエナジードリンク味。モンスターの緑とほぼ同じ味。
赤(pomegranate):正直美味しくない。表現しがたい後味の悪さ。pomegranateの意味は果物のザクロです。
青(ice):あるのは見かけましたが飲んでいないので味は不明。
ピンク(pepper):味はまさにドクターペッパー。ドクターペッパーを飲める人は美味しく飲める。
70-90rub



adrenalin
Adrenalin
柑橘系のエナジードリンク。色は黄色で味は爽やか。
70rub



effect
effect
70rub
普通のエナジードリンクの味。特に目立った味はしていない。



jaguar0_c
Jaguar
エナジードリンクとお酒のカクテル。7.2%。
レッドブルやモンスター(緑)のエナジードリンクとアルコールを混ぜました、という感じです。
URL:http://jaguar-drink.com/
美味しい…んでしょうか?
80rub



[adsense2]


アルコール

kozel
Kozel
チェコのブランドのビールです。ですが、ロシアで売られているこれはロシア製造でロシアで消費されるとのこと。
なのでロシア製…でいいんでしょうか?
飲んだ中で最も美味しかったビールの一つです。
Kozelは日本人が飲みやすいものでしょう。
65rub
他にpremiumという種類もありましたが、それらは輸入しているらしく230rubしました。



zhigulevskoye
ジグレフスコエЖигулёвское
これは完全にロシア純正のビールで、伝統あるもの。ソビエト時代(もしくはもっと前)から作られ続けられています。
現地の方も”хорошо“と言っていました。
実際に飲んでみても美味しい。Kozelよりも苦みがあります。
1つの会社が出しているだけの物ではありません。何というんでしょうか。
何かの製法に従ったものをジグレフスコエ、と呼ぶのだと思いますが良く分かりません。
2.8%以上のアルコールと15%までの麦芽穀物?を添加することが許されているもの、のようですZhigulevskoye -wikipedia 英語



Russian Standard(Русский Стандарт)
russian_standard_bottle

“ロシア スタンダード” とでも呼べばいいのでしょうか。
ロシア色に染まりたいならこれです。ロシア人が少し贅沢したい時にこれを飲むんだそう。
ロシアのウォッカとしては高級な部類に入ります。全く違うロシア人3人に聞いて、3人からお勧めされました。
ベルーガの半分程度の値段で買うことが出来ます。実際飲みましたが、とても飲みやすかったです。
日本のAmazonでありました。

飲む前に冷凍庫に入れてキンキンに冷やしておきましょう。凍ることはありません。
スーパー500rub, 空港1400rub前後
Русский Стандарт (водка) -ロシア語wikipedia
Russian Standard (vodka) -英語wikipedia


(c)Achim Raschka(2012)


ウォッカ
beluga
ベルーガBeluga
ベルーガ(チョウザメの意味)と呼ばれるウォッカが高品質の物のようです。
現地のお酒大好き人に勧められました。
確かに美味しいです。そして高い(スーパーで1000rub, 空港で2000rub)。
空港の出発ロビーで買うことができます…というより、地元のスーパーよりも空港で買えと言われました。
本当かわかりませんが、地元のスーパーだと粗悪品が混ざっていたり、何があるかわからないそうです。
ただ、空港では高いです。
ロシアの有名なウォッカですが、かなり有名なため、日本でもAmazon等で買うことが出来ます。

ベルーガ ウォッカ クラシック 1L 1000ml ベルーガ -Amazon


vodkas

(左)バリザム
薬草を使ったロシアのアルコール。45-50%位。
イェーガーマイスターに近い?
ロシアの知人曰く、「45%と、アルコールが強いからコーヒーや紅茶に飲むことをお勧めする。ウォッカは40%だから、ストレートで飲んで大丈夫。」だそうです。
冷凍庫に入れてから飲むと40%であることを感じさせません。
味は…イェーガーマイスターを思わせますが、はっきりと違うものだと分かります。甘くないコーラやドクターペッパーの味が近いんでしょうか…?(右)のベルーガhunting HURBALよりは薬草感は強いです。

(中央)ウォッカ
種類は調べていませんが、ウォッカです。

(右)ベルーガ hunting HURBAL, URL:http://www.beluga-hunting.com/
BERRYとHURBALがあります。私が買ったのはHURBALで、1100rub位だった気がします。
ホームページによれば、お勧めは冷凍庫(-15℃)で冷やしてストレートで飲むこと。
また、お勧めのカクテルは50mlのHURBALとトニックビターレモン100mlを混ぜてローズマリーを添える事
だそうです。
日本語のページではほぼ紹介されているページが無く、ましてや通販などされていません。
味は、炭酸と砂糖を減らしたコーラに薬草を足した味です。100%のレモンジュースとよく合いました。美味しいです!

ウォッカは地方のスーパーで700mlくらいの瓶で300-500rubで購入できます。これはアル中になれる。