diff --git a/task2/go.mod b/task2/go.mod new file mode 100644 index 0000000..eb09af4 --- /dev/null +++ b/task2/go.mod @@ -0,0 +1,10 @@ +module workshop/task2 + +go 1.22.1 + +require ( + github.com/olekukonko/tablewriter v0.0.5 + gonum.org/v1/gonum v0.15.0 +) + +require github.com/mattn/go-runewidth v0.0.9 // indirect diff --git a/task2/go.sum b/task2/go.sum new file mode 100644 index 0000000..43f9b48 --- /dev/null +++ b/task2/go.sum @@ -0,0 +1,8 @@ +github.com/mattn/go-runewidth v0.0.9 h1:Lm995f3rfxdpd6TSmuVCHVb/QhupuXlYr8sCI/QdE+0= +github.com/mattn/go-runewidth v0.0.9/go.mod h1:H031xJmbD/WCDINGzjvQ9THkh0rPKHF+m2gUSrubnMI= +github.com/olekukonko/tablewriter v0.0.5 h1:P2Ga83D34wi1o9J6Wh1mRuqd4mF/x/lgBS7N7AbDhec= +github.com/olekukonko/tablewriter v0.0.5/go.mod h1:hPp6KlRPjbx+hW8ykQs1w3UBbZlj6HuIJcUGPhkA7kY= +golang.org/x/exp v0.0.0-20231110203233-9a3e6036ecaa h1:FRnLl4eNAQl8hwxVVC17teOw8kdjVDVAiFMtgUdTSRQ= +golang.org/x/exp v0.0.0-20231110203233-9a3e6036ecaa/go.mod h1:zk2irFbV9DP96SEBUUAy67IdHUaZuSnrz1n472HUCLE= +gonum.org/v1/gonum v0.15.0 h1:2lYxjRbTYyxkJxlhC+LvJIx3SsANPdRybu1tGj9/OrQ= +gonum.org/v1/gonum v0.15.0/go.mod h1:xzZVBJBtS+Mz4q0Yl2LJTk+OxOg4jiXZ7qBoM0uISGo= diff --git a/task2/main.go b/task2/main.go new file mode 100644 index 0000000..11fa373 --- /dev/null +++ b/task2/main.go @@ -0,0 +1,58 @@ +package main + +import ( + "fmt" + "os" + + "github.com/olekukonko/tablewriter" + "gonum.org/v1/gonum/mat" +) + +// подсчет критериев и вывод таблицы в консоль +func printTable(A, L, LT mat.Dense) { + t := tablewriter.NewWriter(os.Stdout) + t.SetHeader([]string{"Матрица", "Спектральный", "Объемный", "Угловой"}) + data := [][]string{ + {"A", fmt.Sprint(spectralCriterion(A)), fmt.Sprint(volumeCriterion(A)), fmt.Sprint(angularCriterion(A))}, + {"L", fmt.Sprint(spectralCriterion(L)), fmt.Sprint(volumeCriterion(L)), fmt.Sprint(angularCriterion(L))}, + {"L^T", fmt.Sprint(spectralCriterion(LT)), fmt.Sprint(volumeCriterion(LT)), fmt.Sprint(angularCriterion(LT))}, + } + for _, v := range data { + t.Append(v) + } + t.SetRowLine(true) + t.Render() +} + +func main() { + A := mat.NewDense(4, 4, []float64{5, 2, 3, 4, 2, 4, 2, 2, 3, 2, 8, 2, 4, 2, 2, 9}) + B := mat.NewDense(4, 1, []float64{4, 2, 0, 5}) + L := calculateL(A) + fmt.Println("Матрица L:") + matPrint(L) + // matPrint((*TDense(*L)).T().T()) + fmt.Println("") + fmt.Println("Исходная матрица:") + matPrint(A) + fmt.Println("") + ok, newA := check(*A, *L) + fmt.Println("Вычисленная L*L^T:") + matPrint(newA) + fmt.Println("") + if ok { + fmt.Println("A = L*L^T!") + } else { + fmt.Println("A != L*L^T :(") + } + printTable(*A, *L, *TDense(*L)) + fmt.Println("Решение системы с помощью библиотеки:") + x, xH, eps := matSolve(A, B), sqrtSolve(*L, *B), 1e-12 + matPrint(&x) + fmt.Println("\nРешение системы c помощью разложения Холецкого") + matPrint(&xH) + if mat.EqualApprox(&x, &xH, eps) { + fmt.Printf("\nОтветы совпадают с точностью %s", fmt.Sprint(eps)) + } else { + fmt.Println("Ответы не совпадают") + } +} diff --git a/task2/square_method.go b/task2/square_method.go new file mode 100644 index 0000000..aa57750 --- /dev/null +++ b/task2/square_method.go @@ -0,0 +1,122 @@ +package main + +import ( + "fmt" + "math" + + "gonum.org/v1/gonum/mat" +) + +// вывод матрицы в консоль +func matPrint(X mat.Matrix) { + fa := mat.Formatted(X, mat.Prefix(""), mat.Squeeze()) + fmt.Printf("%v\n", fa) +} + +// вычисление матрицы L по исходной матрице А +func calculateL(A *mat.Dense) (L *mat.Dense) { + n, _ := A.Dims() + res := mat.NewDense(n, n, nil) + res.Set(0, 0, math.Sqrt(A.At(0, 0))) + var el, sum float64 + for i := 0; i < n; i++ { + for j := 0; j <= i; j++ { + sum = 0 + if j != 0 { + for k := 0; k < j; k++ { + sum += res.At(i, k) * res.At(j, k) + } + } + if i == j { + el = math.Sqrt(A.At(i, j) - sum) + } else { + el = (A.At(i, j) - sum) / res.At(j, j) + } + res.Set(i, j, el) + } + } + return res +} + +// проверка разложения: L*L^T == A +func check(A, L mat.Dense) (ok bool, res *mat.Dense) { + n, _ := A.Dims() + res = mat.NewDense(n, n, nil) + L1 := L.T() + res.Mul(&L, L1) + return mat.EqualApprox(&A, res, 1e-12), res +} + +// решение системы методом Холецкого +func sqrtSolve(L, B mat.Dense) (x mat.Dense) { + y := matSolve(&L, &B) + x = matSolve(L.T(), &y) + return x +} + +// ____________________________________________________ +// вычисление спектрального критерия по матрице +func spectralCriterion(A mat.Dense) (res float64) { + var A1 mat.Dense + err := A1.Inverse(&A) + if err != nil { + fmt.Printf("A is not invertible: %v", err) + return + } + res = A.Norm(2) * A1.Norm(2) + return res +} + +// вычисление объемного критерия по матрице +func volumeCriterion(A mat.Dense) (res float64) { + rows, cols := A.Dims() + res = 1.0 + cur := 0.0 + for i := 0; i < rows; i++ { + for j := 0; j < cols; j++ { + el := A.At(i, j) + cur += el * el + } + res *= math.Sqrt(cur) + cur = 0 + } + res = res / math.Abs(mat.Det(&A)) + return res +} + +// вычисление углового критерия по матрице +func angularCriterion(A mat.Dense) (res float64) { + var C mat.Dense + err := C.Inverse(&A) + if err != nil { + fmt.Printf("A is not invertible: %v", err) + return + } + rows, _ := A.Dims() + for i := 0; i < rows; i++ { + res = math.Max(res, math.Abs(mat.Norm(A.RowView(i), 2))*math.Abs(mat.Norm(C.ColView(i), 2))) + } + return res +} + +// Транспонирование матрицы в формате Dense +func TDense(A mat.Dense) (res *mat.Dense) { + n, _ := A.Dims() + re := mat.NewDense(n, n, nil) + for i := 0; i < n; i++ { + for j := 0; j < n; j++ { + re.Set(i, j, A.At(j, i)) + } + } + return re +} + +// нахождение решения СЛАУ вида Ax=b, возвращает вектор X +func matSolve(A mat.Matrix, b mat.Matrix) (X mat.Dense) { + err := X.Solve(A, b) + if err != nil { + fmt.Printf("Ошибка при решении СЛАУ: %v", err) + return + } + return X +} diff --git a/task2/square_test.go b/task2/square_test.go new file mode 100644 index 0000000..a988b70 --- /dev/null +++ b/task2/square_test.go @@ -0,0 +1,52 @@ +package main + +import ( + "fmt" + "math" + "testing" + + "gonum.org/v1/gonum/mat" +) + +func TestDifferentConditionedMatrix(t *testing.T) { + cases := []struct { + inA []float64 + inB []float64 + want bool + }{ + {[]float64{1.0001, 0.9999, 0.9998, 0.9999, 1.0002, 0.9997, 0.9998, 0.9997, 1.0003}, []float64{1, 2, 3}, true}, + {[]float64{2.001, 0.999, 0.998, 0.999, 3.002, 0.997, 0.998, 0.997, 4.003}, []float64{4.1, 0.5, 2.6}, true}, + {[]float64{10}, []float64{2}, true}, + {[]float64{2, 0, 0, 0, 0, 9, 0, 0, 0, 0, 5, 0, 0, 0, 0, 4}, []float64{4, 4.5, 5, 1}, true}, + {[]float64{10, 1, 2, 3, 4, 5, 1, 20, 6, 7, 8, 9, 2, 6, 30, 10, 11, 12, 3, 7, 10, 40, 13, 14, 4, 8, 11, 13, 50, 15, 5, 9, 12, 14, 15, 60}, []float64{9, 1, 3, 8, 12, 6}, true}, + } + for _, c := range cases { + n := math.Sqrt(float64(len(c.inA))) + A := mat.NewDense(int(n), int(n), c.inA) + B := mat.NewDense(int(n), 1, c.inB) + fmt.Println("Матрица:") + matPrint(A) + fmt.Println("") + fmt.Printf("Спектральный критерий обусловленности матрицы: %f\n\n", spectralCriterion(*A)) + L := calculateL(A) + got, LLT := check(*A, *L) + if got != c.want { + t.Errorf("Полученные матрицы не равны") + } else { + fmt.Println("L*L^T =") + matPrint(LLT) + } + fmt.Println("\nРешение системы с помощью библиотеки:") + x, xH, eps := matSolve(A, B), sqrtSolve(*L, *B), 1e-10 + matPrint(&x) + fmt.Println("\nРешение системы c помощью разложения Холецкого") + matPrint(&xH) + if mat.EqualApprox(&x, &xH, eps) { + fmt.Printf("\nОтветы совпадают с точностью %s\n", fmt.Sprint(eps)) + } else { + fmt.Println("Ответы не совпадают") + t.Errorf("Ответы не совпадают") + } + fmt.Println("_______________________________________________") + } +}