Решение систем линейных уравнений

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

а₁₁ x₁ + a₁₂ х₂ +... +а_1n х_n =с₁

а₂₁ x₁ + a₂₂ x₂ +... +а_2n х_n =с₂

. . . . . . . . .

а_n1 х₁ + а_n2 х₂ +... +а_nn x_n =с_n.

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

МЕТОД ИСКЛЮЧЕНИЯ ГАУССА

Наиболее распространенные прямые методы основаны на приведении системы уравнений к "треугольному" виду. При этом одно из уравнений содержит только одну неизвестную, а в каждом следующем добавляется еще по одной неизвестной. При счете вручную приведение к треугольному виду достигается сложением и вычитанием уравнений после умножения их на соответствующие постоянные множители. Выполняя эту процедуру вручную, нетрудно ошибиться, однако она позволяет построить удобный алгоритм численного решения на ЭВМ. Одним из используемых для этого методов является метод Гаусса. Применяя его, сначала нормируют первое уравнение, деля его коэффициенты на а , Затем первое уравнение умножают на первые коэффициенты, а всех других уравнений и последовательно вычитают из остальных уравнений. В результате первая переменная будет исключена из всех уравнений, кроме первого. На следующем этапе решения такая же процедура применяется к остальным n - 1 уравнениям. В результате из оставшихся n - 2 уравнений исключается вторая неизвестная. Всю процедуру повторяют до тех пор, пока после n шагов вся система не будет приведена к треугольному виду. Математически эту процедуру можно описать следующим образом. На k-м шаге процесса исключения новые нормированные коэффициенты k-го уравнения имеют вид

,

а новые коэффициенты в следующих уравнениях будут иметь вид

b_i,j =a_ij -a_ik b_kj>k.

Выполняя эту процедуру, следует помнить, что коэффициенты нижестоящих уравнений а меняются на каждом шаге. Например, коэффициенты b­_ij на следующем шаге превращаются в коэффициенты а_ij . Проиллюстрируем применение описанной процедуры на следующем простом примере. Пусть методом исключения Гаусса требуется решить систему уравнений

Х₁ + Х₂ + Х₃ - Х₄ = 2

Х₁ - Х₂ - Х₃ + Х₄ = 0

2X₁ + X₂ - X₃ +2X₄ = 9

ЗX₁ + X₂ +2Х₃ - X₄ = 7

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

А₁ 1 1 1 -1 2

А₂ 1 -1 -1 1 0

А₃ 2 1 -1 2 9

А₄ 3 1 2 -1 7

Исключая члены, содержащие х₁ , получим

В₁ =А₁ /1 1 1 1 -1 2

В₂ =А₂ -В₁ 0 -2 -2 2 -2

В₃ =А₃ -2B₁ 0 -1 -3 4 5

В₄ =А₄ -ЗВ₁ 0 -2 -1 2 1

После исключения членов с х₂ имеем

В₁ 1 1 1 -1 2

С₂ =В₂ /(-2) 0 1 1 -1 1

С₃ =В₃ +С₂ 0 0 -2 3 6

С₄ =В₄ +2С₂ 0 0 1 0 3

Исключение членов с х₃ дает

В₁ 1 1 1 -1 2

С₂ 0 1 1 -1 1

D₃ =С₃ /(-2) 0 0 1 -3/2 -3

D₄ =C₄ -D₃ 0 0 0 3/2 6

Дойдя до последнего ряда, получим

В₁ 1 1 1 -1 2

С₂ 0 1 1 -1 1

В₃ 0 0 1 -3/2 -3

Е₄ =2D₄ /3 0 0 0 1 4

Возвращаясь к форме уравнений, получим

х₁ + х₂ + х₃ - х₃ = 2,

x₂ + x₃ - x₄ = 1,

x₃ - (3/2) x₄ =-3,

x₄ = 4,

откуда, подставляя значение х в 3-е уравнение, х - во 2-е и т.д., находим решение системы уравнений

х₁ =1, х₂ =2, х₃ =З, х₄ =4.

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

МЕТОД ПРОСТОЙ ИТЕРАЦИИ

Пусть система линейных уравнений

Ах=b

каким-либо образом приведена к виду

х= Cx+f,

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

Исходя из произвольного вектора х ,

х⁽⁰⁾ =

строим итерационный процесс

х^(k+1) = Сх^(k) +f (k=0,1,2,...),

или в развернутой форме

x₁^(k+1) =c₁₁ x₁^(k) + c₁₂ x₂^(k) +...+ c₁₃ x₃^(k) + f₁

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

x_n^(k+1) =c_1n x₁^(k) + c_1n x₂^(k) +...+ c_1n x₃^(k) + f_n

Производя итерации, получим последовательность векторов х⁽¹⁾ ,

x⁽²⁾ ,..., x^(k) ,...

Доказано , что если элементы матрицы С удовлетворяют одному из условий

(i=1,2,...,n)

или

(i=1,2,...,n)

то процесс итерации сходится к точному решению системы х при любом начальном векторе х , т. е.

х = .

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

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

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

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

0 (i=1,2,...,n),

то систему (9.1) можно записать в виде

х₁ =(b₁ - a₁₂ x₂ -...- a_1n x_n ),

х₂ =(b₂ - a₂₁ x₁ -...- a_2n x_n ),

. . . . . .

х_n =(b₁ - a₁₂ x₂ -...- a_1n x_n ),

В этом случае элементы матрицы С определяются следующим образом:

(ij), c_ii =0,

и тогда условия (9.4) и (9.5) соответственно приобретают вид

(i=1,2, ...,n),

(i=1,2, ...,n),

Неравенства (9.7), (9.8) будут выполнены, если диагональные элементы матрицы А удовлетворяют условию

|а_ii |> (i=1,2, ..., n),

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

В т о р о й с п о с о б (на примере оценки погрешности)

Решить систему

1,02 х₁ - 0,05 х₂ - 0,10 х₃ = 0,795,

-0,11 х₁ + 1,03 х₂ - 0,05 х₃ =0,849,

-0,11 х₁ - 0,12 х₂ + 1,04 х₃ =1,398,

произведя три итерации. Указать погрешность полученного результата.

Р е ш е н и е. Матрица данной системы такова, что диагональные элементы близки к единице, а все остальные - значительно меньше единицы. Поэтому для применения метода итераций естественно записать систему (9.11) в виде

х₁ = 0,795 - 0,02 х₁ + 0,05 х₂ + 0,10 х₃

х₂ = 0,849 + 0,11 х₁ - 0,03 х₂ + 0,05 х₃

х₃ = 1,398 + 0,11 х₁ + 0,12 х₂ - 0,04 х₃

условия сходимости (9.4) для полученной системы выполнены. Действительно,

=0,02+0,05+0,10=0,17 < 1,

=0,11+0,03+0,05=0,19 < 1,

=0,11+0,12+0,04=0,27 < 1

Берем в качестве начального вектора х⁽⁰⁾ столбец свободных членов, округлев его элементы до двух знаков после запятой:

х⁽⁰⁾ =

Далее последовательно находим:

при k=1

х₁⁽¹⁾ =0,795-0,016+0,0425+0,140=0,96150,962,

х₂⁽¹⁾ =0,849+0,088-0,0255+0,070=0,98150,982,

x₃⁽¹⁾ =1,398+0,088+0,1020-0,056=1,532;

при k=2

х₁⁽²⁾ =0,978060,978, х₁⁽²⁾ =1,001961,002,

х₃⁽²⁾ =1,560381,560;

при k=3

х₁⁽³⁾ =0,980, х₂⁽³⁾ =1,004, х₃⁽³⁾ =1,563.

Значения неизвестных при k = 2 и k = -3 отличаются не более чем на 310^-3 ; поэтому, если в качестве приближенных значений неизвестных взять

х₁ =0,980, х₂ =1,004, x₃ = 1,563,

то погрешность этих приближенных значений не превзойдет

МЕТОД ЗЕЙДЕЛЯ

Метод Зейделя является модификацией метода простой итерации. Он заключается в том, что при вычислении (k+1)-го приближения неизвестного х при i>1 используются уже вычисленные ранее (k+1)-е приближения неизвестных х х, ..., х . Таким образом, для системы вычисления по методу Зейделя ведутся по формулам

x₁^(k+1) =c₁₁x₁^(k) + c₁₂ x₂^(k) +...+ c_1n x_n^(k) +f₁

x₂^(k+1) =c₂₁x₁^(k+1) + c₂₂ x₂^(k) +... +c_{2 n-1} x_n-1^(k)+ c_2n x_n^(k) +f₂

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

x_n^(k+1) =c_n1x₁^(k+1) + c_n2 x₂^(k+1) +... +c_{n n-1} x_n-1^(k+1)+ c_{n n} x_n^(k) +f₂

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

ДОСТОИНСТВА ИТЕРАЦИОННЫХ МЕТОДОВ

Если метод итераций сходится, он дает следующие преимущества по сравнению с методами, рассмотренными выше.

1) Если итерации сходятся достаточно быстро, т. е. Если для решения системы требуется менее n итераций, то получаем выигрыш во времени, так как число арифметических действий, необходимых для одной итерации, пропорционально n а общее число арифметических действий в методе Гаусса, например, пропорционально n .

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

Последнее обстоятельство часто используется для уточнения значений неизвестных, полученных методом Гаусса.

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

4) процесс итераций приводит к выполнению однообразных операций и сравнительно легко программируется на ЭВМ.

ПРИМЕР ПРОГРАММЫ РЕШЕНИЯ СИСТЕМЫ УРАВНЕНИЙ МЕТОД ГАУССА

Метод Гаусса реализован в виде подпрограммы. Вызывается из основной программы, в которой и задается исходное уравнение.

Обратите особое внимание! Использовалась процедура с параметрами переменными (в заголовке процедуры стоят слова VAR). Это нужно для того чтобы не тратить времени и памяти на размещение внутри процедуры копии массива, а также чтобы не переполнить стек при их передаче в подпрограмму.

type ttarr=array [1..3,1..3] of real; (* объявление типов *)

type tarr=array [1..3] of real;

var

n,i,j:integer; (*обявление переменных*)

a1:ttarr;

x1,b1:tarr;

procedure Gauss(var A:ttarr;var X,B:tarr); (*заголовок процедуры*)

var

h:real;

I,J,K:integer;

begin

for I:=1 to N-1 do

begin

for J:=I+1 to N do

begin

A[J,I]:=-A[J,I]/A[I,I]; (*нормированные кооф.*)

for k:=I+1 to N do

A[J,K]:=A[J,K]+A[J,I]*A[I,K]; (*выч. новых кооф.*)

B[J]:=B[J]+A[J,I]*B[I]

end;

end;

X[N]:=B[N]/A[N,N]; (*обратное решение. Найдена одна*)

for i:=n-1 downto 1 do (*переменная*)

begin

h:=b[i]; (*определение остальных

for J:=I+1 to N do h:=h-X[J]*A[I,J]; переменных*)

X[I]:=h/A[i,i]

end;

end;

begin (*начало основной программы*)

n:=3; (*задание размерности и системы уравнений*)

a1[1,1]:=4; A1[1,2]:=0.24; A1[1,3]:=-0.08; b1[1]:=8; x1[1]:=1.8;

a1[2,1]:=0.09; A1[2,2]:=3; A1[2,3]:=-0.15; b1[2]:=9; x1[2]:=3;

a1[3,1]:=0.04; A1[3,2]:=-0.08; A1[3,3]:=4; b1[3]:=20;x1[3]:=5;

Writeln;

for i:=1 to n do (*Вывод системы на экран*)

begin

for j:=1 to n do write(a1[i,j]:2:2,' ');

writeln(b1[i]:2:2,' ',x1[i]:2:2);

end;

writeln;writeln;

Gauss(A1,x1,b1); (*решение системы методом Гаусса*)

for i:=1 to n do

writeln(' ',x1[i]:2:7) (*Вывод решения на экран*)

end.

ПРИМЕР ПРОГРАММЫ РЕШЕНИЯ СИСТЕМЫ УРАВНЕНИЙ МЕТОДОМ ИТТЕРАЦИЙ

Используем метод итераций для решения того же уравнения. Структура программы такая же. Процедура решения уравнения указанным методом называется Iter. При начальных значениях заданных в переменных х1[1], х1[2], х1[3]

type ttarr=array [1..3,1..3] of real;

type tarr=array [1..3] of real;

var

s,n,i,j:integer; (*s-специально объявлена за пределами процедуры*)

e:real;

a1:ttarr;

z,x1,b1:tarr;

procedure Itter(var A:ttarr;var X,B:tarr;e:real); (*Заголовок процедуры*)

var

I,J,K:integer;

begin

s:=0; (*количество итераций*)

repeat

k:=0; (*Сброс флага точности результата*)

for i:=1 to n do

begin

x[i]:=-b[i];

for j:=1 to n do x[i]:=x[i]+a[i,j]*z[j]; (*вычисление итерации*)

if abs(x[i]/a[i,i])>=e then k:=1; (*проверка точности и сброс*)

x[i]:=z[i]-x[i]/a[i,i] (*флага*)

end;

for i:=1 to n do z[i]:=x[i]; (*замена на новое приближение*)

s:=s+1;

until k=0; (*Конец итераций*)

end;

begin (*основная программа*)

n:=3;e:=1e-4;

a1[1,1]:=4; A1[1,2]:=0.24; A1[1,3]:=-0.08; b1[1]:=8; x1[1]:=2;

a1[2,1]:=0.09; A1[2,2]:=3; A1[2,3]:=-0.15; b1[2]:=9; x1[2]:=3;

a1[3,1]:=0.04; A1[3,2]:=-0.08; A1[3,3]:=4; b1[3]:=20;x1[3]:=5;

Writeln;

for i:=1 to n do

begin

for j:=1 to n do write(a1[i,j]:2:2,' ');

writeln(b1[i]:2:2,' ',x1[i]:2:2);

end;

writeln;writeln;

itter(A1,x1,b1,e);

for i:=1 to n do

writeln(' ',x1[i]:2:7) (*вывод решения*)

Writeln('Кол-во итераций=',s:3)

end.

ПРИМЕР ПРОГРАММЫ РЕШЕНИЯ СИСТЕМЫ УРАВНЕНИЙ МЕТОДОМ ЗЕЙДЕЛЯ

Модифицируем предыдущую программу для решения системы методом Зейделя.

type ttarr=array [1..3,1..3] of real;

type tarr=array [1..3] of real;

var

s,n,i,j:integer;

e:real;

a1:ttarr;

z,x1,b1:tarr;

procedure Itter(var A:ttarr;var z,B:tarr;e:real);

var

I,J,K:integer;

x:tarr;

begin

s:=0;

repeat

k:=0;

for i:=1 to n do

begin

x[i]:=-b[i];

for j:=1 to n do x[i]:=x[i]+a[i,j]*z[j];

if abs(x[i]/a[i,i])>=e then k:=1;

x[i]:=z[i]-x[i]/a[i,i];

z[i]:=x[i] (*отличия здесь*)

end;

s:=s+1;

until k=0;

end;

begin

n:=3;e:=1e-4;

a1[1,1]:=4; A1[1,2]:=0.24; A1[1,3]:=-0.08; b1[1]:=8; x1[1]:=2;

a1[2,1]:=0.09; A1[2,2]:=3; A1[2,3]:=-0.15; b1[2]:=9; x1[2]:=3;

a1[3,1]:=0.04; A1[3,2]:=-0.08; A1[3,3]:=4; b1[3]:=20;x1[3]:=5;

Writeln;

for i:=1 to n do

begin

for j:=1 to n do write(a1[i,j]:2:2,' ');

writeln(b1[i]:2:2,' ',x1[i]:2:2);

end;

writeln;writeln;

itter(A1,x1,b1,e);

for i:=1 to n do

writeln(' ',x1[i]:2:7);

writeln('КОЛ-ВО ИТЕРАЦИЙ=',s:3)

end.