Матрицы и векторы: вычисление обратной матрицы / forpes.ru

Главная
Матрицы и векторы: вычисление обратной матрицы

Матрицы и векторы: вычисление обратной матрицы +7

30.12.2025 12:04

MishIvan 7 9500 Источник

Вычисление обратной матрицы квадратной матрицы высокого порядка, а именно, вычисление алгебраических дополнений и определителя матрицы займёт большое количество машинных ресурсов. Поэтому обращение матрицы при помощи алгебраических дополнений нерационально. В этой статье для вычисления обратной матрицы были использованы нижеследующие методы:

Решение системы AX = E, где A, X, E - квадратные матрицы порядка n, X - обратная A матрица, E - единичная матрица, $E_{ij} = \begin{cases}1 & i = j\\0 & i \neq j\end{cases}$
Использование LU разложения матрицы A = LU [1,2,3]

Обращение матриц было реализовано в виде шаблонов функций шаблона класса MATRIX<T>. В данной статье обсуждается решение вышеуказанными методами обращения квадратной матрицы порядка n и полученные результаты.

Решение системы AX = E

Решение матричного уравнения AX = E для матрицы порядка n A сводится к решению n систем линейных уравнений:

$\begin{matrix} A\cdot \bar{X_{j}} = \bar{e_{j}},\quad j = 1...n \\ \bar{X_{j}} = \begin{bmatrix}X_{1j} \\X_{2j} \\ ... \\ X_{nj} \end{bmatrix},\quad j-ый\;столбец \;матрицы\;X_{ij} \\ \bar{e_{j}},\quad j-ый\;столбец\;матрицы\; E \end{matrix}$

Код реализации решения приводится ниже:

/// <summary>
/// Формирование матриц alpha и gamma в компактной схеме исключения
/// Матрица alpha верхняя треугольная, gamma - нижняя треугольная
/// </summary>
/// <param name="alpha">формируемые матрицы</param>
/// <returns>значение определителя матрицы</returns>
template <typename T>
T MATRIX<T>::FormMatrixCompactScheme(MATRIX& alpha)
{
	if (this->m_columns != this->m_rows) return NAN;
	int n = this->m_rows;

	for (int i = 0; i < n; i++)
	{
		*(alpha.m_data + i * n) = *(this->m_data + i * n);
		if (i > 0) *(alpha.m_data + i) = *(this->m_data + i) / *this->m_data;
	}

	T sum = 0;
	int k = 1, i = 1;
	while (i < n)
	{
		if (k >= n)
		{
			k = 1; i++;
			if (i >= n) break;
		}
		if (i >= k)
		{
			sum = 0.0;
			for (int j = 0; j <= k - 1; j++)
			{
				sum += *(alpha.m_data + i * n + j) * *(alpha.m_data + j * n + k);
			}
			*(alpha.m_data + i * n + k) = *(this->m_data + i * n + k) - sum;
		}
		else
		{
			sum = 0.0;
			for (int j = 0; j <= i - 1; j++)
			{
				sum += *(alpha.m_data + i * n + j) * *(alpha.m_data + j * n + k);
			}
			if (abs(*(alpha.m_data + i * n + i)) < 1.0e-18)
			{
				return 0.0;
			}
			*(alpha.m_data + i * n + k) = (*(this->m_data + i * n + k) - sum) / *(alpha.m_data + i * n + i);
		}
		k++;
	}

	// вычисление определителя
	T det = 1;
	for (int i = 0; i < n; i++)
		det *= *(alpha.m_data + i * n + i);
	return det;
}

/// <summary>
/// Решение системы линейных алгебраических уравнений A*x = b с 
/// верхней треугольной матрицей
/// </summary>
/// <param name="A">верхняя треугольная матрица системы</param>
/// <param name="b">вектор правой части системы</param>
/// <param name="x">решение системы</param>
template <typename T>
void TriangleSolve(MATRIX<T>& A, VECTOR<T>& b, VECTOR<T>& x)
{
	int n = b.m_size;
	VECTOR<T> beta(n);
	T sum = 0;
	*beta.m_data = *b.m_data / *A.m_data;
	for (int i = 1; i < n; i++)
	{
		sum = 0.0;
		for (int j = 0; j <= i - 1; j++)
			sum += *(A.m_data + i * n + j) * *(beta.m_data + j);
		*(beta.m_data + i) = (*(b.m_data + i) - sum) / *(A.m_data + i * n + i);
	}

	// решение системы уравнений с труегольной матрицей		   
	*(x.m_data + n - 1) = *(beta.m_data + n - 1);
	for (int i = n - 2; i >= 0; i--)
	{
		sum = 0.0;
		for (int j = n - 1; j > i; j--)
			sum += *(A.m_data + i * n + j) * *(x.m_data + j);
		*(x.m_data + i) = *(beta.m_data + i) - sum;
	}

}
/// <summary>
/// Копировать вектор в  j-ю колонку матрицы
/// </summary>
/// <param name="v">вектор</param>
/// <param name="j">номер колонки</param>
template <typename T>
void MATRIX<T>::CopyColumn(VECTOR<T>& v, int j)
{
	if (m_columns < j)
	{
		return;
	}
	if (m_rows != v.size())
	{
		return;
	}
	for (int i = 0; i < m_rows; i++)
		*(m_data + i * m_columns + j) = v[i];
}

/// <summary>
/// Вычисляет первую норму матрицы
/// </summary>
/// <returns>значение первой нормы матрицы</returns>
template <typename T>
T MATRIX<T>::normI()
{
	T norm = 0;
	for (int i = 0; i < m_rows; i++)
	{
		T sum = 0;
		for (int j = 0; j < m_columns; j++)
			sum += abs(*(m_data + i * m_columns + j));
		if (sum > norm) norm = sum;
	}
	return norm;
}

/// <summary>
/// Обращение марицы при помощи решения систем уравнений A*X = E,
/// столбцы матрицы X - столбцы обратной к A матрицы,
/// E - единичная матрица
/// </summary>
/// <returns>обратную матрицу</returns>
template <typename T>
MATRIX<T> MATRIX<T>::Invert()
{
	MATRIX A(m_rows, m_columns);
	if (m_rows != m_columns)
		return A;

	MATRIX alpha(m_rows, m_columns);
	T det = FormMatrixCompactScheme(alpha);

	if (abs(det) < 1.0e-36)
		return A;

	for (int i = 0; i < m_rows; i++)
	{
		VECTOR<T> e(m_rows), x(m_rows);
		e[i] = 1;
		TriangleSolve(alpha, e, x);
		A.CopyColumn(x, i);
	}

	// итерационное уточнение обратной матрицы
	if (m_rows >= MIN_SIZE_INVERSION)
	{
		MATRIX E(m_rows, m_columns);
		for (int i = 0; i < m_rows; i++)
			*(E.m_data + i * m_columns + i) = 1;
		int iter = 0;
		T norm;
		do       
        {
			MATRIX R(E - *this * A);
			A = A * (E + R);
			norm = R.normI();
			if (++iter > MAX_ITER_NUMBER) break;
		} while (norm < 1.0e-17);
	}

	return A;
}

Перегрузка операторов *, =, +, - реализована в шаблоне класса MATRIX<T>. Решение системы уравнений осуществлялось с использованием компактной схемы исключения [4].

Обращение матрицы с использованием LU разложения матрицы

Разложение матрицы A = LU на произведение нижней треугольной L и верхней треугольной U может быть использовано не только для решения системы линейных алгебраических уравнений, но и для обращения матрицы. Описание метода обращения матрицы при помощи LU разложения приведён в источниках [2,3] и заключается в следующем:

$\begin{matrix} A = L\cdot U \Leftrightarrow A^{-1} = U^{-1}\cdot L^{-1} \\ U\cdot A^{-1} = L,\;A^{-1}\cdot L=U^{-1} \end{matrix}$

Обозначая искомые элементы матрицы $A^{-1}$ через $x_{ij}$ и учитывая, что треугольные матрицы при обращении сохраняют свою структуру, последние соотношения можно записать в виде:

$\begin{matrix} \begin{bmatrix}u_{11} & u_{12} & ... & u_{1n} \\ 0 & u_{22} & ... & u_{2n} \\ ... & .... & ... & ... \\ 0 & 0 & 0 & u_{nn} \end{bmatrix}\cdot \begin{bmatrix}x_{11} & x_{12} & ... & x_{1n} \\ x_{21} & x_{22} & ... & x_{2n} \\ ... & .... & ... & ... \\ x_{n1} & x_{n2} & ... & x_{nn} \end{bmatrix} = \begin{bmatrix}1 & 0 & ... & 0 \\ * & 1 & ... & 0 \\ ... & .... & ... & ... \\ * & * & ... & 1 \end{bmatrix} \\ \begin{bmatrix}x_{11} & x_{12} & ... & x_{1n} \\ x_{21} & x_{22} & ... & x_{2n} \\ ... & .... & ... & ... \\ x_{n1} & x_{n2} & ... & x_{nn} \end{bmatrix}\cdot \begin{bmatrix}1 & 0 & ... & 0 \\ l_{21} & 1 & ... & 0 \\ ... & .... & ... & ... \\ l_{n1} & l_{n2} & ... & 1 \end{bmatrix} = \begin{bmatrix}* & * & ... & * \\ 0 & * & ... & * \\ ... & .... & ... & ... \\ 0 & 0 & ... & * \end{bmatrix} \end{matrix}$

где символом '*' обозначены элементы, знание значений которых для нахождения обратной матрицы не требуется.

Из этих соотношений можно получить выражения для коэффициентов обратной матрицы $x_{ij}$ :

$\begin{matrix} x_{ij} = -\frac{1}{u_{ii}}\sum\limits_{k=i+1}^n u_{ik}x_{kj}\;(i<j) \\ x_{ji}= \frac{1}{u_{jj}}\sum\limits_{k=j+1}^n (1 - u_{jk}x_{kj})\;(i=j) \\ x_{ij} = -\sum\limits_{k=j+1}^n x_{ik}l_{kj}\;(i > j) \end{matrix}$

Код реализации решения приводится ниже:

/// <summary>
/// Декомпозиция матрицы A = LU
/// L - нижняя треугольная матрица, на главной диагонали которой расположены единицы
/// U - верхняя треугольная матрица
/// </summary>
/// <param name="alfa">Матрицы L и U, ниже главной диагонали которой 
/// расположены внедиагональные элементы L, 
/// на главной диагонали и выше расположены элементы матрицы U</param>
/// <returns>Определитель матрицы</returns>
template <typename T>
T MATRIX<T>::LUDecomposition(MATRIX& alfa)
{
	T det = 0;
	if (m_rows != m_columns) return det;

	// декомпозиция матрицы
	for (int i = 0; i < m_rows; i++)
		for (int j = 0; j < m_columns; j++)
		{
			T val = 0;
			if (i <= j)
			{
				for (int k = 0; k <= i - 1; k++)
					val += *(alfa.m_data + i * alfa.m_columns + k) * *(alfa.m_data + k * alfa.m_columns + j);
				*(alfa.m_data + i * alfa.m_columns + j) = *(m_data + i * m_columns + j) - val;
			}
			else
			{
				for (int k = 0; k <= j - 1; k++)
					val += *(alfa.m_data + i * alfa.m_columns + k) * *(alfa.m_data + k * alfa.m_columns + j);
				if (abs(*(alfa.m_data + j * alfa.m_columns + j)) < 1.0e-36) return 0;
				*(alfa.m_data + i * alfa.m_columns + j) = (*(m_data + i * m_columns + j) - val) / *(alfa.m_data + j * alfa.m_columns + j);

			}
		}
	det = 1;
	for(int i = 0; i < alfa.m_rows; i++)
		det *= *(alfa.m_data + i * alfa.m_columns + i);
	return det;
}

/// <summary>
/// Вычисляет первую норму матрицы
/// </summary>
/// <returns>значение первой нормы матрицы</returns>
template <typename T>
T MATRIX<T>::normI()
{
	T norm = 0;
	for (int i = 0; i < m_rows; i++)
	{
		T sum = 0;
		for (int j = 0; j < m_columns; j++)
			sum += abs(*(m_data + i * m_columns + j));
		if (sum > norm) norm = sum;
	}
	return norm;
}
/// <summary>
/// Обращение матриы с пощью LU разложения
/// </summary>
/// <returns>обратную матрицу</returns>
template <typename T>
MATRIX <T> MATRIX <T>::InvertLT()
{
	MATRIX D(m_rows, m_columns), alpha(m_rows, m_columns);
	if (m_rows != m_columns) // матрица должна быть квадратной
		return D;
	T det = LUDecomposition(alpha); // декомпозиция A = LU
	if (abs(det) < 1.0e-36)
		return D;
	// формирование элементов обратной матрицы
	int n = m_rows;
	for (int i = n-1; i >= 0; i--)
		for(int j = n-1;j >= 0; j--)
		{
			T sum = 0;
			if (i < j)
			{
				for (int k = i + 1; k < n; k++)
					sum -= *(alpha.m_data + i * alpha.m_columns + k) * *(D.m_data + k * D.m_columns + j);
				*(D.m_data + i * D.m_columns + j) =  sum / *(alpha.m_data + i * alpha.m_columns + i);

			}
			else if (i == j)
			{
				for (int k = j + 1; k < n; k++)
					sum += *(alpha.m_data + j * alpha.m_columns + k) * *(D.m_data + k * D.m_columns + j);
				*(D.m_data + j * D.m_columns + i) = (1-sum) / *(alpha.m_data + j * alpha.m_columns + j);

			}
			else
			{
				for (int k = j + 1; k < n; k++)
					sum -= *(alpha.m_data + k * alpha.m_columns + j) * *(D.m_data + i * D.m_columns + k);
				*(D.m_data + i * D.m_columns + j) = sum;

			}
		}
		// итерационное уточнение обратной матрицы
		if (m_rows >= MIN_SIZE_INVERSION)
		{
			MATRIX E(m_rows, m_columns);
			for (int i = 0; i < m_rows; i++)
				*(E.m_data + i * m_columns + i) = 1;
			int iter = 0;
			T norm;
			do
			{
				MATRIX R(E - *this * D);
				D = D * (E + R);
				norm = R.normI();
				if (++iter > MAX_ITER_NUMBER) break;
			} while (norm < 1.0e-17);
		}

	return D;
}

Результаты расчётов

Расчёты по обращению квадратной матрицы порядка n проводились обоими методами для разных значений n. Исходная матрица A порядка n заполнялась при помощи генератора случайных чисел в интервале от -1000 до 1000. В качестве меры погрешности вычислений была использована I-норма

$\parallel A\parallel_{I} = \begin{matrix}max & \sum\limits_{j=1}^n \mid a_{ij}\mid \\ 0\leq i \leq n\end{matrix}$

матрицы $E - A^{-1}A$ , где E - единичная матрица.

Расчёты проводились без итерационного уточнения и с итерационным уточнением обратной матрицы.

Итерационное уточнение осуществлялось следующим образом [2]. Пусть $D_{0}$ первое приближение матрицы $A^{-1}$ . Если при этом $\parallel{R_{0}}\parallel_{I}\leq k < 1,\quad где\; R_{0} = E-AD_{0}$ , то можно построить итерационный процесс:

$\begin{matrix} D_{1} = D_{0}(E-R_{0}),\quad R_{1} = E- AD_{1} \\ D_{2} = D_{1}(E-R_{1}),\quad R_{2} = E- AD_{2} \\ ...................................\\ D_{m} = D_{m-1}(E-R_{m-1}),\quad R_{m} = E- AD_{m} \end{matrix}$

Непосредственной проверкой можно убедиться в том, что [2]:

$R_{m}=R_{m-1}^{2}=R_{m-2}^{4}=...=R_{0}^{2^{m}}$

Откуда далее следует оценка погрешности [2]:

$\parallel{D_{m} - A^{-1}}\parallel_{I}=\parallel{D_{0}}\parallel_{I}\cdot \frac{k^{2^{m}}}{1-k}$

Расчёт без итерационного уточнения

Результаты расчётов без итерационного уточнения матрицы $A^{-1}$ приведены в нижеследующей таблице:

Порядок матрицы A	Решение уравнения AX=E		Декомпозиция A=LU
Порядок матрицы A	Время, сек	Погрешность	Время, сек	Погрешность
50	0,0007115	$1,82559\cdot 10^{-11}$	0,0003418	$1,20276\cdot 10^{-9}$
100	0,0030217	$1,38811\cdot 10^{-10}$	0,0024037	$2,69221\cdot 10^{-8}$
200	0,0238938	$8,36936\cdot 10^{-10}$	0,0193017	$2,20405\cdot 10^{-6}$
500	0,349852	$5,84951\cdot 10^{-9}$	0,288407	$6,42369\cdot 10^{-6}$
1000	2,83838	$4,17686\cdot 10^{-7}$	2,36379	0,0509241
2000	27,9798	$2,53471\cdot 10^{-7}$	33,5439	0,00136141

Расчёт с итерационным уточнением

По приведённым результатам расчёта без итерационного уточнения было сделано заключение о нецелесообразности итерационного уточнения обращения квадратной матрицы порядка ниже 200. Результаты расчётов с итерационным уточнением матрицы $A^{-1}$ приведены в нижеследующей таблице:

Порядок матрицы A	Решение уравнения AX=E		Декомпозиция A=LU
Порядок матрицы A	Время, сек	Погрешность	Время, сек	Погрешность
200	0,0981653	$5,85991\cdot 10^{-12}$	0,0836408	$5,76079\cdot 10^{-12}$
500	1,44216	$1,13022\cdot 10^{-10}$	1,31034	$1,13075\cdot 10^{-10}$
1000	11,4696	$1,70632\cdot 10^{-10}$	10,9775	$2,17188\cdot 10^{-10}$
2000	106,481	$6,39272\cdot 10^{-10}$	109,584	$6,52071\cdot 10^{-10}$

Выводы

Обращение матрицы без итерационного уточнения метод LU декомпозиции показывает лучшие результаты по времени, а метод решения системы AX = E демонстрирует лучшие результаты по точности. При этом, с ростом порядка квадратной матрицы разница в точности расчёта методом LU декомпозиции и методом решения системы AX = E возрастает на несколько порядков, а при обращении матрицы порядка 2000 метод решения системы AX = E демонстрирует лучшие результаты как по времени, так и по точности.

Сравнительный анализ обращения матрицы с итерационным уточнением методом LU декомпозиции и методом решения системы AX = E показывает, что разница как по точности расчёта, так по времени несущественна, а итерационное уточнение вычисления обратной матрицы даёт хорошие результаты по точности методом LU в сравнении с решением без итерационного уточнения.

Список источников

Википедия. LU - разложение/ URL: https://ru.wikipedia.org/wiki/LU-разложение (дата посещения 18.12.2022)
Фаддеев, Д.К. Вычислительные методы линейной алгебры: Учебник/ Фаддеев Д.К, Фаддеева В.Н. — 4‑е изд.— СПб : Лань, 2009. — 736 с.
Численные методы на Python: LU разложение/ URL: https://orion1401.gitbooks.io/numerical_analysys_python/content/lu_razlozhenie.html (дата посещения 22.12.2025)
Корн, Г. Справочник по математике для научных работников и инженеров : определения, теоремы, формулы / Г. Корн, Т. Корн ; пер. с англ. под общ. ред. И. Г. Арамановича. — 3‑е изд. — Москва : Наука, 1973. — 832 с.

Комментарии (7)

Gay_Lussak
30.12.2025 13:17
#29323308
У вас шаблонная функция, а епсилон фиксированная? А вообще для этого есть std::numeric_limits<T>::epsilon
```
abs(*(alfa.m_data + j * alfa.m_columns + j)) < 1.0e-36
```
Что за магическое число?
```
while (norm < 1.0e-17);
```
Уверены, что тут никогда не будет деления на 0?
```
sum / *(alpha.m_data + i * alpha.m_columns + i);
```
А что произойдет, если я передам T = std::string?
```
template <typename T>
T MATRIX<T>::normI()
{
	T norm = 0;
```
Для копирования столбцов по-хорошему нужно перегружать оператор присваивания, а не писать отдельную функцию CopyColumn. Еще и VECTOR по неконстантной ссылке.
```
template <typename T>
void MATRIX<T>::CopyColumn(VECTOR<T>& v, int j)
{
```
Да и зачем вам свой matrix класс, если кроме обращения элементов тут почти ничего нет, не все операторы перегружены, безопасности исключений тут нет. Использовали бы тогда уж std::vector<std::vector>> там хотя-бы более безопасно работать.
Или на питоне, на порядки же проще да и не напишите вы на С++ быстрее numpy.

Ну и вы проверяли на численную стабильность? Обычно для этого используют LUP разложение. А у вас вроде такого нет.
А вообще смысл статьи? Нового ничего нет. Если учебный проект, то сделан он халтурно, как в плане кода, так и в плане математики. Особенности вычислений с плавающей точкой не учтены от слова совсем.

TimurZhoraev
30.12.2025 13:17
#29323680
Если порассуждать про какие-либо принципиально новые методы вычислений то может обратить внимание на численно-символические. Составим простой скрипт на Maxima или любой другой CAS:
```
Q:makelist(makelist(a[j,i],i,1,4),j,1,4);
A: apply ('matrix,Q);
B:A^^-1;
horner(num(B[1,1]));
horner(denom(B[1,1]));
```
Сформировать например, матрицу 4х4, взять ей обратную и посмотреть что получается с числителем и знаменателем, например, для первого элемента применив схему Горнера.

Исходная матрица
Числитель первого элемента обратной матрицы
Знаменатель первого элемента обратной матрицы
Разумеется, что это следует из определителей матриц-миноров, однако, можно выявить индексную арифметику для некоего цифрового фильтра с заданными коэффициентами и операцией умножения с накоплением. В этом случае знаменатель считается один раз (он общий для всех элементов обратной матрицы). А вот числители вычисляются для каждого элемента примерно так:

Вычисление числителя
Умножение с накоплением идёт по внутренним скобкам как коэффициент умноженный на отсчёт. Иными словами мы представляем матрицу в виде массивов отсчётов как это бы делали для цифрового фильтра - умножение на коэффициент. То есть фактически всё распадается на работу с индексами и формирование промежуточных массивов с результатами. И если расположить элементы как показано то получается что-то вроде "бабочки" для БПФ. Наверняка может быть матрицы с размерностью кратной степени 2 имеют такой сверхбыстрый алгоритм.

Nemoumbra
30.12.2025 13:17
#29325036
Ага, сиквел к прошлой статье подъехал. Зачем это на Хабре? Линал запрогали и выложили... Никаких новых идей, как будто лаба во вычматам.
Так, теперь код-ревью. Иду сверху вниз по файлам.
MatrixVectorTemplate.hpp - тут вижу using namespace std;. Вас за такое на части разорвать могут. Кстати, инклюдить надо не math.h, а cmath.
Дальше, у вас умерла кодировка в файле => русский язык выглядит, как тарабарщина.
Зачем метод класса inline int size()? Инлайнить компилятор будет и без ваших подсказок. А почему не помечено словом const? Хотите, чтобы у константных векторов нельзя было спросить размер?

Почему указатель m_data инициализируете числовой константой 0? У вас отняли nullptr? Ну хотя бы NULL, но всё равно плохо.
А почему не используете списки инициализации в конструкторах?
Из вашей реализации копи-конструктора я вижу, что у вас отвратительные знания основ ООП в C++, т.к. вы зачем-то проверяете, нет ли у вас там аллоцированного массива и предварительно удаляете, если есть.
В реализации оператора присваивания отсутствует защита от присваивания себе.
Почему-то вы считаете, что умножать векторы разной длины - норм (возвращаем ноль), а вот складывать - не норм (делаем throw "Âåêòîðû èìåþò ðàçíóþ ðàçìåðíîñòü";).
Почему бросаете не что-то, наследующеся от std::exception? Почему в блоке после выброса исключения какой-то ещё код?
Кстати, почему вы понапихали френдов в код? Обычно бинарные операторы, тип результата которых совпадает с типами обоих аргументов, реализуют через копию + составное присваивание...
А ещё странно, что вы где-то используете this для доступа к полям, а где-то - нет. Нужен общий стиль. C-style касты понапиханы... В целом, очень много непонятных названий, особенно всякие однобуквенные переменные. Дефайны в коде есть для констант - тоже плохо... Магические числа в рандомных местах (всякие эпсилоны в сравнениях).
std::rand используете вместо качественных генераторов, а в функции randomDouble при неудачном рнг делите на ноль. Кстати, непонятно, почему дабл в названии есть - это можно инстанцировать и для целых типов.

Уфф... Язык знаете плохо, общая задумка слабая, совсем непонятно, зачем работаете с голой памятью вместо уже имеющихся в языке классов.

BorisU
30.12.2025 13:17
#29325900
а теперь сравните скорость вашего велосипеда с каким-нибудь настоящим BLAS

Матрицы и векторы: вычисление обратной матрицы +7

Решение системы AX = E

Обращение матрицы с использованием LU разложения матрицы

Результаты расчётов

Расчёт без итерационного уточнения

Расчёт с итерационным уточнением

Выводы

Список источников

Комментарии (7)

Gay_Lussak

TimurZhoraev

Nemoumbra

BorisU