Точное решение, метод наименьших квадратов и сопряженных градиентов

· lsqr(A, В)—возвращает точное решение X СЛУ А*Х=В, если матрица последовательная, в противном случае — возвращает решение, полученное итерационным методом наименьших квадратов. Матрица коэффициентов А должна быть прямоугольной размера тхя, а вектор-столбец правых частей уравнений В должен иметь размер т. Условие m>=n может быть и необязательным. Функция 1 sqr начинает итерации от начальной оценки, по умолчанию представляющей собой вектор размером п, состоящий из нулей. Итерации производятся или до сходимости к решению, или до появления ошибки, или до достижения максимального числа итераций (по умолчанию равного min(20, m, n) — либо 20, либо числу уравнений, либо числу неизвестных). Сходимость достигается, когда отношение вторых норм векторов norm(B-Ax)/norm(B) меньше или равно погрешности метода tol (по умолчанию 1е-б);

· lsqr(A.B,tol) — возвращает решение с заданной погрешностью (порогом отбора) tol;

· lsqr(A,b.tol .maxlt) — возвращает решение при заданном максимальном числе итераций maxit вместо, возможно, чересчур малого числа, заданного по умолчанию;

· lsqr(A,b.tol .maxit,M) и lsqr(A,b,tol .maxit.Ml.M2) — при решении используются матрица предусловий М или М=М1*М2, так что производится решение системы inv(M)*A*x=inv(M)*b относительно х. Если Ml или М2 — пустые матрицы, то они рассматривается как единичные матрицы, что эквивалентно отсутствию входных условий вообще;

· lsqr(A.B,tol .maxit.Ml.M2.X0) — точно задается начальное приближение Х0. Если Х0 — пустая матрица, то по умолчанию используется вектор, состоящий из нулей;

· X = lsqr(A,B.tol .maxit,Ml.M2.X0) — при наличии единственного выходного параметра возвращает решение X. Если метод 1 sqr сходится, выводится соответствующее сообщение. Если метод не сходится после максимального числа итераций или по другой причине, на экран выдается относительный остаток попп(В-А*Х)/ norm(B) и номер итерации, на которой метод остановлен;

· [X.flag] = lsqr(A.X.tol.maxit.Ml.M2.X0) — возвращает решение X и флаг flag. описывающий сходимость метода;

· [X.flag.relres] = lsqr(A,X,tol.maxit,Ml.M2.X0) — также возвращает относительную вторую норму вектора остатков rel res=norm(B-A*X)/norm(B). Если флаг flag равен 0, то relres<tol;

· [X.flag.relres.iter] = bicg(A,B.tol,maxit,Ml,M2.X0) — также возвращает номер итерации, на которой был вычислен X. Значение iter всегда удовлетворяет условию 0<iter<maxit;

· [X.flag.relres,iter,resvec]= lsqr(A.B.tol.maxit.Ml.M2.X0) — также возвращает вектор вторых норм остатков 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 предоставляют следующие данные о сходимости решения:

o flag=0 - решение сходится при заданной точности tol и числе итераций не более заданного maxit;

o flag=l - число итераций равно заданному maxit, но сходимость не достигнута;

o flag=2 - матрица предусловий М плохо обусловлена;

o flag=3 - процедура решения остановлена, поскольку две последовательные оценки решения оказались одинаковыми;

o fl ag=4 - одна из величин в процессе решения вышла за пределы допустимых величин чисел (разрядной сетки компьютера).

Если значение flag больше нуля, то возвращается не последнее решение, а то решение, которое имеет минимальное значение отношения вторых норм векторов norm(B-A*x)/norm(B).

Пример:

» А=[0 012; 1300; 0101; 1010];

» В=[11; 7; 6; 4];

Введенные в этом примере матрица А и вектор В будут использованы и в других примерах данного раздела. В примере процесс итераций сходится на пятом шаге с относительным остатком (отношением вторых норм векторов невязки и свободных членов) 1.9 10- ¹³ .

Пример:

» lsqr(A,B.1e-6.5)

Isqr converged at iteration 5 to a solution

with relative residual

1.9e-013

ans =

1.0000

2.0000

3.0000

4.0000

Метод минимизации обобщенной невязки

Итерационный метод минимизации обобщенной невязки также реализован в системе MATLAB. Для этого используется функция gmres:

· gmres (А, В. restart) — возвращает решение X СЛУ А*Х=В.

А —квадратная матрица. Функция gmres начинает итерации от начальной оценки, представляющей собой вектор размера n , состоящий из нулей. Итерации производятся либо до сходимости к решению, либо до появления ошибки, либо до достижения максимального числа итераций. Сходимость достигается, когда относительный остаток norm(B-A*X)/norm(B) меньше или равен заданной погрешности (по умолчанию 1е-6). Максимальное число итераций — минимум из n/restart и 10. Функция gmres (...) имеет и ряд других форм записи, аналогичных описанным для функции bieg(...). Пример:

» А=[0 0 1 2; 1 3 0 0; 0 1 0 1; 1 0 1 0];

» В=[11; 7; 6; 4];

» gmres(A,B)

GMRES(4) converged at Iteration 4 to a solution with relative residual

1e-016

ans =

1.0000

2.0000

3.0000

4.0000

Вычисление нулей функции одной переменной

Ряд функций системы MATLAB предназначен для работы с функциями. По аналогии с дескрипторами графических объектов могут использоваться объекты класса дескрипторов функций, задаваемых с помощью символа @, например:

» fe=@exp.

Примечание

Под функциями понимаются как встроенные функции, например sin(x) или ехр(х),так и функции пользователя, например f(x), задаваемые как M-файлы-функции.

Численные значения таких функций, заданных дескрипторами, вычисляются с помощью функции feval:

» feval(fe.1.0)

ans =

2.7183

Для совместимости с прежними версиями можно записывать функции в символьном виде в апострофах, использование функции eval для их вычисления может быть более наглядно, не нужно создавать m-файл, но в учебном курсе мы будем стараться использовать новую нотацию, с использованием дескрипторов функций и feval, так как при этом программирование становится «более объектно-ориентированным», повышается скорость, точность и надежность численных методов. Поэтому, хотя везде в нижеследующем тексте вместо @fun можно подставить и символьное значение функции в апострофах, мы будем использовать нотацию @fun в дидактических целях. Все же иногда в интерактивном режиме можно использовать старую запись, чтобы не создавать m-файл функции.

Довольно часто возникает задача решения нелинейного уравнения вида f(x) = 0 или f ₁ (х) = f ₂ (х). Последнее, однако, можно свести к виду f(x) =f ₁ (х) - f ₂ (х) = 0. Таким образом, данная задача сводится к нахождению значений аргумента х функции f(x) одной переменной, при котором значение функции равно нулю. Соответствующая функция MATLAB, решающая данную задачу, приведена ниже:

· fzero(‘fun’,x) (или равнозначная запись fzero(@fun,x) )— возвращает уточненное значение х, при котором достигается нуль функции fun, представленной в символьном виде, при начальном значении аргумента х. Возвращенное значение близко к точке, где функция меняет знак, или равно NaN, если такая точка не найдена;

· fzero(@fun,[xl x2]) — возвращает значение х, при котором fun(x)=0 с заданием интервала поиска с помощью вектора x=[xl х2], такого, что знак fun(x(1) отличается от знака fun(x(2)). Если это не так, выдается сообщение об ошибке. Вызов функции fzero с интервалом гарантирует, что fzero возвратит значение, близкое к точке, где fun изменяет знак;

· fzero(@fun,x.tol) — возвращает результат с заданной погрешностью tol;

· fzero(@fun,x.tol .trace) — выдает на экран информацию о каждой итерации;

· fzero(@fun,х.tol .trace,Р1.Р2,...) — предусматривает дополнительные аргументы, передаваемые в функцию fun(x.Pl,P2,...). При задании пустой матрицы для tol или trace используются значения по умолчанию. Пример:

fzero(fun,x,[ ],[ ],Р1).

Для функции fzero ноль рассматривается как точка, где график функции fun пересекает ось х, а не касается ее. В зависимости от формы задания функции fzero реализуются следующие хорошо известные численные методы поиска нуля функции: деления отрезка пополам, секущей и обратной квадратичной интерполяции. Приведенный ниже пример показывает приближенное вычисление р/2 из решения уравнения cos(x)=0 с представлением косинуса дескриптором:

» х= fzero(@cos.[1 3]) %аналог х= fzero(‘cos’.[1 3])

x =

1.5708

В более сложных случаях настоятельно рекомендуется строить график функции f(x) для приближенного определения корней и интервалов, в пределах которых они находятся. Ниже дан пример такого рода (следующий листинг представляет собой содержимое m-файла fun1.m):

%Функция, корни которой ищутся

function f=funl(x) %сохранить как М-файл с именем funl

f=0.25*x+sin(x)-1;

Построим график функции:

» х=0:0.1:10;

» plot(x,funl(x));grid on;

Из рисунка нетрудно заметить, что значения корней заключены в интервалах [0.5 1], [2 3] и [5 6]. Найдем их, используя функцию fzero:

» xl=fzero(@funl.[0.5 1])

xl =

0.8905

» x2=fzero(@funl.[2 3])

x2 =

2.8500

» x3=fzero(@funl,[5.6])

x3 =

5.8128

» x3=fzero(@funl,5,0.001)

x3 =

5.8111

Обратите внимание на то, что корень хЗ найден двумя способами и что его значения в третьем знаке после десятичной точки отличаются в пределах заданной погрешности tol =0.001. К сожалению, сразу найти все корни функция fzero не в состоянии.

Решим эту же систему при помощи функции fsolve из пакета Optimization Toolbox, которая решает систему нелинейных уравнений вида f(x)=0 методом наименьших квадратов, ищет не только точки пересечения, но и точки касания, fsolve имеет почти те же параметры (дополнительный параметр — задание якобиана) и почти ту же запись, что и функция lsqnonneg, подробно рассмотренная ранее. Пример:

»fsolve(@funl,0:10 )

ans =

Columns 1 through 7

0.8905 0.8905 2.8500 2.8500 2.8500 5.8128 5.8128

Columns 8 through 11

5.8128 2.8500 2.8500 10.7429

!!!!!!!!!!!!!!

Для решения систем нелинейных уравнений следует также использовать функцию solve из пакета Symbolic Math Toolbox. Эта функция способна выдавать результат в символьной форме, а если такого нет, то она позволяет получить решение в численном виде. Пример:

» solve('0.25*x + sin(x) -1)

ans =

.89048708074438001001103173059554

Не нашли, что искали? Воспользуйтесь поиском по сайту: