Итерации

6. Итерации

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

x_k+1=F(x_k), k=0, 1, 2, ... , x₀ задано,

где F – заданное преобразование y=F(x), x₀ – как-то выбранное начальное приближение, x_k – значение переменной x на k-й итерации, а сама переменная x может быть любой – числом, вектором или матрицей. Предел итераций, если он существует, будем обозначать через X, и тогда уравнение

X=F(X) (1)

должно обязательно выполняться для итерационного преобразования F. Это обстоятельство помогает выбирать различные варианты для F. Все решения X этого уравнения называются неподвижными точками преобразования F(x). Все такие x₀, для которых последовательность x_k сходится, образуют область сходимости итерационного преобразования F. Cкорость сходимости итераций x_k характеризуется поведением числовой последовательности

v_k=norm(x_k+1-x_k)/norm(x_k-x_k-1),

где норма может выбираться по-разному. Для сходимости x_k уже не обязательно, чтобы существовал предел V последовательности v_k, но очень часто для сходящихся итераций он существует, и если это так, то пусть a=abs(V). Тогда при a=1 сходимость x_k будет крайне медленной, при 0<a<1 это будет сходимость по геометрической прогрессии со знаменателем q=V, а при a=0 последовательность x_k будет сходиться быстрее любой геометрической прогрессии. Покажем теперь на простейших примерах, как все это выглядит на деле. Чтобы иметь возможность разнообразия вариантов, задачу возьмем нелинейной. Рвссмотрим уравнение

(x-2)(x-3)=0 или f(x)ºx²-5x+6=0 c корнями x₁=2, x₂=3, (2)

построим для нахождения его решений x₁=2, x₂=3 несколько итерационных преобразований или схем F и проанализируем их работу.

1.Пусть для уравнения (2)

x=F(x), где y=F(x)=(x²+6)/5.

Обязательное условие (1) для преобразования F выполнено, однако при этом переходе появилось дополнительное решение x=inf (µ). Потеря некоторых или приобретение новых решений часто случается при переходе от исходного уравнения к его итерационной форме. Переходя к нужной нам индексной записи, будем иметь

x(k+1)=F(x(k)), k=1:n,

где начальное приближение x(1) и число итераций n должны быть заданы. Будем накапливать значения x(k) в переменной x, а текущее значение x(k) обозначим через xt. Нужные вычисления реализует строка

1;xt=0; n=100; x=xt; TF='(xt^2+6)/5'; for k=1:n,xt=eval(TF); x=[x,xt]; end, plot(x), grid

Здесь записаны текстовая переменная TF и команда eval (она интерпретирует TF как выражение xt^2+6)/5 и выполняет его). После выполнения строки 1 на графике отобразится 101 значение x(k), включая начальное приближение. Итерации сошлись к px=2. Строка

2;xt=0; n=100; x=xt; TF='xt=(xt^2+6)/5;'; for k=1:n,eval(TF),x=[x,xt]; end, plot(x), grid

произведет те же вычисления, но обратите внимание на то, как в ней записаны TF и eval.

Хотя внешне сходимость x(k) к x₁=2 не вызывает сомнений, это в действительности всегда нужно проверять более тщательно. Так как y'=F'(x)=2x/5, то F'(2)=0.8<1, а F'(3)=1.2>1, но итерации, как известно, могут сойтись только к такой неподвижной точке X преобразования F, для которой |F'(X)|£1. Отсюда ясно, почему X=2. Однако далеко не всегда можно найти F'(x) аналитически, и поэтому в общем случае приходится использовать вычислительный подход. При известных уже x(k) его реализует строка

3;w=2:n; v=(x(w+1)-x(w))./(x(w)-x(w-1)); plot(v), grid

Из графика видно, что v(k) действительно сходятся к a=0.8, как и должно быть согласно теории, которая здесь выглядит очевидной.

Теперь будем варьировать начальное приближение, выполняя строки 1 и 3.

(a) xt=1.5, 1.9, 1.99, 1.9999. При последних двух значениях xt на графике из строки 3 появятся осцилляции справа, так как здесь разности x(k+1)-x(k) и тем самым значения v(k) теряют точность: x(k) уже сошлись со многими знаками.

(a') xt=2 . График из строки 2 вообще пуст, потому что тогда все v(k)=0/0=NaN.

(b) xt=2.01, 2.5, 2.9, 2.99, 2.9999. В последнем случае x(k) вначале довольно долго (до k=20) задерживаются в районе неподвижной точки x₂=3 (они все время уходят от нее, но на графике этого не видно), а затем примерно за 60 итераций монотонно движутся к x₁=2.

(b') xt=3 . Все x(k)=3, а график из строки 3 пуст. Это произошло потому, что при xt=3 F(xt)=15/5=3 получается без ошибок округления.

(b'') xt=3.01 . Пределами для x(k) и v(k) будет inf, так что x₂=3 является неустойчивой неподвижной точкой преобразования F: при малейшем сдвиге x₀ с x₂ в пределе итераций получится либо устойчивая неподвижная точка x₁, либо inf – приобретенная неподвижная точка. Проварьируем этот сдвиг:

(с) xt=3+1e-15, 3+1e-14, 3+1e-8 . В 1-м варианте ухода xt не происходит, во 2-м тенденция ухода уже зародилась и он неизбежно произойдет с увеличением числа итераций, в 3-м уход проявился в полной мере уже на 100 итерациях.

(с') xt=-2.5, -2.9, -2.99, -3, -3.01, -100, 100 . При xt=-3 снова получим x₂ за одну итерацию, а при xt= -3.01 и далее получим inf. Таким образом, при вещественном x₀ итерации сходятся к устойчивой неподвижной точке x₁=2 при -3<x₀<3 и к inf при |x₀|>3. Просчитаем тот случай, когда |x₀|=3. Этот расчет выполняется строкой (редактируем строку 1)

4;n=100; fi=-pi:pi/20:pi; xt=3*exp(i*fi); TF='(xt.^2+6)/5'; for k=1:n,xt=eval(TF);end, plot(xt,'.')

На графике (он комплекснозначный) видны 4 точки: точка z=3, точки z=3±s*i с малым s>0 и точка z=2. Чтобы разобраться в результате, выполним строку 4 с n=1000 и, сделав окно MATLAB'а полным, выдадим

5;imag(xt')

Образы первой и последней точки начальной окружности слегка отличаются от z=2, а ее 21-я точка (она соответствует fi=0) есть z=3. Это получилось потому, что sin(pi) и sin(-pi) не равны нулю в точности. Снова сделаем окно MATLAB'а небольшим и выполним строку 4 с радиусом 3.01, а затем строки

6;sum(isnan(xt))

7;find(isnan(xt))

и получим, что точки первоначальной xt с номерами 1, 21, 41 обращаются в inf (здесь nan=inf/inf). Для радиуса 5 их будет 21, для радиуса 5.6 их будет 41, т.е. все.

Границу области сходимости обычно трудно исследовать аналитически, даже в таком простом примере, как этот, и сама она не определяется из условия |F'(x)|=1 или |x|=2.5. Условие |F'(X)|<1 гарантирует устойчивость неподвижной точки X преобразования F(x), условие |F'(X)|>1 является достаточным для ее неустойчивости, а в случае |F'(X)|=1 она может быть как устойчивой, так и неустойчивой. Неустойчивые неподвижные точки сравнительно редки в вычислительных задачах, но их полезно иметь в виду при исследовании численных алгоритмов. В нашем случае мы установили только, что множество |x|£3 принадлежит области сходимости итераций F, но вовсе не описали границу этой области. Мы установили также, что при |x|³5.6 итерации сходятся к inf, и поэтому inf является устойчивой неподвижной точкой F. Чтобы получить представление о границе области сходимости, выполним строки

8;r=3:.1:5.6; z=[];n=100; for kr=r,xt=kr*exp(i*fi); for k=1:n,xt=eval(TF); end,z=[z;xt]; end

9;zn=isnan(z); Z=r'*exp(i*fi); plot(Z(zn)), axis equal

и увидим приближенный график границы – это вовсе не окружность. Комагнда axis equal выбирает по осям одинаковый масштаб (это действует только на текущий plot; у команды axis много других опций). Чтобы построить границу аккуратнее, выполним строку 4 с шагом по fi, равным pi/180, затем строку 8 с r=3:.02:5.6 и строку 9. Получим график границы, выполнив строки

10;zn=isnan([z;z(end,:)]); zn=diff(zn)~=0; plot(Z(zn)), hold on

11;y=3*exp(i*fi); plot(y,'m'), axis equal, hold off

и увидим отличие области сходимости от круга |z|£3.

Найдем те точки, которые перейдут в z=3 на первых 10 итерациях. Для этого придется рассмотреть обратное к F преобразование x=±(5y-6)^.5 и сделать с ним 10 итераций, задав начальное значение y=3. Строка

12;n=10; yt=3; for k=1:n,yt=(5*yt-6).^.5; yt=[yt,-yt];end, plot(yt,'.')

показывает, что при этом получится практически та же линия, что и в строке 10. Удивительно, что это верно и для других точек z из области сходимости преобразования F, но с некоторыми вкраплениями внутрь области сходимости. Брать n большим нельзя, так как в конце длина вектора yt равна 2ⁿ.

2.Теперь представим (2) в виде

x=F(x), где y=F(x)=5-6/x, так что F'(x)=6/x² , F'(3)=2/3<1, |F'(x)|<1 при |x|>2.45.

Следовательно, устойчивая неподвижная точка – это x₀=3. Получим резултат сразу для "всех" вещественных x₀=xt:

1;n=100; xt=-5:.5:5; x=xt; TF='5-6./xt'; for k=1:n,xt=eval(TF);end, plot(x,xt,'.')

Все пределы равны 3, выпадает только неустойчивая точка xt=2, для которой итерации опять идут без ошибок округления. Возникает предположение, что вся комплексная плоскость с выколотой точкой x₁=2 стягивается итерациями в точку x₂=3, хотя не везде |F'(x)|<1. Посмотрим, как это можно увидеть на графике, выполнив программу

2;n=4; fi=-pi:pi/20:pi; xt=4*exp(i*fi); TF='5-6./xt'; z=xt;

3;for k=1:n, xt=eval(TF); z=[z;xt];end, plot(z','.'), axis equal

Начальное приближение (окружность радиуса 4 с центром в нуле) итерируется 4 раза и все результаты выдаются на график. Окружности (на графике их 5) довольно быстро стягиваются в точку x₂=3. Применим zoom, чтобы посмотреть малые окружности.

Сделаем разрезы на начальной окружности:

4;n=4; fi=-pi:pi/20:pi*.75; fi(end-4)=[]; xt=4*exp(i*fi); TF='5-6./xt'; z=xt;

и снова выполним строку 3. Мы увидим, что 1-я итерация меняет направление обхода окружности на противоположное, остальные – нет. Выполним строку 4 с n=20 и затем строку 3. С помощью zoom дойдем до самой малой окружности (она белого цвета) и убедимся, что она лежит правее точки z=3 примерно в полосе от 3.0001 до 3.0004.

Потерь и приобретений других неподвижных точек здесь нет.

3. Решим нашу задачу методом Ньютона. Если x'' удовлетворяет уравнению f(x'')=0, а x' находится вблизи x'' и используется для приближенной аппроксимации f(x''), то

0=f(x'')=f (x')+f'(x')(x''-x') или x''=x'-f(x')/f'(x') при условии, что f'(x'')¹0 и тем самым f'(x')¹0.

Это и есть итерации по Ньютону для уравнения f(x)=0. При переходе к индексной форме записи для f(x)=x²-5x+6 и f'(x)=2x-5 из (2) получим

x(k+1)=x(k)-f(x(k))/f'(x(k)) или y=F(x), где F(x)=x-f(x)/(f'(x), F'(x)=2f(x)/(f'(x))² .

Так как f'(x_1,2)¹0, можно использовать эти итерации по Ньютону (часто их называют методом касательной), а поскольку F'(x_1,2)=0, следует ожидать их быстрой сходимости и того, что теперь оба решения x_1,2 будут устойчивыми неподвижными точками итерационного преобразования F. Неподвижной точкой будет и inf, но она "не наша", и природа ее, как мы увидим ниже, гораздо сложнее.

Выделим TF в отдельную строку и проведем наши стандартные вычисления

1;TF='xt-(xt.^2-5*xt+6)./(2*xt-5)';

2;xt=0; n=100; x=xt; for k=1:n,xt=eval(TF); x=[x,xt];end, plot(x), grid

3;w=2:n; v=(x(w+1)-x(w))./(x(w)-x(w-1)); plot(v), grid

Предел итераций при x₀=0 равен 2, а сходимость имеет порядок выше первого, поскольку v(k) до потери ими точности успели подоойти к нулю. Чтобы определить порядок сходимости, отредактируем строку 2 на предмет оценки квадратичной сходимости:

4;w=2:n;v=(x(w+1)-x(w))./(x(w)-x(w-1)).^2; plot(v(1:6)), grid

Теперь до потери точности v(k) успевают подойти к 1 (у v(7) точность уже потеряна), так что сходимость итераций здесь квадратичная. Напомним, что точность теряется у v(k), но не у x(k). Чтобы увидеть потерю точности у v(k), выполним строку 4, выдавая в plot 7 точек.

При x₀=4 предел итераций равен 3, а v(k) успевают подойти к -1 (у v(6) точность уже потеряна). Таким образом, в зависимости от начального приближения x₀=xt итерации могут сходиться к обоим решениям задачи.

Так как теперь преобразование y=F(x) уже не является дробно-линейным, рассмотрим трансформацию комплексной прямоугольной области в процессе итераций. Создадим такую область xt из 41²=1681 комплексных точек и проитерируем ее (при этом в командном окне будет много сообщений о делении на нуль):

5;x=-10:.5:10; [X,Y]=meshgrid(x); xt=X+i*Y; n=100; for k=1:n,xt=eval(TF);end, plot(xt,'.')

На графике будут оба решения задачи и некоторое множество точек на прямой Re(z)=2.5. Сделаем только 10 итераций и вставим выдачу графика на каждой итерации:

6;x=-10:.5:10; [X,Y]=meshgrid(x); xt=X+i*Y;n=10; for k=1:n,xt=eval(TF); plot(xt,'.'), pause(0), end

Конечный результат тот же, что и при 100 итерациях, но видна динамика его формирования.

Полуплоскость Re(z)<2.5 стягивается итерациями в точку x₁=2, а полуплоскость Re(z)>2.5 - в точку x₂=3: прямая Re(z)=2.5 переходит сама в себя. Точка x₁=2 имеет красный цвет потому, что в нее перешли первые 25 столбцов матрицы xt, а остаток rem(25,7)=4, но 4-й цвет – красный. Далее идут зеленые точки (26-й столбец) и синяя x₂=3, поскольку rem(41,7)=6 и 6-й цвет – синий. Но всегда лучше проверить такие выводы численно: найдем

7;z=xt(:); [sum(abs(z-2)<.1), sum(abs(z-3)<.1), sum(abs(real(z)-2.5)<.1)]

 (=1025=41*25), (=615=41*15), (=38<41 на 3).

Так как потерь в матрице-таблице xt MATLAB не допускает, выясним, во что перешли эти 3 точки - в inf или NaN:

8;z=xt(:,26); [sum(isinf(z)), sum(isnan(z))] =0, =3 .

Индексы этих трех точек из 26-го столбца находятся командой

9;find(isnan(z))' =20 21 22 ,

и это есть точки

z₁=2.5-0.5i, z₂=2.5, z₃=2.5+0.5i.

Теперь разберемся, как преобразуется итерациями F прямая Re(z)=2.5. Так как прямая Re(z)=2.5 переходит в себя, ее можно параметризовать с помощью переменной y=Im(z), и пусть y_k=Im(F^k(Re(z)=2.5), k=1, 2,... , т.е. после k итераций y переходит в y_k(y). Ясно, что на k-й итерации в inf перейдут только те точки из всего множества y_k-1, которые равны нулю. Поскольку -inf<y<inf, на 1-й итерации это случится только с одной точкой z₁=2.5 (для нее y=0). Выдадим на график y₁ после 1-й итерации:

10;xt=2.5+(-10:.1:10)'*i; y=imag(xt); n=1; for k=1:n,xt=eval(TF);end, plot(y,imag(xt)), grid

и увидим, что функция y₁(y) пересекает ось абсцисс только два раза (при этом y₁ дважды непрерывно пробегает от -inf до inf), так что на 2-й итерации в inf перейдут только две точки (это уже известные нам z₁=2.5-0.5i, и z₃=2.5+0.5i). Выполним строку 10 с n=2 и снимем с графика строкой

11;q=ginput;q(:,1)'

четыре точки пересечения кривой y₂(y) с осью абсцисс: приближенно это -1.2079 -0.2056 0.2055 1.2079,

тогда как для y₁(y) это были значения -0.5 и 0.5. Каждая переходящая в inf на k-й итерации точка порождает слева и справа от себя две такие точки, которые перейдут в inf на (k+1)-й итерации. Поэтому с ростом k точки y_k(y)=0, с одной стороны, уходят в обе стороны от нуля, а с другой - все чаще появляются в тех местах, где они однажды уже появились. Всего на k-й итерации в inf перейдет 2^k-1 точек. В действительности с ростом k они всюду плотно заполняют прямую Re(z)=2.5, а между ними y_k(y) непрерывно пробегает все значения от -inf до inf. Чтобы быть в этом уверенными, выполним последнюю программу

12;hy=1.0033e-4; xt=2.5+(-1:hy:1)'*i; w=2:length(xt); n=100; for k=1:n,xt=eval(TF); end

13;pzn=sum(sign(xt(w-1)).*sign(xt(w))<0)

которая зафиксирует 674 перемены знака (pzn) у функции y₁₀₀(y) для -1£y£1 с hy=1.0033e-4 (из-за недостаточной малости hy далеко не все они учтены). Для -10£y£-8 с hy=1.0033e-4 pzn=675. Для симметричного относительно точки y=0 отрезка 8£y£10 с тем же hy получим pzn=702. Ошибки при последовательном вычислении F^k(z) на прямой Re(z)=2.5 будут ничтожными, так что значения y_k(y) вычисляются с высокой точностью, но от шага hy число найденных перемен знака в y_k(y), конечно, зависит: при k=100 число всех нулей функции y_k(y) равно 2^100-1=6.3383e29.

Сделаем краткое резюме относительно неподвижной точки inf преобразования F. Ее следует рассматривать как неустойчивую, поскольку любая ее окрестность с разрезом Re(z)=2.5 в комплексной плоскости с ростом k "уходит" от нее, а из этого разреза все больше точек попадает в нее, и эти ее дискретные прообразы все плотнее (в пределе всюду плотно) заполняют этот разрез, но все их множество, поскольку оно счетно, имеет меру нуль. Отсюда следует, что на прямой Re(z)=2.5 почти всюду не существует теоретического предела итераций. Из-за ошибок округления, хотя они и малы, прообразы inf почти никогда не попадают в inf при вычислениях и потому ведут себя так же (т.е. хаотически), как и не прообразы.

Заметим теперь, что порядок действий в F можно переставить:

F(x)=x/2+1.25+0.25/(2x-5) и F(x)=(x²-6)/(2x-5)– это две другие математически эквивалентные записи F. Выполним строку

14;TF='xt/2+1.25+.25./(2*xt-5)';

а затем строки 5 и 7 и получим тот же результат, что и выше. Но строка

15;TF='(xt.^2-6)./(2*xt-5)';

и строки 5, 8 дадут только устойчивые неподвижные точки преобразования F и, как и раньше, три точки NaN, так что с такой TF мы не увидели бы рассмотренной выше картины. Выполним строку 6 с n=100 и увидим, как все происходит до конца.

Пример 3 показывает, что, пользуясь довольно простыми командами MATLAB'а и руководствуясь здравым смыслом, можно проанализировать достаточно тонкие математические детали задачи.