Численные методы (пакет MATLAB). Пантелеев И.Н

относительно векторной функции F векторного аргумента x. Следовательно, исходную задачу можно рассматривать как

задачу о нулях нелинейного отображения F : Rn ® Rn . В этой

постановке данная задача является прямым обобщением зада-

чи о нахождении решения нелинейного уравнения для случая пространств большей размерности. Это означает, что можно строить методы ее решения как на основе обсужденных в пре- дыдущей главе подходов, так и осуществлять формальный пе- ренос выведенных для скалярного случая расчетных формул. Однако не все результаты и не все методы оказывается воз- можным перенести формально (например, метод половинного деления). В любом случае следует позаботиться о пра- вомерности тех или иных операций над векторными перемен- ными и векторными функциями, а также о сходимости полу- чаемых таким способом итерационных процессах. Отметим, что переход от n= 1 к n ≥2 вносит в задачу нахождения нулей нелинейного отображения свою специфику, учет которой при-

вел к появлению новых методов и различных модификаций уже имеющихся методов. В частности, большая вариативность

методов решения нелинейных систем связана с разнообразием способов, которыми можно решать линейные алгебраические задачи, возникающие при пошаговой линеаризации данной не- линейной вектор-функции F(x) .

Начнем изучение методов решения нелинейных систем с метода простых итераций.

Пусть система (IV.1) преобразована к следующей эквива-

или в компактной записи:

Для задачи о неподвижной точке нелинейного отображения F = Rn ® Rn запишем формальное рекуррентное равенство:

x (k +1) = F(x (k ) ), (IV.5)

где k определяет метод простых итераций, для задачи (IV.3) и k= 0,1,2,...,n .

Если начать процесс построения последовательности (x (k ) ) с некоторого вектора x (0) = (x1(0), x2(0) ,K, xn(0) )T и продолжить вы-

числительный процесс по формуле IV.5), то при определенных условиях данная последовательность со скоростью геометриче-

ской прогрессии будет приближаться к вектору x* — неподвиж- ной точке отображения Ф(х).

Справедлива следующая теорема, которую мы приводим без доказательства.

Теорема IV.1. Пусть функция Ф(х) и замкнутое множество M Í D(F) Í Rn таковы, что:

1.F(x)Î M, "x Î M.

Тогда Ф(х) имеет в М единственную неподвижную точку x* ; последовательность (x (k ) ), определяемая (IV.5), сходится

Отметим низкую практическую ценность данной теоремы из-за не конструктивности ее условий. В случаях, когда вы-

брано хорошее начальное приближение x (0) решению x* , больший практический интерес представляет следующая тео- рема.

Теорема IV.2. Пусть F(x) дифференцируема в замкну-

IV.1 с М = S .

Запишем метод последовательных приближений (IV.5) в развернутом виде:

Сравнение (IV.6) с вычислительной формулой метода про- стой итерации решения систем линейных уравнений (III.13) обнаруживает их сходство. Учитывая, что в линейном случае, как правило, более эффективен метод Зейделя, а данном слу- чае также может оказаться более эффективным его многомер- ный аналог, называемый методом покоординатных итераций:

Заметим, что, как и для линейных систем, отдельные урав- нения в (IV.7) неравноправны, т, е. перемена местами уравне- ний системы (IV.3) может изменить в некоторых пределах чис- ло итераций и вообще ситуацию со сходимостью последова- тельности итераций. Для того чтобы применить метод простых итераций (IV.6) или его "зейделеву" модификацию (IV.7) к ис- ходной системе (IV.1), необходимо сначала тем или иным спо- собом привести эту систему к виду (IV.3). Это можно сделать, например, умножив (IV.2) на неособенную nхn матрицу А и

эквивалентна исходной и имеет вид, аналогичный уравнению в методе итераций в одномерном случае. Проблема состоит лишь

вправильном подборе матричного параметра.

4.2.Метод Ньютона решения систем нелинейных

уравнений

Для решения системы (IV.3) будем пользоваться методом

Тогда точный корень уравнения (IV.2) можно представить в виде:

x* = x (k ) + Dx (k ), (IV.9)

где Dx (k ) = (Dx1(k ), Dx2(k ),K, Dxnk ) — поправка (погрешность кор- ня).

Подставляя выражение (IV.9) в (IV.2), имеем

Предполагая, что функция F(x) непрерывно дифференцируе-

ма в некоторой выпуклой области, содержащей x и x (k ) , раз- ложим левую часть уравнения (IV.10) по степеням малого век-

тора Dx (k ) , ограничиваясь линейными членами:

Из формул (IV.11) и (IV.12) видно, что под производной F ′(x)следует понимать матрицу Якоби системы функций

f1 , f2 ,K, fn относительно переменных x1 , x2 ,Kxn , т. е.

поэтому формула (IV.12) может быть записана в следующем виде:

Отсюда видно, что метод Ньютона для решения системы (IV.1) состоит в построении итерационной последовательности:

x (k +1) = x (k ) -W −1 (x (k ) )F(x (k ) ), (IV.15)

где k =0, 1,2,

Если все поправки становятся достаточно малыми, счет прекращается. Иначе новые значения хi, используются как при- ближенные значения корней, и процесс повторяется до тех пор, пока не будет найдено решение или не станет ясно, что полу- чить его не удастся.

Аналогично находятся дальнейшие приближения. Результаты вычислений приведены в табл. 4.1.

Таблица 4.1.

Последовательные приближения корней

Останавливаясь на приближении .x(3), будем иметь: x= 0.7852;

у = 0.4966; г = 0.3699.

4.3. Решение нелинейных систем методами спуска

Общий недостаток всех рассмотренных ранее методов ре- шения систем нелинейных уравнений состоит в локальном ха- рактере сходимости, затрудняющем их применение в случаях (достаточно типичных), когда существуют проблемы с выбо- ром начального приближения, обеспечивающего сходимость итерационной вычислительной процедуры. В этих случаях можно использовать численные методы оптимизации — дан- ный раздел вычислительной математики выделяется в само-

стоятельную дисциплину. Для использования наглядной гео-

метрической интерпретации приводимых ниже рассуждений и их результатов ограничимся, как и в предыдущем пункте, рас- смотрением системы, состоящей из двух уравнений с двумя не-

известными

Из функций f(x,y), g(x,y) системы (IV.16) образуем новую функцию.

F(x, y) = f (x, y)2 + g(x, y)2 . (IV.17)

Так как эта функция не отрицательная, то найдется точка (вообще говоря, не единственная) (х*, у*) такая, что

F(x, y) ³ F(x*, y *) ³ 0 "(x, y)Î R2 .

Рис. IV.1. Пространственная интерпретация метода наиско- рейшего спуска для функции (IV.17)

Следовательно, если тем или иным способом удается полу- чить точку (x*,y*), минимизирующую функцию Ф(х,у), и если

при этом окажется, что min F(x, y) = F(x*, y *) = 0 , то точка

(x, y ) R2

(х*,у*) — истинное решение системы (IV.16), поскольку

Последовательность точек (xk,yk) — приближений к точке (х*, у*) минимума функции Ф(x,у) — обычно получают по ре-

где k =0,1,2,..., (pk,qk)r - вектор, определяющий направление ми- нимизации, a αk— скалярная величина, характеризующая ве- личину шага минимизации (шаговый множитель). Учитывая геометрический смысл задачи минимизации функций двух пе- ременных Ф(х, у) — "спуск на дно" поверхности z = Φ(x:,у) (рис. IV.1), итерационный метод (IV.18) можно назвать мето- дом спуска, если вектор (pk,qk)T при каждом k является направ- лением спуска (т. е. существует такое а >0, что Ф(хk + apk,yk +αqk)<Φ(xk,yk)), и если множитель ak подбирается так, чтобы

выполнялось условие релаксации Ф(xk+1,yk+1)<Ф(xk,yk) озна- чающее переход на каждой итерации в точку с меньшим значе-

нием минимизируемой функции.

Таким образом, при построении численного метода вида (IV.18) минимизации функции Ф(х,у) следует ответить на два главных вопроса: как выбирать направление спуска (pk,qk)T и как регулировать длину шага в выбранном направлении с по- мощью скалярного параметра — шагового множителя ak. При- ведем простые соображения по этому поводу.

При выборе направления спуска естественным является выбор такого направления, в котором минимизируемая функция убы-

вает наиболее быстро. Как известно из математического ана- лиза функций нескольких переменных, направление наиболь- шего возрастания функции в данной точке показывает ее гра- диент в этой точке. Поэтому примем за направление спуска

вектор

— антиградиент функции Ф(х,у). Таким образом, из семейств (IV.18) выделяем градиентный метод.

Оптимальный шаг в направлении антиградиента — это та- кой шаг, при котором значение Ф(xk+1,yk+1 ) наименьшее среди всех других значений Ф(x,у) в этом фиксированном направле- нии, т. е. когда точка (xk +1, yk+1 ) является точкой условного

минимума. Следовательно, можно рассчитывать на наиболее быструю сходимость метода (IV.19), если полагать в нем

Такой выбор шагового множителя, называемый исчерпы- вающим спуском, вместе с формулой (IV.19) определяет метод наискорейшего спуска.

Геометрическая интерпретация этого метода хорошо вид- на на рис. IV.1, IV.2. Характерны девяностоградусные изломы траектории наискорейшего спуска, что объясняется исчерпы- ваемостью спуска и свойством градиента (а значит, и антигра- диента) быть перпендикулярным к линии уровня в соответст- вующей точке.

Наиболее типичной является ситуация, когда найти точно (аналитическими методами) оптимальное значение ak не уда- ется. Следовательно, приходится делать ставку на применение каких-либо численных методов одномерной минимизации и находить ak в (IV.18) лишь приближенно.

Несмотря на то, что задача нахождения минимума функции одной переменой ϕk (a) = Φ(xk − aΦ′x (xk , yk ), yk − aΦ′y (xk , yk )). намного проще, чем решаемая задача, применение тех или

иных численных методов нахождения значений ak = argmin ϕk (a) с той или иной точностью требует вычисле-

ния нескольких значений минимизируемой функции.

Так как это нужно делать на каждом итерационном шаге, то

при большом числе шагов реализация метода наискорейшего спуска в чистом виде является достаточно высокозатратной.

Существуют эффективные схемы приближенного вычисления квазиоптимальных ak, в которых учитывается специфика мини- мизируемых функций (типа сумм квадратов функций).

Зачастую успешной является такая стратегия градиентного метода, при которой шаговый множитель ak в (IV.18) берется либо сразу достаточно малым постоянным, либо предусматри- вается его уменьшение, например, делением пополам для удов- летворения условию релаксации на очередном шаге. Хотя каж- дый отдельный шаг градиентного метода при этом, вообще го- воря, далек от оптимального, такой процесс по числу вычисле- ний функции может оказаться более эффективным, чем метод наискорейшего спуска. Главное достоинство градиентных ме- тодов решения нелинейных систем — глобальная сходимость. Нетрудно доказать, что процесс градиентного спуска приведет к какой-либо точке минимума функции из любой начальной точки. При определенных условиях найденная точка минимума будет искомым решением исходной нелинейной системы.

Главный недостаток — медленная сходимость. Доказано, что сходимость этих методов только линейная, причем, если для многих методов, таких как метод Ньютона, характерно ус- корение сходимости при приближении к решению, то здесь имеет место скорее обратное. Поэтому есть необходимость по- строения гибридных алгоритмов, которые начинали бы поиск искомой точки — решения данной нелинейной системы — гло- бально сходящимся градиентным методом, а затем производи-

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

Примечание.

Порядком сходимости последовательности (xk) к х* назы-

при всех k>k0.

Разработан ряд методов решения экстремальных задач, ко-

торые соединяют в себе низкую требовательность к выбору начальной точки и высокую скорость сходимости. К таким ме- тодам, называемым квазиньютоновскими, можно отнести, на- пример, метод переменной метрики (Дэвидона-Флетчера- Пауэлла), симметричный и положительно определенный мето- ды секущих (на основе формулы пересчета Бройдена).

При наличии негладких функций в решаемой задаче следу- ет отказаться от использования производных или их аппрок- симаций и прибегнуть к так называемым методам прямого по- иска (циклического покоординатного спуска, Хука и Дживса, Роленброка и т. п.). Описание упомянутых и многих других

методов такого типа можно найти в учебной и в специальной литературе, посвященной решению экстремальных задач.

Для разных семейств численных методов минимизации мо- гут быть рекомендованы свои критерии останова итерацион- ного процесса. Например, учитывая, что в точке минимума дифференцируемой функции должно выполняться необходи- мое условие экстремума, на конец счета градиентным методом можно выходить, когда достаточно малой становится норма градиента. Если принять во внимание, что минимизация при- меняется к решению нелинейной системы, то целесообразнее отслеживать близость к нулю минимизируемой неотрицатель- ной функции, т. е. судить о точности получаемого приближе- ния по квадрату его евклидовой метрике.

Для решения n-мерной системы (IV.1) следует свести зада- чу к решению экстремальной задачи:

F(x) = å fi 2 (x) ® min .

i=1n

Рассмотрим далее примеры реализации некоторых алгорит- мов поиска экстремумов функций, зависящих от нескольких переменных, в пакете MATLAB.

Пример IV.1. Алгоритм поиска экстремума функции с шагом, не зависящим от свойств минимизируемой функции.

Простейший вариант метода наискорейшего спуска рассмотрим на примере поиска минимума квадратической функции

f (x, y) = x2 + μy2 двух переменных с оврагом, пологость ко-

торого определяется параметром μ. Решение данной задачи в пакете MATLAB находится выполнением следующей по- следовательности команд:

1.Создайте файл F_L4.m (листинг IV.1), содержащий описание функции, возвращающей значения функции f(x,y) в узлах коор- динатной сетки.

Листинг 4.1. Файл F_L4.m function z=F_L4 (x,y,mu) N=length(x);

z=zeros(N); for i=l:N for j=l:N

z(i,j)=x(i).^2+mu*y(j).^2; end;

end;

2.Постройте графики исследуемой функции при различных значениях параметра μ.

» N=23;

» Xmin=-5; Xmax=5; » Ymin=-5;Ymax=5; » i=l:N;j=l:N;

» x(i)=Xmin+i*(Xmax-Xmin)/N; » y(j)=Ymin+j*(Ymах-Ymin)/N;

»Ml=F_L4(x,y,0.5); M2=F_L4(x,y,l) ;M3=F_L4(x,y, 1.5) ;

»[X Y]=meshgrid(x,y);

»surfc(X,Y,Ml); colormap gray

»surfc(X,Y,M2); colorap gray

»surfc(X,Y,M3); colormap gray

Из рис. IV.3-IV.5 видно, что при μ = l функция f(x,y,μ) пред- ставляет собой параболоид вращения, при μ>1 параболоид становится эллиптическим, "вытягиваясь" вдоль оси оХ (при μ < 1 — вдоль оси оY).

Рис. IV.3. Поверхность и карта линий уровня функции f (x, y) = x2 + 0.5 × y2