4.2.2. Метод Крамера

Значения неизвестных x_i (i=1,2,...,n) могут быть получены по формулам Крамера:

x_i=detA_i/detA ,

где матрица А_i определяется из матрицы А заменой ее i-того столбца столбцом свободных членов. Такой способ решения линейных систем с n-неизвестными приводит к вычислению n+1 детерминантов порядка n, что очень трудоемко при достаточно большом числе n.

4.2.3. Метод Гаусса

Наиболее распространенным методом решения систем алгебраических уравнений является метод Гаусса, в основе которого лежит идея последовательного исключения неизвестных. Существуют разные вычислительные схемы, которые реализуют этот метод. Рассмотрим одну из них - схему единичного деления.

Пусть задана система линейных уравнений n-го порядка, детерминант которой отличный от нуля. Предположим, что матрица коэффициентов не имеет нулевых диагональных элементов. Если такие имеются, то соответствующей перестановкой строки их всегда можно сделать ненулевыми.

Метод Гаусса заключается в следующем:

1. Из всех уравнений, кроме первого, исключаются члены, которые содержат x₁. Для этого из второго, третьего,....., n-го уравнений системы почленно, включая правые части, отнимается первое уравнение, разделенное на а₁₁ и умноженное соответственно на а₂₁, а₂₂, а₂₃,.....,а_n1. В результате этой операции порядок всех уравнений, за исключением первого, понижается на единицу.

2. Вновь полученное второе уравнение делится на а₂₂ и аналогичным способом, начиная с третьего уравнения, исключаются все элементы, которые содержат х₂.

3. Повторяем эту процедуру n-1 раз, то есть каждый раз исключаем неизвестные из ниже расположенных уравнений. В результате получаем ступенчатую треугольную систему уравнений, эквивалентную начальной, в которой последнее уравнение получает только одну неизвестную.

4. Решение ступенчатой системы уравнений осуществляется путем последовательного расчета неизвестных, начиная с последнего уравнение.

Пример

Найти корни следующей системы уравнений:

2.0*x₁+1.0*x₂-0.1*x₃= 3.7

0.4*x₁+0.5*x₂-4.0*x₃=13.4

0.3*x₁-1.0*x₂+1.0*x₃= 1.3

Для решения системы используем метод Гаусса. Из второго и третьего уравнений исключим члены, которые содержат х₁. Для этого сначала поделим первое уравнение на а₁₁=2.0. Получим:

х₁+0.5*x₂-0.05*x₃=1.85

Потом умножим полученное уравнение на а₂₁ и а₃₁, отнимем его из второго и третьего уравнений соответственно. Таким образом, получим систему с двумя неизвестными:

0.3*x₂+4.02*x₃=12.66

1.15*x₂+1.015*x₃=0.745

Поделив первое уравнение полученной системы на а₂₂ и умножив его на а₃₂, отнимем его из второго уравнения этой системы.

16.425*x₃=49.275

Таким образом, эквивалентная система имеет вид:

x₁+0.5*x₂-0.05*x₃=1.85

x₂+13.4*x₃=42.2

16.425*x₃=49.275

Из полученной эквивалентной системы последовательно найдем:

x₃=3.0

x₂=42.2-13.4*3.0=2.0

x₁=1.85-0.5*2.0+0.05*0.3=1.0.

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

4.2.4. Блок-схема программы для решения систем линейных уравнений методом Гаусса

Блок-схему программы можно поделить условно на трех части.

Первая - реализует алгоритм исключения неизвестных; вторая - обратную подстановку; третья - определяет наибольший коэффициент при х и осуществляет перестановку уравнений в случае необходимости. В блок-схеме все отношения обозначаются одним символом m. Если вычисления в программе организованы верно, то каждый этап вычислений требует не более одного отношения. При анализе блок-схемы надо четко различать смысл индексов i, j, k: k - определяет номер того уравнения, которое отнимается от последних, а также номер того неизвестного, которое исключается из последних n-k уравнений, i - означает номер уравнения, из которого в данный момент исключается неизвестная, j - означает номер столбца.

Блок-схема последовательного исключения неизвестных показана на рис.4.8. Обратная подстановка для поиска значений неизвестных задается следующими формулами:

x_N=b_N^(N-1)/a_NN^(N-1)

x_N-1=[b_N-1^(N-2)-a_{N-1, N}^(N-2)*x_N]/a_{n-1, n-1}^(N-2)

..................................

x_j=(b_j^(j-1)-a_jn^(j-1)*x_n-...-a_jj^(j-1)*x_j+1]/a_jj ^(j-1)

для j=n-2, n-3, ..., 1.

Блок-схема обратной подстановки показана на рис.4.9. Заметим, что в блок-схеме (рис.4.8) все индексы в процессе вычислений увеличиваются, а в блок-схеме (рис.4.9) один из индексов, а именно i, уменьшается. Третья блок-схема перестановки уравнений (рис.4.10) входит в блок-схему (рис.4.8). Задача нахождения наибольшего коэффициента при x_k и перестановки уравнений при необходимости связана с уменьшением ошибки округления, которая возникает при расчетах с числами с плавающей запятой. Во-первых, перестановка может быть нужна, чтобы получить a_kk0 и уклониться от деления на 0. А также при использовании перестановки возможно уменьшение ошибки округления.

Рис.4.8 Рис.4.9

Сначала рассмотрим числовой пример.

Решим систему методом последовательного исключения:

3.241*10⁰*x₁+1.600*10²*x₂=1.632*10²

1.020*10⁴*x₁+1.540*10³*x₂=1.174*10⁴

Точное решение этой системы следующее:

x₁=1.000*10⁰;

x₂=1.000*10⁰;

Будем решать эти уравнения методом исключение в том порядке, в котором они записаны, используя при вычислении числа с 4 значащими цифрами в мантиссе. Так как а₁₁0, то первое и единственное отношение:

m=1.020*10⁴/(3.241*10⁰⁾)=3.147*10³

Превращенное второе уравнение имеет вид:

5.730*10^(-1)*x₁-5.020*10⁵*x₂=-5.019*10⁵

Коэффициент при х₁ должен был равняться 0, но ошибка округления не разрешает получить точный результат. Так как этот коэффициент не участвует в дальнейших вычислениях, то принимаем его равным 0 и образовываем обратную подстановку:

x ₂=-5.019*10⁵/(-5.020*10⁵)=9.998*10^(-1)

Из первого получаем: x₁=9.873*10^(-1)

Подставим уравнение так, чтобы максимальный по модулю элемент попал

на главную диагональ:

1.020*10⁴*x₁+1.540*10³*x₂=1.174*10⁴

3.241*10⁰*x₁+1.600*10²*x₂=1.632*10²

Теперь отношение

m=3.241*10⁰/(1.020*10⁴)=3.177*10^(-4)

Преобразованное второе уравнение:

0.000*x₂+1.595*10²*x₂=1.595*10²

Откуда x₂ = 1.000*10⁰ ; x₁=1.000*10⁰

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

ответы могут не содержать ни одной верной Рис.4.10

цифры. Поэтому при решении СЛР на ЭВМ

методом исключения должно предусмотреть возможность перестановки уравнений.

Итак, задача сводится к перестановке n-k+1 последних уравнений так, чтобы наибольший по модулю коэффициент при x_k попал на главную диагональ. Описанный способ решения СЛР называют методом главного элемента. При рассмотрении блок-схемы необходимо помнить, что при переходе к данному этапу расчетов индекс k имеет некоторое определенное значение, а i = k+1.

1. Сначала в программу вводится вспомогательный индекс l, которому предоставляется значение k.

2. Первое сравнение выполняется между а_k,l-м элементом, который лежит на главной диагонали, и следующим за ним a_i,k. Если |a_ik|>|a_lk|, то индексу l дают значение i, и дальнейшее сравнение выполняется со следующим элементом. Поэтому индекс l является номером максимального по модулю элемента. Индекс i пробегает по ходу программы значения от k+1 по n включительно, и в конце этого цикла индекс l определяет номер наибольшего по модулю элемента в k-ом столбце. При перестановке уравнения этот коэффициент должен попасть на главную диагональ. Может случится так, что значение а_ik уже является наибольшим по модулю элементом. Поэтому проверка выполняется сразу, и в таком случае перестановка не проводится. Фактически процесс состоит в перестановке коэффициентов, один из которых берется в уравнении k, а второй -в уравнении l, каким бы ни было значения l. Перестановку можно провести с помощью трехступенчатого процесса. Эту перестановку необходимо произвести для этих паров коэффициентов двух уравнений:

DO j=k TO N; TEMP=A(k,j); A(k,j)=A(l,j); A(l,j)=TEMP; END;

Аналогично для свободных членов b. Перед возвращением к основному процессу вычислений необходимо восстановить значение индекса i, которое было сначала и использовалось при сравнении. Так как перед началом процесса индекс i имел значения k+1, то в конце блока ему присваивается предшествующее значение i=k+1.

4.2.5 Метод простой итерации

При большом числе неизвестных линейной системы схема метода Гаусса, которая дает верное решение, становится очень громоздкой. В таком случае для поиска корней системы иногда лучше пользоваться приближенными методами вычислений - итерационными. Итерационные методы притягивают своей самонастраиваемостью и простотой реализации на ЭВМ. Пусть дана система линейных алгебраических уравнений (4.7) с неособенной матрицей. Чтобы применить к ней метод итераций, ее необходимо привести к виду:

X=CX+f,

(4.8)

где С - некоторая матрица, а f - вектор-столбец.

Для системы уравнений вида:

АX=В

положим, что диагональные элементы матрица А а_ii0.

В системе выразим Х₁ через 1-е уравнение, Х₂ - через 2-е и т.д.:

Обозначим b_i/a_ii = f, -a_ij/a_ii = C_ij, тогда

X₁=f+C₁₂X₂+C₁₃X₃ +C₁₄X₄ +…+C_1nX_n. или:

X=CX+f,

то есть уравнение (4.8).

Итерационный процесс для начала вычислений требует задание начальных условий. Итерационный процесс строится от начального вектора

.

X(k+1) =СХ(k)+f (k=0,1,2,....)

(4.9)

Сделав итерацию, получим последовательность векторов x⁽¹⁾, x⁽²⁾ ,x⁽³⁾ ,... x^(k). Если последовательность приближений имеет границу, то эта граница есть решением системы (4.7):

(4.10)

Итерационный процесс заканчивается, если |x^(k+1)-x^k|<, где  - заданная точность решения. Процесс итераций (4.9) хорошо сходится (то есть количество приближений, необходимых для получения корней системы (4.7) с заданной точностью, небольшое), если модули диагональных коэффициентов системы (4.7) больше суммы модулей оставшихся коэффициентов уравнения (свободные члены не учитываются). Это условие является условием сходимости. Отметим одну из важных особенностей метода итераций. Сходящийся процесс итераций владеет свойствами самоисправления, то есть отдельная ошибка в вычислениях не отражается на окончательном результате, так как ошибочное приближение можно рассматривать как новое начальное приближение. Начальный вектор x⁰ может быть выбран свободно. Иногда берут x⁰=f, чаще в качестве компонент вектора x⁰ берут приближенные значения неизвестных, полученные грубой прикидкой. Чем более близко начальное приближение к значениям корней системы, тем скорее сойдется итерационный процесс. Если матрица А неособенна, то систему (4.7) с помощью совокупности элементарных преобразований всегда можно привести к виду (4.8).

К элементарным преобразованиям матриц относятся:

а) замена строк (столбцов);

б) умножение всех элементов какой-нибудь строки (столбца) на одну и ту же величину, отличную от нуля;

в) прибавление к элементам какой-нибудь строки (столбца), соответствующих элементов второй строки, столбца, умноженных на одну и ту же величину.

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

Пример

Для нахождения коэффициентов эмпирической зависимости необходимо решить следующую систему уравнений:

2*x₁-0.5*x₂+4*x₃=22 (A)

3*x₁+3.5*x₂+3*x₃=35 (Б)

x₁+0.5*x₂-x₃=-1 (B)

В уравнении (A) коэффициент при х₃ по модулю больший суммы модулей остальных коэффициентов, поэтому можно принять это уравнение как третье уравнение новой системы. Чтобы получить второе уравнение с максимальным по модулю коэффициентом при х₂, достаточно составить разность (Б)-(A):

х₁+4*x₂-x₃=13.

При преобразовании начальной системы к виду (4.8) были использованы уравнения (A) и (Б). Тогда, в первое уравнение обязательно должно войти уравнение (В), например, используя следующее преобразование (А)+4*(В):

6*x₁+1.5*x₂+0*x₃=18.

Таким образом, получена эквивалентная система уравнений, которая удовлетворяет условиям сходимости:

6*x₁+1.5*x₂=18 (А1)

x₁+4*x2-x₃=13 (Б1)

2*x₁-0.5*x₂+4*x₃=22 (В1)

Решая первое уравнение относительно х₁, второе - относительно х₂, третье - относительно х₃, получим:

x₁=3-0.25*x₂

x₂=3.25-0.25*x₁+0.25*x₃ (4.11)

x₃=5.5-0.5*x₁+0.125*x₂

В качестве начального вектора x⁰ возьмем элементы столбца свободных членов, округляя их значения до одного знака после запятой: x₁⁰=3; x₂⁰=3.3; x₃⁰=5.5. Точность вычисления =0.01. Последовательно вычисляем:

при k=1

x₁⁽¹⁾=3-0.25*3.3=2.175

x₂⁽¹⁾=3.25-0.25*3+0.25*5.5=3.875

x₃⁽¹⁾=5.5-0.5*3+0.125*3.3=4.413

при k=2

x₁⁽²⁾=3-0.25*3.875=2.031

x₂⁽²⁾=3.25-0.25*2.175+0.25*4.413=3.810

x₃⁽²⁾=5.5-0.5*2.175+0.125*3.875=4.897 и т.д.

Результаты дальнейших расчетов сведем в таблицу.

K	0	1	2	3	4	5	6
X1 X2 X3	3 3,3 5,5	2,175 2,875 4,413	2,031 3,875 4,413	2,031 3,810 4,961	2,009 3,978 4,961	2,006 3,991 4,993	2,002 3,997 4,996

Так как модули различий значения x^(k) при k=5 и k=6 составляют |x₁⁶-x₁⁵|=0.004, |x₂⁶-x₂⁵|=0.006, |x₃⁶-x₃⁵|=0.003, что меньше заданной точности , то в качестве решения принимаем: х₁=2.002 ; х₂=3.997 ; х₃=4.9967.

Для сравнения точные значения неизвестных: х₁=2; х₂2=4; х₃=5.