Bluepostech

まじめなことを書きます

ガウス積分チートシート

数学

理系学生としてあってはならないことだが、僕は頻繁にガウス積分を忘れる。導出の方法は覚えているのだが、毎回導出するのはかなりあほらしいしまちがいのもとにもなりかねないので、ここにまとめておこうと思う。

1.一番基本的なケース(e^{-x^2})

\displaystyle \int^{\infty}_{-\infty} e^{-x^2}dx

積分の積に関する可換性と変数変換を用いる。

I=\displaystyle \int^{\infty}_{-\infty} e^{-x^2}dx =\int^{\infty}_{-\infty} e^{-y^2}dyとおくと、

\begin{eqnarray} I^2 &=& \displaystyle \left( \int^{\infty}_{-\infty} e^{-x^2}dx\right) \cdot \left( \int^{\infty}_{-\infty} e^{-y^2}dy \right) \\ &=& \displaystyle \int^{\infty}_{x=-\infty}\int^{\infty}_{y=-\infty}e^{-(x^2+y^2)}dxdy \end{eqnarray}

ここで2次元極座標に変換する、すなわち

\left\{ \begin{eqnarray} r &=& \sqrt{x^2+y^2} \\ \theta &=& \mathrm{arctan}\left( \frac{y}{x} \right) \end{eqnarray} \right.

なる変数変換を行う。ただし、変域はxもyも負の無限大から正の無限大だったので、点対称となり\thetaの項は2\piとなる。すなわち

\begin{eqnarray} I^2 &=& \displaystyle \int^{\infty}_{0} e^{-r^2}\cdot 2\pi rdr \\ &=&-\pi \displaystyle\int^{\infty}_{0} -2r\cdot e^{-r^2} dr \\ &=& \pi \left[e^{-r^2}\right]^{0}_{\infty} \\ &=& \pi \end{eqnarray}

となる。ここの積分がわからない場合は4を参照されたい。e^xの正値性より積分も正なので、

I = \displaystyle \int^{\infty}_{-\infty} e^{-x^2}dx = \sqrt{\pi}
を得る。

--------------

2.指数の肩が定数倍されているケース(e^{-ax^2})

\displaystyle \int^{\infty}_{-\infty} e^{-ax^2}dx \hspace{5mm}(a>0)

これは単なる置換積分の問題で、y=\sqrt{a}xと置換すれば、
\displaystyle \int^{\infty}_{-\infty} e^{-ax^2}dx=\sqrt{\frac{\pi}{a}}
となる。

--------------

3.積分区間が0以上のみのケース

\displaystyle \int^{\infty}_{0} e^{-ax^2}dx \hspace{5mm}(a>0)

これも、e^{-ax^2}が偶関数であることを考えれば単純で、2の結果を半分にすればよい。すなわち
\displaystyle \int^{\infty}_{0} e^{-ax^2}dx =\frac{1}{2}\sqrt{\frac{\pi}{a}}
また、次のように考えてもよい。1の議論で「積分は全空間だから角度積分は2\piとなる」とした。それに対して積分区間が0から無限大の場合、{I^2}の積分領域は座標平面でx、yがともに正である第一象限のみとなり、I^2は1/4、結局Iは1/2となるのである。

--------------

4. xがかかっているケース(xe^{-x^2})

\displaystyle \int^{\infty}_{0} xe^{-ax^2}dx

1の最後の部分の積分そのものなのだが、置換積分である。詳しくやっておく。

\displaystyle \int^{\infty}_{0} xe^{-x^2}dx = -\frac{1}{2}\int^{\infty}_{0} -2xe^{-x^2}dx
u=-x^2とおくと、du=-2xdx、u:0\to-\inftyとなり、
\begin{eqnarray}\displaystyle \int^{\infty}_{0} xe^{-x^2}dx &=& -\frac{1}{2}\int^{-\infty}_{0} e^{u}du \\ &=& \frac{1}{2}[e^u]^{0}_{-\infty} \\ &=& \frac{1}{2} \end{eqnarray}

ちなみに一応言っておくが、被積分関数は奇関数と偶関数の積だから、
\displaystyle \int^{\infty}_{-\infty} xe^{-x^2}dx=0
である。

係数がついて\displaystyle \int^{\infty}_{0} xe^{-ax^2}dxとなると、
\displaystyle \int^{\infty}_{0} xe^{-ax^2}dx = \frac{1}{2a}
となる。

--------------

5.x^nがかかっているケース(x^n e^{-ax^2})

I_n=\displaystyle \int^{\infty}_{0} x^n e^{-ax^2}dx\hspace{5mm}(n\ge 0)

これをaの関数と見て微分すると、
\displaystyle\frac{dI_n}{da}=-\int^{\infty}_0 x^{n+2} e^{-ax^2}dx=-I_{n+2}
となる。

3,4を見るとわかるとおり、この積分はnが偶数か奇数かで振る舞いが変わる。3がI_0、4がI_1に相当するので、その結果をaで微分することによって具体的なnの値に対して答えが得られる。

低次の場合についてまとめると、以下のようになる。（というかおそらくここが最も求められているところだと思う）。

n I_n 結果
0 \displaystyle\int^\infty_0 e^{-ax^2}dx \displaystyle\frac{1}{2}\sqrt{\frac{\pi}{a}}
1 \displaystyle\int^\infty_0 xe^{-ax^2}dx \displaystyle\frac{1}{2a}
2 \displaystyle\int^\infty_0 x^2e^{-ax^2}dx \displaystyle\frac{1}{4a}\sqrt{\frac{\pi}{a}}
3 \displaystyle\int^\infty_0 x^3e^{-ax^2}dx \displaystyle\frac{1}{2a^2}
4 \displaystyle\int^\infty_0 x^4 e^{-ax^2}dx \displaystyle\frac{3}{8a^2}\sqrt{\frac{\pi}{a}}
5 \displaystyle\int^\infty_0 x^5e^{-ax^2}dx \displaystyle\frac{1}{a^3}

積分が全領域の場合は、偶数次の場合については2倍すればよい。
奇数次のものは、先述のとおり0である。

計算がまちがっていたらコメントしてください。

【python】fortranのnamelistで書かれた設定ファイルを読み込む（f90nml）

python fortran

以前fortranのnamelistという機能を紹介した。

数字を羅列しただけの平文ファイルとちがってnamelistはfortran独自のものなので、ほかの言語で読み込む場合は何らかの工夫が必要になる。

「fortranで計算した結果をpythonのmatplotlibでグラフ化したい」というような場合は多々あるだろう。pythonではそのために「f90nml」というライブラリが用意されている。

インストールについてはpipが使える環境であれば簡単にできる。またgit cloneしてインストールすることもできる。インストールについては以下を参照されたい。

f90nml 0.2.1 : Python Package Index

ただしこのドキュメントには1ヶ所間違いがある。管理者権限のないユーザでも--userをつければインストールできるとあるが、オプションをつける位置が間違っており、正しくは以下のようにする必要がある。

 python setup.py install --user

f90nmlではreadというメソッドでnamelistをリストとして読み込み、writeというメソッドで出力できる。

具体的な使い方として、前回の記事で使ったparams.inのhogeというリストのyyという要素を取得するには、次のようにする。

 import f90nml
params=f90nml.read('params.in')
yy=params['params']['yy']

 

【fortran】設定値を別ファイルで与える（namelist）

fortran

大規模計算プログラムには大量の設定すべきパラメータがある。物理定数のように内部で与えておけばよいものもあるが、頻繁に変えて実行するものもあり、内部で与えるとその都度コンパイルしなければならず効率が悪い。

簡単に思いつくのは、「平文のまま設定値を順番に書き下してread文で読ませる」という方法だ。しかしこれだと設定ファイルが数値（時々文字列や論理値）の羅列になってしまい、可読性に欠ける。

そこで、fortranにはnamelistという機能がある。これはソースファイルとは別の設定ファイルに決まった書き方で書いた設定値を読み書きできる機能で、これを使えば比較的わかりやすい書き方で設定値を与えることができる。

使うには、まずは実行コード側でnamelistを定義する必要がある。これにはnamelist文を使って次のようにする。

 namelist/hoge/ xx,yy,zz

ただしこの文の前で変数xx、yy、zzがそれぞれ定義されている必要がある。*1設定ファイルとはいえ別ファイルなのでopenで開いてreadで読む必要があり、writeで書き出すこともできる。そのためread、writeにはnmlというパラメータが準備されている。

一方で、設定ファイルのほうは次のように記述する。

&hoge
xx=1d0,
yy=2,
zz='hoge',
/
!

注意しなければならないのは、「最後のスラッシュのあとに改行が必要」ということである。どういうわけかgfortranではこれがないと「End of file」というエラーが出るようだ。

ここではhogeというリストひとつしか記述していないが、ひとつの設定ファイルに複数個のリストを書くこともできる。この場合リスト単位で読み書きができる。

まとめて、上記の設定ファイルを「params.in」としてそれを読みとって「out.in」という名前で保存するプログラムを書くと、次のようになる。

program main
  implicit none
  double precision :: xx
  integer :: yy
  character :: zz
  namelist/hoge/xx, yy, zz

  open(10,file='params.in')
  read(10,nml=hoge)

  open(11,file='out.in')
  write(11,nml=hoge)

  close(10)
  close(11)
end program main

 

*1:implicit noneが適用されている場合。fortranプログラミングでは変数の書き間違いによるバグを防ぐためにimplicit noneは強く推奨される。