В 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, построив графики

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=[], так что все матрицы C_k=[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, которые определяются как все решения l₁, ..., l_m характеристического уравнения 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) третье из чисел l_k, дальше всего (=898) – 41-е, а модуль их отношения c(A)ºc=6.6040e+3 называется числом обусловленности (condition number) матрицы A и отражает гораздо более глубокие ее свойства, чем величина ее определителя (которая по существу вообще ничего не отражает ни в практическом, ни в теоретическом плане). При обращении [V,D]=eig(A) в диагональной матрице D расположены l_k (раньше они были в векторе-столбце sp), а в матрице V в виде столбцов – соответствующие им собственные векторы v_k с единичной среднеквадратичной нормой

sum ( v_k ² ))^1/2=1, k =1: m .

Команда eig (и целый ряд основанных на ней) имеет неплохую точность при решении спектральных задач средней сложности и является очень информативной. В нашем примере есть относительно близкие друг к другу l_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=f ¹ 0, A*du=df ¹ 0, re(f)=max(abs(df))/max(abs(f)) .

Тогда u¹0 и du¹0 ввиду невырожденности 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). Пусть ошибка