Так как теперь преобразование 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; fork=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.0033e-4; xt=2.5+(-1:hy:1)'*i; w=2:length(xt); n=100; fork=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'а и руководствуясь здравым смыслом, можно проанализировать достаточно тонкие математические детали задачи.
7. Системы линейных алгебраических уравнений
В MATLAB'e имеется несколько команд для решения таких задач. Рассмотрим здесь наиболее употребительную из них – оператор \ (обратный слэш или деление слева). Для чисел a\b=a-1b в отличие от того, что a/b=ab-1. Для невырожденной квадратной матрицы A и вектора-столбца f вектор-стобец u=A\f есть решение системы линейных уравнений Au=f, т.е. в записи по аналогии с числами u=A-1f, но матрица A-1 в действительности не вычисляется, а система Au=f решается методом исключения Гаусса (это значительно дешевле), с которым каждый из нас как-то знаком еще со школы. Правая часть f может быть и матрицей, и тогда каждый столбец матрицы u=A\f есть решение системы для соответствующего столбца из f, т.е. командой A\f можно сразу решить много систем, и поступать так всегда выгодно, ибо тогда матрица A будет разлагаться на множители по методу Гаусса только один раз. Если матрица A не квадратная, вектор u=A\f есть решение системы Au=f в смысле метода наименьших квадратов, но мы не будем здесь разбирать этот вариант оператора \ , поскольку не должны сколько-нибудь заметно вдаваться в детали численных методов. Таким оброазом, оператор \ – одна из самых мощных команд системы. А поскольку решение систем Au=f с квадратными матрицами – наиболее часто встречающаяся в вычислительной практике задача, необходимо иметь представление хотя бы об одном из методов ее реализации.
1.Решим и проанализируем на точность систему Au=f, выполнив следующую строку:
1;x=1:.1:5; m=length(x); A=toeplitz(exp(x)); ut=sin(x)'; f=A*ut; u=A\f;
Здесь задается вектор-строка x=1:.1:5 из 41 элемента, затем по вектору exp(x) командой toeplitz вычисляется квадратная матрица A размеров 41*41, далее в виде вектора-столбца задается точное решение ut, по нему строится правая часть f и, наконец, находится численное решение системы Au=f. Выясним вид A, построив графики
2;mesh(A) plot(A(1,:)) plot(A(20,:)) plot(A(41,:))
Теперь выполним команды
3;plot(ut) plot(f) plot([ut,u])
и увидим, что f сильно отличается от u, а точное и численное решения графически совпадают (изображающая u фиолетовая линия накрыла желтую линию ut). Иногда нужно сравнивать сильно разномасштабные кривые (у нас это u и f). Для такого сравнения следует провести нормировку каждой из них на свой абсолютный максимум. В нашем случае график
4;plot([u/max(abs(u)),f/max(abs(f))])
позволяет изобразить динамику изменения кривых u и f относительно друг друга.
2.Команда det вычисляет определитель. Найдем
1;max(abs(u-ut)) (=9.7167e-13) det(A) (=4.1058e-9)
откуда видно, что u и ut совпадают примерно в 12 знаках, хотя определитель det(A) системы на первый взгляд очень мал. Тем не менее попробуем решить нашу систему по правилу Крамера, обозначив такое решение через uc:
2;d=det(A); uc=zeros(m,1); for k=1:m,uc(k)=det([A(:,1:k-1),f,A(:,k+1:m)])/d; end, plot([ut,uc])
Здесь при k=1 A(:,1:k-1)=[], а при k=m A(:,k+1:m), поскольку 1:0=[] и m+1:m=[], так что все матрицы Ck=[A(:,1:k-1),f,A(:,k+1:m)], k=1:m, сформированы правильно. Из графика видно, что так найденное uc снова совпадает с точным решением ut. Теперь
3;max(abs(uc-ut)) (=9.6934e-13) т.е. погрешность практически не изменилась. Длительное время считали, что прималых d систему решать численно нельзя.
3.В действительности трудности численного решения системы связаны с возможной близостью к нулю собственных значений матрицы A, которые определяются как все решения l1, ..., lm характеристического уравнения det(A-lE)=0, где E=eye(m) – единичная (с нулями вне главной диагонали) матрица порядка m. За последние 10 лет были созданы достаточно хорошие программы для нахождения собственных значений. В MATLAB'е это делает команда eig, которую мы и применим сейчас:
1;sp=eig(A); plot(sp), [y,yi]=min(abs(sp)), [z,zi]=max(abs(sp)), c=z/y
Из графика видно, что все собственные значения вещественны (по оси абсцисс отложены значения индекса) и никак не упорядочены. Ближе всего к нулю (=0.136) третье из чисел lk, дальше всего (=898) – 41-е, а модуль их отношения c(A)ºc=6.6040e+3 называется числом обусловленности (condition number) матрицы A и отражает гораздо более глубокие ее свойства, чем величина ее определителя (которая по существу вообще ничего не отражает ни в практическом, ни в теоретическом плане). При обращении [V,D]=eig(A) в диагональной матрице D расположены lk (раньше они были в векторе-столбце sp), а в матрице V в виде столбцов – соответствующие им собственные векторы vk с единичной среднеквадратичной нормой