最近遇到的一个求解雅可比迭代的问题，这个计算方法在 python 中有现成的库，但是在 golang 中没找到相应的实现。
于是根据雅可比行列式的推导实现了一个 golang 版本的雅可比迭代。
?
雅可比迭代
推导
一个 \(N \times N\) 的线性方程组 。
\[Ax = b
\]
其中：
\[A = \begin{bmatrix}
a_{11} & a_{12} & \cdots & a_{1n} \a_{21} & a_{22} & \cdots & a_{2n} \\vdots & \vdots & \ddots & \vdots \a_{n1} & a_{n2} & \cdots & a_{nn}
\end{bmatrix},
x = \begin{bmatrix}
x_{1} \x_{2} \\vdots \x_{n}
\end{bmatrix},
b = \begin{bmatrix}
b_{1} \b_{2} \\vdots \b_{n}
\end{bmatrix}
\]
将 \(A\) 分解成对角矩阵 \(D\) 和 其余部分 \(R\)
?
\[A = D + R，其中 D = \begin{bmatrix}
a_{11} & 0 & \cdots & 0 \0 & a_{22} & \cdots & 0 \\vdots & \vdots & \ddots & \vdots \0 & 0 & \cdots & a_{nn}
\end{bmatrix},
R = \begin{bmatrix}
0 & a_{12} & \cdots & a_{1n} \a_{21} & 0 & \cdots & a_{2n} \\vdots & \vdots & \ddots & \vdots \a_{n1} & a_{n2} & \cdots & 0
\end{bmatrix}
\]
线性方程组可以改写为：
\[Dx = b - Rx
\]
雅可比迭代就是在每一次的迭代中，上一次算出的 x 被用在右侧，用来计算左侧新的x。
这个过程用公式表示如下：
\[x^{(k+1)} = D^{-1}(b - Rx^{(k)})
\]
每个元素就是：
\[x^{(k+1)}_{i} = \frac{1}{a_{ii}}(b_{i} - \sum_{j \neq i}a_{ij}x^{(k)}_{j}), i = 1,2,\cdots,n.
\]
也就是：
\[x^{(k+1)}_{i} = \frac{b_{i}}{a_{ii}} - \sum_{j \neq i}\frac{a_{ij}}{a_{ii}}x^{(k)}_{j}, i = 1,2,\cdots,n.
\]
实现函数
主要根据最后一步推导出的公式，实现如下：
package mainimport ( "fmt" "math")// Jacobi 迭代法 输入系数矩阵 mx、值矩阵 mr、迭代次数 n、误差 c(以 list 模拟矩阵 行优先)func Jacobi(mx [][]float64, mr []float64, n int, c float64) ([]float64, error) { if len(mx) != len(mr) { return nil, fmt.Errorf("系数矩阵和值矩阵长度不匹配") } x := initX(len(mr)) // 迭代初始值 //迭代次数控制 for count := 0; count < n; count++ { nx := initX(len(x)) for i := 0; i < len(x); i++ { nxi := mr[i] for j := 0; j < len(mx[i]); j++ { if j != i { nxi = nxi + (-mx[i][j])*x[j] } } nxi = nxi / mx[i][i] // 迭代计算得到的下一个 xi 值 nx[i] = nxi } lc := make([]float64, 0) // 存储两次迭代结果之间的误差的集合 for i := 0; i < len(x); i++ { lc = append(lc, math.Abs(x[i]-nx[i])) } if max(lc) < c { fmt.Printf("一共迭代 %d 次\n", count+1) return nx, nil } x = nx } return nil, fmt.Errorf("超过最大迭代次数，未找到解")}func initX(n int) []float64 { x := make([]float64, n) // 迭代初始值 return x}func max(x []float64) float64 { var m float64 for _, a := range x { if a > m { m = a } } return m}
测试代码：
package mainimport ( "fmt" "log" "testing")func TestJacobi(t *testing.T) { mx := [][]float64{{8.0, -3.0, 2.0}, {4.0, 11.0, -1.0}, {6.0, 3.0, 12.0}} mr := []float64{20.0, 33.0, 36.0} ret, err := Jacobi(mx, mr, 100, 1e-20) if err != nil { log.Fatalln(err) } fmt.Printf("result: %v\n", ret)}
测试结果如下：
$ go test -v=== RUN TestJacobi一共迭代 39 次result: [3 2 1]
雅可比迭代的改进
推导
上述的公式有一个可以改进的地方，每次迭代时，**i **之前的元素可以用本次迭代已经算出的值来替代。
也就是最终的公式变成：
\[x^{(k+1)}_{i} = \frac{b_{i}}{a_{ii}} - \sum_{j=1}^{i-1}\frac{a_{ij}}{a_{ii}}x^{(k+1)}_{j} - \sum_{j=i+1}^{n}\frac{a_{ij}}{a_{ii}}x^{(k)}_{j}, i = 1,2,\cdots,n.
\]
实现方式
func Jacobi2(mx [][]float64, mr []float64, n int, c float64) ([]float64, error) { if len(mx) != len(mr) { return nil, fmt.Errorf("系数矩阵和值矩阵长度不匹配") } x := initX(len(mr)) // 迭代初始值 //迭代次数控制 for count := 0; count < n; count++ { nx := initX(len(x)) for i := 0; i < len(x); i++ { nxi := mr[i] for j := 0; j <= i-1; j++ { nxi = nxi + (-mx[i][j])*nx[j] } for j := i + 1; j < len(mx[i]); j++ { nxi = nxi + (-mx[i][j])*x[j] } nxi = nxi / mx[i][i] // 迭代计算得到的下一个 xi 值 nx[i] = nxi } lc := make([]float64, 0) // 存储两次迭代结果之间的误差的集合 for i := 0; i < len(x); i++ { lc = append(lc, math.Abs(x[i]-nx[i])) } if max(lc) < c { fmt.Printf("一共迭代 %d 次\n", count+1) return nx, nil } x = nx } return nil, fmt.Errorf("超过最大迭代次数，未找到解")}
测试函数：
func TestJacobi2(t *testing.T) { mx := [][]float64{{8.0, -3.0, 2.0}, {4.0, 11.0, -1.0}, {6.0, 3.0, 12.0}} mr := []float64{20.0, 33.0, 36.0} ret, err := Jacobi2(mx, mr, 100, 1e-20) if err != nil { log.Fatalln(err) } fmt.Printf("result: %v\n", ret)}
?
测试结果如下：
$ go test -v=== RUN TestJacobi一共迭代 39 次result: [3 2 1]--- PASS: TestJacobi (0.00s)=== RUN TestJacobi2一共迭代 19 次result: [3 2 1]--- PASS: TestJacobi2 (0.00s)PASSok sample 0.009s
?
计算结果一样，但是迭代次数减少了很多。
经验首页