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=3s*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 равна 2n.
2.Теперь представим (2) в виде
x=F(x), где y=F(x)=5-6/x, так что F'(x)=6/x2 , F'(3)=2/3<1, |F'(x)|<1 при |x|>2.45.
Следовательно, устойчивая неподвижная точка - это x0=3. Получим резултат сразу для "всех" вещественных x0=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, для которой итерации опять идут без ошибок округления. Возникает предположение, что вся комплексная плоскость с выколотой точкой x1=2 стягивается итерациями в точку x2=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) довольно быстро стягиваются в точку x2=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)=x2-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))2 .
Так как f'(x1,2)0, можно использовать эти итерации по Ньютону (часто их называют методом касательной), а поскольку F'(x1,2)=0, следует ожидать их быстрой сходимости и того, что теперь оба решения x1,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
Предел итераций при x0=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 точек.
При x0=4 предел итераций равен 3, а v(k) успевают подойти к -1 (у v(6) точность уже потеряна). Таким образом, в зависимости от начального приближения x0=xt итерации могут сходиться к обоим решениям задачи.
Так как теперь преобразование y=F(x) уже не является дробно-линейным, рассмотрим трансформацию комплексной прямоугольной области в процессе итераций. Создадим такую область xt из 412=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 стягивается итерациями в точку x1=2, а полуплоскость Re(z)>2.5 - в точку x2=3: прямая Re(z)=2.5 переходит сама в себя. Точка x1=2 имеет красный цвет потому, что в нее перешли первые 25 столбцов матрицы xt, а остаток rem(25,7)=4, но 4-й цвет - красный. Далее идут зеленые точки (26-й столбец) и синяя x2=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 ,
и это есть точки
z1=2.5-0.5i, z2=2.5, z3=2.5+0.5i.
Теперь разберемся, как преобразуется итерациями F прямая Re(z)=2.5. Так как прямая Re(z)=2.5 переходит в себя, ее можно параметризовать с помощью переменной y=Im(z), и пусть yk=Im(Fk(Re(z)=2.5), k=1, 2,... , т.е. после k итераций y переходит в yk(y). Ясно, что на k-й итерации в inf перейдут только те точки из всего множества yk-1, которые равны нулю. Поскольку -inf<y<inf, на 1-й итерации это случится только с одной точкой z1=2.5 (для нее y=0). Выдадим на график y1 после 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
и увидим, что функция y1(y) пересекает ось абсцисс только два раза (при этом y1 дважды непрерывно пробегает от -inf до inf), так что на 2-й итерации в inf перейдут только две точки (это уже известные нам z1=2.5-0.5i, и z3=2.5+0.5i). Выполним строку 10 с n=2 и снимем с графика строкой
Страницы: 1, 2, 3, 4, 5, 6, 7
