面向对象有限元编程|自定义求解器之共轭梯度法

共轭梯度法是方程组求解的一种迭代方法。这种方法特别适合有限元求解，因为该方法要求系数矩阵为对称正定矩阵，而有限元平衡方程的系数矩阵正好是对称正定矩阵(考虑边界条件)。同时，共轭梯度法也适合并行计算。

• 算法原理

对于方程组    ，假定    是对称正定矩阵，采用共轭梯度法算法步骤如下：取初始值    

这里    迭代持续进行，直到向量    的模达到一个较小的值，也就是误差允许范围之内。

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就是整体节点力向量。