■フォークト関数の計算（その７）
前回のコラムでは，複素誤差関数を用いたVoigt関数の計算方法に対して，収束半径拡大を目指した漸近展開・連分数展開を施しました．精度については，ベキ級数展開する計算法が優れているのですが，収束域が狭いという欠点があります．収束半径外に対しては，Voigt関数に対する諸値を強制的に0にする対処法も考えられるのですが，それでは面白くないので，漸近展開・連分数展開を用いることにしたわけです．
その場合，複素誤差関数ではうまくいくのですが，超幾何関数に漸近展開を用いると複雑かつ時間がかかりすぎてうまく計算できません．また，計算がかえって複雑になるというのでは，精度の点でも向上が望めません．そこで，今回のコラムでは，プログラムの
1)速度
2)精度
3)収束半径
の性能比較を行ってみます．これらは互いに密接にリンクした問題なので，こっちを立てればあっちが立たず，二兎を追う者一兎も得ずという現象が起こってきます．
＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝
【１】フォークト関数値
ｘ ベキ級数展開 漸近展開 連分数展開 数値積分
0 .208709 × .208830 .208328
1 .165796 × .165761 .165640
2 .090715 × .090701 .090720
3 .043386 × .043388 .043390
4 .022814 × .022814 .022814
5 .013885 × .013885 .013885
6 9.33768D-3 × 9.33768D-3 9.37768D-3
7 6.77860D-3 6.77860D-3 6.77860D-3 6.77860D-3
8 5.22479D-3 5.13668D-3 5.13668D-3 5.13664D-3
9 × 4.03053D-3 4.03053D-3 4.03053D-3
10 × 3.24874D-3 3.24873D-3 3.24873D-3
11 × 2.67525D-3 2.67525D-3 2.67525D-3
12 × 2.24184D-3 2.24184D-3 2.24184D-3
13 × 1.90618D-3 1.90618D-3 1.90618D-3
14 × 1.64086D-3 1.64086D-3 1.64086D-3
15 × 1.42746D-3 1.42746D-3 1.42746D-3
20 × 7.99783D-4 7.99784D-4 7.99784D-4
30 × 3.54466D-4 3.54466D-4 3.54466D-4
40 × 1.99193D-4 1.99193D-4 1.99193D-4
50 × 1.27426D-4 1.27426D-4 1.27426D-4
100 × 3.18374D-5 3.18374D-5 3.18374D-5
＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝
ベキ級数展開（テイラー展開）に基づく方法は，ｘの小さいところで正確ですが，ｘが少し大きくなると発散してしまいます．f(z)=exp(x^2)erfc(z)という式を使っているので，ｚが大きくなると桁あふれしてしまうのです．逆に，漸近展開（ローラン展開）はｘの大きいときに正確な値がでますが，ｘが小さいときは発散しました（×印）．この収束域はプログラムの収束判定基準（この場合最適打ち切り項数の設定）によって多少変わってきますが，基本的には両者はまったく正反対の結果となることが示されます．
その点，ラプラス・ヤコビの式：
exp(z^2)Erfc(z)=1/2z 1/2z^2 1/z^2 3/2z^2 2/z^2
1+ 1+ 1+ 1+ 1+
を利用した漸近展開（連分数変換）は広い収束域をもち，連分数展開によって収束半径の制限は解消されました．このように，連分数展開は元になる漸近展開（ローラン展開）と比べて，収束半径が大きいことが特長です．
一般にデジタル数値計算では，連分数展開は誤差の累積がひどく，実用にならないと評価されてしまったようですが，この例のように，漸近展開など収束半径が０である形式的ベキ級数を連分数変換すると，広い収束域をもつので実用上大変有用と考えられます．
もちろん，ラプラス・ヤコビの式は漸近展開を連分数に直したものであって，本質的に漸近展開であることには変わりありませんから，ｘの小さいところでの精度はベキ級数展開より劣ります．しかし，収束半径が大きいというメリットは何物にも代え難いし，捨てがたいものがあります．
一方，ガウス・エルミート則の数値積分を１６点法で行ってみたところ，ｘの小さいところでの精度はベキ級数展開よりも劣りますが，収束半径が大きく，漸近展開は不要と考えられました．というよりも，元々，この方法には収束半径の問題が存在しないのです．そのため，複数のバンドがあるケースのあてはめに都合がよく，多少の精度の点について目をつぶることにすれば，及第の値が算出されます．また，ガウス・エルミート則には，速度が圧倒的に速いというアドバンテージもあります．
＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝
【２】ヤコビアンの値（df/dσ）
ｘ ベキ級数展開 漸近展開 連分数展開 数値積分
0 -.099109 × -.099351 -.101466
1 -.031254 × -.031359 -.033547
2 .034482 × .034512 .033373
3 .027155 × .027175 .027103
4 .010177 × .010177 .010177
5 4.07235D-3 × 4.07233D-3 4.07224D-3
6 1.80006D-3 × 1.80006D-3 1.80005D-3
7 9.21548D-4 9.19110D-4 9.19108D-4 9.19108D-4
8 × 5.20027D-4 5.20047D-4 5.20047D-4
9 × 3.17027D-4 3.17046D-4 3.17046D-4
10 × 2.04562D-4 2.04567D-4 2.04567D-4
11 × 1.38015D-4 1.38025D-4 1.38025D-4
12 × 9.65579D-5 9.65616D-5 9.65616D-5
13 × 6.95901D-5 6.96092D-5 6.96092D-5
14 × 5.14531D-5 5.14621D-5 5.14621D-5
15 × 3.88709D-5 3.88755D-5 3.88755D-5
20 × 1.21296D-5 1.21387D-5 1.21387D-5
30 × 2.34264D-6 2.35744D-6 2.35744D-6
40 × 7.43326D-7 7.49159D-7 7.49159D-7
50 × 3.04866D-7 3.06394D-7 3.06394D-7
100 × × 1.91113D-8 1.91113D-8
＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝
漸近展開はローラン級数形式になっているのですが，ｘが大きくなるとやがて計算がおかしくなってしまいました．その点，連分数変換したものは収束半径に問題はなさそうです．
また，ヤコビアンも連分数展開を用いて計算することにしたいのですが，一般的にはこれはうまくいきません．しかし，複素誤差関数の場合は，うまい手が利用できます．
{exp(z^2)erfc(z)}'=2zexp(z^2)erfc(z)-2/π
これを使えば，複素誤差関数を用いるVoigt関数の漸近計算では，ヤコビアン計算の速度を速くすることが可能であって，速度についてはこれ以上速いアルゴリズムは望めないほどです．しかも，収束半径が大きく，数値微分よりもずっと精度の良い値が得られました．
＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝
【３】ヤコビアンの値（df/dβ）
ｘ ベキ級数展開 漸近展開 連分数展開 数値積分
0 -.109600 × -.109479 -.106862
1 -.060630 × -.060659 -.059685
2 5.03448D-3 × 4.99098D-3 4.88441D-3
3 .022504 × .022502 .022460
4 .017968 × .017969 .017966
5 .012370 × .012370 .012370
6 8.74356D-3 × 8.74331D-3 8.74331D-3
7 6.46103D-3 6.46113D-3 6.46113D-3 6.46113D-3
8 6.16704D-3 4.95890D-3 4.95890D-3 4.95890D-3
9 × 3.92281D-3 3.92281D-3 3.92281D-3
10 × 3.17952D-3 3.17952D-3 3.17952D-3
11 × 2.62868D-3 2.62868D-3 2.62868D-3
12 × 2.22093D-3 2.22093D-3 2.22093D-3
13 × 1.88278D-3 1.88278D-3 1.88278D-3
14 × 1.62358D-3 1.62358D-3 1.62358D-3
15 × 1.41442D-3 1.41442D-3 1.41442D-3
20 × 7.95722D-4 7.95724D-4 7.95724D-4
30 × 3.53667D-4 3.53673D-4 3.53673D-4
40 × 1.98942D-4 1.98943D-4 1.98943D-4
50 × 1.27323D-4 1.27324D-4 1.27324D-4
100 × 3.18214D-5 3.18310D-5 3.18310D-5
＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝
【４】積分値（−∞〜ｘ）
ｘ 超幾何関数 クンマー変換 複素誤差関数 漸近展開
0 .5 .5 .5 ×
1 .693560 .693560 .693560 ×
2 .820959 .820959 .820959 ×
3 .885114 .885114 .885114 ×
4 .916717 .916717 .916717 ×
5 .934508 .934508 .934508 ×
6 .945914 .945914 .945914 .945914
7 .953885 .953885 .953895 .953885
8 .959791 .959785 .959776 .959785
9 .982600 .964334 × .964334
10 × .967953 × .967953
11 × .970901 × .970900
12 × .973350 × .973350
13 × .975417 × .975417
14 × × × .9771855
15 × × × .9787160
20 × × × .9840578
30 × × × .9893818
40 × × × .9920389
50 × × × .9936321
100 × × × .9968167
＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝
積分計算についても，収束半径の大きい計算法が望まれます．複素関数：f(z)=exp(z^2)erfc(z)の不定積分は
∫f(z)dz=-z^2/√π2F2(1,1,3/2,2,z^2)+z1F1(1/2,3/2,z^2)
で与えられますが，この式は実際の数値計算には有用ではありません．それよりも
∫f(z)dz=Σ(-1)^n*z^(n+1)/Γ(n/2+1)/(n+1)
の方がはるかに有用なのです．
複素関数：f(z)=exp(z^2)erfc(z)の漸近展開の連分数展開はラプラス・ヤコビの式：
exp(z^2)Erfc(z)=1/2z 1/2z^2 1/z^2 3/2z^2 2/z^2
1+ 1+ 1+ 1+ 1+
があるため解決できたのですが，しかし，複素誤差関数の積分
∫1/z2F0(a,b:-1/z^2)dz=logz+ab*(1/z^2/2)+a(a+1)b(b+1)/2!*(-1/z^4/4)+...
の連分数展開は簡単な式には収まりそうもありません．
そのため，積分計算の収束域を拡大させることができず，小さいところで精確な値がでる算法と大きいところで精確な値がでる算法を使い分けることが必要となります．結局，ｘの小さいところと大きいところに分けて計算せざるを得ないのですが，その境界をどうやって決めたらよいかという問題が残ることになり，どこから漸近展開に切り替えるべきかという検討も必要となります．
＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝
【５】まとめ
積分表示のままで要求される数値積分（この場合はガウス・エルミート法がベスト）を行い，その精度でうまくいく場合もあるのですが，このプログラムを作った所期の目的は，Voigt関数を解析的に表現することによって，Voigt関数の関数値，ヤコビアンともに高精度で計算することです．
複素誤差関数経由でVoigt関数の値を計算する方法は，連分数展開を用いることによって，収束域が広くなり収束半径の問題をクリアすることができました．しかし，この連分数展開にはｘが小さいところの精度が数値積分並みという欠点があります．
たった１つの方法であらゆる場合の計算が可能であればそれに越したことはないのですが，それはかなり難しそうです．したがって，Voigt関数の高精度計算では２つの方法を組み合わせること，すなわち，ｘの小さいところで正確な方法（ベキ級数）とｘの大きいところで正確な方法（漸近展開＋連分数変換）を切り替えることが現実的であり，かつ，正当ではないかと考えられました．どこをその境にするかは問題があるのですが，・・・
実は，このような分割計算法を採用しているものに，
The collected algorithms from ACM
http://www.acm.org/calgo
の680があります．ACM/680の方法では，境界を調整するアルゴリズムが記載されていますので，ぜひご参照下さい．
ここで作成したプログラムは，ACM/680と比較してもヤコビアンの精度が高いので，精度面では圧倒的に有利です．速度の点でも，ACM/680と同程度の時間でした．これまでは，near Gaussian では，シンプレックス法でパラメータの最適化を行っていたのですが，ヤコビアンが正確に計算できるようになると，ニュートン法に切り換えることが可能になりますから，劇的な改善が望めそうです．
また，数値積分であるガウス・エルミート法は予想していたよりも精確であって，収束半径の問題もないようです．最初にこれで計算してから解析的な方法に切り替えるのも１つの手かもしれません．
＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝
