最近遇到的一个求解雅可比迭代的问题，这个计算方法在 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\)
?
线性方程组可以改写为：
\[Dx = b - Rx \]雅可比迭代就是在每一次的迭代中，上一次算出的 x 被用在右侧，用来计算左侧新的x，
这个程序用公式表示如下：
每个元素就是：
\[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 **之前的元素可以用本次迭代已经算出的值来替代，
也就是最终的公式变成：
实作方式
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
?
计算结果一样，但是迭代次数减少了很多，
0 评论