2.0000 2.0000
3.0000 3.0000 .
Появившиеся в первом столбце R цифры 0 как бы "наведены" значениями из второго столбца, и таких мелких шероховатостей у команды format немало. А цифры 0 во втором столбце R говорят о том, что уже появилась погрешность в определении корней. Она, конечно, еще очень мала:
3;((R(:,2)-R(:,1))./R(:,1))'
дает относительную ошибку для корней
ans = 1e-14 *( 0.1110 -0.0444 -0.1184)
- это означает, что в численном результате верны примерно 14 знаков.
Для n=10 выполнение строки (ее можно создать из строк 2 и 3)
4;n=10; vn=1:n; pn=poly(vn'); wn=roots(pn); R=[vn',sort(wn)]; R1=(R(:,2)-R(:,1))./R(:,1)
говорит о том, что точность результата постепенно падает. Строка
5;me=max(abs(R1))
дает me=4.2009e-10, т.е. теперь для всех корней верны только 9 знаков (me - max. error). Но корни еще остаются вещественными, поскольку
6;iwn=sum(abs(imag(wn)))
дает iwn=0.
Выполним строку 4 для n=20 и затем строкой 5 найдем максимальную относительную ошибку me=0.0073, так что теперь для всех корней верны только 2 знака, но результат пока еще остается вещественным, поскольку после строки 6 получим iwn=0. Сравним точные и вычисленные корни графически: график
7;plot(R), grid
показывает, что погрешность для некоторых корней уже видна на глаз - для корней 10:17 желтая и фиолетовая линии слегка различаются.
Теперь приступим к самому интересному в данном примере. При n=20 pn(2)= -210 (это коэффициент при x19). Прибавим к нему 1e-7, т.е. внесем в него малое возмущение примерно в 10-м знаке, и повторим расчет с отредактированной строкой 4:
8;n=20; vn=1:n; pn=poly(vn'); pn(2)=pn(2)+1e-7; wn=roots(pn); R=[vn',wn], plot(R), grid
Несмотря на такое малое возмущение в коэффициенте pn(2), некоторые корни стали комплексными. Это видно, во-первых, из выдачи на экране (их мнимые части достигают по модулю 2.7), во-вторых, из того, что теперь строкой 6 получим iwn=18.67, и из графика. Если построить график
9;plot(R,'.'), grid
т.е. убрать соединения между точками, результат будет выглядеть более рельефно. На таких графиках режим zoom работает.
Выясним, почему так сильно изменились результаты при внесении столь малого возмущения. Обозначим через p(x) наш невозмущенный полином pn при n=20 и через a его второй коэффициент:
p(x)=prod(x-k),k=1:20, или p(x)=x20+ax19+ ... +20! , a=-210.
Тогда для корней x=1:20 по теореме о производной неявной функции будем иметь
p/x*dx/da + p/a = 0,
откуда
dx/da = -p/a / p/x.
У нас p/a =x19, а полином p/x находится как polyder(pn) (см. начало этой темы). Поэтому для вычисления dxda на множестве vn наших корней сначала выполним отредактированную строку 4 с n=20
10;n=20; vn=1:n; pn=poly(vn');
а затем строку (на графике значения функции определены только для x=1:20)
11;dpn=polyder(pn); dxda=-(vn.^19)./polyval(dpn,vn); plot(log10(abs(dxda)),'.'), grid
и увидим, что уже при x=8 dx/da=105 и будет еще больше с ростом x. Поэтому внесение в коэффициент a возмущения 10-7 должно в обязательном порядке заметно сказаться на значениях некоторых корней, каким бы методом они ни находились. Более того, если эти необходимые изменения никак не проявились, метод следует забраковать.
6. Итерации
Итерации являются с точки зрения программирования одним из самых эффективных способов организации вычислений. Простые итерации в общем случае представляются в виде
xk+1=F(xk), k=0, 1, 2, ... , x0 задано,
где F - заданное преобразование y=F(x), x0 - как-то выбранное начальное приближение, xk - значение переменной x на k-й итерации, а сама переменная x может быть любой - числом, вектором или матрицей. Предел итераций, если он существует, будем обозначать через X, и тогда уравнение
X=F(X) (1)
должно обязательно выполняться для итерационного преобразования F. Это обстоятельство помогает выбирать различные варианты для F. Все решения X этого уравнения называются неподвижными точками преобразования F(x). Все такие x0, для которых последовательность xk сходится, образуют область сходимости итерационного преобразования F. Cкорость сходимости итераций xk характеризуется поведением числовой последовательности
vk=norm(xk+1-xk)/norm(xk-xk-1),
где норма может выбираться по-разному. Для сходимости xk уже не обязательно, чтобы существовал предел V последовательности vk, но очень часто для сходящихся итераций он существует, и если это так, то пусть a=abs(V). Тогда при a=1 сходимость xk будет крайне медленной, при 0<a<1 это будет сходимость по геометрической прогрессии со знаменателем q=V, а при a=0 последовательность xk будет сходиться быстрее любой геометрической прогрессии. Покажем теперь на простейших примерах, как все это выглядит на деле. Чтобы иметь возможность разнообразия вариантов, задачу возьмем нелинейной. Рвссмотрим уравнение
(x-2)(x-3)=0 или f(x)x2-5x+6=0 c корнями x1=2, x2=3, (2)
построим для нахождения его решений x1=2, x2=3 несколько итерационных преобразований или схем F и проанализируем их работу.
1.Пусть для уравнения (2)
x=F(x), где y=F(x)=(x2+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) к 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)
Страницы: 1, 2, 3, 4, 5, 6, 7
