Решение систем алгебраических уравнений

 

Решение систем линейных уравнений – одна из наиболее часто встречающихся задач в научных и инженерных расчетах. Многие приложения приводят к этой задаче непосредственно, кроме того системы линейных уравнений могут возникать как один из этапов в решении более сложных проблем. При этом обычно число уравнений равно числу неизвестных. Такую систему можно записать в матричном виде A x=b, где А – заданная квадратная матрица порядка n, b – заданный вектор-столбец с n компонентами и x – неизвестный вектор-столбец с n компонентами. В свернутом виде система может быть записана как , i=1,2,…,n. Тогда в развернутом виде

a₁₁ x₁+ a₁₂ x₂ +…+ a_1n x_n =b₁,

a₂₁ x₁+ a₂₂ x₂ +…+ a_2n x_n =b₂, (78)

……………………………………

a_n1 x₁+ a_n2 x₂ +…+ a_nn x_n =b_n .

 

Решением системы является вектор X^*={x₁*, x₂*,…, x_n*}, при подстановке которого в каждое уравнение системы получаем истинные равенства. Для решения СЛАУ применяют прямые (точные) и итерационные (приближенные) методы. Метод решения задачи относится к классу точных, если предполагается отсутствие округлений и с его помощью можно найти решение в результате конечного числа арифметических и логических операций.

Среди прямых методов (используют конечное число преобразований системы) можно выделить методы Крамера, Гаусса, Гаусса–Жордана и др.

В частности, метод Гаусса (гауссовых исключений) состоит в приведении исходной системы путем последовательных исключений неизвестных (при помощи элементарных операций над строками: перестановка строк; умножение строки на число отличное от нуля; сложение строк) к эквивалентной системе с треугольной (ступенчатой) матрицей (прямой ход метода Гаусса):

 

x₁+ с₁₂ x₂ +…+ с_1n x_n =d₁,

c₂₂ x₂ +…+ c_2n x_n =d₂,

…………………

x_n =d_n,

 

далее ее решают по рекуррентным формулам

 

x_n =d_n, , i=n – 1,n – 2,…, 1.

 

Итерационные методы реализуют стратегию постепенного приближения к решению X^* путем построения соответствующей последовательности векторов X⁰, X¹,…, X^k. Вектор X⁰ называется начальным приближением, а переход от текущего k -го приближения (вектор X^k) к k+1 -му приближению – итерацией. Данный переход реализуется при помощи функционального оператора X^k+1 = F(X^k), вид которого определяется заданной системой и конкретным методом.

Для итерационных методов (метод простой итерации, метод Зейделя и др.) решение получают всегда приближенным, в то же время необходимо отметить ряд преимуществ, которыми обладают последние по сравнению с прямыми методами:

1) меньшие затраты времени при получении решения;

2) значительно меньшие погрешности округления;

3) самоисправление случайных ошибок на очередном шаге;

4) особенная эффективность при работе с разряженными матрицами (большое количество коэффициентов которых равно нулю), а также для систем большой размерности (10³–10⁶).

Метод простой итерации реализует переход к k+1 -му приближению, используя при определении каждой координаты нового вектора координат только предыдущего k- го приближения X^k+1 = F(X^k). Пусть дана СЛАУ вида

 

a₁₁ x₁+ a₁₂ x₂ +…+ a_1n x_n + a_1n+1 =0,

a₂₁ x₁+ a₂₂ x₂ +…+ a_2n x_n + a_2n+1 =0, (79)

……………………………………

a_n1 x₁+ a_n2 x₂ +…+ a_nn x_n + a_nn+1 =0.

 

Выразив из каждого уравнения соответствующую координату, получим

 

x₁=с₁₁ x₁+ с₁₂ x₂ +…+ с_1n x_n + с_1n+1,

x₂=с₂₁ x₁+ с₂₂ x₂ +…+ с_2n x_n + с_2n+1, (80)

……………………………………

x_n=с_n1 x₁+ с_n2 x₂ +…+ с_nn x_n + с_nn+1 .

Для обеспечения условий сходимости (сходимость, или устойчивость, метода обеспечивается при любом начальном приближении X⁰, если норма матрицы, составленной из коэффициентов только при неизвестных в системе (2) меньше единицы, т.е., например, при i= 1,2,…,n) при переходе от системы (1) к системе (2) можно воспользоваться следующими преобразованиями:

 

, , (81)

 

где .

Если после данных преобразований сходимость не обеспечивается (единичные случаи), исходную систему необходимо преобразовать в эквивалентную. Это достигается путем попарного сложения или вычитания уравнений.

Далее текущее приближение X^k подставляют в правую часть системы (2), тогда в левой части получают координаты нового уточненного приближения X^k+1. В качестве начального приближения X⁰ можно выбрать столбец свободных членов b или нулевой вектор. Окончанием итерационного процесса может служить выполнение условия

 

,

где e – заданная погрешность приближенного решения.

 

Метод Зейделя является модификацией метода простой итерации. В целях ускорения сходимости он реализует переход к k+1 -му приближению, используя при определении координат нового вектора не только предыдущее k- е приближение, но и вычисленные на данный момент координаты k+1- го приближения X^k+1 = F(X^k, X^k+1):

 

x₁^k+1=с₁₁ x₁^k+ с₁₂ x₂^k+…+ с_1n-1 x_n-1^k +с_1n x_n^k+ с_1n+1,

x₂^k+1=с₂₁ x₁^k+1+ с₂₂ x₂^k+…+ с_2n x_n^k + с_2n+1, (82)

……………………………………

x_n^k+1=с_n1 x₁^k+1+ с_n2 x₂^k+1+…+ с_nn x_n^k + с_nn+1 .

 

Дальнейшее ускорение сходимости последовательности приближений достигается в методе последовательной верхней релаксации. Для специального вида матриц, например симметричных, положительно определенных, имеется также семейство методов сопряженных градиентов или сопряженных направлений [9].

В целом для обеспечения сходимости итерационного метода достаточно, чтобы выполнялось условие || A ||<1 или || С ||<1 по какой-либо норме матрицы, согласованной с нормой векторов. Если для векторов x =(x₁,x₂,…,x_n) введена норма ||x||, то согласованной с ней нормой матриц называют величину . Чаще всего используют нормы векторов и соответствующих согласованных норм матриц, представленные в табл. 2.

 

Таблица 2

№ п/п	Норма вектора	Согласованная норма матрицы	Функция MathCAD
1.			norm1(A)
2.			normi(A)
3.		или с учетом ||A||2£||A|| e	norm2(A)     norme(A)

 

Примечание: максимальное собственное значение матрицы .

В общем случае система может включать m уравнений относительно n неизвестных

 

a₁₁ x₁+ a₁₂ x₂ +…+ a_1n x_n =b₁,

a₂₁ x₁+ a₂₂ x₂ +…+ a_2n x_n =b₂,

……………………………………

a_m1 x₁+ a_m2 x₂ +…+ a_{m n} x_n =b_m .

 

Тогда данная система совместна, если она имеет хотя бы одно решение, и несовместна, если не имеет ни одного.

 

Во многих случаях применение линейных моделей даже для приближенного описания реальных явлений уже недостаточно. Одной из наиболее известных классических нелинейных задач является нахождение нулей полинома. Среди источников систем нелинейных уравнений можно выделить задачи вычисления интегралов (квадратурные формулы Гаусса для вычисления площадей), задачи минимизации функции, задачи интерполяции данных в нелинейных моделях, усложненные модели популяций и т.д. [9]. При этом между линейными и нелинейными уравнениями можно выделить следующие фундаментальные отличия:

1. Любая невырожденная СЛАУ имеет единственное решение. Для нелинейных уравнений вещественных решений может не существовать вообще (как, например, для уравнения f(x)=sinx+2=0), или их может быть много, как в уравнении f(x)=sinx=0.

2. Любая СЛАУ может быть точно решена за фиксированное, конечное число арифметических операций (примерно n³/3 для метода исключений Гаусса). Только некоторые нелинейные уравнения специального вида могут быть решены точно за конечное число шагов (полиномы до 4-й степени).

3. Значительные трудности обеспечения гарантий существования решения или сходимости к нему для СНЛУ.

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

 

В общем случае систему нелинейных уравнений можно представить в виде F(X)=0, или

 

f₁ (x₁, x₂,…,x_n)=0,

f₂ (x₁, x₂,…,x_n)=0, (83)

………………………

f_n (x₁, x₂,…,x_n)=0.

 

Тогда решением системы будет вектор X^*={x₁*, x₂*,…, x_n*}, при подстановке которого в каждое уравнение системы получаем истинные равенства F(X^*)=0. С учетом использования приближенных методов понятие «решение» расширяют. Говорят, что приближенно удовлетворяет («решает») системе(у) F(X)=0, если , или приближенно равен решению X^*, если .

Решение отдельного нелинейного уравнения f(x)=0 (f(x) – непрерывная функция на интервале [a,b]) с учетом возможности множества решений включает два этапа. На первом, который называется отделением корней, находятся интервалы [a,b] _k, длина которых не больше единицы, а функция на концах имеет противоположные знаки. В качестве алгоритма можно использовать простой перебор значений аргумента с некоторым шагом h (x_i+1=x_i+h), вплоть до наступления моментов f(x_i+1)×f(x_i)<0. На втором для каждого выделенного k -го интервала уточняются значения корня одним из рассматриваемых далее методов.

Метод половинного деления (бисекции, дихотомии, проб) на каждом шаге итерации (сужения) делит текущий интервал [a,b] пополам: с=(a+b)/2. Если выполняется условие |f(c)|<e, то c является искомым корнем. В противном случае выбирается та часть интервала ([a,c] или [c,b]), на границах которого функция имеет различные знаки. Далее процедура повторяется. Для устранения ситуаций «зацикливание», а также ускорения процесса нахождения корня желательно ограничить число приближений к корню 30 итерациями.

Метод простой итерации (последовательных приближений) предполагает замену исходного уравнения вида f(x)=0 на равносильное ему уравнение x=g(x). Далее текущее k -е приближение x^k подставляют в правую часть преобразованного уравнения, получая в левой части более точное k+1 -е x^k+1=g(x^k). Процесс начинается с произвольного выбора начального приближения x⁰ и заканчивается при выполнении условия | x^k+1 – x^k | <e или | f(x^k) | <e. Контроль сходимости метода может быть осуществлен при помощи условия |g¢(x^k)|=|(g(x^k)–g(x^k-1))/(x^k – x^k-1)|<1, если оно не выполняется, необходимо прекращать вычисления.

Кроме того, могут быть использованы метод Ньютона (касательных), метод хорд (ложного положения, пропорциональных частей), их комбинация, а также метод секущих [9].

При решении систем нелинейных уравнений можно выделить метод простой итерации и метод Ньютона.

Метод простой итерации является развитием аналогичного метода для одного уравнения (аналогичен методу простой итерации для решения СЛАУ) и основан на допущении о возможности преобразования системы (47) к виду (48)

 

x₁=g₁ (x₁, x₂,…,x_n),

x₂=g₂ (x₁, x₂,…,x_n), (84)

………………………

x_n=g_n (x₁, x₂,…,x_n).

 

Успех процесса приближения во многом определяется удачным выбором начального вектора X^k=1 (он должен быть достаточно близким к истинному решению). Метод гарантированно сходится, если выполняется условие || J(x₁, x₂,…,x_n)||<1, где J(x₁, x₂,…,x_n) – матрица Якоби для функций g₁,…, g_n, а || J || – ее норма. Тогда для системы из двух уравнений

 

,

 

а начальное приближение должно удовлетворять условиям и .

 

Метод Ньютона реализует идею предварительной линеаризации исходной СНЛУ, он обеспечивает значительно более быструю сходимость и является наиболее распространенным методом. В основе его лежит разложение функции в ряд Тейлора. Тогда для системы (47), состоящей из двух уравнений для двух неизвестных, получим

 

(85)

 

где R₁ и R₂ – члены ряда более высоких порядков.

Если приращения Dx₁ и Dx₂ выбрать так, чтобы величины x₁ + Dx₁ и x₂ + Dx₂ были близки к корням, то левые части соответствующих уравнений системы (49) будут примерно равны нулю. Отбросив члены более высоких порядков, сведем задачу к решению системы линейных алгебраических уравнений

 

или ,

где решением является вектор DХ={Dx₁, Dx₂ }. Найденные значения Dx₁, Dx₂ используют в виде поправок к текущему приближению X^k+1 = X^k + DХ^k. Расчет прекращается при выполнении условия DХ^k<e.