连分数近似算法
(2019年2月19日注:这篇文章原先发在自己github那边的博客,时间是2017年2月5日)
这道题源自数学实验上面的一组实验,当时困扰了我特别久,题目的内容是用matlab求出π的连分数展开及每层迭代的值。
因为matlab的数值精度的问题,当你运行3+16450/16421时,就算你设置了format rat,也没有办法显示准确的有理分数的值。原书是用mathmatica做的,直接一句命令就出来了。
其实matlab里面也是有相应的连分数展开的命令的,详细可以查看以下:
但是显然精度高时,就会遇到那天叙森写的学习笔记的另外一面,精度要求太高,函数无法满足要求,举个例子。
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了。
好了,说了这么多,该介绍下连分数是什么了,其实也不难,据说拉马努金对这个非常的了解。Sn=1a1+1a2+1a3+…
为了算法的简便,我们设定的是这组数列全部都是正整数。对于有理数,我们知道,这组数列一定是有限数列,对于无理数,这组数列大部分情况下是无限数列。我们可以从某一个数开始,设定它为0,按照这样子的顺序,加出来一个有理分数,这个有理分数就是我们需要的值,举个例子。π=13+17+115+…
算法的大致步骤如下:
(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重复上述过程
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秒钟左右。
【推荐】编程新体验,更懂你的AI,立即体验豆包MarsCode编程助手
【推荐】凌霞软件回馈社区,博客园 & 1Panel & Halo 联合会员上线
【推荐】抖音旗下AI助手豆包,你的智能百科全书,全免费不限次数
【推荐】博客园社区专享云产品让利特惠,阿里云新客6.5折上折
【推荐】轻量又高性能的 SSH 工具 IShell:AI 加持,快人一步