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'; fork=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;'; fork=1:n,eval(TF),x=[x,xt]; end, plot(x), grid
произведет те же вычисления, но обратите внимание на то, как в ней записаны TF и eval.
Хотя внешне сходимость x(k) к x1=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) задерживаются в районе неподвижной точки x2=3 (они все время уходят от нее, но на графике этого не видно), а затем примерно за 60 итераций монотонно движутся к x1=2.
(b') xt=3 . Все x(k)=3, а график из строки 3 пуст. Это произошло потому, что при xt=3 F(xt)=15/5=3 получается без ошибок округления.
(b'') xt=3.01 . Пределами для x(k) и v(k) будет inf, так что x2=3 является неустойчивой неподвижной точкой преобразования F: при малейшем сдвиге x0 с x2 в пределе итераций получится либо устойчивая неподвижная точка x1, либо 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 снова получим x2 за одну итерацию, а при xt= -3.01 и далее получим inf. Таким образом, при вещественном x0 итерации сходятся к устойчивой неподвижной точке x1=2 при -3<x0<3 и к inf при |x0|>3. Просчитаем тот случай, когда |x0|=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)), holdon
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 итерации могут сходиться к обоим решениям задачи.