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, 3, 4, 5, 6, 7
