2.3 Итерационные методы решения слау
С помощью итерационных
методов строится последовательность
векторов {x(k)}
такая, что
,
где x
– точное решение
системы. Эффективность таких методов
определяется скоростью сходимости
последовательных приближений x(k)
к x,
причем немаловажно, чтобы вычисления
на каждой итерации были не слишком
затратными.
Итерационные способы построения решения систем линейных алгебраических уравнений являются серьёзной альтернативой прямым методам, по крайней мере, в случаях, когда размерность этих систем велика. Особенно широкого применения методы итераций получили для систем уравнений, к которым приводят численные методы решения дифференциальных уравнений в частных производных. Матрицы таких систем имеют большое число нулевых элементов, и в отличие от прямых методов итерационные не увеличивают число ненулевых элементов матрицы в процессе вычислений.
Отметим еще одну важную особенность итерационных методов: как правило, они менее чувствительны к погрешностям округления, возникающим в реальных компьютерных расчетах из-за конечности разрядной сетки. Эти методы даже используют для уточнения решения, полученного прямыми методами. С другой стороны, если матрица системы плохо обусловлена, их скорость сходимости может быть очень медленной.
2.2.2.1 Методы Якоби и Гаусса–Зейделя
Довольно широкий класс итерационных методов в своей основе использует следующую идею: пусть ищется решение невырожденной системы
Ax=b, (2.9)
которая может быть преобразована к эквивалентной системе
x=Bx+c (2.10)
с тем же самым вектором неизвестных x и некоторыми новыми матрицей B и вектором c. Далее строится итерационный процесс по правилу
x(k+1)=Bx(k)+c, k = 0, 1, 2,…, (2.11)
где
–
заданный вектор неизвестных – начальное
приближение. В качестве условия завершения
счёта обычно используют малость величин
||x(k+1)–x(k)||
и/или ||Аx(k)–b||.
Вычисления приближений по формуле
(2.11) иногда называют обобщённым методом
простой итерации.
Конкретные методы отличаются способом приведения исходной системы к виду (2.10), а также условиями и параметрами сходимости итерационного процесса (2.11). К примеру, представим матрицу системы в виде
A=L+D+R
где A – диагональная, L и R – соответственно левая и правая строго треугольные (т.е. с нулевой диагональю) матрицы. Тогда система (2.9) может быть записана в виде
Lx+Dx+Rx=b
или, если на главной диагонали исходной матрицы нет нулей,
x = –D–1(L+R)x + D–1b,
т.е. в (2.10) и (2.11) следует положить
B = –D–1(L+R); c = D–1b.
Построенный таким образом итерационный процесс
x(k+1)= –D–1(L+R)x(k)+D–1b, k = 0, 1, 2,…,
называется методом Якоби. Его легко записать в развернутом виде, если учесть, что D–1, как и D, – матрица диагональная, но с элементами dii = 1/aii:
(2.12)
Под методом Гаусса–Зейделя (или просто Зейделя) понимают такую модификацию простых итераций (2.11), когда найденные при вычислении (k+1)-го приближения компоненты вектора x(k+1) сразу же используются в вычислениях. В координатной форме итерационный процесс Зейделя имеет вид
Рассматривая данную модификацию применительно к итерациям метода Якоби (2.12), получим
(2.13)
где k = 0, 1, 2,…, а задается.
Заметим, что итерационный процесс (2.13) укладывается в общее представление (2.11), при этом
B = – (L+D)–1R; c = (L+D)–1b.
Точные формулировки условий сходимости методов Якоби и Зейделя можно найти в литературе по численным методам (см., например, [2, 5]). Отметим здесь лишь, что для определенного класса матриц эти методы непременно сходятся, а именно для матриц A с диагональным преобладанием, т.е. таких, для которых справедливо
для всех i
= 1, 2, …, n.
Кроме того, метод Зейделя сходится для систем с симметричными положительно определенными матрицами, а это значит, что им можно решать системы общего вида (2.12), но применяя к видоизмененной форме ATAx=ATb (сходимость при этом может сильно замедляться).
Ниже приводится текст программы по реализации метода Зейделя типа (2.13). Для счёта используется система
> restart;
> with(linalg);
> A:=matrix([[4.,2,-1,0],[2,4,2,1],[-1,2,8,4],
[0,1,4,10]]);
> b:=vector([9.,11,7,28]); # 1)
> x:=vector(4); x0:=vector(4); x:=[0.,0.,0.,0.];# 2)
> eps:=1.0e-5: n:=4: k:=0:
> r:=evalm(1.*A&*x-b);
> while sqrt(dotprod(r,r))>eps do
for i from 1 to n do
s:=0.; x0[i]:=x[i];
for j from 1 to n do
if i<>j then s:=(s+A[i,j]*x[j]);fi; od;
x[i]:=(b[i]-s)/A[i,i];
r[i]:=x[i]-x0[i];
od;
k:=k+1;
if k>1000 then
print(`метод не сходитсЯ`); `stop`(1); fi;
od:
> evalm(x); k;
Указание. Если в процессе счета наблюдается зацикливание (т.е. метод не сходится), рекомендуется провести аналогичный расчет для эквивалентной системы Wx=g, где W=ATA, g=ATb. Для этого между строками 1 и 2 в данной программе нужно вставить
> W:= evalm(transpose(A)&*A);
> g:= evalm(transpose(A)&*b);
а в остальном тексте программы вместо матрицы A и вектора b вставить соответственно W и g.
2.2.2.2 Метод сопряженных градиентов
Еще один класс методов итерационного решения СЛАУ – это так называемые методы вариационного типа. К ним относятся методы минимальных невязок, минимальных ошибок, наискорейшего спуска, сопряженных градиентов и др. Их главная особенность в том, что непосредственное решение системы подменяется решением эквивалентной экстремальной задачи.
Например, один из наиболее популярных и хорошо разработанных методов подобного типа – метод сопряженных градиентов (МСГ) – минимизирует квадратичный функционал
.
Можно показать, что если матрица A симметричная и положительно определенная, то вектор x, для которого F(x) достигает минимума, одновременно является решением системы Ax=b.
Приведем без вывода алгоритм МСГ.
1. Задать x(0) (начальный вектор) и число (уровень допустимых погрешностей).
2. Вычислить вектор r(0) = Ax(0)–b (невязка начального приближения – разность между левой и правой частью системы).
3. Положить v(0) = r(0), k = 0 (номер итерации).
4. Вычислить число
.
5. Вычислить вектор x(k+1) = x(k)+kv(k) (новое приближение).
6. Вычислить вектор r(k+1) = Ax(k+1) – b (невязка (k+1)-го приближения).
7. Вычислить число
.
8. Вычислить вектор v(k+1) = r(k+1) + kv(k).
9. Положить kk+1.
10. Проверить условие ||r(k+1)||<; если “да”, остановить счёт и вывести результаты, если “нет”, перейти к шагу 4.
Интересно отметить, что МСГ занимает особое положение среди методов решения СЛАУ, ибо доказано, что векторы невязок {r(k), k = 0, 1, 2, …} образуют ортогональную систему векторов, следовательно, нулевая невязка r(k)=0 (что равносильно точному решению) должна появиться не позднее n-й итерации, где n – порядок системы уравнений. Таким образом, МСГ должен быть отнесен к прямым методам. Те не менее на практике он применяется именно как итерационный метод, поскольку реальный вычислительный процесс не идеален из-за неизбежных ошибок округления, и на n-м шаге может быть не достигнута нужная точность. Кроме того, для систем большой размерности для получения решения с приемлемой точностью может понадобиться значительно меньшее, чем n, число итераций. Поэтому в силу указанных обстоятельств к МСГ иногда употребляют термин «полуитерационный метод».
Запрограммируем алгоритм МСГ в системе Maple (на примере той же системы, что была решена методом Зейделя).
> restart;
> with(linalg);
> A:=matrix([[4.,2,-1,0],[2,4,2,1],[-1,2,8,4], [0,1,4,10]]);
> b:=vector([9.,11,7,28]); r:=vector(4):
> x:=eval(b): eps:=0.0001: k:=0:
> r:=evalm(A&*x-b); v:=eval(r); t1:=dotprod(r,r);
> while sqrt(t1)>eps do
alpha:=-evalf(dotprod(r,v)/dotprod(A&*v,v));
x:=evalm(x+alpha*v): print(eval(x));
t2:=t1; r:=evalm(A&*x-b); t1:=dotprod(r,r);
beta:=t1/t2;
v:=evalm(r+beta*v);
k:=k+1;
if k>1000 then print(`метод не сходитсЯ`);
`stop`(1); fi;
od:
> eval(x); k;
4
Замечание. Для успешной работы программы необходимо отслеживать, чтобы все вычисления проводились в числовом формате (а не в символьном). Прямым указанием на такой режим счёта, как обычно, является явная запись чисел с точкой. Обратите внимание, что здесь так записаны по одному элементу матрицы A и вектора b. В качестве дополнительной меры можно рекомендовать видоизменить команды вычисления r:
r:=evalm(1.0*A&*x-b);
а также все вычисления производить посредством evalf.
Если матрица системы не симметричная и не положительно определенная, сходимость МСГ не гарантируется. В этом случае можно рекомендовать применить алгоритм метода сопряженных градиентов к эквивалентной системе ATAx=ATb, у которой матрица обладает необходимыми свойствами. Скорость сходимости при этом может быть низкой.