连分数近似算法

Maths

      这道题源自数学实验上面的一组实验，当时困扰了我特别久，题目的内容是用matlab求出π的连分数展开及每层迭代的值。
      因为matlab的数值精度的问题，当你运行3+16450/16421时，就算你设置了format rat，也没有办法显示准确的有理分数的值。原书是用mathmatica做的，直接一句命令就出来了。
      其实matlab里面也是有相应的连分数展开的命令的，详细可以查看以下：

Continued fraction expansion

      但是显然精度高时，就会遇到那天叙森写的学习笔记的另外一面，精度要求太高，函数无法满足要求，举个例子。

1
2
3
4
5
6
>> rat(pi,1e-31)
ans =
3 + 1/(7 + 1/(16 + 1/(-294 + 1/(3 + 1/(-4 + 1/(5 + 1/(-15)))))))
>> rat(pi,1e-18)
ans =
3 + 1/(7 + 1/(16 + 1/(-294 + 1/(3 + 1/(-4 + 1/(5 + 1/(-15)))))))

      其实后来我测试自己的算法也是有问题的，但是至少比之前迭代不到10次就直接误差为0了要好得多，究其原因还是因为前面运算的时候用的是double型数，这里我尝试过，如果在一开始就使用符号运算的话，不仅运算时间无法忍受，而且很容易迭代到20次以后就直接程序堆栈溢出或者NaN或者是Inf了。
      好了，说了这么多，该介绍下连分数是什么了，其实也不难，据说拉马努金对这个非常的了解。

$$
S_{n}=
\begin{equation}
\cfrac{1}{a_1 +
\cfrac{1}{a_2 +
\cfrac{1}{a_3 + \dotsb}}}
\end{equation}
$$

      为了算法的简便，我们设定的是这组数列全部都是正整数。对于有理数，我们知道，这组数列一定是有限数列，对于无理数，这组数列大部分情况下是无限数列。我们可以从某一个数开始，设定它为0，按照这样子的顺序，加出来一个有理分数，这个有理分数就是我们需要的值，举个例子。

$$
π=3+
\begin{equation}
\cfrac{1}{7 +
\cfrac{1}{15 +
\cfrac{1}{292 + \dotsb}}}
\end{equation}
$$

      算法的步骤大概如下:
            (1) 根据输入的Number值，取不足近似整数(matlab里面就是fix函数)，比如π就是取3，同时将3存入数组P的第一个值P[1]=fix(Number)(matlab的数组从1开始的)
            (2) 做运算π-3，作为数组A的第一个值A[1]
            (3) P[2]=1/A[1]
            (4) A[2]=fix(P[2])
      下一个P和A重复上述过程。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
function result=test(Number,n)
%Number=(sqrt(5)-1)/2;	%测试用
%n=100;					%测试用
result=0;
P=zeros(1,n);A=zeros(1,n);
P(1)=fix(Number);
A(1)=Number-fix(Number);
for i=2:n-1
    P(i)=1/A(i-1);
    A(i)=P(i)-fix(P(i));
end
P=sym(fix(P));    %转成符号计算
for j=n-1:-1:2
    result=1/(P(j)+result);
end
result=result+P(1);

      测试以后发现，这组算法还是有问题的，前面也说过，最大问题就是数组A是一个double型的数组，把精度截断了。我尝试过用sym型来写，但是会在某个数值的时候因为大数相乘(分子分母比较大的时候，通分造成的原因)直接超过1e308变成inf，然后就NaN了。综合上面，只能放弃掉一些精度。经过测试，直接sym型的时候连分数能写到20层，如果将A进行double化，目前没发现层数的问题，但是有效的层数，100层的话大约在70层左右，vpa以后精度在1e-15左右，还是可以的。计算处理速度上基本没有大的问题，100层的时间大约是2秒钟左右，300层的时间大约是3秒钟左右。

文章目录
|