流れる空の中で数学を。

とある数学好きの「手作りすうがく」と「気ままな雑記」。

exp(x)を多項式で近似したときの零点の話

プチ研究数学プログラミング追記あり、一部未完成 twitter オリジナル

Amazon.co.jp: Yoshiki Ueoka:作品一覧、著者略歴

数学関連の絶版本・品切れ本をコチラから購入できます！

twitterで見かけた指数関数の話

　twitterでtsujimotterさんが、指数関数について以下のようなツイートをしていた。

ふと思った疑問。指数関数のテイラー展開って、有限項で打ち切った多項式の零点と指数関数の零点ってどういう関係になっているだろう。
— tsujimotter (@tsujimotter) September 30, 2018

　これが少し気になったので、Wolfram alphaを使って10次の多項式まで計算してみたところ零点が放物線っぽく並んでいるのを見つけた。もっと高い次数まで計算してみたかったのだが、Wolfram alpha(無料版)は10次までしか相手してくれなかった。

　そこで、勉強中のPythonを使って計算してみることにした、というのが今回するお話です。

指数関数の多項式近似とスケーリング

　指数関数を $N$ 次の多項式で近似すると、次のようになる。

$exp(x) \simeq \sum_{n=0}^{N} \frac{x^n}{n!}$ ……(1)

　PythonのNumpyには、多項式の零点を計算してくれる便利な関数numpy.rootsがあるので今回はそれを使うことにした。ところが、このまま計算すると、 $N$ 次の項が小さくなり過ぎて桁落ちする恐れがあるので、スターリングの公式(スターリングの近似 - Wikipedia)を参考にして、

$x=\frac{Ny}{e}$ ……(2)

と置き換えて、変数のスケールを事前に変えておく。つまり、

$\sum_{n=0}^{N} \frac{(N/e)^n}{n!} y^n$ ……(3)

という多項式の零点を数値計算する。こうすることで、 $N$ 次の項が極端に小さくなることを回避できる。このスケールした多項式の零点が $y_0$ と求まったら、元の多項式(1)の零点は、

$x_0=\frac{N}{e} y_0$ ……(4)

と求まることになる。

　大量の零点を計算することになるので、ちゃんと計算できているかどうかのチェックが必要になる。そこで、求まった零点からそれらの零点を持つ多項式を逆に再計算して、定数項を1になるように規格化し、元の多項式と比較を行った。今回の記事に載せる計算結果は、すべて、元の多項式と零点から復元した多項式の最大の係数が、有効数字で8桁以上で一致することを確認してある。

零点の分布

　まずは、次数 $N$ が小さい方から多項式(1)の零点の分布をみていく。まずは、Wolfram Alphaでやったやつの再計算結果から。

f:id:FoxQ:20181005111524p:plain

f:id:FoxQ:20181005111547p:plain

f:id:FoxQ:20181005111555p:plain

f:id:FoxQ:20181005111600p:plain

このように、無事再現できた。

　次は、次数を20から60まで一気に増やしていく。 $N=60$ に近づくにつれて、放物線というよりは、馬蹄形のようになっていることがわかった。

f:id:FoxQ:20181005111740p:plain

f:id:FoxQ:20181005111747p:plain

f:id:FoxQ:20181005111752p:plain

　さらに、次数を増やす。すると妙なことが起こった。
追記(2018/10/07)：後日、これ以降の次数の数値計算に数値誤差があることが判明しました。きちんと計算できているのは、 $N=60$ くらいまでのようです。

f:id:FoxQ:20181005111908p:plain

　謎の分岐が発生した。さらに、次数を増やしていく。

f:id:FoxQ:20181005111952p:plain

f:id:FoxQ:20181005112006p:plain

f:id:FoxQ:20181005112032p:plain

f:id:FoxQ:20181005112040p:plain

f:id:FoxQ:20181005112051p:plain

　 $N=170$ まで計算したところで、多項式の復元時の誤差が有効数字8桁程度になったり、 $N=180$ でプログラムがエラーを吐いたのでここで計算を打ち切った。

　最終的に、零点の分布が描いたのは、二つの弧とその間にぽつぽつとした点があり、それらの弧に二本のツノを引っ付けたような図形となった。
最初の放物線っぽい形から予想できるような、単純な形にはならなかったのである。

零点の絶対値の範囲

　また、twitterでtaketo1024さんが(1)の多項式の零点の絶対値の範囲の見積もりがあることを教えてくれた。

気になって調べたら、このノートに指数関数の n 次近似式の零点は n/e^2 < |z| < n なる範囲にあることの証明が載ってました😆https://t.co/JLUzrBQe2M
— さのたけと (@taketo1024) September 30, 2018

　ノートによると、この不等式は $N$ が十分大きいところで成立するらしい。

　そこで、今回私が計算した零点の絶対値をとり、その最大値と最小値をプロットしてみた。青色の線が零点の絶対値の最小値で、赤色が最大値である。また、ピンクと水色の線は、上記のツイートのノートに書いてあった最大値と最小値である。

$N=100$ までの図をひとまず見ると、不等式の範囲内にあることがわかる。

f:id:FoxQ:20181005113053p:plain

　ところが、 $N=170$ までを見ると、不等式の下限を下回ってしまった。

f:id:FoxQ:20181005113106p:plain

　これが、numpyの限界なのか、 $N$ の大きさが十分でないためなのかは、はっきりとしない。現状、 $N=60 \sim 70$ 付近で、なぜか零点の絶対値の最小値が頭打ちになっているようにも見える。これはちょうど、零点の分布の図で２個目の『弧』が現れ始めたところである。

追記(2018/10/07):この最小値が頭打ちになっているあたりから数値誤差が強く出ているようです。正しい零点の極限分布について、以下のサイトにまとめを見つけたので、そちらを見てください。

数学の散歩道の「３．指数関数が定める多項式のゼロ点集合の極限形」

http://www.eonet.ne.jp/~kotani/sanpo.html

　このような現象がなぜおきるのか？数学的にありうることなのか？単に数値計算による誤差などの限界から来るものなのか？

　いずれもよくわからない、というのが今の感想だ。みなさんは、どう思いますか？