共軛梯度法是方程組求解的一種迭代方法。這種方法特別適合有限元求解，因?yàn)樵摲椒ㄒ笙禂?shù)矩陣為對(duì)稱正定矩陣，而有限元平衡方程的系數(shù)矩陣正好是對(duì)稱正定矩陣(考慮邊界條件)。同時(shí)，共軛梯度法也適合并行計(jì)算。

• 算法原理

對(duì)于方程組 ，假定 是對(duì)稱正定矩陣，采用共軛梯度法算法步驟如下：取初始值

這里 迭代持續(xù)進(jìn)行，直到向量 的模達(dá)到一個(gè)較小的值，也就是誤差允許范圍之內(nèi)。

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



    #向量點(diǎn)乘
    @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) )

以上是求解器類(lèi)，保存在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) #創(chuàng)建一個(gè)求解器的實(shí)例cls

x = cls.CGsolver() #調(diào)用共軛梯度法求解
print(x) 

以上是測(cè)試文件。

實(shí)際在有限元分析使用時(shí)，A就是剛度矩陣，b就是整體節(jié)點(diǎn)力向量。