Matlab

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=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 и снимем с графика строкой

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; for k=1:n,xt=eval(TF); end

13;pzn=sum(sign(xt(w-1)).*sign(xt(w))<0)

которая зафиксирует 674 перемены знака (pzn) у функции y100(y) для -1y1 с hy=1.0033e-4 (из-за недостаточной малости hy далеко не все они учтены). Для -10y-8 с hy=1.0033e-4 pzn=675. Для симметричного относительно точки y=0 отрезка 8y10 с тем же 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, которые определяются как все решения 1, ..., m характеристического уравнения det(A-E)=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) третье из чисел k, дальше всего (=898) - 41-е, а модуль их отношения c(A)c=6.6040e+3 называется числом обусловленности (condition number) матрицы A и отражает гораздо более глубокие ее свойства, чем величина ее определителя (которая по существу вообще ничего не отражает ни в практическом, ни в теоретическом плане). При обращении [V,D]=eig(A) в диагональной матрице D расположены k (раньше они были в векторе-столбце sp), а в матрице V в виде столбцов - соответствующие им собственные векторы vk с единичной среднеквадратичной нормой

sum(vk2))1/2=1, k=1:m.

Команда eig (и целый ряд основанных на ней) имеет неплохую точность при решении спектральных задач средней сложности и является очень информативной. В нашем примере есть относительно близкие друг к другу k. Чтобы в этом убедиться, выполним строку

2;ssp=sort(sp); v=1:m-1; di=diff(ssp); min(abs([di./ssp(v);di./ssp(v+1)]))

(нормировка разностей делается на каждый член пары) и получим 0.0044, т.е. есть такая пара, у которой совпадает более двух знаков, так что наш пример в спектральном плане не совсем тривиальный. Теперь сделаем матрицу A вырожденной путем дублирования ее строк и посмотрим, как это отразится на результатах eig:

3;r=zeros(1,m-1); for k=1:m-1, v=[1:k,k,k+2:m]; C=A(v,:); r(k)=min(abs(eig(C))); end, plot(r)

Результаты следует признать неплохими - max(r)=1e-13 и есть даже несколько чистых нулей.

Чтобы получить представление о собственных векторах преобразования A, выполним строку

4;[V,D]=eig(A); D=V'*V; D(1:m+1:m^2)=0; mcv=max(abs(D(:)))

и получим mcv=7.5460e-16. Это означает, что собственные векторы с высокой степенью точности ортогональны, как и должно быть для симметричной матрицы A.

4.Теоретическая оценка погрешности при решении линейных систем. Объясним смысл числа обусловленности c(A), который имеет фундаментальное значение для вычислительных методов линейной алгебры. Пусть

A*u=f0, A*du=df0, re(f)=max(abs(df))/max(abs(f)) .

Тогда u0 и du0 ввиду невырожденности A. Будем рассматривать df как ошибку в f (ввиду линейности задачи это не приводит к потере общности - считать вектор df малым не нужно). Оценим величину

re(u)=max(abs(du))/max(abs(u)) для всех допустимых f и df,

т.е. возникающую при решении нашей системы относительную ошибку, которая является весьма общей характеристикой матрицы A. Если y, z и c - числа, определенные в строке

1;sp=eig(A); plot(sp), [y,yi]=min(abs(sp)), [z,zi]=max(abs(sp)), c=z/y

(это строка 1 из предыдущего примера), то

max(abs(du))z -1max(abs(df)),

поскольку du=A-1u и z -1 - максимальное по модулю собственное значение для A-1, причем равенство для max(abs(du)) теоретически возможно. Аналогично

max(abs(u))y-1max(abs(f)),

поскольку y -1 - минимальное по модулю собственное значение для A-1, и знак равенства здесь также возможен. Объединяя оценки для числителя и знаменателя re(u), будем иметь

re(u) z -1/y -1(max(abs(df))/max(abs(f))) или re(u) (y/z) re(f), т.е. re(u) c*re(f),

причем равенство обязательно достигается хотя бы для одной пары f и df. Другими словами, число обусловленности c(A) есть точная верхняя граница роста относительной ошибки при решении нашей системы: если относительная ошибка правой части равна re(f), то относительная ошибка решения re(u) не превзойдет ее более чем в c(A) раз.

Ввиду важности этого результата полезно понимать его и на пальцах: df переходит в du с возрастанием компонент не более чем в z -1 раз, а f переходит в u с уменьшением компонент не более чем в y -1 раз, так что

re(u) (z -1/y -1) re(f).

Мы установили смысл числа c(A), исходя из максимум-нормы векторов f, когда norm(f)=max(abs(f)), но более точное рассмотрение показывает, что это верно и для среднеквадратичной нормы

norm(f)=(sum(abs(f.^2))^0.5.

Команда cond(A) вычисляет значение c(A) для квадратной невырожденной матрицы A, но в cond спектр sp ищется не для A, а для B=A'A, ибо тогда вся спектральная задача для B становится заведомо вещественной и, более того, все ее собственные векторы взаимно ортогональны (последнее обстоятельство существенно уменьшает среднеквадратичную ошибку при этих сложных вычислениях). Для получения наших y и z нужно лишь извлечь квадратный корень из y(B) и z(B) - это также делается в команде cond(A). Для нашего примера c(A)=

2;cond(A) (=6.6040e+3)

и совпадает со значением c из строки 1, поскольку A - симметричная матрица.

Несмотря на внешнюю простоту, понятие числа обусловленности было четко сформулировано только в середине 1960-х гг., т.е. примерно через 15 лет после появления первых ЭВМ, а широко использоваться стало еще позже. Для вычислительных методов оно оказалось важнее понятия линейной зависимости, которое теперь некоторым образом выражается через него, но мы уже не будем здесь этого уточнять.

5.Практическая оценка погрешности. Число обусловленности может быть не всегда подходящим для оценки погрешности. Это так при большом n - тогда решение спектральной задачи в команде cond будет слишком дорогим. Это так и в том случае, когда ошибка df специфически неравномерна и имеет характерный профиль - тогда желательно иметь и профиль для поточечных ошибок du(k), чего, конечно, никак не может дать cond(A). В таких случаях приходится прибегать к непосредственному моделированию ошибок, которое состоит в задании некоторого множества случайных возмущений правой части и решении всех таких систем с последующим поточечным описанием границ получившихся решений. Проиллюстрируем этот способ оценки ошибки на нашем примере.

Восстановим этот пример:

1;hx=.1; x=1:hx:5; m=length(x); A=toeplitz(exp(x)); ut=sin(x)'; f=A*ut; u=A\f;

(здесь введен шаг hx). Пусть ошибка

df(x) в f(x) есть равномерно распределенная случайная величина,

по модулю не превосходящая значения g(x)=1e-4*f(t)dt в пределах от 1 до x.

В строке

2;g=1e-4*hx*abs(cumsum(f)); rand('state',0); n=30; V=(1:m)'*ones(1,n); plot([f/max(abs(f)),g/max(g)]), grid

задается профиль g=g(x) максимально допустимой ошибки в точке x, приводится в исходное состояние счетчик случайных чисел rand, задается число n=30 возмущений правой части f и матрица V размеров m*n из продублированного n раз столбца (1:m)'. График показывает, что самые большие ошибки там, где abs(f) мало. В строке

3;F=f(V)+g(V).*(2*rand(m,n)-1); U=A\F;plot(U), pause, plot([g(V).*(2*rand(m,n)-1),g])

создается матрица F возмущенных правых частей из m строк и n столбцов путем прибавления к размноженному n раз вектору f возмущений в нужных границах, решаются все системы AU=F и выводятся на график все решения U и все возмущения вместе с их границей. Из последнего графика видно, что возмущения заданы правильно. В строке

4;ua=max(U,[],2); ui=min(U,[],2); plot([ui,ut,ua])

путем обработки U вдоль строк находятся поточечные границы решений (ua - верхняя, ui - нижняя) и строится график точного решения и этих границ.

Результат кажется нам плохим потому, что на графике

5;plot(F)

возмущений просто не видно (масштаб f забил их), а в U они отразились слишком сильно. Но формально он не противоречит оценке, связанной с числом обусловленности c(A). Действительно, без учета профиля df максимальная ошибка me вычисляется как

6;sp=eig(A); me=min(abs(sp))^(-1)*max(g)

и равна 0.6064, тогда как

7;max([ua-ut;ut-ui]) (=0.5180)

и лишь немного подрастет с увеличением n (при n=100 это 0.5540), т.е. не превзойдет значения me, что и свидетельствует о формальной согласованности обоих методов оценки погрешности. Чтобы окончательно избавиться от ощущения какой-то несогласованности наших методов, применим критерий cond(A) для среднеквадратичной нормы. Для этогог нам нужно пройтись по всем столбцам k=1:n (сейчас n=30) матриц U и F и вычислить

max( (norm(du)/norm(u)) / (norm(df)/norm(f)) ) .

Запишем это в виде (максимально используйте карман при наборе строк)

8;max((sum((U-ut(V)).^2).^.5./sum(U.^2).^.5)./(sum((F-f(V)).^2).^.5./sum(F.^2).^.5))

и получим 3.6484e+3, что меньше c(A)=6.6040e+3, так что и здесь мы не обнаружили противоречия.

Если подойти к оценке погрешности упрощенно, построив график

9;gm=max(g); plot(A\[f-gm,f,f+gm])

то все три линии на нем практически совпадут. Таким способом можно моделировать ошибку, если преобразование A-1 монотонно, т.е. если при f1<f2 обязательно u1<u2 или, наоборот, обязательно u1>u2. Но у нас это не так: на графике (эта строка получается из строки 9)

10;gm=100*max(g);plot(A\[f-gm,f,f+gm])

желтая линия (она соответствует решению для правой части f-gm<f) то выше, то ниже фиолетовой. Именно из-за немонотонности преобразования A-1 получается такой заметный разброс в U. С помощью этого приема можно быстро выяснить монотонность преобразования y=Q(x) (не обязательно линейного), что далеко не всегда удается определить теоретически: если все отклики Y=Q(X), x-a<X<x+a, a>0, лежат между Q(x-a) и Q(x+a), то преобразование Q монотонно, и нужно лишь взять значение a таким, чтобы результат был виден на графике (для этого нам пришлось увеличить max(g) в 100 раз).

6.Посмотрим, как изменятся результаты нашего примера при увеличении m - порядка матрицы A. Выполним, не меняя смысла задачи, отредактированную строку 1 предыдущего примера

1;clear all, hx=.01; x=1:hx:5; A=toeplitz(exp(x)); ut=sin(x)'; f=A*ut;u=A\f; plot([ut,u]), c=cond(A)

Здесь шаг hx уменьшен в 10 раз, так что теперь порядок m=401 - довольно высокий; c(A)= 6.0804e5 возросло почти в 100 раз, т.е. обусловленность A заметно ухудшилась (примерно в 102 раз), но

2;max(abs(u-ut)) (=1.3153e-9)

еще достаточно мал, хотя и возрос примерно в 103 раз, т.е. больше, чем c(A). Такое расхождение с теорией как бы предупреждает о том, что даже при сохранении смысла задачи увеличение ее размерности не позволяет автоматически применять критерий числа обусловленности к оценке ошибок округления. К выбору числа m нужно всегда относиться с повышенным вниманием.

Чтобы получить представление о собственных векторах преобразования A, выполним строку

3;[V,D]=eig(A); D=V'*V; m=length(x); D(1:m+1:m^2)=0; mcv=max(abs(D(:)))

и получим mcv=2.2985e-15, т.е. степень ортогональности остается удивительно высокой. Жордановы клетки порядка выше первого могут быть тогда, когда mcv(A)>0.99.

Мы рассмотрели этот пример так подробно, чтобы показать исключительно высокие возможности MATLAB'а в том, что касается анализа результатов.

Заключение

MATLAB - высокоуровневая система программирования, позволяющая резко сократить затраты труда при проверке алгоритмов и проведении прикидочных расчетов. Возможность проведения больших расчетов на MATLAB'е определяется в основном теми затратами времени, на которые может пойти пользователь: здесь приходится выбирать между легкостью и наглядностью программирования и представления результатов, с одной стороны, и затратами времени на счет - с другой. Система очень удобна для освоения и апробации численных методов, что мы и хотим показать здесь прежде всего. Именно поэтому она рекомендуется как одна из основных для физиков и многих других естественно-научных специальностей в ведущих американских университетах. Детальное освоение любой большой программной системы - это достаточно длительный процесс, основу которого составляют индивидуальная работа, и наши занятия призваны дать лишь первоначальный импульс этому процессу в отношении MATLAB'а. Темы 2 - 4 представляют сравнительно элементарное введение, а в остальных рассматриваются более сложные примеры, показывающие, как можно использовать программные и графические возможности системы для исследования численных алгоритмов.

Страницы: 1, 2