√2の肩に無限に√2を乗せると2になる振る舞いを matplotlib でリアルタイム描画する

:::note info
これは [ZOZO Advent Calendar 2022](https://qiita.com/advent-calendar/2022/zozo) カレンダー Vol.5の17日目の記事です。昨日の投稿は @cease2xist さんの 「[AWS Lambda SnapStartを有効にする時にハマったポイント](https://zenn.dev/papadino/articles/282fc66d5c123332824e)」でした。
:::

# 概要
√2の肩に無限に√2を乗せると2に収束することが数学的に証明されています[1, 2]。数式で書くと下記の式(1)になり Infinite Power Tower などと呼ばれています。本記事では、 Python の描画ライブラリ matplotlib のリアルタイム描画機能を使って Infinite Power Tower の収束の振る舞いを視覚的に示します。

```math
\def\rddots#1{\cdot^{\cdot^{#1}}}
{{{\sqrt{2}}^{\sqrt{2}}}^{\sqrt{2}}}^{\rddots \cdot} = 2  \tag{1}
```

# なぜ√2の肩に無限に√2を乗せると2になるのか？
一旦、数学的な厳密性を考えないで式(1)の直感的な"式変形"を行います。それは
```math
\def\rddots#1{\cdot^{\cdot^{#1}}}
s = {{{\sqrt{2}}^{\sqrt{2}}}^{\sqrt{2}}}^{\rddots \cdot} \tag{2}
```
と定義して関係式
```math
s = {\sqrt{2}}^{\large s} \tag{3}
```
を得るというものです。この"式変形"は 「$s$ では ${\sqrt{2}}$ を無限に肩に乗せる操作をしているのだから(無限 - １)回肩に乗せる操作も同じ $s$ として扱えるだろう」という意図で行っています。この式(3)は容易に解くことができて $s = 2, 4$ という解が得られます。しかし、そもそも $s$ は式(2)で定義される何かしらの値なので $s=2=4$ となり、明らかに間違った結果になります。この原因は無限の操作を厳密に取り扱わないで式(3)を導入したことにあります。

そこで、無限の操作を厳密に取り扱うために以下の数列を導入します。
```math
a_1 = {\sqrt{2}}, a_{n+1} = {\sqrt{2}}^{\large a_n}  \tag{4}
```
つまり、この数列 $a_n$ の極限が $\lim_{n{\rightarrow}{\infty}}a_n = 2$ であることを示すことが式(1)を証明することになります。詳細な証明は省くのですが、式(3)に対して数学的帰納法のアプローチを取ることで「$a_n$ が単調増加数列」「$a_n {\leq} 2$」の条件を導くことができます。つまり、この2つの条件から数列 $a_n$ の極限が2であることが証明され、式(1)を証明したことになります。詳しい証明は参考文献[1, 2]を参照して下さい。

# 一般化した Infinite Power Tower
Infinite Power Tower は以下のように一般化して考えることもできます。
```math
\def\rddots#1{\cdot^{\cdot^{#1}}}
{\large x^{x^x}}^{\rddots \cdot} = y  \tag{5}
```
実はこの式は $x$ が $e^{-e} {\leq} x {\leq} e^{1/e}$ （$e$ はネイピア数)の範囲で収束することが Leonhard Euler によって証明されています[2]。この収束条件はおよそ $0.06 {\lesssim} x {\lesssim} 1.44$ となります。そこで、本記事では、式(1)だけでなく、 $x=0.5$ の場合（式(6)）の収束についても数値計算の力を借りて見ていきます。
```math
\def\rddots#1{\cdot^{\cdot^{#1}}}
{{0.5}^{{0.5}^{0.5}}}^{\rddots \cdot} = y   \tag{6}
```

# 数列の収束の振る舞いをリアルタイム描画
## リアルタイム描画を利用する理由
もし式(4)のように事前に収束値がわかっており、単調増加数列であることも手計算でわかっている場合はその振る舞いを描画することは簡単です。例えば、収束値に近い値になった時点で数値計算を止めるようプログラムし、その結果を描画することで全体の振る舞いを確認することができます。一方で、一般化した式(5)の場合は事前に収束値を求めることが難しく、単調増加数列であるかについても自明ではないために数値計算を自動で止めるようプログラムすることは難しいです。
そこで、数列の振る舞いをリアルタイムに描画し、十分に振る舞いを確認できたと思う任意のタイミングで描画処理を手動で止められる実装を行うことにしました。

## matplotlib を使ってリアルタイム描画
Python の matplotlib を使うことで計算結果をリアルタイムに描画することができます[3]。そこで、今回はこの機能を使って数列の繰り返し計算結果をグラフに描画しました。具体的な実装は以下です。
```python
import numpy
import matplotlib.pyplot as plt


# 初期値の設定
start_idx = 1                 # インデックスの最初の数字
a_1 = 2 ** 0.5                # 初項: 本記事では √2 or 0.5 で実験
x = numpy.array([start_idx])  # 描画するx座標リストの初期値
y = numpy.array([a_1])        # 描画するy座標リストの初期値

# リアルタイム描画用設定
fig, ax = plt.subplots()
(lines,) = ax.plot(x, y)

n = start_idx
while True:
    if n == 1:
        a_n = a_1
    else:
        x = numpy.append(x, n)
        y = numpy.append(y, a_n)

    # 数列の更新
    a_n = a_1 ** a_n

    # matplotlib のリアルタイム描画
    ax.set_xlim((x.min(), x.max()))  # x座標の描画範囲
    ax.set_ylim((y.min(), y.max()))  # y座標の描画範囲
    lines.set_data(x, y)
    # どの値に収束しているかわかり易いように数列の最終インデックスの値を描画
    ax.set_title(f"a_{n}: {str(y[n-1])}")
    plt.pause(0.2)  # 描画間隔[sec]

    # インデックス更新
    n += 1

    # 安全のためにn=300で強制終了
    if n == 300:
        break
```

描画結果は以下です。縦軸が $a_n$, 横軸が $n$ に対応していて $n$ が増加するに従って $a_n$ が 2 に収束していることがわかります。加えて、収束の早さについても視覚的に簡単に把握することができます。
![sq2.gif](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2953210/adbcac48-a9c7-c801-af8b-6817c4cf4c4c.gif)


また上記のコードの初項 $a_1$ を 0.5 にすることで式(6)の振る舞いも確認しました。結果は以下で振動しながら 0.64118 あたりに収束していることがわかります。
![05.gif](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2953210/c66e9234-a404-46d3-9b81-083b4a4affce.gif)
実際に $0.5^{0.64118} {\sim} 0.64118$ となり数値計算の結果が妥当であることがわかります。

# まとめ
Infinite Power Tower という手計算だけでは理解することが難しい極限的な数列の振る舞いを matplotlib のリアルタイム描画によって視覚的に示すことができました。ぜひ皆さんも上記のコードの初項 $a_1$ を自由に変えて気になる Infinite Power Tower を描画して遊んでみて下さい。ただし、数値誤差が含まれるので収束値が極端に小さい場合などは正しい解を得られない可能性がありますのでご注意下さい。


# 参考資料
[1] [YouTubeチャンネル 予備校のノリで学ぶ「大学の数学・物理」  『√2の肩に無限に√2を乗せたらなぜ2になるのか』](https://youtu.be/LABUpfm-rF0)
[2] [Luca Moroni 『The strange properties of the infinite power tower』 arXiv:1908.05559](https://arxiv.org/abs/1908.05559)
[3] [matplotlibでリアルタイム描画](https://qiita.com/hausen6/items/b1b54f7325745ae43e47)

※ 上記、[2]はあくまで参考資料であって原著論文という訳ではありません。