连分数近似算法

  (2019年2月19日注:这篇文章原先发在自己github那边的博客,时间是2017年2月5日)

  这道题源自数学实验上面的一组实验,当时困扰了我特别久,题目的内容是用matlab求出π的连分数展开及每层迭代的值。

  因为matlab的数值精度的问题,当你运行3+16450/16421时,就算你设置了format rat,也没有办法显示准确的有理分数的值。原书是用mathmatica做的,直接一句命令就出来了。

  其实matlab里面也是有相应的连分数展开的命令的,详细可以查看以下:

Continued fraction expansion

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

>> 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} = \frac{1}{a_{1}+\frac{1}{a_{2}+\frac{1}{a_{3}+…}}}$$

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

  算法的大致步骤如下:

    (1) 根据输入的Number值,取不足近似整数(matlab里面就是fix函数),比如π就是取3,同时将3存入数组P的第一个值$P[1]=fix(Number)$(matlab的数组从1开始的)

    (2) 做运算$π-3$,作为数组A的第一个值$A[1]$

    (3) $P[2]=\frac{1}{A[1]}$

    (4) $A[2]=fix(P[2])$

  下一个P和A重复上述过程

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秒钟左右。

posted @ 2019-02-19 17:08  JCChan  阅读(2268)  评论(1编辑  收藏  举报