共轭梯度法是方程组求解的一种迭代方法。这种方法特别适合有限元求解,因为该方法要求系数矩阵为对称正定矩阵,而有限元平衡方程的系数矩阵正好是对称正定矩阵(考虑边界条件)。同时,共轭梯度法也适合并行计算。
对于方程组 ,假定 是对称正定矩阵,采用共轭梯度法算法步骤如下:取初始值
这里 迭代持续进行,直到向量 的模达到一个较小的值,也就是误差允许范围之内。
python代码如下:
import math
import numpy as np
class solver:
def __init__(self, A, b, eps, max_iter):
self.A = A
self.b = b
self.eps = eps
self.max_iter = max_iter
def CGsolver (self):
iter = 0
x0 = 0.
g0 = -self.b
r0 = self.b
err = 1e6
while ( err > self.eps and iter < self.max_iter ):
tmp1 = self.vec_dot(g0, g0)
err = math.fabs(tmp1)
tmp_row = np.dot(self.A, r0)
tmp2 = self.vec_dot(tmp_row, r0)
alpha = tmp1 / tmp2
x1 = x0 + alpha * r0
g1 = g0 + alpha * tmp_row
tmp3 = self.vec_dot(g1, g1)
beta = tmp3 / tmp1
r1 = -g1 + beta * r0
x0 = x1
r0 = r1
g0 = g1
iter += 1
return x1
#向量点乘
@staticmethod
def vec_dot(v1, v2):
assert len(v1) == len(v2), 'Length of vectors must be same'
return sum( a * b for a, b in zip(v1, v2) )
以上是求解器类,保存在Modsolver.py文件中。
import numpy as np
import Modsolver
A = np.array( [ [5, 7, 6, 5],
[7, 10, 8, 7],
[6, 8, 10, 9],
[5, 7, 9, 10] ] )
b = np.array( [62, 87, 91, 90] )
cls = Modsolver.solver(A, b, 1e-4, 500) #创建一个求解器的实例cls
x = cls.CGsolver() #调用共轭梯度法求解
print(x)
以上是测试文件。
实际在有限元分析使用时,A就是刚度矩阵,b就是整体节点力向量。