连分数近似算法
(2019年2月19日注:这篇文章原先发在自己github那边的博客,时间是2017年2月5日)
这道题源自数学实验上面的一组实验,当时困扰了我特别久,题目的内容是用matlab求出π的连分数展开及每层迭代的值。
因为matlab的数值精度的问题,当你运行3+16450/16421时,就算你设置了format rat,也没有办法显示准确的有理分数的值。原书是用mathmatica做的,直接一句命令就出来了。
其实matlab里面也是有相应的连分数展开的命令的,详细可以查看以下:
但是显然精度高时,就会遇到那天叙森写的学习笔记的另外一面,精度要求太高,函数无法满足要求,举个例子。
>> 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秒钟左右。