Материалы представленные на сайте исключительно с целью ознакомления пользователям Интернета и не преследуют коммерческих целей или нарушение авторских прав.
© 2018 - 2024
1. Простой пример: найти корни уравнения x*sin(x^2)=0 на отрезке [0,3]. Программа:
1; x =0:.01:3; f = x .* sin ( x .^2); plot ( x ,[ f ;0* f ] ), grid
Ginput
В команде ginput точка снимается нажатием левой клавиши мыши, Enter – выход из ginput.
Проверим это другим способом:
3; nx = length ( x ); w =1: nx -1; x ( find ( f ( w ).* f ( w +1)<0| f ( w )==0)) Отв: 0, 1.77, 2.5.
Эту строку можно упростить:
4; nx = length ( x ); w =1: nx -1; x ( f ( w ).* f ( w +1)<0| f ( w )==0)
Матрицы и векторы с элементами 0-1.
2. Сложный пример – неявные функции. Построим график неявной функции f(x,y)ºx3y-2xy2+y-0.2=0, x,y=[0, 1]. Это выполнит программа
1; h =.02; x =0: h :1; [ X , Y ]= meshgrid ( x ); f = X .^3.* Y -2* X .* Y .^2+ Y -.2;
2;v=[0,0]; contour(x,x,f,v), grid
На графике зеленая линия (справа она двузначная) представляет искомый результат. Область в первом квадранте между этими кривыми обозначим через G. Эту задачу совсем непросто сделать в других системах программирования прежде всего потому, что вычисление образующих линии уровней точек – в общем случае очень сложная процедура.
Выясним, какой знак имеет f в области G, для чего выполним
3; mesh ( x , x , f .*( f >0))
Это пример трехмерной, т.е. xyz-графики. В ней цвет используется для изображения амплитуды (значения z),
изменяясь с ростом z от темносинего через голубой, зеленый и желтый до темнокрасного.
Вычислим площадь S этой области:
4; S = h ^2* sum ( f (:)>=0) (S=0.7296).
Для h=0.01 выполним строку 1, затем строку 4 и получим S=0.7204, а для h=0.005 найдем S=0.7152. При интегрировании всегда естественно делать такие проверки.
Выясним, какой объем заключен между поверхностью f(x,y) и областью G, где f(x,y)>=0. Для этого снова возьмем в строке 1 h=0.02 и вычислим
5; V = h ^2* sum ( f ( f >=0)) (V=0.1268)
Для h=0.01 V=0.1235, а для h=0.005 V=0.1219. Теперь не нужно писать f(:), поскольку f(f>=0) есть вектор.
Конечно, эти результаты приближенные (с точностью до 1 - 2%), но отметьте, как быстро и просто они были получены. Такие приемы можно применять для решения достаточно широкого круга задач.
Выполним строку
6; C = contour ( x , x , f ); clabel ( C )
которая зашлет числовую информацию о графике в матрицу C и построит график, выбрав значения уровней автоматически. Из матрицы C можно последовательно выбирать все кривые.
Обобщения. Графическим способом можно решать системы уравнений и уравнения в комплексной плоскости. Команда contour3 строит линии уровней для функций f(x,y,z), при этом сетки по аргументам всегда должны быть прямоугольными.
Полиномы
По степени применимости, по разнообразию и качеству соответствующих команд скалярные полиномы – следующие за матрицами математические объекты в MATLAB'е. Полином
p(x)=anxn+an-1xn-1+...+a0 задается вектором-строкой p из чисел an, an-1, ... , a0,
т.е. коэффициентами, расположенными в порядке убывания показателя степени. Его степень n задавать не надо, поскольку n=length(p)-1; полином может быть и константой – тогда n=0; коэффициенты ak – любые комплексные числа. Вектор p интерпретируется системой как полином только тогда, когда он задается в качестве параметра для одной из команд, производящих вычисления с полиномами. Так как в этих командах не проверяется условие an¹0, надо стараться самим соблюдать его, поскольку иногда это может служить источником ошибок.
Основные команды для действий с полиномами таковы:
conv(p,q) – произведение полиномов p и q. Название команды происходит от слова convolution (свертка), поскольку коэффициенты произведения действительно получаются как компоненты свертки векторов p и q.
[q,r]=deconv(b,a) – частное (q) и остаток (r) от деления b на a, так что conv(a,q)+r=b.
residue(b,a) – разложение рациональной функции b(x)/a(x) на элементарные дроби над полем комплексных чисел с выделением целой части. Если a(x) имеет кратные или близкие друг к другу корни, результаты могут быть неверными, поскольку такая задача плохо обусловлена. Плохая обусловленность, т.е. крайне сильная зависимость результата от коэффициентов, иллюстрируется заключительным примером из этой темы.
p=poly(r) – построение полинома по корням, заданным в векторе-столбце r. Для квадратной матрицы r полином p будет ее характеристическим многочленом.
polyval(p,x) – поэлементное вычисление значений полинома p на множестве x, где x может быть как вектором, так и матрицей. Размеры результата совпадают с size(x).
polyder(p) – производная от p.
roots(p) - вектор-столбец, содержащий все корни полинома. Их порядок не определен. Сказанное по поводу неустойчивости результатов для команды residue точно так же относится и к команде roots. Корни полинома вычисляются как собственные значения некоторой матрицы A(p) того же порядка.
Приведем несколько примеров по применению этих команд.
1. Построим график полинома p(x)=x3-x+2 на отрезке -1<=x<=1. Это выглядит так:
p =[1,0,-1,2]; x =-1:.01:1; f = polyval ( p , x ); plot ( x , f ), grid
Полином не имеет корней на заданном отрезке. Это подтверждает и команда
roots ( p )'
которая даст
ans = -1.5214 0.7607 - 0.8579i 0.7607 + 0.8579i.
2.Разделим предыдущий полином на x-3:
[q,r]=deconv(p,[1,-3])
Тогда
q = 1 3 8, r = 0 0 0 26.
Другими словами, частное q(x)=x2+3x+8, а остаток r=26.
3. Разложим функцию (x-3)/p(x) на элементарные дроби:
1;[r,s,k]=residue([1,-3],p); r', s', k'
Для r' получим вектор из трех компонент r1, r2, r3 :
-0.7607 0.3803 - 0.4289i 0.3803 + 0.4289i,
для s' - также вектор из трех компонент s1, s2, s3 :
-1.5214 0.7607 - 0.8579i 0.7607 + 0.8579i
и k=[] (это означает, что целой части в разложении нет – действительно, у числителя первая, а у знааменателя третья степени). Компоненты векторов r и s означают, что
(x-3)/p(x)=sum(ri/(x-si), i=1:3.
Команда residue работает и в обратную сторону:
2;[q,p]=residue(r,s,k)
восстановит исходные числитель и знаменатель:
q = 0 1 -3 (он получился точно),
p = 1.0000 -0.0000 -1.0000 2.0000 (здесь уже сказались ошибки округления и старший коэффициент не равен 1 автоматически).
4. В заключение приведем сложный пример (Уилкинсон, 1963), показывающий, что иногда, несмотря на хорошую разделенность корней полинома, их вычисленные значения могут очень сильно зависеть от значений некоторых коэффициентов просто потому, что производные корней по этим коэффициентам – очень большие по модулю числа. Такие задачи называются плохо обусловленными и всегда требуют повышенного внимания независимо от того, каким методом они решаются. В то же время они в наибольшей степени стимулируют теоретические исследования по оценке точности машинных вычислений, и пример Уилкинсона – одна из первых классических задач такого рода.
Перейдем теперь к описанию примера. Пусть vn=1:n, где n>1 - целочисленный параметр, и pn=poly(v') - полином с корнями 1:n, которые хорошо отделены друг от друга, а wn=roots(pn) - вектор-столбец с вычисленными корнями полинома pn. Проведем сравнение vn' и wn для различных n. Начнем с n=2:
1;n=2; vn=1:n; pn=poly(vn'); wn=roots(pn); [vn',wn]
и получим ans=1 2 2 1
откуда видно, что элементы в wn нужно упорядочить. Выполняя при n=3 отредактированную строку 1
2;n=3; vn=1:n; pn=poly(vn'); wn=roots(pn); R=[vn',sort(wn)]
найдем R = 1.0000 1.0000
2.0000 2.0000
3.0000 3.0000 .
Появившиеся в первом столбце R цифры 0 как бы "наведены" значениями из второго столбца, и таких мелких шероховатостей у команды format немало. А цифры 0 во втором столбце R говорят о том, что уже появилась погрешность в определении корней. Она, конечно, еще очень мала:
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. Сравним точные и вычисленные корни графически: график
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, и из графика. Если построить график
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 должно в обязательном порядке заметно сказаться на значениях некоторых корней, каким бы методом они ни находились. Более того, если эти необходимые изменения никак не проявились, метод следует забраковать.
Итерации
Итерации являются с точки зрения программирования одним из самых эффективных способов организации вычислений. Простые итерации в общем случае представляются в виде
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)
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'а полным, выдадим
Imag ( xt ')
Образы первой и последней точки начальной окружности слегка отличаются от z=2, а ее 21-я точка (она соответствует fi=0) есть z=3. Это получилось потому, что sin(pi) и sin(-pi) не равны нулю в точности. Снова сделаем окно MATLAB'а небольшим и выполним строку 4 с радиусом 3.01, а затем строки
Sum(isnan(xt))
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 и снимем с графика строкой
11; q = ginput ; q (:,1)'
четыре точки пересечения кривой y2(y) с осью абсцисс: приближенно это -1.2079 -0.2056 0.2055 1.2079,
тогда как для y1(y) это были значения -0.5 и 0.5. Каждая переходящая в inf на k-й итерации точка порождает слева и справа от себя две такие точки, которые перейдут в inf на (k+1)-й итерации. Поэтому с ростом k точки yk(y)=0, с одной стороны, уходят в обе стороны от нуля, а с другой - все чаще появляются в тех местах, где они однажды уже появились. Всего на k-й итерации в inf перейдет 2k-1 точек. В действительности с ростом k они всюду плотно заполняют прямую Re(z)=2.5, а между ними yk(y) непрерывно пробегает все значения от -inf до inf. Чтобы быть в этом уверенными, выполним последнюю программу
12; hy =1.0033 e -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) у функции y100(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. Ошибки при последовательном вычислении Fk(z) на прямой Re(z)=2.5 будут ничтожными, так что значения yk(y) вычисляются с высокой точностью, но от шага hy число найденных перемен знака в yk(y), конечно, зависит: при k=100 число всех нулей функции yk(y) равно 2100-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)=(x2-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'а и руководствуясь здравым смыслом, можно проанализировать достаточно тонкие математические детали задачи.
Дата: 2019-05-28, просмотров: 138.