2. Точные методы решения СЛАУ (LU-разложение, метод Холецкого).

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

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=

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("Ответы не совпадают")
+	}
+}

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
+}

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("_______________________________________________")
+	}
+}