Для решения систем линейных алгебраических уравнений (СЛАУ) большой размерности, а также систем, имеющих разреженные матрицы, применение точных методов (например, метод Гаусса) не является целесообразным, так как сказывается ограниченность разрядной сетки и накапливается погрешность округления.
Система с разреженной матрицей получается, например, при составлении системы уравнений теплового баланса технологической схемы или же при численном расчете потенциалов из дифференциального уравнения Гельмгольца на равномерной сетке. Есть и другие примеры, которые приводят нас к необходимости решать СЛАУ с разряженной матрицей.
Рассмотрим здесь простейший метод решения линейных систем - метод простых итераций.
Так как я не люблю сложные теоретические усложнения, то начну с простого примера системы с тремя неизвестными.
Допустим, у вас есть следующая система уравнений:
![Здесь x1, x2, x3 - неизвестные, а a[i][j] и b[i] - известные постоянные](https://avatars.dzeninfra.ru/get-zen_doc/1584893/pub_6027115cc219c97e32ed9e8f_60271afb0454f6146aa766db/scale_1200)
В матричном виде эту систему можно было бы записать так:
Выразим из 1-го уравнения системы первую неизвестную x1
Выразим из 2-го уравнения системы вторую неизвестную x2
Выразим из 3-го уравнения системы третью неизвестную x3
Заметим, что элементы, лежащие на главное диагонали матрицы не должны быть нулевые. Ну а если такое случается, то надо переставлять местами уравнения, так, чтобы выполнялось условие неравенства нулю диагональных элементов.
Итак, все эти три уравнения мы можем записать общей формулой, если перейдем к индексам: первый индекс будет пробегать по строкам, а второй по столбцам. В нашей теории нумерация начинается с 1.
Обратите внимание, что в скобках, где находится сумма произведений, не встречается элемента с совпадающими индексами (тот а[k][k] который лежит на диагонали), поэтому при суммировании мы должны его пропустить.
В принципе, этот пропуск можно реализовать конструкцией continue при программировании алгоритма. Но если вам это не нравится, то общую формулу суммы можно переписать в другом виде:
Постараюсь сделать схему решения еще более понятной для программирования:
Важные уточнения:
columns - количество столбцов
k - текущий корень системы, совпадает с номером строки матрицы
i - номер столбца, индекс, по которому проводится суммирование
x[i] - старые корни, полученные на предыдущей итерации
x[k] - новый корень, которые считается на базе корней, полученных на предыдущей итерации
Ещё хочется обратить ваше внимание на то, что во многих языках программирования принято начинать нумерацию с нуля (!). Поэтому в наших формулах произойдет сдвиг индекса, это нужно не забыть исправить, чтобы не выйти за диапазон изменения допустимых значений индекса при записи решения в массив / список / вектор.
Иногда, в книгах по численным методам еще используется верхний индекс (чаще всего в скобках, чтоб не спутали со степенью), чтобы разграничить где решение на новой итерации, а где решение на предыдущей итерации.
У нас получается рекуррентная формула:
На мой взгляд, не нужно особо обращать внимание на верхние индексы, показывающие где новая итерация, а где старая. Так как начинающих это только сбивает с толку. В реальной программе это будет перезапись одной и той же переменной, БЕЗ как либо верхних индексов.
Приведу очень простой пример, почему на практике не нужны загромождения верхними индексами для рекуррентной формулы:
Теперь опять вернемся к нашей задаче. Итерационный процесс можно запустить до бесконечности и далее... Шучу, нельзя. Компилятор быстро спустит вас с небес на землю! Но когда останавливаться? Теория говорит нам, что можно остановиться при достижении необходимой точности.
Точность часто обозначается буквой экспилон ε, ну а в программировании взяли за правило обозначать eps ( сокращение от epsilon ), так как греческими буквами нам нельзя пользоваться при написании названий переменных.
Для остановки и проверки точности, а также завершения вычислений на основе вычисления расстояния между соседними элементами в последовательности итераций можно воспользоваться Евклидовой метрикой.
В этом случае, условие завершения вычислений следующее:
По сути, это что-то вроде относительной погрешности, учитывающий все решения, а именно среднеквадратичные отклонения, чтобы исключить влияние знака. Потому что при сложении абсолютных разниц у нас могла произойти иллюзия достижения необходимой точности и иллюзия маленького разброса, если бы одна разница была бы большой со знаком "+", а другая разница примерно такой же большой, но со знаком "-". Они бы просто друг друга аннигилировали и дали бы нам ложный результат. Именно поэтому используются квадраты отклонений, а потом уже накладывается корень.
Для того, чтобы мы могли считать Евклидову метрику на каждой итерации, нам необходимо сохранять вектор решения на предыдущей итерации. То есть считая новый набор корней на новой итерации, мы должны сохранить старый набор корней на предыдущей итерации, чтобы было потом с чем сравнивать.
На практике же, я бы посоветовал помимо погрешности ограничивать выполнения цикла каким-либо предельным количеством итераций. Например, чтобы максимальное число итераций было 100~500. Для адекватной сходимости этого достаточно за глаза. Так как, если метод сходится, то система достигает нужной точности примерно за 20-50 итераций (чаще всего, если точность eps = 0.0001).
Уже совсем другое дело, если решение расходится. Именно тогда вас спасает ограничение по числу итераций, чтобы не войти в бесконечный цикл.
Следует отметить, что скорость сходимости итерационного процесса выше для матриц, у которых элементы главной диагонали велики по сравнению с внедиагональными элементами. В связи с этим, перед началом численного решения задачи желательно преобразовать систему уравнений так, чтобы преобладали диагональные элементы.
Также, чтобы запустить рекуррентную формулу и начать любой метод итераций, нам нужно начальное приближение корней. Можно взять любое, но в разумных пределах. К примеру, можно использовать столбец свободных коэффициентов в матричном представлении системы. Или просто обнулить начальное приближение. Или везде поставить единицы. Если метод будет сходиться, то он в любом случае придет к корням с нужным приближением.
Порядок решения СЛАУ методом Якоби такой:
● Приведение системы уравнений к виду, в котором на каждой строчке выражено какое-либо неизвестное значение системы.
● Произвольный выбор нулевого решения, в качестве него можно взять вектор-столбец свободных членов.
● Производим подстановку произвольного нулевого решения в систему уравнений, полученную под пунктом 1.
● Осуществление дополнительных итераций, для каждой из которых используется решение, полученное на предыдущем этапе.
Итак, теперь всё сказано и можно приступать к реализации.
Я набросаю алгоритм на Python:
Библиотека с книгами для физиков, математиков и программистов
Репетитор IT mentor в VK
Репетитор IT mentor в Instagram
Репетитор IT mentor в telegram
