Add calculation of qualitive criteria for a SLAE with matrix A and te…

…st on various matrices

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,47 @@
+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})
+	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(":(")
+	}
+	printTable(*A, *L, *TDense(*L))
+}

task2/square_method.go

@@ -0,0 +1,105 @@
+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 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
+}

task2/square_test.go

@@ -0,0 +1,39 @@
+package main
+
+import (
+	"fmt"
+	"math"
+	"testing"
+
+	"gonum.org/v1/gonum/mat"
+)
+
+func TestIllConditionedMatrix(t *testing.T) {
+	cases := []struct {
+		in   []float64
+		want bool
+	}{
+		{[]float64{1.0001, 0.9999, 0.9998, 0.9999, 1.0002, 0.9997, 0.9998, 0.9997, 1.0003}, true},
+		{[]float64{2.001, 0.999, 0.998, 0.999, 3.002, 0.997, 0.998, 0.997, 4.003}, true},
+		{[]float64{10}, true},
+		{[]float64{2, 0, 0, 0, 0, 9, 0, 0, 0, 0, 5, 0, 0, 0, 0, 4}, 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}, true},
+	}
+	for _, c := range cases {
+		n := math.Sqrt(float64(len(c.in)))
+		A := mat.NewDense(int(n), int(n), c.in)
+		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("get %q, want %q", fmt.Sprint(c.in), fmt.Sprint(c.want))
+		} else {
+			fmt.Println("L*L^T =")
+			matPrint(LLT)
+		}
+		fmt.Println("_______________________________________________")
+	}
+}