民科吧 关注:312,831贴子:4,718,478

祖冲之割圆术的简易改进

只看楼主收藏回复

一楼祭天


IP属地:澳大利亚来自iPhone客户端1楼2024-03-18 02:47回复
    在复现了割圆术及对应数学原理后,可以通过简易的一阶差分观察对已有的结果序列进行微扰,从而提高割圆术的收敛速度。泰勒展开为微扰方法提供数学原理支持。


    IP属地:澳大利亚来自iPhone客户端2楼2024-03-18 03:47
    回复
      会有人对这个感兴趣吗?昨天无聊的时候研究了一下,算是一个原创不过没什么实际价值的研究。


      IP属地:澳大利亚来自iPhone客户端3楼2024-03-18 03:49
      回复
        众所周知,割圆术的原理是通过圆内接正多边形的周长与圆半径之比来近似描述π。
        由于三角不等式的存在,正m*n边形的周长严格大于正n边形的周长。基于此,我们可以构造一个数列,得到π越来越紧的下界。


        IP属地:澳大利亚来自iPhone客户端4楼2024-03-20 02:10
        回复
          割圆术可以以任何正多边形为起点,但有些对于计算而言显得更友好。三角形是一个好的起点:单位圆内接正三角形的边长是可以简易表示的,即3^0.5。当然,也可以从正六边形开始,它的边长甚至是不带根号的1。那么我们从正三角形开始,以2为公比开始割圆,割圆数列p0的前两个值是3^1.5/2和3。
          让数列开始的方法是简单的,需要的工具有两个:余弦定理和余弦半角公式。我们希望计算第n项时可以利用第n-1项的信息,这不仅有利于简化计算,在之后的改进中也会起到关键的作用。
          我们有余弦公式的特殊形式cos(2x)=2cos(x)^2-1,换元可得cos(2/x)=((1/2)*(1+cos(x)))^0.5。构建一个有关余弦的数列a_n,有递推公式
          a_(n+1)=((1/2)*(1+a_n))^0.5。
          然后我们需要余弦公式。对于我们的情况,它写作c=(2*(1-cos(x))^0.5,这样我们就描述了正n边形的边长,由此,结合之前的递推公式和初始条件,我们有:
          p_n=3*2^(n-2)*(2*(1-a_n))^0.5。
          代n=1 a_1=-1/2入得p_1=3^1.5/2,p_2=3。验算正确。
          根据之前的分析,我们得到π的精确表达式:π=lim(n–>∞)p_n。这是由于我们已经知晓了π的存在性与数列的单调性得出的结论:单调有界必收敛。


          IP属地:澳大利亚来自iPhone客户端5楼2024-03-20 02:29
          回复
            每一项与π的精确值之差如第一张图所示。第二张图为对数图,良好的线性证明它的收敛速度为指数,即计算N位精度的π的时间复杂度为O(N)。
            这里必须强调:我们默认计算每一项需要的时间相等,这也是为什么之前给出了递推公式。另外涉及开根的方法并不适合计算机进行高精度运算,这里仅为学术讨论。
            图中给出的最后一项是p0_13的误差,除表示误差的蓝线之外,还给出了用橙线表示的通项差。如果对古中国数学史有了解,应该知道这里就是祖冲之计算的终点:他把圆周率的精度提高到了3.1415926…。虽然我并不清楚他是否严格证明了3.1415927是圆周率的上界,但后来的故事说明了这是正确的。
            以上是祖冲之的工作。
            而接下来的分析,则会说明他本有机会使用并不大太多的计算量去得到更高精度的圆周率——当然建立在割圆术上,前提是它愿意把根开到他想要的位数上。
            我画出了通项差,它似乎也和数列的误差一样,随着n的线性增大而呈指数缩小,甚至好像有一样的底数——它们是平行的。我们能否利用这个性质在已有的数列上改进割圆术的结果?


            IP属地:澳大利亚来自iPhone客户端6楼2024-03-20 02:54
            收起回复
              最近看量子化学的微扰理论,想到了这个,就捣鼓了一下。不过感觉这种没什么实用价值只是比较有趣的东西还是发到民科吧比较合适


              IP属地:澳大利亚来自iPhone客户端7楼2024-03-20 03:05
              回复
                如果有人感兴趣的话,也可以证明一下这个方法的收敛速度确实是指数。不过真的有人喜欢听数学证明吗


                IP属地:澳大利亚来自iPhone客户端8楼2024-03-20 03:13
                回复
                  自己顶一下


                  IP属地:澳大利亚来自iPhone客户端9楼2024-03-20 06:47
                  回复
                    那么,把差分值与误差值作商,会得到什么呢?运算的结果是,它会越来越接近1/3,而且收敛速度非常快。
                    那么我们自然而然地可以想到一个微扰:把p0_(n+1)加上它与p0_n的差的1/3,是否可以提高数列的收敛速度?
                    依照这个想法,我们构建一个新的数列p1,它的通项是:
                    p1_n=p0_n+(p0_n - p0_(n-1))/3
                    这个定义会存在问题,因为不存在p0_0,导致p1_n没有定义。所以,修正定义为:
                    p1_n=p0_n+(p0_n - p0_(n-1))/3 if n>1
                    else: p1_n=p0_n。
                    我们可以对数列进行收敛速度的分析,图中的绿线就是p1的表现。可以看出,它的收敛速度比p0(也就是祖冲之的原始割圆术)快得多,是原本的2倍。这意味着p1_7的精度就已经达到了3.1415926,相较于p0_13可以少开一半次数的根号。这意味着只需要割到正216边形(p0_6)就够了,而不是至少正12288边形(p0_12)(实际割到了24576边形(p0_13))。


                    IP属地:澳大利亚来自iPhone客户端10楼2024-03-20 18:52
                    收起回复
                      那这里就是终点了吗?我们能不能如法炮制一下?继续对p1序列的误差与每步增长作商,发现得到的值越来越接近1/15(这里不包括p1_1与p1_2的变化率,因为p1_1与p0_1是一样的,一阶微扰没有改良p0_1)。于是我们继续写:
                      p2_n=p1_n + (p1_n - p1_(n-1))/15 if n > 2
                      else: p2_n=p1_n
                      不出所料地,p2比p1的收敛速度更快,已经是p0的三倍。p2_5就已经达到了3.1415926(5)的精度,这意味着只需要割到正32边形!须知割圆术的计算量来源于开根,开根不是基础运算,需要迭代求解。而这个方法可以把开根转化成加减乘除,客观上会大幅度提高运算效率。


                      IP属地:澳大利亚来自iPhone客户端11楼2024-03-20 19:09
                      回复
                        没有人理本民科吗


                        IP属地:澳大利亚来自iPhone客户端12楼2024-03-20 19:35
                        回复
                          沉了啊


                          IP属地:澳大利亚来自iPhone客户端13楼2024-03-20 21:45
                          回复
                            挺好,但是我看不懂


                            IP属地:北京来自Android客户端14楼2024-03-21 00:40
                            收起回复
                              显然,这样的微扰可以继续下去。对于计算到第n项的p0数列,最多可以微扰到n-1阶,即我们在计算了前n项后,可以期望获得的π的有效位数是O(n^2/2),其中需要进行2n-1次开根计算(除a1外的余弦值迭代和使用余弦值计算边长)和n(n-1)/2组简单的四则运算(微扰方法的作差,乘上系数和求和)。
                              显然,改良方法的性质太好了:它不需要进行额外的开根运算,可以通过简单的加减乘除去获得π越来越紧的下界,效率甚至有拉马努金的影子(当然,拉马努金不需要开2n-1次根,而且效率应该还是高一些)。于是我们在为祖冲之没有发现此改良而惋惜的同时不免要问:为什么我们可以如此简单粗暴地对原始割圆术进行优化?它背后有着怎样的数学原理?


                              IP属地:澳大利亚来自iPhone客户端15楼2024-03-21 03:33
                              回复