高斯消元法python编程_割圆术计算圆周率与矩阵高斯消元法（Python）

课堂笔记整理：割圆术计算圆周率的矩阵方法，高斯消元法及其Python求解。

割圆术

历史上阿基米德使用了接多边形利用周长计算圆周率的割圆术，刘徽和祖冲之使用了接多边形利用面积计算圆周率的割圆术。这里以周长法为例。

设圆直径

$d=1$ ，圆的周长介于内接正

$k$ 边形与外接正

$k$ 边形之间：

$\pi_{k}=k l_{k}$

$k l_{k}^{\mathrm{in}}<\pi<k l_{k}^{\mathrm{cir}}$

其中

$\pi _k$ 为正

$k$ 边形的周长。

阿基米德通过割正96边形得到圆周率介于

$3.140845$ 和

$3.142857$ 之间。

矩阵加速

根据“正多边形周长逼近圆周率”的思想，可以考虑把正多边形的周长按边数

$k$ 的倒数进行展开：

$\pi_{k}=\pi_{\infty}+\frac{c_{1}}{k}+\frac{c_{2}}{k^{2}}+\frac{c_{3}}{k^{3}}+\cdots$

其中

$\pi _k$ 为正多边形周长，展开式的第一项

$\pi _\infty$ 为圆周率精确值，与

$k$ 无关，可以理解为分母上为

$k$ 的

$0$ 阶，因此即为

$\pi _\infty$ 。第二项为

$\frac{1}{k}$ 的一阶项，相应展开系数为

$c_1$ ......以此类推。

在实际求解中，不可能展开到无穷多阶，因此需要取截断，例如此处只取4项，即展开到

$\frac{1}{k}$ 的

$3$ 次方阶。

对圆内接正

$k_1$ 边形，可得展开式：

$\pi_{\infty}+c_{1} \cdot \frac{1}{k_{1}}+c_{2} \cdot \frac{1}{k_{1}^{2}}+c_{3} \cdot \frac{1}{k_{1}^{3}}=\pi_{k_{1}}$

同理，对圆内接正

$k_2$ 、

$k_3$ 、

$k_4$ 边形，可得展开式

$\pi_{\infty}+c_{1} \cdot \frac{1}{k_{2}}+c_{2} \cdot \frac{1}{k_{2}^{2}}+c_{3} \cdot \frac{1}{k_{2}^{3}}=\pi_{k_{2}}$

$\pi_{\infty}+c_{1} \cdot \frac{1}{k_{3}}+c_{2} \cdot \frac{1}{k_{3}^{2}}+c_{3} \cdot \frac{1}{k_{3}^{3}}=\pi_{k_{3}}$

$\pi_{\infty}+c_{1} \cdot \frac{1}{k_{4}}+c_{2} \cdot \frac{1}{k_{4}^{2}}+c_{3} \cdot \frac{1}{k_{4}^{3}}=\pi_{k_{4}}$

将该方程组写为矩阵形式：

$\left[\begin{array}{cccc} \mathbf{1} & \frac{\mathbf{1}}{\mathbf{k}_{\mathbf{1}}} & \frac{\mathbf{1}}{\mathbf{k}_{\mathbf{1}}^{2}} & \frac{\mathbf{1}}{\mathbf{k}_{\mathbf{1}}^{3}} \\ \mathbf{1} & \frac{\mathbf{1}}{\mathbf{k}_{\mathbf{2}}} & \frac{\mathbf{1}}{\mathbf{k}_{\mathbf{2}}^{2}} & \frac{\mathbf{1}}{\mathbf{k}_{\mathbf{2}}^{3}} \\ \mathbf{1} & \frac{\mathbf{1}}{\mathbf{k}_{3}} & \frac{\mathbf{1}}{\mathbf{k}_{3}^{2}} & \frac{\mathbf{1}}{\mathbf{k}_{3}^{3}} \\ \mathbf{1} & \frac{\mathbf{1}}{\mathbf{k}_{\mathbf{4}}} & \frac{\mathbf{1}}{\mathbf{k}_{\mathbf{4}}^{2}} & \frac{\mathbf{1}}{\mathbf{k}_{4}^{3}} \end{array}\right] \left[\begin{array}{l} \pi_{\infty} \\ c_{1} \\ c_{2} \\ c_{3} \end{array}\right]=\left[\begin{array}{l} \pi_{k_{1}} \\ \pi_{k_{2}} \\ \pi_{k_{3}} \\ \pi_{k_{4}} \end{array}\right]$

4个方程对应4个未知数：

$\pi_\infty$ 、

$c_1$ 、

$c_2$ 、

$c_3$ ，求解即可得到圆周率的值。

下面使用高斯消元法进行求解。

高斯消元法

高斯消元法的思想是通过初等变换进行逐次消元，把方程组转化为等价的上三角方程组，之后从最后一个未知数开始逐次回带进行求解。

对于

$n\times n$ 矩阵方程：

$\left[\begin{array}{ccccc} A_{11} & A_{12} & A_{13} & \dots & A_{1 n} \\ A_{21} & A_{22} & A_{23} & \dots & A_{2 n} \\ A_{31} & A_{32} & A_{33} & \dots & A_{3 n} \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ A_{n 1} & A_{n 2} & A_{n 3} & \dots & A_{n n} \end{array}\right]\left[\begin{array}{c} x_{1} \\ x_{2} \\ x_{3} \\ \vdots \\ x_{n} \end{array}\right]=\left[\begin{array}{c} b_{1} \\ b_{2} \\ b_{3} \\ \vdots \\ b_{n} \end{array}\right]$

将第1行乘消元系数

$-\frac{A_{i1}}{A_{11}}$ 之后加到第

$i$ 行，可以消去

$A_{i1}$ 项。对第2行到第n行进行该操作，可以将第1个元消去，使矩阵第1列中第1行以下元素全为0。

之后开始消第2个元：将第2行乘

$-\frac{A_{i2}}{A_{22}}$ 加到第i行，可以消去

$A_{i2}$ 项。对第3行到第n行进行该操作，可以将第2个元消去，使矩阵第2列中第2行以下元素全为0。

......

重复操作

$(n-1)$ 次，可以把系数矩阵化为上三角矩阵。

为更清楚地展示高斯消元法的原理，此处使用for循环进行描述。一般而言，使用numpy中的矩阵写法比for循环高效很多。

在逐次消元中，对第

$i（i=1,2...n-1）$ 个元进行消元操作，为第一层for循环。对第

$i$ 个元的操作过程中，将第

$i$ 列中

$A_{ii}$ 以下的元素

$A_{ji}$ 都消为0，为第二层for循环。在初等行变换中，

$A_{ji}$ 所在第

$j$ 行的每一个元素

$A_{jk},k=i,i+1...n$ 都需要进行相应的变换，为第三层for循环。列矢量

$\mathbf{b}$ 中每一行元素

$b_{j}$ 也相应进行变换，

$b_j$ 每个元素单独为1行，因此跟随第

$j$ 的变换，位于第2层循环内。

化为上三角矩阵后，可以由第

$n$ 行解得

$x_n$ ，而

$x_n$ 代入第

$(n-1)$ 行，可以解得

$x_{n-1}$ ，将

$x_{n}$ 和

$x_{n-1}$ 代入第

$(n-2)$ 行......以此类推，得到方程的解。这一过程同样可以用for循环写，回带求解是逐次消元之后的步骤，因此该循环位于逐次消元的3重循环之外。

Python计算圆周率

def Gaussian(A, B): # 定义高斯消元函数

dimen = len(A)

for i in range(1, dimen): # 处理第i个元

for j in range(i, dimen):

ratio = A[j][i-1] / A[i-1][i-1] # 消元系数

for k in range(i-1, dimen):

A[j][k] = A[j][k] - A[i-1][k] * ratio # 消元

B[j] = B[j] - B[i-1] * ratio

B[dimen-1] = B[dimen-1] / A[dimen-1][dimen-1]

for i in range(dimen-2, -1, -1):

for j in range(dimen-1, i, -1):

B[i] = B[i] - A[i][j] * B[j]

B[i] = B[i] / A[i][i]

pi = B[i]

print('pi = ' + str(pi))

return B

编写函数时，有几个地方容易踩坑：range(a, b,1) 函数产生从

$a$ 开始逐个递增1直到

$b-1$ 的列表，range(b, a, -1) 产生从

$b$ 开始逐个递减1直到

$a+1$ 的列表，而矩阵的行（列）角标从0开始。

回到计算圆周率的矩阵方程：

取

$k_1=8,k_2=16,k_3=32,k_4=64$ ，即内接正8、16、32、64边形，根据割圆术可得相应多边形周长为：

$3.061467, 3.121445, 3.136548, 3.140331$ ，使用高斯消元法求解为：

A = [[1, 1 / 8 ** 2, 1 / 8 ** 3, 1 / 8 ** 4], [1, 1 / 16 ** 2, 1 / 16 ** 3, 1 / 16 ** 4], [1, 1 / 32 ** 2, 1 / 32 ** 3, 1 / 32 ** 4], [1, 1 / 64 ** 2, 1 / 64 ** 3, 1 / 64 ** 4]]

B = [3.061467, 3.121445, 3.136548, 3.140331]

Gaussian(A, B)

输出：

pi = 3.14159273968254

同理，取外接正8、16、32、64边形，根据割圆术可得相应多边形周长为：

$3.313708, 3.182598, 3.151725, 3.144118$ ，计算得：

pi = 3.1415906412698416

相比之下，阿基米德割圆周长为正96边形，得到圆周率位于

$3.140 845$ ~

$3.142 857$ 之间（与

$\pi$ 前3位数字相符），刘徽割圆面积为正192边形，得到圆周率位于

$3.141 0$ ~

$3.1427$ 之间（与

$\pi$ 前3位数字相符）。

使用矩阵改进的割圆术，在最多割圆为正64边形的情况下，可以求得圆周率位于

$3.14159273968254$ ~

$3.1415906412698416$ 之间（与

$\pi$ 前6位数字相符），而祖冲之得到类似精度的圆周率

$3.141 596$ ~

$3.141 597$ （与

$\pi$ 前6位数字相符），则割圆为正24576边形。

Reference：

2、周善贵，《计算物理》课程讲义

高斯消元法python编程_割圆术计算圆周率与矩阵高斯消元法（Python）相关推荐

割圆术c语言程序设计,c语言实现割圆术计算圆周率.pdf
割圆术计算圆周率 "割圆术"是我国数学家刘徽创立的一种求圆周率的方法.思想是当圆的内接正多边形的边数无限大时内接正多边形的面积就无限趋近于圆的面积,即所谓 "割之弥细 ...
高等数学Mathematica实验题——2.1 - 15 用割圆术计算圆周率（Calcaluation of π with cyclotomic method ）
题目: 2.1-15. 用割圆术计算圆周率(题干描述详见教材) an=Sqrt[2-2Sqrt[1-(an-1/2)^2]] Sn=3*2^(n-1)*Sn
使用割圆术计算圆周率
使用割圆术计算圆周率理论基础圆周长计算公式: C = pi * d (C是圆的周长,pi是圆周率,d是直径) 推导可得 pi = C/d 所以我们只需要测量一定直径对应的周长即可算出圆周率. 思考 ...
大学生计算机python_人人都能学计算机：计算机科学入门与Python编程_学堂在线章节测试答案...
查看答案人人都能学计算机:计算机科学入门与Python编程_学堂在线章节测试答案单击图层调板下方的新图层按钮可以产生新图层.A:错B:对在图示的薄壁杆件截面图形中,形心与弯曲中心重合的截面有() ...
python 科学计算基础教程电子版-自学Python 编程基础、科学计算及数据分析
自学Python 编程基础.科学计算及数据分析 epub pdf mobi txt 下载自学Python 编程基础.科学计算及数据分析 epub pdf mobi txt 下载 ☆☆☆☆☆ 李金著 ...
python编程入门与案例详解-自学Python 编程基础、科学计算及数据分析
自学Python 编程基础.科学计算及数据分析 epub pdf mobi txt 下载自学Python 编程基础.科学计算及数据分析 epub pdf mobi txt 下载 ☆☆☆☆☆ 李金著 ...
《读九章学python》如何用Python编程实现少广术？
卷四少广以御积幂方圆随着人工智能概念的大火,其重要的支持语言Python也一路高歌猛,Python的设计哲学是"优雅.明确.简单". <九章算术>我国现存的最古老 ...
Python编程新手看过来，如何求素数（Python学习教程）
本期的Python学习教程是针对新入门Python编程的新手来写的:关于怎么求素数! 一.什么是素数? 素数就是质数,通俗点说就是只能被1和其本身整数的数就是素数(1除外) 举个例子: 2,3,4,5 ...
python 项目学编程_《从问题到程序：用Python学编程和计算》——3.5 练习-阿里云开发者社区...
复习下面概念:数值积分,区间分割法,舍入误差,简单重复,累积,累积变量,生成和筛选,递推,递推变量,素数(质数),因子和真因子,哥德巴赫猜想,输入循环,输入控制的循环,递归定义,递归函数,循环定义,无 ...
高斯消元法python编程_高斯消元法的Python实现
高斯消元法节约内存的算法: 例1:用高斯消元法求解线性方程组: # -*- coding: utf-8 -*- """ 求解线性方程组: [12 -3 3] [x1] ...

高斯消元法python编程_割圆术计算圆周率与矩阵高斯消元法（Python）

高斯消元法python编程_割圆术计算圆周率与矩阵高斯消元法（Python）相关推荐

最新文章

热门文章