收敛太慢了
看看马青公式，计算前5项，得到祖冲之的精度。
\
Clear["Global`*"];(*清除所有变量*)
$MaxExtraPrecision=500;(*最大额外精度*)
aa=Sum[(-1)^(k+1)*x^(2k-1)/(2k-1),{k,1,4}];(*计算前n项和*)
pi=(16*aa/.x->1/5)-(4*aa/.x->1/239)(*用马青公式计算圆周率*)
N(*计算差值*)
N(*计算近似值*)
得到
3.141592682404399517240259836073575860489757996720135506301754703591490348774810995008612800369424553
如果计算72项，大概就能得到圆周率的100位
\
差值
-4.7347325070038906601393332841512003488000614441363296637426861454146\
16030131177859665455168151268452*10^-103
估计71项达不到精度，但是我没仔细研究！
nyy 发表于 2023-12-25 12:12
看看马青公式，计算前5项，得到祖冲之的精度。
\
如果要计算到小数点后100位，那么需要2^168=3.741444191*10^(50)正多边形，代码如下
Clear["Global`*"];(*清除所有变量*)
a=Sqrt;(*正方形的初始边长*)
n=4;(*正方形边数*)
Do];(*由n边形的边长,计算出2n边形的边长*)
n=n*2;(*边数从n变成2n*)
pi=a*n/2,(*周长除以直径=圆周率*)
{k,1,166}]
N
aa=N(*鲁道夫的圆周率*)
bb=N(*真实的圆周率*)
结果
3.141592653589793238462643383279502884197169399375105820974944592307816406286208998628034825342117067945
3.141592653589793238462643383279502884197169399375105820974944592307816406286208998628034825342117067982
需要连续开167次方！
nyy 发表于 2023-12-25 12:12
看看马青公式，计算前5项，得到祖冲之的精度。
\
利用马青公式\(48\arctan\frac{1}{18}+32\arctan\frac{1}{57}-20\arctan\frac{1}{239}\)
只需要计算前40项就可以了！
代码如下：
Clear["Global`*"];(*清除所有变量*)
$MaxExtraPrecision=100;(*最大额外精度*)
aa=Sum[(-1)^(k+1)*x^(2k-1)/(2k-1),{k,1,50}];(*计算前n项和*)
pi=(48*aa/.x->1/18)+(32*aa/.x->1/57)-(20*aa/.x->1/239)(*用马青公式计算圆周率*)
N(*计算差值*)
N(*计算近似值*)
N(*计算近似值*)
结果如下
3.1415926535897932384626433832795028841971693993751058209749445923078164062862089986280348253421170679809
3.1415926535897932384626433832795028841971693993751058209749445923078164062862089986280348253421170679821
可以看到，实际上精确到小数点后101位。
nyy 发表于 2023-12-29 11:47
利用马青公式\(48\arctan\frac{1}{18}+32\arctan\frac{1}{57}-20\arctan\frac{1}{239}\)
只需要计算前40 ...
利用马青公式
\(\frac{\pi}{4}=44\arctan\frac{1}{57}+7\arctan\frac{1}{239}-12\arctan\frac{1}{682}+24\arctan\frac{1}{12943}\)
只需要计算29项
Clear["Global`*"];(*清除所有变量*)
$MaxExtraPrecision=250;(*最大额外精度*)
aa=Sum[(-1)^(k+1)*x^(2k-1)/(2k-1),{k,1,29}];(*计算前n项和*)
(*用马青公式计算圆周率*)
pi=4*((44*aa/.x->1/57)+(7*aa/.x->1/239)-(12*aa/.x->1/682)+(24*aa/.x->1/12943));
N(*计算差值*)
N(*计算近似值*)
N(*计算近似值*)
计算结果
3.1415926535897932384626433832795028841971693993751058209749445923078164062862089986280348253421170679822
3.1415926535897932384626433832795028841971693993751058209749445923078164062862089986280348253421170679821
我偶然发现的一个圆周率公式,大家看下：一个简洁的圆周率算法 - happytdw的文章 - 知乎
https://zhuanlan.zhihu.com/p/675031866 tdwhh 发表于 2023-12-31 11:58
我偶然发现的一个圆周率公式,大家看下：一个简洁的圆周率算法 - happytdw的文章 - 知乎
https://zhuanlan.z ...
(*一个简洁的圆周率算法 https://zhuanlan.zhihu.com/p/675031866*)
Clear["Global`*"];(*清除所有变量*)
a0=1;a1=2;(*a数列的初始值*)
b0=1;b1=2;(*b数列的初始值*)
(*利用循环进行递推*)
Do[a2=a1+n*a0;(*计算递推结果*)
b2=b1+n*b0+1;
{a0,a1}={a1,a2};(*重新赋值*)
{b0,b1}={b1,b2},
{n,2,400}]
pi=2*(b2/a2)^2//Simplify;
N
pi的计算结果
62025000939678909516687735193323735679065222430471033763381665839374761995927718792304418729658751210049860488860911551680939456432144359841211072335895952147921464697513968242090392315718012516015590713627465244700797183151760707771499894364547214710591597600846408933165356599304275471466095871602326246357731878903070418372922554993470258688672327713486649550161027062315402388427258430020421815659675995218714419996458444204554859904759312408758535748607611738041254470374876593240980407603750246924235349756903590684663913440602428379474402975747442120724740808012952797639745546189266638840413443786459171826980252856226982197660944851286153199683302766825452008610849213435634382955651685398647960622185880437762598492223473579124892773167361353256715117973875293103124963822411256705801871887131620488818568592134867440992370758611285265299055526977458774244746074556800428542081/19743170989658702829708591285591630717562365062427886961099403631133476119547148422881762950355106763464658566139692506973077016062198092073282093446192327692126784166349645485356214315502798552923766493857569813524669593818230628367479232839474331330314390642177373411824165311786333767561565008246688177567595388974795002389384579994311747141634363559263064669613868071103416913625011554865450732523170394255632143364262004632226712388029603885440999693919866053401215899783057464554443190561413222910839475667718942779276999809824558685043390141945559782178717329802237035412583201025565795456311920839240762285416938939318009154045190203091851759480430913872416213859324276124729267616043236366753593124787498572461727221865110910684771119362709648153306017044421356246752838000108858440613003052598961859241908248055209419156532740031605893034256340075591729806631243285637929893888
误差
4.932603769219693390126654336135238877664753636994259082389243155687944028218791267849224234479444015*10^-17
根据我的计算结果，收敛速度好像并不快，不知道原理是什么 本帖最后由 nyy 于 2024-1-3 13:27 编辑
tdwhh 发表于 2023-12-31 11:58
我偶然发现的一个圆周率公式,大家看下：一个简洁的圆周率算法 - happytdw的文章 - 知乎
https://zhuanlan.z ...
利用这儿的代码
https://bbs.emath.ac.cn/forum.php?mod=redirect&goto=findpost&ptid=17586&pid=98458
(*拉马努金计算圆周率*)
Clear["Global`*"];(*清除所有变量*)
aa=Sum[(4k)!*(1103+26390k)/((k!)^4*396^(4k)),{k,0,2}]//Simplify
pi=9801/(2*Sqrt)/aa
N
圆周率的近似表达式=
\[\frac{2286635172367940241408 \sqrt{2}}{1029347477390786609545}\]
与精确值的误差
5.68242325601395950808108722896*10^-24
这个表达式，比上面的表达式计算要小很多，且精度要搞几个有效数字。
所以我推测你的算法不仅仅计算量大且精度小。
这个是我的论据。
nyy 发表于 2024-1-3 13:26
利用这儿的代码
https://bbs.emath.ac.cn/forum.php?mod=redirect&goto=findpost&ptid=17586&pid=98458
计算前五项的值
(*拉马努金计算圆周率*)
Clear["Global`*"];(*清除所有变量*)
aa=Table)/Sum[(4k)!*(1103+26390k)/((k!)^4*396^(4k)),{k,0,n}],{n,0,5}]
计算结果
\[\begin{array}{c}
\frac{9801}{2206 \sqrt{2}} \\
\frac{2510613731736 \sqrt{2}}{1130173253125} \\
\frac{2286635172367940241408 \sqrt{2}}{1029347477390786609545} \\
\frac{17252765328978109815564789153792 \sqrt{2}}{7766473062254307011793347201855} \\
\frac{1131379202490552979877435552947122965839872 \sqrt{2}}{509299577881529611662930757403081523769055} \\
\frac{128805730098892711723125911845114081418091536842752 \sqrt{2}}{57982950211280781944919792648021104999982386829481} \\
\end{array}\]
页:
1
[2]