Данная статья посвящена нахождению корней систем линейных алгебраических уравнений. Методы численного решения таких систем подразделяются на два типа: прямые (конечные) и итерационные (бесконечные). Оба метода полезны и удобны для практических вычислений, и каждый из них имеет свои преимущества и недостатки.
Данный метод является наиболее известным и широко применяемым.
Рассмотрим систему из 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" |
| _[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] |.
Метод исключений имеет то преимущество, что он конечен и, теоретически, с его помощью можно решить любую невырожденную систему. Итерационный метод сходится только для специальных систем уравнений, однако когда итерационные методы сходятся, они, как правило, предпочтительнее:
А вот и исходник по теме:
{
Writeln('Labinskiy Nikolay aka e-moe (c) 2005');
Writeln('See you later on forum.sources.ru ;)');
}
{$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.