本文介绍如何通过数值方式求解圆周率π.
如图, 假设圆的半径为 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
来源: https://www.qcloud.com/developer/article/1345702