少しは研究のマニアックな話があっても良いかと思い、物理の数式が出てくる話題を書きたいと思います。実験研究者ですが、実験結果が物理法則の予言通りになっているか確認することが必要です。ブラウン運動の研究では、基本的に確率過程を考えていくこととなり、そんな実験において、実験結果をFokker-Planck方程式で考えるケースが出てきました。
$$\frac{\partial P(z,t)}{\partial t}=\frac{\partial}{\partial z} (\frac{mg}{\gamma}P(z,t)+\frac{D}{2\gamma^2}\frac{\partial P(z,t)}{\partial z})$$
Fokker-Planck方程式の細かな説明は置いておき、簡単に言えば、確率微分方程式で、ブラウン運動などの確率過程の時間変化を追いかける方程式になります。この方程式は解くことが出来て、解析解は次のようになります。（初期条件などの説明は省略）
$$P(z,t)=\frac{2\alpha}{1-e^{-2\alpha H}}e^{-2\alpha z}+\sum_{n=1}^{\infty} \frac{(e^{\alpha H} cos(n\pi)-1)4\alpha n^2 \pi ^2}{(\alpha ^2 H^2 + n^2 \pi ^2)^2} [cos(\frac{n\pi z}{H})-\frac{\alpha H}{n\pi}sin(\frac{n\pi z}{H})]e^{-(\gamma_n t+\alpha z)}$$
無限級数の形式であるが、解析解があるので、実験条件の数値を代入しチマチマと計算していると、うまく計算できない場合が発生。最初はエクセルで計算していたから、エクセルのせいにしていたのだけどPythonを使って、ちゃんと計算しても結果がスムーズにならない。（物理的には確率の計算なので、滑らかな曲線にならないと、どこかがオカシイ）
いろいろ試行錯誤しているとプログラムで扱っている変数に問題があることが分かった。Pythonで小数を取り扱うと倍精度浮動小数(double)という精度で計算することになる。コンピュータの中は２進法で数値を取り扱うのだが、１０進法で考えると、約１５～１６桁の桁数で小数を取り扱うことになる。解析解をみると総和の記号Σが出ているが、そこで、足し算している最初の数字が１７桁近くになることが判明し、倍精度浮動小数では、四捨五入が行われ、十分な精度で計算できないことが分かったのである。
ここからは、完全にエクセルでの計算をあきらめ、Pythonでの数値の取扱いについて勉強です。ここで、Pythonの良いところを一つ知りました。整数の取扱いに関して、そのサイズに制限がないのです。例えば、１から１００まで掛け合わせると（１００！＝）math.factorial(100)
93326215443944152681699238856266700490715968264381621468592963895217599993229915608941463976156518286253697920827223758251185210916864000000000000000000000000
と求めることも簡単にできるのです。
次に、小数・指数の取扱いについて、Pythonで手法を検索して、数値の桁数を増やす方法として
from decimal import getcontext, Decimal
getcontext().prec = 20
とdecimalモジュールを使って、２０桁の計算を実行させることにしました。ただ、このdecimalモジュールで計算の有効桁数を増やすことはできるのだけど、sin関数やcos関数をnumpyのモジュールで計算することができないようである。（numpyの計算では、倍精度浮動小数に固定されているみたい。）解決策として、decimalのドキュメントから、decimalで動作する独自関数を作成した・・・・
しかし、ソースコードをよく見ると、単にテーラー展開を計算しているだけであった。
$$sin x=x-\frac{1}{3!}x^3+\frac{1}{5!}x^5-\dotsb$$
$$cos x=1-\frac{1}{2!}x^2+\frac{1}{4!}x^4-\dotsb$$
この展開を２０桁（decialで規定した桁数）まで数値が変化しなくなるまで続けているだけだった。普段から学部１～２年生にテーラー展開は大事だよと言い続けているのだけど、それを再確認できたプロフラムコードでした。
