Задача:
Для заданной матрицы [latex]A(n, n)[/latex] найти обратную [latex]A^{-1}[/latex], используя итерационную формулу:
[latex]A_{k}^{-1} = A_{k-1}^{-1}(2E-A A_{k}^{-1}),[/latex]
где [latex]E[/latex] — единичная матрица, [latex]A_{0}^{-1} = E[/latex]. Итерационный процесс заканчивается, если для заданной погрешности [latex]\varepsilon[/latex] справедливо:
[latex]\left| det(A A_{k}^{-1})-1 \right| \le \varepsilon[/latex]
Анализ задачи:
Прежде чем приступать к решению средствами языка C++, я создал прототип в системе компьютерной алгебры Wolfram Mathematica, с которым планировал сверяться при тестировании программы. Тем не менее, внимательный анализ показал, что с таким выбором начального приближения процесс уже на пятом шаге в значительной мере расходится даже для матриц размера 2*2. После уточнения условий и анализа дополнительного материала, посвященного методу Ньютона-Шульца, исходное приближение было изменено (по результатам исследования «An Improved Newton Iteration for the Generalized Inverse of a Matrix, with Applications», Victor Pan, Robert Schreiber, стр. 8):
[latex]A_{0}^{-1} =\frac { A }{ \left\| A \right\|_{1} \left\| A \right\| _{\infty } }[/latex], где [latex]{ \left\| A \right\| }_{ 1 }=\max _{ i }{ \sum _{ j=0 }^{ n-1 }{ \left| { a }_{ ij } \right| } } [/latex], [latex]{ \left\| A \right\| }_{ \infty }=\max _{ j }{ \sum _{ i=0 }^{ n-1 }{ \left| { a }_{ ij } \right| } }[/latex].
Эффективность предложенного подхода иллюстрируют результаты работы прототипа:
Следовательно, из пространства задачи можно переместиться в пространство решения и составить алгоритм реализации предложенного метода на языке C++.
Тесты:
[latex]n[/latex] | [latex]A[/latex] | [latex]A^{-1}[/latex] | Результат |
3 | 1 2 35 5 711 13 7 | -0.964771 0.430661 -0.0171830.723533 -0.447973 0.137884 0.172358 0.155200 -0.086211 | Пройден |
2 | 1 22 1 | -0.333067 0.666400 0.666400 -0.333067 | Пройден |
5 | 1 0 0 0 00 1 0 0 00 0 1 0 00 0 0 1 00 0 0 0 1 | 1 0 0 0 00 1 0 0 00 0 1 0 00 0 0 1 00 0 0 0 1 | Пройден |
3 | 1 2 34 5 67 8 9 | Матрица вырождена | Пройден |
Программный код:
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 |
#include <cstdio> #include <algorithm> #include <cmath> //очистить выделенную память void clear(double **arr, int n) { for(int i = 0; i < n; i++) delete[] arr[i]; delete[] arr; } //создать копию массива double** clone(double **arr, int n) { double **newArr = new double*[n]; for(int row = 0; row < n; row++) { newArr[row] = new double[n]; for(int col = 0; col < n; col++) newArr[row][col] = arr[row][col]; } return newArr; } //print the matrix void show(double **matrix, int n) { for(int row = 0; row < n; row++){ for(int col = 0; col < n; col++) printf("%lf\t", matrix[row][col]); printf("\n"); } printf("\n"); } //матричное умножение матриц double** matrix_multi(double **A, double **B, int n) { double **result = new double*[n]; //заполнение нулями for(int row = 0; row < n; row++){ result[row] = new double[n]; for(int col = 0; col < n; col++){ result[row][col] = 0; } } for(int row = 0; row < n; row++){ for(int col = 0; col < n; col++){ for(int j = 0; j < n; j++){ result[row][col] += A[row][j] * B[j][col]; } } } return result; } //умножение матрицы на число void scalar_multi(double **m, int n, double a){ for(int row = 0; row < n; row++) for(int col = 0; col < n; col++){ m[row][col] *= a; } } //вычисление суммы двух квадратных матриц void sum(double **A, double **B, int n) { for(int row = 0; row < n; row++) for(int col = 0; col < n; col++) A[row][col] += B[row][col]; } //вычисление определителя double det(double **matrix, int n) //квадратная матрица размера n*n { double **B = clone(matrix, n); //приведение матрицы к верхнетреугольному виду for(int step = 0; step < n - 1; step++) for(int row = step + 1; row < n; row++) { double coeff = -B[row][step] / B[step][step]; //метод Гаусса for(int col = step; col < n; col++) B[row][col] += B[step][col] * coeff; } //Рассчитать определитель как произведение элементов главной диагонали double Det = 1; for(int i = 0; i < n; i++) Det *= B[i][i]; //Очистить память clear(B, n); return Det; } int main() { //Исходная матрица, динамический двухмерный массив int n; scanf("%d", &n); double **A = new double*[n]; for(int row = 0; row < n; row++) { A[row] = new double[n]; for(int col = 0; col < n; col++) scanf("%lf", &A[row][col]); } /* Численное вычисление обратной матрицы по методу Ньютона-Шульца 1. Записать начальное приближение [Pan, Schreiber]: 1) Транспонировать данную матрицу 2) Нормировать по столбцам и строкам 2. Повторять процесс до достижения заданной точности. */ double N1 = 0, Ninf = 0; //норма матрицы по столбцам и по строкам double **A0 = clone(A, n); //инициализация начального приближения for(size_t row = 0; row < n; row++){ double colsum = 0, rowsum = 0; for(size_t col = 0; col < n; col++){ rowsum += fabs(A0[row][col]); colsum += fabs(A0[col][row]); } N1 = std::max(colsum, N1); Ninf = std::max(rowsum, Ninf); } //транспонирование for(size_t row = 0; row < n - 1; row++){ for(size_t col = row + 1; col < n; col++) std::swap(A0[col][row], A0[row][col]); } scalar_multi(A0, n, (1 / (N1 * Ninf))); //нормирование матрицы //инициализация удвоенной единичной матрицы нужного размера double **E2 = new double*[n]; for(int row = 0; row < n; row++) { E2[row] = new double[n]; for(int col = 0; col < n; col++){ if(row == col) E2[row][col] = 2; else E2[row][col] = 0; } } double **inv = clone(A0, n); //A_{0} double EPS = 0.001; //погрешность if(det(A, n) != 0){ //если матрица не вырождена while(fabs(det(matrix_multi(A, inv, n), n) - 1) >= EPS) //пока |det(A * A[k](^-1)) - 1| >= EPS { double **prev = clone(inv, n); //A[k-1] inv = matrix_multi(A, prev, n); //A.(A[k-1]^(-1)) scalar_multi(inv, n, -1); //-A.(A[k-1]^(-1)) sum(inv, E2, n); //2E - A.(A[k-1]^(-1)) inv = matrix_multi(prev, inv, n); //(A[k-1]^(-1)).(2E - A.(A[k-1]^(-1))) clear(prev, n); } //вывод матрицы на экран show(inv, n); } else printf("Impossible\n"); clear(A, n); clear(E2, n); return 0; } |
Программа доступна для тестирования по ссылке: http://ideone.com/7YphoX.
Алгоритм решения
При решении данной задачи оправдывает себя подход «разделяй и властвуй»: вычисление обратной матрицы подразумевает промежуточные этапы, корректной реализации которых следует уделить особое внимание. Наверняка стоило бы написать класс matrix и реализовать перегрузку операций, но задачу можно решить и без применения объектно-ориентированных средств, пусть от этого решение и потеряет в изящности.
Можно выделить следующие подзадачи:
- Инициализация динамического двухмерного массива и передача его в функцию.
- Вычисление определителя матрицы (с применением метода Гаусса).
- Вычисление нормы матрицы.
- Транспонирование матрицы.
- Умножение матрицы на скаляр.
- Матричное умножение матриц.
- Сложение двух матриц.
- Непосредственно итерационный процесс Ньютона-Шульца
Ниже приведено пояснение к подходу к реализации некоторых подзадач:
- Выделение памяти для хранения массива происходит не на этапе запуска программы, после компиляции. Для инициализации использован конструктор new
- Вычисление определителя можно разбить на два последовательных шага:
- Приведение матрицы к верхнетреугольному виду (методом Гаусса).
- Вычисление определителя как произведения элементов главной диагонали.
Если матрица вырождена, то дальнейшие вычисления не производятся.
- Нормы матрицы вычисляются как максимальные значения суммы элементов по столбцам и строкам соответственно.
- При транспонировании обмениваются местами элементы, симметричные главной диагонали.
- При умножении матрицы на скаляр каждый элемент матрицы умножается на соответствующее вещественное число.
- При перемножении двух квадратных матриц используется промежуточный массив для хранения результата вычислений.
- Сложение двух матриц аналогично попарному сложению элементов, расположенных на соответствующих позициях.
- Максимально допустимая погрешность для метода Ньютона-Шульца [latex]\varepsilon = 0.001[/latex]. Программа использует локальный массив для хранения [latex]A_{k}^{-1}[/latex], инициализация которого происходит в теле цикла.
Технические детали реализации:
При выполнении подзадач часто приходится использовать локальные массивы, так что для очистки выделенной под них памяти создана отдельная функция clear().
Отличная работа, но метки и ссылка на ideone всё же нужны.
Как раз вспомнил и зашел их добавить, когда увидел ваш комментарий. Исправлено.
Проверил на нескольких подходящих примерах. Работает отлично.
Вы разобрались с довольно сложным материалом. Молодец!
Зачтено