最近遇到的一个求解雅可比迭代的问题，这个计算方法在 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 main
import (
"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 main
import (
"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)
PASS
ok sample 0.009s
计算结果一样，但是迭代次数减少了很多。