课堂笔记整理：割圆术计算圆周率的矩阵方法，高斯消元法及其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、周善贵，《计算物理》课程讲义
