数值方式求解圆周率π

本文介绍如何通过数值方式求解圆周率π。

如图,假设圆的半径为1,可知圆的周长为2π,我们现在只需要用积分的方法求出 1/4 周长,即为π/2。

现在讲解如何用积分的方式求出 1/4 周长。

根据勾股定理,可知

Δy=\sqrt{1-x^2}-\sqrt{1-(x+Δx)^2}

Δs=\sqrt{(Δx)^2+(Δy)^2}\\\ \ \ \ \ =\sqrt{(Δx)^2+(1-x^2)+[1-(x+Δx)^2]-2\sqrt{1-x^2}\sqrt{1-(x+Δx)^2}}\\\ \ \ \ \ =\sqrt{2-2x^2-2xΔx-2\sqrt{1-x^2}\sqrt{1-(x+Δx)^2}}\\\ \ \ \ \ =\sqrt{2}\sqrt{1-x^2-xΔx-\sqrt{1-x^2}\sqrt{1-(x+Δx)^2}}

为了处理式子中根号下的Δx,我们对式子做泰勒展开

f(Δx)=\sqrt{1-(x+Δx)^2)}

分别求出其一阶、二阶导数

f^{(1)}(Δx)=\frac{1}{2}\frac{-2(x+Δx)}{\sqrt{1-(x+Δx)^2}}=-\frac{x+Δx}{\sqrt{1-(x+Δx)^2}}

f^{2}(Δx)=-\frac{(x+Δx)^2}{\sqrt{[1-(x+Δx)^2]^3}}-\frac{1}{\sqrt{1-(x+Δx)^2}}

忽略高阶无穷小量后

f(Δx)≈f(0)+f^{(1)}(0)Δx+\frac{1}{2}f^{(2)}(0)(Δx)^2\\\ \ \ \ \ \ \ \ \ \ \ ≈\sqrt{1-x^2}-\frac{xΔx}{\sqrt{1-x^2}}-\frac{x^2(Δx)^2}{2\sqrt{[1-x^2]^3}}-\frac{(Δx)^2}{2\sqrt{1-x^2}}

把此式代回Δs

Δs=\sqrt{2}\sqrt{1-x^2-xΔx-\sqrt{1-x^2}(\sqrt{1-x^2}-\frac{xΔx}{\sqrt{1-x^2}}-\frac{x^2(Δx)^2}{2\sqrt{[1-x^2]^3}}-\frac{(Δx)^2}{2\sqrt{1-x^2}})}\\\ \ \ \ \ =\sqrt{2}\sqrt{1-x^2-xΔx-(1-x^2)+xΔx+\frac{x^2(Δx)^2}{2(1-x^2)}+\frac{(Δx)^2}{2}}\\\ \ \ \ \ =\sqrt{\frac{x^2}{1-x^2}+1}Δx\\\ \ \ \ \ =\sqrt{\frac{1}{1-x^2}}Δx

故根据积分公式, 1/4 周长为

\frac{π}{2}=\int_0^1{ds}\\\ \ \ =\int_0^1{\sqrt{\frac{1}{1-x^2}}dx}

采用最简单的定积分数值计算解法,直接取1e-11为迭代步长进行求和,解出定积分

#include <math.h>
#include <stdio.h>

int main()
{
    double delta_x = 1e-11;
    double s = 0;
    unsigned long long i = 0;
    for (double x = 0; x <= 1; x+=delta_x) {
        double delta_s = sqrt(1.0/(1.0-x*x)) * delta_x;
        s += delta_s;
        
        if (0 == i % 10000000000) {
            printf("current x:%lf\n", x);
        }
        i ++;
    }
    
    return 0;
}

运行代码,经过超级漫长的等待,可以看到2*s已经有圆周率3.14159的影子了。

7149DE4C9B432C2D7F32E7DAB1CA8966.png

原创声明,本文系作者授权云+社区发表,未经许可,不得转载。

如有侵权,请联系 yunjia_community@tencent.com 删除。

编辑于

我来说两句

0 条评论
登录 后参与评论

相关文章

来自专栏Aloys的开发之路

C语言中随机数相关问题

用C语言产生随机数重要用到rand函数、srand函数、及宏RAND_MAX(32767),它们均在stdlib.h中进行了声明。 int rand(void)...

87880
来自专栏封碎

当今世界最为经典的十大算法 博客分类: 经典文章转载 算法数据结构网络应用数据挖掘J#

本文转载自July CSDN博客:http://blog.csdn.net/v_JULY_v/archive/2011/03/07/6228235.aspx

26120
来自专栏程序员宝库

码农也要学算法

当“人工智能”、“AlphaGo”、“无人驾驶”、“智能投顾”等词语不断在人们视野中出现的时候,意味着我们正步入一个算法的时代。计算机通过提供给人类每天要面临的...

516100
来自专栏生信技能树

比较不同的对单细胞转录组数据聚类的方法

背景介绍 聚类之前必须要对表达矩阵进行normalization,而且要去除一些批次效应等外部因素。通过对表达矩阵的聚类,可以把细胞群体分成不同的状态,解释为什...

973120
来自专栏小樱的经验随笔

【机器学习笔记之一】深入浅出学习K-Means算法

摘要:在数据挖掘中,K-Means算法是一种 cluster analysis 的算法,其主要是来计算数据聚集的算法,主要通过不断地取离种子点最近均值的算法。 ...

30990
来自专栏机器之心

教程 | 如何从TensorFlow转入PyTorch

选自Medium 作者:Illarion Khlestov 机器之心编译 参与:李泽南、黄小天 当我第一次尝试学习 PyTorch 时,没几天就放弃了。和 Te...

797150
来自专栏机器学习算法与Python学习

Encoder-Decoder自动生成对联,要试试么?

另外,点击阅读原文尝试微软的自动对联系统(http://duilian.msra.cn/app/couplet.aspx)

14400
来自专栏编程心路

算法妙应用-算法的复杂度

对于任何一个程序来说,都可以从三个方面进行分析,分别是 输入、处理、输出,也即 IPO(Input、Process、Output),这种分析方法对硬件和软件程序...

12930
来自专栏机器之心

EMNLP 2018 | 结合通用和专用NMT的优势,CMU为NMT引入「语境参数生成器」

神经机器翻译(NMT)无需单独训练或调整系统的任何部分就可以直接建模源语言到目标语言的映射。这使得 NMT 快速发展,并在许多大规模环境中成功应用 (Wu et...

14210
来自专栏PPV课数据科学社区

强大的PyTorch:10分钟让你了解深度学习领域新流行的框架

摘要: 今年一月份开源的PyTorch,因为它强大的功能,它现在已经成为深度学习领域新流行框架,它的强大源于它内部有很多内置的库。本文就着重介绍了其中几种有特色...

38590

扫码关注云+社区

领取腾讯云代金券