数值方式求解圆周率π

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

如图，假设圆的半径为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