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

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

割圆术

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

设圆直径

equation?tex=d%3D1 ,圆的周长介于内接正

equation?tex=k 边形与外接正

equation?tex=k 边形之间:

equation?tex=%5Cpi_%7Bk%7D%3Dk+l_%7Bk%7D+

equation?tex=k+l_%7Bk%7D%5E%7B%5Cmathrm%7Bin%7D%7D%3C%5Cpi%3Ck+l_%7Bk%7D%5E%7B%5Cmathrm%7Bcir%7D%7D

其中

equation?tex=%5Cpi+_k 为正

equation?tex=k 边形的周长。

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

equation?tex=3.140845

equation?tex=3.142857 之间。

矩阵加速

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

equation?tex=k 的倒数进行展开:

equation?tex=%5Cpi_%7Bk%7D%3D%5Cpi_%7B%5Cinfty%7D%2B%5Cfrac%7Bc_%7B1%7D%7D%7Bk%7D%2B%5Cfrac%7Bc_%7B2%7D%7D%7Bk%5E%7B2%7D%7D%2B%5Cfrac%7Bc_%7B3%7D%7D%7Bk%5E%7B3%7D%7D%2B%5Ccdots

其中

equation?tex=%5Cpi+_k 为正多边形周长,展开式的第一项

equation?tex=%5Cpi+_%5Cinfty 为圆周率精确值,与

equation?tex=k 无关,可以理解为分母上为

equation?tex=k

equation?tex=0 阶,因此即为

equation?tex=%5Cpi+_%5Cinfty。第二项为

equation?tex=%5Cfrac%7B1%7D%7Bk%7D 的一阶项,相应展开系数为

equation?tex=c_1......以此类推。

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

equation?tex=%5Cfrac%7B1%7D%7Bk%7D

equation?tex=3 次方阶。

对圆内接正

equation?tex=k_1 边形,可得展开式:

equation?tex=%5Cpi_%7B%5Cinfty%7D%2Bc_%7B1%7D+%5Ccdot+%5Cfrac%7B1%7D%7Bk_%7B1%7D%7D%2Bc_%7B2%7D+%5Ccdot+%5Cfrac%7B1%7D%7Bk_%7B1%7D%5E%7B2%7D%7D%2Bc_%7B3%7D+%5Ccdot+%5Cfrac%7B1%7D%7Bk_%7B1%7D%5E%7B3%7D%7D%3D%5Cpi_%7Bk_%7B1%7D%7D

同理,对圆内接正

equation?tex=k_2

equation?tex=k_3

equation?tex=k_4 边形,可得展开式

equation?tex=%5Cpi_%7B%5Cinfty%7D%2Bc_%7B1%7D+%5Ccdot+%5Cfrac%7B1%7D%7Bk_%7B2%7D%7D%2Bc_%7B2%7D+%5Ccdot+%5Cfrac%7B1%7D%7Bk_%7B2%7D%5E%7B2%7D%7D%2Bc_%7B3%7D+%5Ccdot+%5Cfrac%7B1%7D%7Bk_%7B2%7D%5E%7B3%7D%7D%3D%5Cpi_%7Bk_%7B2%7D%7D

equation?tex=%5Cpi_%7B%5Cinfty%7D%2Bc_%7B1%7D+%5Ccdot+%5Cfrac%7B1%7D%7Bk_%7B3%7D%7D%2Bc_%7B2%7D+%5Ccdot+%5Cfrac%7B1%7D%7Bk_%7B3%7D%5E%7B2%7D%7D%2Bc_%7B3%7D+%5Ccdot+%5Cfrac%7B1%7D%7Bk_%7B3%7D%5E%7B3%7D%7D%3D%5Cpi_%7Bk_%7B3%7D%7D

equation?tex=%5Cpi_%7B%5Cinfty%7D%2Bc_%7B1%7D+%5Ccdot+%5Cfrac%7B1%7D%7Bk_%7B4%7D%7D%2Bc_%7B2%7D+%5Ccdot+%5Cfrac%7B1%7D%7Bk_%7B4%7D%5E%7B2%7D%7D%2Bc_%7B3%7D+%5Ccdot+%5Cfrac%7B1%7D%7Bk_%7B4%7D%5E%7B3%7D%7D%3D%5Cpi_%7Bk_%7B4%7D%7D

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

equation?tex=%5Cleft%5B%5Cbegin%7Barray%7D%7Bcccc%7D+%5Cmathbf%7B1%7D+%26+%5Cfrac%7B%5Cmathbf%7B1%7D%7D%7B%5Cmathbf%7Bk%7D_%7B%5Cmathbf%7B1%7D%7D%7D+%26+%5Cfrac%7B%5Cmathbf%7B1%7D%7D%7B%5Cmathbf%7Bk%7D_%7B%5Cmathbf%7B1%7D%7D%5E%7B2%7D%7D+%26+%5Cfrac%7B%5Cmathbf%7B1%7D%7D%7B%5Cmathbf%7Bk%7D_%7B%5Cmathbf%7B1%7D%7D%5E%7B3%7D%7D+%5C%5C+%5Cmathbf%7B1%7D+%26+%5Cfrac%7B%5Cmathbf%7B1%7D%7D%7B%5Cmathbf%7Bk%7D_%7B%5Cmathbf%7B2%7D%7D%7D+%26+%5Cfrac%7B%5Cmathbf%7B1%7D%7D%7B%5Cmathbf%7Bk%7D_%7B%5Cmathbf%7B2%7D%7D%5E%7B2%7D%7D+%26+%5Cfrac%7B%5Cmathbf%7B1%7D%7D%7B%5Cmathbf%7Bk%7D_%7B%5Cmathbf%7B2%7D%7D%5E%7B3%7D%7D+%5C%5C+%5Cmathbf%7B1%7D+%26+%5Cfrac%7B%5Cmathbf%7B1%7D%7D%7B%5Cmathbf%7Bk%7D_%7B3%7D%7D+%26+%5Cfrac%7B%5Cmathbf%7B1%7D%7D%7B%5Cmathbf%7Bk%7D_%7B3%7D%5E%7B2%7D%7D+%26+%5Cfrac%7B%5Cmathbf%7B1%7D%7D%7B%5Cmathbf%7Bk%7D_%7B3%7D%5E%7B3%7D%7D+%5C%5C+%5Cmathbf%7B1%7D+%26+%5Cfrac%7B%5Cmathbf%7B1%7D%7D%7B%5Cmathbf%7Bk%7D_%7B%5Cmathbf%7B4%7D%7D%7D+%26+%5Cfrac%7B%5Cmathbf%7B1%7D%7D%7B%5Cmathbf%7Bk%7D_%7B%5Cmathbf%7B4%7D%7D%5E%7B2%7D%7D+%26+%5Cfrac%7B%5Cmathbf%7B1%7D%7D%7B%5Cmathbf%7Bk%7D_%7B4%7D%5E%7B3%7D%7D+%5Cend%7Barray%7D%5Cright%5D++%5Cleft%5B%5Cbegin%7Barray%7D%7Bl%7D+%5Cpi_%7B%5Cinfty%7D+%5C%5C+c_%7B1%7D+%5C%5C+c_%7B2%7D+%5C%5C+c_%7B3%7D+%5Cend%7Barray%7D%5Cright%5D%3D%5Cleft%5B%5Cbegin%7Barray%7D%7Bl%7D+%5Cpi_%7Bk_%7B1%7D%7D+%5C%5C+%5Cpi_%7Bk_%7B2%7D%7D+%5C%5C+%5Cpi_%7Bk_%7B3%7D%7D+%5C%5C+%5Cpi_%7Bk_%7B4%7D%7D+%5Cend%7Barray%7D%5Cright%5D

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

equation?tex=%5Cpi_%5Cinfty

equation?tex=c_1

equation?tex=c_2

equation?tex=c_3 ,求解即可得到圆周率的值。

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

高斯消元法

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

对于

equation?tex=n%5Ctimes+n 矩阵方程:

equation?tex=%5Cleft%5B%5Cbegin%7Barray%7D%7Bccccc%7D+A_%7B11%7D+%26+A_%7B12%7D+%26+A_%7B13%7D+%26+%5Cdots+%26+A_%7B1+n%7D+%5C%5C+A_%7B21%7D+%26+A_%7B22%7D+%26+A_%7B23%7D+%26+%5Cdots+%26+A_%7B2+n%7D+%5C%5C+A_%7B31%7D+%26+A_%7B32%7D+%26+A_%7B33%7D+%26+%5Cdots+%26+A_%7B3+n%7D+%5C%5C+%5Cvdots+%26+%5Cvdots+%26+%5Cvdots+%26+%5Cddots+%26+%5Cvdots+%5C%5C+A_%7Bn+1%7D+%26+A_%7Bn+2%7D+%26+A_%7Bn+3%7D+%26+%5Cdots+%26+A_%7Bn+n%7D+%5Cend%7Barray%7D%5Cright%5D%5Cleft%5B%5Cbegin%7Barray%7D%7Bc%7D+x_%7B1%7D+%5C%5C+x_%7B2%7D+%5C%5C+x_%7B3%7D+%5C%5C+%5Cvdots+%5C%5C+x_%7Bn%7D+%5Cend%7Barray%7D%5Cright%5D%3D%5Cleft%5B%5Cbegin%7Barray%7D%7Bc%7D+b_%7B1%7D+%5C%5C+b_%7B2%7D+%5C%5C+b_%7B3%7D+%5C%5C+%5Cvdots+%5C%5C+b_%7Bn%7D+%5Cend%7Barray%7D%5Cright%5D

将第1行乘消元系数

equation?tex=-%5Cfrac%7BA_%7Bi1%7D%7D%7BA_%7B11%7D%7D 之后加到第

equation?tex=i 行,可以消去

equation?tex=A_%7Bi1%7D 项。对第2行到第n行进行该操作,可以将第1个元消去,使矩阵第1列中第1行以下元素全为0。

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

equation?tex=-%5Cfrac%7BA_%7Bi2%7D%7D%7BA_%7B22%7D%7D 加到第i行,可以消去

equation?tex=A_%7Bi2%7D 项。对第3行到第n行进行该操作,可以将第2个元消去,使矩阵第2列中第2行以下元素全为0。

......

重复操作

equation?tex=%28n-1%29 次,可以把系数矩阵化为上三角矩阵。

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

在逐次消元中,对第

equation?tex=i%EF%BC%88i%3D1%2C2...n-1%EF%BC%89 个元进行消元操作,为第一层for循环。对第

equation?tex=i 个元的操作过程中,将第

equation?tex=i 列中

equation?tex=A_%7Bii%7D 以下的元素

equation?tex=A_%7Bji%7D 都消为0,为第二层for循环。在初等行变换中,

equation?tex=A_%7Bji%7D 所在第

equation?tex=j 行的每一个元素

equation?tex=A_%7Bjk%7D%2Ck%3Di%2Ci%2B1...n 都需要进行相应的变换,为第三层for循环。列矢量

equation?tex=%5Cmathbf%7Bb%7D 中每一行元素

equation?tex=b_%7Bj%7D 也相应进行变换,

equation?tex=b_j 每个元素单独为1行,因此跟随第

equation?tex=j 的变换,位于第2层循环内。

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

equation?tex=n 行解得

equation?tex=x_n ,而

equation?tex=x_n 代入第

equation?tex=%28n-1%29 行,可以解得

equation?tex=x_%7Bn-1%7D ,将

equation?tex=x_%7Bn%7D

equation?tex=x_%7Bn-1%7D 代入第

equation?tex=%28n-2%29 行......以此类推,得到方程的解。这一过程同样可以用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) 函数产生从

equation?tex=a 开始逐个递增1直到

equation?tex=b-1 的列表,range(b, a, -1) 产生从

equation?tex=b 开始逐个递减1直到

equation?tex=a%2B1 的列表,而矩阵的行(列)角标从0开始。

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

equation?tex=%5Cleft%5B%5Cbegin%7Barray%7D%7Bcccc%7D+%5Cmathbf%7B1%7D+%26+%5Cfrac%7B%5Cmathbf%7B1%7D%7D%7B%5Cmathbf%7Bk%7D_%7B%5Cmathbf%7B1%7D%7D%7D+%26+%5Cfrac%7B%5Cmathbf%7B1%7D%7D%7B%5Cmathbf%7Bk%7D_%7B%5Cmathbf%7B1%7D%7D%5E%7B2%7D%7D+%26+%5Cfrac%7B%5Cmathbf%7B1%7D%7D%7B%5Cmathbf%7Bk%7D_%7B%5Cmathbf%7B1%7D%7D%5E%7B3%7D%7D+%5C%5C+%5Cmathbf%7B1%7D+%26+%5Cfrac%7B%5Cmathbf%7B1%7D%7D%7B%5Cmathbf%7Bk%7D_%7B%5Cmathbf%7B2%7D%7D%7D+%26+%5Cfrac%7B%5Cmathbf%7B1%7D%7D%7B%5Cmathbf%7Bk%7D_%7B%5Cmathbf%7B2%7D%7D%5E%7B2%7D%7D+%26+%5Cfrac%7B%5Cmathbf%7B1%7D%7D%7B%5Cmathbf%7Bk%7D_%7B%5Cmathbf%7B2%7D%7D%5E%7B3%7D%7D+%5C%5C+%5Cmathbf%7B1%7D+%26+%5Cfrac%7B%5Cmathbf%7B1%7D%7D%7B%5Cmathbf%7Bk%7D_%7B3%7D%7D+%26+%5Cfrac%7B%5Cmathbf%7B1%7D%7D%7B%5Cmathbf%7Bk%7D_%7B3%7D%5E%7B2%7D%7D+%26+%5Cfrac%7B%5Cmathbf%7B1%7D%7D%7B%5Cmathbf%7Bk%7D_%7B3%7D%5E%7B3%7D%7D+%5C%5C+%5Cmathbf%7B1%7D+%26+%5Cfrac%7B%5Cmathbf%7B1%7D%7D%7B%5Cmathbf%7Bk%7D_%7B%5Cmathbf%7B4%7D%7D%7D+%26+%5Cfrac%7B%5Cmathbf%7B1%7D%7D%7B%5Cmathbf%7Bk%7D_%7B%5Cmathbf%7B4%7D%7D%5E%7B2%7D%7D+%26+%5Cfrac%7B%5Cmathbf%7B1%7D%7D%7B%5Cmathbf%7Bk%7D_%7B4%7D%5E%7B3%7D%7D+%5Cend%7Barray%7D%5Cright%5D++%5Cleft%5B%5Cbegin%7Barray%7D%7Bl%7D+%5Cpi_%7B%5Cinfty%7D+%5C%5C+c_%7B1%7D+%5C%5C+c_%7B2%7D+%5C%5C+c_%7B3%7D+%5Cend%7Barray%7D%5Cright%5D%3D%5Cleft%5B%5Cbegin%7Barray%7D%7Bl%7D+%5Cpi_%7Bk_%7B1%7D%7D+%5C%5C+%5Cpi_%7Bk_%7B2%7D%7D+%5C%5C+%5Cpi_%7Bk_%7B3%7D%7D+%5C%5C+%5Cpi_%7Bk_%7B4%7D%7D+%5Cend%7Barray%7D%5Cright%5D

equation?tex=k_1%3D8%2Ck_2%3D16%2Ck_3%3D32%2Ck_4%3D64 ,即内接正8、16、32、64边形,根据割圆术可得相应多边形周长为:

equation?tex=3.061467%2C+3.121445%2C+3.136548%2C+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边形,根据割圆术可得相应多边形周长为:

equation?tex=3.313708%2C+3.182598%2C+3.151725%2C+3.144118 ,计算得:

pi = 3.1415906412698416

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

equation?tex=3.140+845 ~

equation?tex=3.142+857 之间(与

equation?tex=%5Cpi 前3位数字相符),刘徽割圆面积为正192边形,得到圆周率位于

equation?tex=3.141+0 ~

equation?tex=3.1427 之间(与

equation?tex=%5Cpi 前3位数字相符)。

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

equation?tex=3.14159273968254 ~

equation?tex=3.1415906412698416 之间(与

equation?tex=%5Cpi 前6位数字相符),而祖冲之得到类似精度的圆周率

equation?tex=3.141+596 ~

equation?tex=3.141+597 (与

equation?tex=%5Cpi 前6位数字相符),则割圆为正24576边形。

Reference:

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

评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值