面向对象有限元编程|自定义求解器之共轭梯度法
共轭梯度法是方程组求解的一种迭代方法。这种方法特别适合有限元求解，因为该方法要求系数矩阵为对称正定矩阵，而有限元平衡方程的系数矩阵正好是对称正定矩阵(考虑边界条件)。同时，共轭梯度法也适合并行计算。算法原理对于方程组 ，假定 是对称正定矩阵，采用共轭梯度法算法步骤如下：取初始值 这里 迭代持续进行，直到向量 的模达到一个较小的值，也就是误差允许范围之内。python代码如下：import mathimport numpy as npclass 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 npimport ModsolverA = 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) #创建一个求解器的实例clsx = cls.CGsolver() #调用共轭梯度法求解print(x) 以上是测试文件。实际在有限元分析使用时，A就是刚度矩阵，b就是整体节点力向量。来源：数值分析与有限元编程
