最近遇到的一個求解雅可比迭代的問題，這個計算方法在 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
計算結果一樣，但是迭代次數減少了很多。