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

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

割圆术

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

设圆直径

，圆的周长介于内接正

边形与外接正

边形之间：

其中

为正

边形的周长。

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

和

之间。

矩阵加速

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

的倒数进行展开：

其中

为正多边形周长，展开式的第一项

为圆周率精确值，与

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

的

阶，因此即为

。第二项为

的一阶项，相应展开系数为

......以此类推。

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

的

次方阶。

对圆内接正

边形，可得展开式：

同理，对圆内接正

、

、

边形，可得展开式

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

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

、

、

、

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

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

高斯消元法

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

对于

矩阵方程：

将第1行乘消元系数

之后加到第

行，可以消去

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

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

加到第i行，可以消去

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

......

重复操作

次，可以把系数矩阵化为上三角矩阵。

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

在逐次消元中，对第

个元进行消元操作，为第一层for循环。对第

个元的操作过程中，将第

列中

以下的元素

都消为0，为第二层for循环。在初等行变换中，

所在第

行的每一个元素

都需要进行相应的变换，为第三层for循环。列矢量

中每一行元素

也相应进行变换，

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

的变换，位于第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) 函数产生从

开始逐个递增1直到

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

开始逐个递减1直到

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

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

取

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

，使用高斯消元法求解为：

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边形，根据割圆术可得相应多边形周长为：

，计算得：

pi = 3.1415906412698416

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

~

之间（与

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

~

之间（与

前3位数字相符）。

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

~

之间（与

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

~

（与

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

Reference：

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