Двунаправленный метод сопряженных градиентов
Решение СЛУ с разреженной матрицей возможно также известным
двунаправленным методом сопряженных градиентов.
Он реализован указанной ниже функцией.
biсд(А. В) — возвращает решение X СЛУ А*Х=В. Матрица коэффициентов А должна быть квадратной размера
пхп,
а вектор-столбец правых частей уравнений В должен иметь длину
п.
Функция bicg начинает итерации от начальной оценки, по умолчанию представляющей собой вектор размером
п,
состоящий из нулей. Итерации производятся или до сходимости к решению, или до появления ошибки, или до достижения максимального числа итераций (по умолчанию равно min(20,n) — либо 20, либо числу уравнений). Сходимость достигается, когда относительный остаток norm(B-A*x)/norm(B) меньше или равен погрешности метода (по умолчанию le-6). Благодаря использованию двунаправленного метода сопряженных градиентов bicg сходится за меньшее число итераций, чем lsqr (в нашем примере быстрее на одну итерацию), но требует квадратную матрицу А, отбрасывая информацию, содержащуюся в дополнительных уравнениях, в то время как 1 sqr работает и с прямоугольной матрицей;
bicgCA.B.tol) — выполняет и возвращает решение с погрешностью (порогом отбора) tol;
bicg(A,b.tol, max it) — выполняет и возвращает решение при заданном максимальном числе итераций maxit;
bicgCA.b.tol,maxit,М) и bicg(A,b.tol .maxit,Ml,М2) — при решении используются матрица предусловий М или М=М1*М2, так что производится решение системы inv(M)*A*x=inv(M)*b относительно х. Если Ml или М2 — пустые матрицы, то они рассматривается как единичные матрицы, что эквивалентно отсутствию входных условий вообще;
bicgCA.B.tol, maxit. Ml. M2.X0) — точно задается начальное приближение Х0. Если Х0 — пустая матрица, то по умолчанию используется вектор, состоящий из нулей;
X = bi eg (А, В, tol, maxit. Ml, M2.X0) — при наличии единственного выходного параметра возвращает решение X. Если метод bicg сходится, выводится соответствующее сообщение. Если метод не сходится после максимального числа итераций или по другой причине, на экран выдается относительный остаток norm(B-A*X)/ norm(B) и номер итерации, на которой метод остановлен;
[X.flag.relres] = bicg(A,X,tol .maxit.Ml,M2.X0) — также возвращает относительную вторую норму вектора остатков relres=nQnr)(B-A*X)/norm(B). Если флаг flag равен 0, то rel res<tol;
[X, flag, rel res, iter] = bicgCA.B.tol,maxit,Ml,M2.XO) — также возвращает номер итерации, на которой был вычислен X. Значение iter всегда удовлетворяет условию 0<iter<maxit;
[X.flag.relres.iter.resvec] = bicgCA.B.tol,maxit,Ml,M2.XO) — также возвращает вектор вторых норм остатков resvec для каждой итерации начиная с res-vec(l)=norm(B-A*X0). Если флаг flag равен 0, то resvec имеет длину iter+1 и resvec(end)<tol*norm(B). Возможны значения flag, равные 0, 1, 2, 3 и 4. Эти значения предоставляют следующие данные о сходимости решения:
flag=0 - решение сходится при заданной точности tol и числе итераций не более заданного maxit;
flag=l - число итераций равно заданному maxit, но сходимость не достигнута;
f l ag=2 - матрица предусловий М плохо обусловлена;
fl ag=3 - процедура решения остановлена, поскольку две последовательные оценки решения оказались одинаковыми;
fl ag=4 - одна из величин в процессе решения вышла за пределы допустимых величин чисел (разрядной сетки компьютера).
Пример:
» bicg(A.B)
BICG converged at iteration 4
to a solution with relative residual
2.3e-015
ans=
1.0000
2.0000
3.0000
4.0000
[X.flag] = bicg(A,X.tol ,maxit.Ml,M2,X0) — возвращает решение Х и флаг flag, описывающий сходимость метода.
