====== Решение систем линейных алгебраических уравнений ====== Данная статья посвящена //нахождению корней систем линейных алгебраических уравнений//. Методы численного решения таких систем подразделяются на два типа: //прямые (конечные)// и //итерационные (бесконечные)//. Оба типа полезны и удобны для практических вычислений и каждый из них имеет свои преимущества и недостатки. ==== Метод исключений (метод Гаусса) ==== Данный метод является наиболее известным и широко применяемым. Рассмотрим систему из **n **уравнений с **n **неизвестными. Обозначим неизвестные через //x[1], x[2], ... , x[n]// и запишем систему в следующем виде: a[1,1]*x[1] + a[1,2]*x[2] + ... + a[1,n]*x[n] = b[1] a[2,1]*x[1] + a[2,2]*x[2] + ... + a[2,n]*x[n] = b[2] .................................................... a[n,1]*x[1] + a[n,2]*x[2] + ... + a[n,n]*x[n] = b[n] Предполагается, что в силу расположения уравнений, //a[i,i] <> 0//. Если это не так, то, меняя уравнения местами, добиваемся выполнения этого условия. Введем **n-1** множителей: m[i] = a[i,1]/a[1,1] , где i = 2, 3, ... , n и вычтем из каждого **i**-го уравнения первое, умноженное на **m**, обозначая a'[i,j] = a[i,j] - m[i]*a[1,j] b'[i] = - m[i]*b[i] i = 2, 3, ... , n; j = 1, 2, ... , n. Для всех уравнений, начиная со второго, получим: a[i,1] = 0 , где i = 2, 3, ... , n. Преобразованная система уравнений запишется в следующем виде: a[1,1]*x[1] + a[1,2]*x[2] + ... + a[1,n]*x[n] = b[1] 0 + a'[2,2]*x[2] + ... + a'[2,n]*x[n] = b'[2] .................................................... 0 + a'[n,2]*x[2] + ... + a'[n,n]*x[n] = b'[n] Продолжая таким же образом, исключаем **x[2]** из последних **n-2** уравнений, затем **x[3]** из последних **n-3** и т.д. На некотором **k**-м этапе мы исключаем **x[k]** с помощью множителей: m^(k-1)_[i] = a^(k-1)_[i,k]/a^(k-1)_[k,k] , i = k+1, ... , n, , незабывая, что a(k-1)[k,k] <> 0. == Примечание == **^(k-1)** обозначает на **k-1**-й итерации, **_[i,k]** - индекс эл-та. Тогда a^(k)_[i,j] = a^(k-1)_[i,j] - m^(k-1)_[i]*a^(k-1)_[k,j] b^(k)_[i] = - m^(k-1)_[i]*b^(k-1)_[k] i = k+1, k+2, ... , n; j = k, k+1, ... , n. Окончательно, треугольная система уравнений записывается: a[1,1]*x[1] + a[1,2]*x[2] + ... + a[1,n]*x[n] = b[1] a[2,2]*x[2] + ... + a[2,n]*x[n] = b[2] .................................................... a[n,n]*x[n] = b[n] Обратная подстановка для нахождения значений неизвестных задается формулами: x[n] = b^(n-1)_[n] / a^(n-1)_[n,n] x[n-1] = (b^(n-2)_[n-1] - a^(n-2)_[n-1,n]*x[n]) / a^(n-2)_[n-1,n-1] ................................................................... x[j] = ( b^(j-1)_[j] - a^(j-1)_[j,n]*x[n] - ... - - a^(j-1)_[j,j+1]*x[j+1] ) / a^(j-1)_[j,j] j = n-2, n-3, ... , 1. Перед началом процесса исключения очередного неизвестного может понадобиться переставить уравнения в системе, чтобы __избежать деления на нуль__. Кроме этого, при перестановке уравнений необходимо добиваться того, чтобы | a^(k-1)_[k,k] | >= | a^(k-1)_[i,k] | Поэтому перед началом процесса исключения очередного неизвестного необходимо переставить уравнения так, чтобы максимальный по модулю коэффициент при **x[k]** попал на главную диагональ (данный способ решения часто называется //методом главного элемента//). ==== Итерационный метод Гаусса-Зейделя ==== Данный метод отличается простотой и легкостью программирования, малой ошибкой округления, но сходится при выполнении некоторых условий.\\ Рассмотрим все ту же систему из **n **уравнений с **n **неизвестными.\\ По прежнему полагаем, что **a[i,i] <>0** для всех **i**. Тогда **k**-е приближение к решению будет задаваться формулой: x^(k)_[i] = ( b[i] - a[i,1]*x^(k)_[1] - ... - a[i,j-1]*x^(k)_[i-1] - a[i,j+1]*x^(k-1)_[i+1] - a[i,n]*x^(k-1)_[n] ) / a[i,i] i = 1, 2, ... , n. Как известно, любой итерационный процесс, как и настоящий мужик, должен имееть свой конец :) Поэтому вычисления нужно продолжать до тех пор, все x^(k)_[i] не станут достаточно близки к x(k-1)[i], т.е. пока для всех **i** не будет выполнятся неравенство max| x^(k)_[i] - x^(k-1)_[i] | < Eps , где **Eps** - некоторое положительное число, задающее точность вычислений. Итерационный метод //**Гаусса-Зейделя**// сходится, если __удовлетворяется достаточное условие сходимости__: //для всех **i** модуль диагонального коэффициента должен быть не меньше суммы модулей остальных коэффициентов данной строки// | a[i,i] | >= | a[i,1] | + | a[i,2] | + ... + | a[i,j-1] | + | a[i,j+1] | + ... + | a[i,n] | , а хотябы для обной строки это неравенство должно выполнятся строго | a[i,i] | > | a[i,1] | + | a[i,2] | + ... + | a[i,j-1] | + | a[i,j+1] | + ... + | a[i,n] |. ==== Сравнение методов ==== Метод исключений имеет то преимущество, что он конечен и теоретически, с его помощью можно решить любую невырожденную систему. Итерационный метод сходится только для специальных систем уравнений, однако когда итерационные методы сходятся, они, обычно, предпочтительнее: * Время вычислений пропорционально n2 на итерацию, в то время как для метода исключений время пропорционально n3. Если для решения системы требуется меньше чем n итераций, то общие затраты машинного времени будут меньше. * Как правило, ощибки округления в итерационном методе меньше. Часто, это соображение может оказаться достаточно важным, что оправдывает дополнительные затраты машинного времени. Многие систему, возникающие на практике, имеют среди коэффициентов большое кол-во нулей. В этом случае итерационные методы, если они сходятся, в высшей степени предпочтительны т.к. в методе исключений получается треугольная система, которая в качестве коэффициентов нулей обычно не содержит. А вот и исходник по теме: {$N+,G+,X+} const Comments: boolean = false; { нужно ли печатать пошаговые рез-ты расчета } Eps = 0.00001; { Точность вычислений } n_max = 4; { Кол-во ур-й и неизвестных } const_AnB: array[1..n_max,1..n_max+1] of double = (( -8, 2, 17, -5, -119.97), { -8*X[1] + 2*X[2] + 17*X[3] - 5*X[4] = -119.97 и т.д. ... } ( 4, -22, 6, 5, 55.73), ( 15, 3, -5, -5, 18.77), ( -4, -4, 5, 14, -79.42)); { Описание типов матрицы уравнеиний и массива найденных Х } Type TMatr = array[1..n_max,1..n_max+1] of double; Type TXMarr = array[1..n_max] of double; procedure PrintX(const X: TXMarr; n: byte; cnt: word); { Печатает найденные Х } var k, { Позиция заданного корня в массиве } i: byte; { Номер текущего корня для вывода } begin if n = 0 then exit; k:=low(X); if cnt = 0 then begin write(' N x1'); for i:=2 to n do write(' x',i); writeln; end; write(cnt:2); for i:=1 to n do begin Write(' ',X[k]:3:5); inc(k); end; writeln; { При решении систем, с большим кол-вом неизвестных, возможно, прийдется переделать эту процедуру для корректного вывода всех решений. } end; { PrintX } procedure PrintMatr(var AnB: TMatr; n: byte); { Печатает заданную матрицу } var i,j: byte; begin for i:=1 to n do begin for j:=1 to n+1 do write(AnB[i,j]:3:5,' '); writeln; end; writeln; { При решении систем, с большим кол-вом неизвестных, возможно, прийдется переделать эту процедуру для корректного вывода всех решений. } end; { PrintMatr } function Normalize(var AnB: TMatr; n: byte): byte; { Располагает на главной диагонали наибольшие элементы в столбцах Если один из них равен нулю,то систему решить не получится ... Возвращает 0 если систему заданными методами решить не получится, 1 - с-му можно решить только методом Гаусса 2 - с-му можно решать любым методом } var max: double; { Переменная для поиска макс. эл-та } imax: byte; { Переменная для поиска макс. эл-та } tmp: double; { Переменная для обмена элементов } i,j: byte; { Счетчики циклов } begin Normalize:=2; for j:=1 to n do begin max:=Abs(AnB[1,j]); imax:=1; for i:=2 to n do if Abs(AnB[i,j]) > max then begin max:=Abs(AnB[i,j]); imax:=i; end; if max < Eps then begin Normalize:=0; Writeln('Error: a[i,i]=0!'); exit; end else if j <> imax then for i:=1 to n+1 do begin tmp:=AnB[j,i]; AnB[j,i]:=AnB[imax,i]; AnB[imax,i]:=tmp; end; end; for i:=1 to n do begin tmp:=0; for j:=1 to n do if i <> j then tmp:=tmp+Abs(AnB[i,j]); if Abs(AnB[i,i]) < tmp then Normalize:=1; end; end; { Normalize } procedure Metod1(const Matr; n: byte); var k : byte; M_ : TMatr absolute Matr; AnB : TMatr; X,M : TXMarr; function Calc(k: byte): boolean; var i,j : byte; begin Calc:= true; for i:=k+1 to n do begin M[i]:= AnB[i,k]/AnB[k,k]; if M[i] > 1 then begin Calc:= false; exit; end; end; for i:=k+1 to n do for j:=k to n+1 do AnB[i,j]:=AnB[i,j]-M[i]*AnB[k,j]; if Comments then PrintMatr(AnB,n); end; { Calc } function Summ(j: byte): double; var i : byte; res : double; begin res:=AnB[j,n+1]; for i:=n downto j+1 do res:=res - AnB[j,i]*X[i]; Summ:=res; end; { Summ } begin Writeln('Метод Гаусса:'); AnB:=M_; if Normalize(AnB,n) = 0 then begin writeln('Error: Невозможно решить систему уравнений!'); exit; end; for k:=1 to n-1 do if not Calc(k) then begin Writeln('Error: Невыполняется условие | a^(k-1)_[k,k] | >= | a^(k-1)_[i,k] | !'); exit; end; X[n]:=AnB[n,n+1]/AnB[n,n]; for k:=n-1 downto 1 do X[k]:=Summ(k)/AnB[k,k]; PrintX(X,n,0); Writeln; end; { Metod1 } procedure Metod2(const Matr; n: byte); var M_ : TMatr absolute Matr; AnB : TMatr; X : TXMarr; i,j: byte; tmp: double; delta: double; cnt: word; function Summ(i: byte): double; var j: byte; res: double; begin res:=AnB[i,n+1]; for j:=1 to n do if i <> j then res:=res-AnB[i,j]*X[j]; Summ:=res; end; { Summ } begin writeln('Метод Гаусса-Зейделя:'); AnB:=M_; For i:=1 to n do X[i]:=0; if Normalize(AnB,n) <> 2 then begin writeln('Error: !'); exit; end; delta:=1; cnt:=0; while delta > Eps do begin if Comments then PrintX(X,n,cnt); delta:=0; for i:=1 to n do begin tmp:=Summ(i)/AnB[i,i]; if delta < Abs(tmp-X[i]) then delta:=Abs(tmp-X[i]); X[i]:=tmp; end; inc(cnt); end; PrintX(X,n,cnt); writeln('Заданная точность вычислений достигнута на ',cnt,' шаге.'); Writeln; end; { Metod2 } begin Writeln(#10,'Численное решение систем линейных алгебраических уравнений:',#10); Metod1(const_AnB,n_max); Writeln('Press Enter to continue'); readln; Metod2(const_AnB,n_max); Writeln('Press Enter if you wanna DOS ;)'); readln; end. Автор: Labinskiy Nikolay aka e-moe (c) 2005