一种求 π 的迭代构想
zulux
·
·
算法·理论
感谢 @Lexot 对 x_0 范围提出的指正。
有一种思路：
因为
\sin \pi=0
所以求 \pi，可以看成求方程
\sin x=0
的正根。
当然，\sin x=0 有很多根：
0,\pi,2\pi,3\pi,\cdots
如果初值选在 \pi 附近，比如 x_0=3，那么目标就是收敛到 \pi。
这时自然可以想到牛顿法。
对
f(x)=\sin x
使用牛顿法：
x_{n+1}=x_n-\frac{f(x_n)}{f'(x_n)}
因为
f'(x)=\cos x
所以：
x_{n+1}=x_n-\frac{\sin x_n}{\cos x_n}
也就是：
x_{n+1}=x_n-\tan x_n
这个迭代在靠近 \pi 时是二阶收敛的。
所谓二阶收敛，大概意思是：如果当前误差是 e_n，下一步误差大约会变成 e_n^2 量级。也就是说，一旦接近正确答案，正确位数会增长得很快。
但我就想：牛顿法是二阶收敛，那有没有办法构造更高阶的迭代？
一个高阶求 \pi 的想法
我的想法是，不直接对 \sin x 用牛顿法，而是考虑：
g(x)=\frac1{\sin x}
因为当 x\to \pi 时，\sin x\to 0，所以 g(x) 在 x=\pi 附近会趋向无穷。
换句话说，\pi 不是 g(x) 的零点，而是 g(x) 的一个极点。
对这种函数，可以尝试利用高阶导数的比值构造迭代。
我想到的公式是：
x_{n+1}
=
x_n+
(p-1)
\frac{
\left.\dfrac{d^{p-2}}{dx^{p-2}}\left(\dfrac1{\sin x}\right)\right|*{x=x_n}
}{
\left.\dfrac{d^{p-1}}{dx^{p-1}}\left(\dfrac1{\sin x}\right)\right|*{x=x_n}
}
其中 p\ge 2。
如果写得紧凑一点，就是：
x_{n+1}
=
x_n+
(p-1)
\frac{
g^{(p-2)}(x_n)
}{
g^{(p-1)}(x_n)
},
\quad
g(x)=\frac1{\sin x}
并希望有：
\lim_{n\to\infty}x_n=\pi
为什么这个公式看起来有道理？
在 x=\pi 附近，令误差：
e=x-\pi
因为
\sin x=\sin(\pi+e)=-\sin e
而当 e 很小时：
\sin e\approx e
所以：
\sin x\approx -e
于是：
\frac1{\sin x}\approx -\frac1e
也就是说，g(x)=\frac1{\sin x} 在 \pi 附近的主要部分很像：
-\frac1{x-\pi}
而对
h(x)=\frac1{x-\pi}
求导，有：
h^{(k)}(x)=(-1)^k k!(x-\pi)^{-k-1}
于是高阶导数的比值大概会把 x-\pi 提取出来。
这就说明，用
\frac{g^{(p-2)}(x)}{g^{(p-1)}(x)}
这种结构，确实有机会构造出一个把 x 往 \pi 推的高阶迭代。
直观地说，这个公式不是像牛顿法那样看函数的切线，而是利用更多阶导数的信息，试图一次性榨出更多误差项。
特殊情况：当 p=2 时
当 p=2，公式变成：
x_{n+1}
=
x_n+
\frac{
g(x_n)
}{
g'(x_n)
}
其中：
g(x)=\frac1{\sin x}
计算导数：
g'(x)=-\frac{\cos x}{\sin^2 x}
所以：
\frac{g(x)}{g'(x)}
=
\frac{\frac1{\sin x}}{-\frac{\cos x}{\sin^2 x}}
-\frac{\sin x}{\cos x}
-\tan x
因此：
x_{n+1}=x_n-\tan x_n
刚好退化成普通牛顿法。
这点我感觉还挺有意思：也就是说，这个公式可以看作牛顿法的一种高阶推广。
这个公式的优点和问题
这个公式理论上的优点是：
如果 p 越大，它利用的导数阶数越高，局部收敛阶也可能越高。也就是说，在已经离 \pi 比较近的时候，它可能比普通牛顿法收敛得更猛。
但是它也有很现实的问题：
第一，高阶导数不好算。
第二，高阶迭代不等于高效算法。
一个算法快不快，不只看迭代次数，还要看每一步的计算量。如果一步里要算特别复杂的高阶导数，那么总耗时可能反而更大。
第三，初值很重要。
这种迭代应该是局部收敛的，也就是初始值要离 $\pi$ 比较近。比如 $x_0=3$ 这种就比较自然。如果初值乱选，可能跑到别的根，甚至发散。
所以我不敢说它一定能超过 Chudnovsky 这种顶级公式，但从迭代法角度看，它至少是一个挺有意思的高阶构造。
---
## 公式
公式如下：
$$
x_{n+1}
=
x_n+
(p-1)
\frac{
\left.\dfrac{d^{p-2}}{dx^{p-2}}\left(\dfrac1{\sin x}\right)\right|*{x=x_n}
}{
\left.\dfrac{d^{p-1}}{dx^{p-1}}\left(\dfrac1{\sin x}\right)\right|*{x=x_n}
}
$$
并猜想在合适初值下：
$$
\lim_{n\to\infty}x_n=\pi
$$
其中 $p$ 越大，理论上的局部收敛阶越高。
本人初一蒟蒻，公式是闲得无聊乱推出来的。如果有大佬能严格证明它的收敛阶、收敛范围，或者指出它和已有高阶迭代法之间的关系，欢迎指教。
轻喷 qwq。