Реферат: СИНГУЛЯРНОЕ РАЗЛОЖЕНИЕ В ЛИНЕЙНОЙ ЗАДАЧЕ МЕТОДА НАИМЕНЬШИХ КВАДРАТОВ

МИНИСТЕРСТВО ОБРАЗОВАНИЯ РОССИЙСКОЙФЕДЕРАЦИИ

Математическийфакультет

Кафедраприкладной математики

ДИПЛОМНЫЙ ПРОЕКТ

сингулярное разложение в линейной задаче методанаименьших квадратов

Заведующий кафедройприкладной

математики          

Исполнил:                                                           

Научный руководитель

          

 

Владикавказ 2002

/>СОДЕРЖАНИЕ

ВВЕДЕНИЕ… 3

Глава 1. Метод наименьших квадратов… 7

1.1. Задача наименьших квадратов… 7

1.2. Ортогональное вращение Гивенса… 9

1.3. Ортогональное преобразование Хаусхолдера… 10

1.4. Сингулярное разложение матриц… 11

1.5. QR–разложение… 15

1.6. Число обусловленности… 20

глава 2. Реализация сингулярного разложения… 25

2.1. Алгоритмы… 25

2.2. Реализация разложения… 27

2.3. Пример сингулярного разложения… 29

глава 3. Использование сингулярного разложения  в методенаименьших квадратов              33

ЗАКЛЮЧЕНИЕ… 38

ЛИТЕРАТУРА… 39

ПРИЛОЖЕНИЕ 1. Исходные тексты программы… 40

ПРИЛОЖЕНИЕ 2. контрольный пример… 45

/>/>/>ВВЕДЕНИЕ

Метод наименьших квадратов обычно используется как составнаячасть некоторой более общей проблемы. Например, при необходимости проведенияаппроксимации наиболее часто употребляется именно метод наименьших квадратов.На этом подходе основаны: регрессионный анализ в статистике, оцениваниепараметров в технике и т.д.

Большое количество реальных задач сводится к линейной задаченаименьших квадратов, которую можно сформулировать следующим образом.

Пусть даны действительная m´n–матрица A ранга k£min(m,n) и действительный m–векторb. Найти действительный n–вектор x0,минимизирующий евклидову длину вектора невязки Ax–b.

Пусть ynмерный вектор фактических значений, x nмерныйвектор значенийнезависимой переменной, b – коэффициенты в аппроксимации yлинейной комбинацией n заданных базисных функций j:

/>.

Задача состоит в том, чтобы в уравнении подобрать такие b,чтобы минимизировать суммы квадратов отклонений e=y–Xb, где X –есть так называемая матрица плана, в которой строкамиявляются nмерныйвектора с компонентами, зависящими от xj:/> каждая строка соответствуетопределенному значению xj.Коэффициенты можно найти решая нормальные уравнения />, откуда />. Покажем это. Возведем вквадрат выражение для е:

/>

т. к.   />.

Это выражение имеет экстремум в точке, где />=0

Откуда и получаем />.

Следует отметить, что последнее выражение имеет вопределенной степени формальный характер, т. к. решение нормальных уравнений,как правило, проводится без вычисления обратной матрицы (метод Крамера) такимиметодами как метод Гаусса, Холесского и т. д.

Пример. Пусть заданы результатычетырех измерений (рис. 1): y=0приx=0; y=1 приx=1; y=2приx=3; y=5приx=4. Задача заключается втом, чтобы провести через эти точки прямую /> такимобразом, чтобы сумма квадратов отклонений была минимальна. Запишем уравнение,описывающее проведение прямой /> порезультатам измерений. Мы получаем переопределенную систему:

/>

или Xb=y. Нам понадобится матрица XTXи обратная к ней:

/>

Тогда решение b=(XTX)-1XTyпо методу наименьших квадратов будет иметь вид />

Таким образом, оптимальная прямая задается уравнением /> Метод точечной квадратичнойаппроксимации (метод наименьших квадратов) не предполагает, что мы должныприближать экспериментальные данные лишь с помощью прямых линий. Во многихэкспериментах связи могут быть нелинейными, и было бы глупо искать для этихзадач линейные соотношения. Пусть, например, мы работаем с радиоактивнымматериалом. Тогда выходными данными у являются показания счетчикаГейгера в различные моменты времени t. Пусть наш материал представляетсобой смесь двух радиоактивных веществ, и мы знаем период полураспада каждогоиз них, но не знаем, в каких пропорциях эти вещества смешаны. Если обозначитьих количества через С и D, то показания счетчика будут вести себяподобно сумме двух экспонент, а не как прямая:

/>.                                                   (1)

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

/>

Рис. 1. Аппроксимация прямой линией.

соответствовать (1). Вместо этого мы имеем серию показаний счетчика/> в различные моменты времени/>, и (1) выполняется лишьприближенно:

/>

Если мы имеем более двух показаний, m>2, то точноразрешить эту систему относительно C и D практически невозможно.Но мы в состоянии получить приближенное решение в смысле минимальных квадратов.

Ситуация будет совершенно иной, если нам известны количествавеществ C и D и нужно отыскать коэффициенты l и m.Это нелинейная задача наименьших квадратов, и решить еесущественно труднее. Мы по–прежнему будем минимизировать сумму квадратовошибок, но сейчас она уже не будет многочленом второй степени относительно l и m,так что приравнивание нулю производной не будет давать линейных уравнений дляотыскания оптимальных решений.

/>/>/>Глава 1. Метод наименьшихквадратов/>1.1. Задача наименьшихквадратов

Задача наименьших квадратов заключается в минимизация евклидовой длины вектора невязок || Ax-b ||.

Теорема 1. Пусть Аm´nматрица ранга k, представленная в виде 

A=HRKT                                                     (2)

где H ортогональная m´mматрица; Rm´nматрица вида

/>,                                              (3)

где: R11  –  kxkматрица ранга k; K ортогональная  kxkматрица. Определим вектор

/>                                          (4)

и введем новую переменную

/>.                                      (5)

Определим /> как единственное решениесистемы  R11y1=g1. Тогда:

1.      Все решения задачи о минимизации ||Ax-b|| имеют вид />,где y2 произвольно.

2.      Любой такой вектор /> приводит к одному и тому жевектору невязки />.                                             (6)

3.      Для нормыrсправедливо />

4.      Единственным решением минимальной длины является вектор />

Доказательство. В выражении для квадрата нормыневязки заменим A на HRKT всоответствии с (2) и умножая на ортогональную матрицу HT (умножение на ортогональную матрицу не меняет евклидову норму вектора) получим

/>                    (7)

Далее из (3) и (5) следует, что

/>.

Из (4) следует

/>

Подставляя оба последних выражения в (7) получим

/>

Последнее выражение имеет минимальное значение /> при R11y1=g1,а в этом уравнении единственным решением является />,так как ранг матрицы R11 равен к.Общее решение y выражается формулой />,где y2 произвольно. Для вектора /> имеем

/>,

что устанавливает равенство (3). Среди векторов /> наименьшую длину имеет тот,для которого y2=0. Отсюда следует,что решением наименьшей длины будет вектор />.Теорема доказана.

Всякое разложение матрицы А типа (2) мы будемназывать ортогональным разложением А. Заметим, что решениеминимальной длины, множество всех решений и минимальное значение для задачиминимизации ||Ax-b|| определяются единственным образом. Они независят от конкретного ортогонального разложения.

При проведении разложения необходимо приводить матрицы кдиагональному виду. Для этого обычно используются два преобразования: Гивенса иХаусхолдера, оставляющие нормы столбцов и строк матриц неизменными.

/>1.2. Ортогональноевращение Гивенса

Лемма. Пусть дан 2–вектор />, причем /> либо />.Существуетортогональная 2´2 матрица такая,что:

/>                      (8)

Доказательство. Положим:

/>.

Далее прямая проверка.

Матрица преобразования представляет собой матрицу вращений

/>

или отражений

/>

/>1.3. Ортогональноепреобразование Хаусхолдера

Применяется для преобразованияматриц к диагональному виду. Матрица преобразования представляет из себяследующее выражение: />,                                    (9)

или, если вектор v нормирован,т.е. используется вектор единичной длины />, то />. В обоих случаях H –симметричная и ортогональная матрица. Покажем это:

/>.

Отсюда следует: что />, т.е. симметричность иортогональность. В комплексном случае матрица  /> эрмитова[1] и унитарна[2].Предположим, что дан вектор х размерности m, тогда существуетматрица H  такая, что />, где

/>

а s = +1,при положительной первой компоненте вектора х и = –1, при отрицательной.

Доказательство. Положим  /> действительная матрица.Любую действительную матрицу можно привести в треугольному виду

/>

Далее принимаем во внимание то, что /> и получаем следующее:

/>

/>/>/>/>1.4.Сингулярное разложение матриц

Пусть X – матрица данных порядка Nxp, где N>p,и пусть r – ранг матрицы X. Чаще всего r=p, но приводимыйниже результат охватывает общий случай, он справедлив и при условии r<p.

Теорема о сингулярном разложении утверждает, что

/>                                                            (10)

где V – матрица порядка Nxr, столбцы которойортонормированы, т.е. />; U –матрица с ортонормированными столбцами порядка pxr; такимобразом, />; Г – диагональная матрицапорядка rxr, диагональные элементы которой />, называемые сингулярнымичислами матрицы X,положительны. Используя диагональные элементы /> матрицыГ, столбцы /> матрицы V, и столбцы/> матрицы U, сингулярноеразложение матрицы X, определяемое по (10), можно записать в виде:

/>                                      (11)

Имеют место следующие фундаментальные соотношения.

·       Квадратная симметричная матрица XX' порядка NxN,имеет r положительных и N–r нулевых собственных чисел.Положительными собственными числами XX' являются />, а соответствующими собственнымизначениями – />. Таким образом, сингулярныезначения /> – это положительныеквадратные корни из положительных собственных чисел матрицы XX', астолбцы матрицы V – соответствующие собственные векторы.

·       Квадратная симметричная матрица X'X порядка pxp,имеет r положительных и p–r нулевых собственных чисел.Положительными собственными числами X'X являются />, а соответствующими собственнымизначениями – />, таким образом, сингулярныезначения /> – это положительныеквадратные корни из положительных собственных чисел матрицы X'X, астолбцы матрицы U – соответствующие собственные векторы.

Положительные собственные числа матрицы X'X и XX'совпадают и равны />. Более того, еслиum – собственный вектор матрицы X'X, а vm – собственный вектор матрицы XX',соответствующие одному и тому же собственному числу />,то um и vm связаныследующим соотношением

/>                        (12)

Эти соотношения дают возможность вычислять />, зная />, и наоборот. В компактнойформе эти соотношения можно записать следующим образом:

/>.                                                (13)

Исследование матрицы X'X в факторном анализеназывается R-модификацией, а XX'Q–модификацией.Соотношения (12)–(13) показывают, что результаты Q–модификации можнополучить по результатам R–модификации и наоборот.

Практическая последовательность нахождения сингулярногоразложения следующая.

1.  Вычисляется X'X или XX', в зависимости от того, порядоккакой матрицы меньше. Предположим, что в данном случае это X'X.

2.  Вычисляются положительные собственные числа /> матрицыX'X и соответствующие им собственные векторы />.

3.  Находятся сингулярные числа />.

4.  Вычисляются /> по соотношению (11).

Пусть в разложении (11) собственные числа расположены впорядке убывания. Аппроксимационные свойства соотношения (11) являются ещеболее фундаментальными, чем само соотношение. Эти свойства вытекают из решенияследующих двух задач.

Задача 1. Дана симметричная матрица S,порядка pxp и ранга r с неотрицательными собственнымизначениями. Требуется найти симметричную матрицу Т, размерности pxp,с неотрицательными собственным значениями заданного ранга k, k<r,являющуюся наилучшей аппроксимацией матрицы S в смысле наименьшихквадратов.

Задача 2. Дана прямоугольная матрица X,порядка Nxp и ранга r и число k<r. требуетсянайти матрицу W порядка pxp и ранга k, наилучшимобразом аппроксимирующую матрицу X в смысле наименьших квадратов.

Решением этих двух задач являются матрицы:

/>                                  (14)

представляющие собой суммы k первых членов в соответствующемразложении. Матрицы T и W называются наилучшими в смысленаименьших квадратов “матричными аппроксимациями меньшего ранга” для матриц Sи X соответственно. Свойство наилучшей аппроксимации в смысле наименьшихквадратов можно выразить следующим образом: матрица T ближе всего кматрице S в том смысле, что сумма квадратов всех элементов матрицы S–Tминимальна. Аналогично матрица W ближе всего к матрице X в томсмысле, что минимальна сумма квадратов элементов матрицы X–W. Меройблизости или качества аппроксимации считается относительная величина />, т.е. сумма r–kнаименьших собственных чисел матрицы X’X. Иногда мерой качества аппроксимациисчитается относительная величина

/>                                              (15)

или функция от нее.

Рассмотрим наиболее распространенный случай p=r.

Матрица S может быть ковариационной матрицей pлинейно независимых переменных. Матрица T также может представлять собойковариационную матрицу p переменных, но так как ранг матрицы T k<p,то эти p переменных линейно зависят от k переменных. Такимобразом, p исходных переменных, ковариационная матрица которых есть S,могут быть приближенно выражены через k переменных.

Во второй задаче исходную матрицу X порядка Nxpможно выразить как X=VГU’, где V – матрица порядка Nxpc ортонормированными столбцами; Г – диагональная матрица порядка pxpU – квадратная ортогональная матрица порядка pxp.

Матричную аппроксимацию меньшего ранга W можнопредставить в виде

/>

где /> – состоит изпервых k столбцов матрицы V, /> –из первых k строк или столбцов матрицы Г, а /> –из первых k столбцов матрицы U. поскольку W»X, то

/>                                               (16)

При умножении этой матрицы справа на /> получаем

/>                                              (17)

Матрица /> порядкаpxk определяет преобразование строк матрицы X из евклидова p–мерногопространства в евклидово k–мерное пространство; уравнение (16)показывает, что существует преобразование матрицы X порядка Nxpв матрицу /> порядка Nxk.Матрица X содержит N точек в p–мерном евклидовомпространстве, которые приближенно могут быть спроектированы в k–мерноеевклидово пространство. матрица /> определяеткоординаты этих точек в k–мерном евклидовом пространстве.

/>1.5. QR–разложение

Теорема 2. Пусть Аm´nматрица. Существует ортогональная m´mматрица Qтакая, что в матрице QA=R подглавной диагональю стоят только нулевые элементы.

Доказательство. Выберем ортогональную m´mматрицу Qв соответствии с преобразованиемХаусхолдера (9), так, чтобы первый столбец Q1A имел нулевые компоненты со 2–ой по m–ю. Далеевыбираем ортогональную (m-1)´(m–1)–матрицуP2 следующим образом. Будучиприменена к  m–1 вектору,составленному из компонент со 2–ой по m–ю второго столбца матрицы Q1A, онааннулирует компоненты с 3–ей по m–ю этого вектора. Матрицапреобразования

/>

ортогональна, и Q2Q1A имеет впервых двух столбцах нули под главной диагональю. Продолжая таким образом,можно построить произведение, состоящее максимум из n ортогональныхпреобразований, которое трансформирует А к верхней треугольной форме.Формальное доказательство можно получить методом конечной индукции.

Полученное представление матрицы произведением ортогональнойи верхней треугольной матриц называется QR–разложением.

Теорема 3. Пусть Аm´nматрица ранга к, причем k<n£m. Существуют ортогональная m´mматрица Q и m´nматрица перестановки Pтакие, что

/>,                                                    (18)

где R – верхняя треугольная к´кматрица ранга к.

Доказательство. Выберем матрицу перестановки Ртаким образом, чтобы первые к столбцов матрицы AP, были линейнонезависимы. Согласно теореме 2, найдется ортогональная m´m–матрица Q такая, что QAPбудет верхней треугольной. Поскольку первые к столбцов АР линейнонезависимы, это будет верно для первых к столбцов QAP.

Все элементы матрицы QAP, стоящие на пересечениистрок с номерами к+1,...,m и столбцов с номерами к+1,...,n, будутнулями. В противном случае rankQAP>k, что противоречит предположению rankA=k.Итак, QAP имеет форму, указанную правой частью (4). Теорема доказана.

Подматрицу[R:T] из правой части (18) можнотеперь преобразовать к компактной форме, требуемой от матрицы R изтеоремы 2. Это преобразование описывает следующая лемма.

Лемма 1. Пусть [R:T]– к´к–матрица, причем R имеетранг к. Существует ортогональная n´n–матрица W такая, что

/>

где /> – нижняятреугольная матрица ранга к.

Доказательство леммы вытекает из теоремы 3, еслиотождествить величины n, k, [R:T], W из формулировки леммыс соответствующими величинами m, n, AT,QT теоремы 3.

Используя теорему 3 и лемму 1 можно доказать следующую теорему.

Теорема 4. Пусть Аm´n–матрица ранга к .Найдутся ортогональная m´m–матрицаН и ортогональная n´n–матрицаК такие, что

/>                                        (19)

причем R11 –невырожденная треугольная к´к–матрица.

Заметим, что выбором Н и К в уравнении (19)можно добиться, чтобы R11 былаверхней или нижней треугольной.

В (19) матрица А представлена произведением A=HRKT, где R – некоторая прямоугольная матрица,ненулевые компоненты которой сосредоточены в невырожденной треугольнойподматрице. Далее мы покажем, что эту невырожденную подматрицу R можноупростить далее до невырожденной диагональной матрицы. Это разложение тесносвязано со спектральным разложением симметричных неотрицательно определенныхматриц ATA и AAT (см. 11).

Теорема 5. Пусть Аm´n–матрица ранга k. Тогдасуществуют ортогональная m´m–матрицаU, ортогональная n´n–матрицаV и диагональнаяm´n–матрицаS такие, что

UTAV=S,A=USVT                                      (20)

Матрицу S можно выбрать так, чтобы ее диагональныеэлементы составляли невозрастающую последовательность; все эти элементынеотрицательны и ровно к из них строго положительны.

Диагональные элементы S называются сингулярнымичислами А. Докажем сперва лемму для специального случая m=n=rankA.

Лемма 2. Пусть Аn´n–матрица ранга n. Тогдасуществует ортогональная n´n–матрицаU, ортогональная n´n–матрицаV и диагональная n´n–матрицаS такие, что UTAV=S, A=USVT ипоследовательные диагональные элементы S положительны и не возрастают.

Доказательство леммы. Положительно определеннаясимметричная матрица ATA допускает спектральное разложение

ATA=VDVT,                                        (21)

где V – ортогональная n´n–матрица, а D –диагональная матрица, причем диагональные элементы D положительны и невозрастают. Определим S как диагональную n´n–матрицу, диагональные элементыкоторой суть положительные квадратные корни из соответствующих диагональныхэлементов D. Таким образом

D=STS=S2, S-1DS-1=I.                                  (22)

Определим матрицу

U=AVS-1                                                     (23)

Из (21), (22), (23) и ортогональности V следует, что

UTU=S-1VTATAVS-1=S-1DS-1=Iт.е. U ортогональна. Из (23) и ортогональностиV выводим USVT=AVS-1SVT=AVVT=A Лемма доказана.

Доказательство теоремы 5. Пусть A=HRKT, где H, R, KT имеют свойства, указанные в теореме 4. Так как R11 из (19) – невырожденная треугольная к´к–матрица, то согласно лемме 2,можно написать

/>                                                  (24)

Здесь /> и /> – ортогональные к´к–матрицы, а /> – невырожденнаядиагональная матрица, диагональные элементы которой положительны и невозрастают. Из (24) следует, что матрицу R в уравнении (19) можнозаписать в виде

/>                                                  (25)

где:

/>       –ортогональная m´m–матрица;

/>       –ортогональная n´n–матрица;

/>           –ортогональная m´n–матрица;

Теперь, определяя U и V формулами

/>                                       (26)

заключаем из (24) – (26), что A=USVT, где U, S, Vимеют свойства, указанные в формулировке теоремы 5. Это завершаетдоказательство.

  Заметим, что сингулярные числа матрицы А определеныоднозначно, в то время, как в выборе ортогональных матриц U, V естьпроизвол. Пусть s – сингулярноечисло А, имеющее кратность l. Это значит, что для упорядоченныхсингулярных чисел найдется индекс I такой, что

/>

  Положим k=min(m,n), и пусть Q –ортогональная к´к–матрицавида

/>

Здесь Р – ортогональная l´l–матрица Если A=USVT – сингулярное разложение А и si=…=si+l-1, то сингулярным разложением Абудет также и />, где          />.

/>/>/>1.6. Число обусловленности

Некоторые вычислительные задачи поразительно чувствительны кизменению данных. Этот аспект численного анализа не зависит от плавающейарифметики или выбранного алгоритма.

Например:

Найти корни полинома: (x-2)2=10-6

Корни этого уравнения есть 2+10-3 и 2-10-3.Однако изменение свободного члена на 10-6 может вызвать изменение вкорнях, равное 10-3.

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

По определению число обусловленности есть величина />. Для болееподробного описания числа обусловленности нам понадобится понятие нормы в пространствевекторов и матриц.

Нормой вектора x в пространстве векторов /> называется функционал,обозначаемый />, удовлетворяющий следующимусловиям:

1)     положительной определенности – />

2)     положительной однородности – />;

3)     неравенству треугольника – />.

Нормой квадратной матрицы А в пространствематриц, согласованной с нормой вектора /> называетсяфункционал  />, удовлетворяющийусловиям 1 – 3 для нормы вектора:

1)     />;

2)     />

3)     />

4)     мультипликативное неравенство – />

Наиболее употребимы следующие нормы для векторов:

·       норма суммы модулей           />

·       евклидова норма                     />

·       норма максимума модуля      />

Нормы матриц:

·       />

·       />

·       />

Здесь /> являютсясингулярными числами[3]матрицы А; это положительные значения квадратных корней /> из собственных значений /> матрицы АТА(которая при невырожденной матрице А положительно определена[4], впротивном случае положительно полуопределена (неотрицательно определена[5]) и поэтомуимеет только вещественные собственные значения ³0). Для вещественных симметричных матриц сингулярные числа равны абсолютнымвеличинам собственных значений: />.

Умножение вектора х на матрицу А приводит кновому вектору Ах, норма которого может очень сильно отличаться от нормывектора х.

Область изменений может быть задана двумя числами

/>

Максимум и минимум берутся по всем ненулевым векторам.Заметим, что если А вырождена, то m=0. Отношение M/mназывается числом обусловленности матрицы А,

/>                                (7)

Рассмотрим норму обратной[6]матрицы />.

Для матрицы А существует сингулярное разложение />, тогда />, отсюда />. Аналогично для обратнойматрицы /> и />. Отсюда следует, чтособственные числа матрицы /> – 1//> есть величины, обратныесобственным числам матрицы /> – />. При этом очевидно, что />. Из последнего выражениявместе с (7) следует />.Таким образом обусловленность матрицы равна произведению нормы матрицы на нормуобратной матрицы.

Рассмотрим систему уравнений Ax=b, и другую систему,полученную изменением правой части: A(x+Dx)=b+Db.Будем считать Db ошибкой в bDx соответствующей ошибкой в x,хотя нам нет необходимости считать ошибки малыми. Поскольку A(Dx)=Db,то определения M и  m немедленно приводят к неравенствам /> Следовательно, при m¹0,

/>

Величина /> естьотносительное изменение правой части, а величина />  –относительная ошибка, вызванная этим изменением. Аналогичные выкладки можнопровести не только с элементами вектора правой части но и с элементами самойматрицы А и найти зависимость между относительным изменением элементов матрицыи относительной ошибкой вызванной этим изменением. Отсюда следует, что числообусловленности выполняет роль множителя в увеличении относительной ошибки.

Приведем некоторые свойства числа обусловленности. Ясно, чтоM³m и поэтому cond(А)³1. Если Р – матрицаперестановок[7],то компоненты вектора Px лишь порядком отличаются от компонент вектора х.Отсюда следует, что /> и cond(P)=1. В частности cond(I)=1. Если А умножается на скаляр с, тоcond(cА)= cond(А). Если D – диагональная матрица, то />

/>глава 2. Реализациясингулярного разложения/>2.1. Алгоритмы

QR–алгоритм начинается с разложения матрицы поГрамму-Шмидту />, затем меняютсяместами сомножители: /> Эта матрицаподобна первоначальной, /> Этотпроцесс продолжается, причем собственные значения не изменяются:

/>

Эта формула описывает QR–алгоритм без сдвигов. Обычновремя которое тратится на такой процесс пропорционально кубу размерностиматрицы – n3. Необходимо процессускорить, для чего используется предварительное приведение матрицы А к формеХессенберга[8]а также используется алгоритм со сдвигом. Форма Хессенберга представляет изсебя верхнюю треугольную матрицу (верхняя форма Хессенберга) у которойсохранена одна диагональ ниже главной, а элементы ниже этой диагонали равнынулю. Если матрица симметрична, то легко видеть, что матрица Хессенбергапревращается в трехдиагональную матрицу[9].При использовании матрицы Хессенберга время процесса пропорционально n2, а при использовании трехдиагональнойматрицы –n.

Можно использовать другие соотношения

/>

где Qs – унитарная, а Ls– нижняя треугольная матрица. Такой алгоритм носит название QL–алгоритма.

В общем случае, когда все собственные значения матрицыразличны, последовательность матриц As имеет пределом нижнюю треугольную матрицу />,диагональные элементы которой представляют собой собственные значения матрицы А,расположенные в порядке возрастания их модулей. Если матрица А имееткратные собственные значения, то предельная матрица не является треугольной, асодержит диагональные блоки порядка p, соответствующие собственномучислу /> кратности p.

В общем случае, наддиагональный  элемент /> матрицы As на s-ом шаге асимптотически равен />, где kij– постоянная величина. Сходимость QL–алгоритма вообще говорянедостаточна. Сходимость можно улучшить, если на каждом шаге вместо матрицы As использовать матрицу As-ksI (QL–алгоритм со сдвигом).Последовательность вычислений в этом случае описывается следующимисоотношениями:

/>

которые определяют матрицу />.При этом асимптотическое поведение элемента /> определеносоотношением />, а не />, как прежде. Если сдвиг ks выбрать близко к величине /> (наименьшее собственноезначение), то в пределе внедиагональные элементы первой строки будут оченьбыстро стремиться к нулю. Когда ими можно пренебречь, элемент /> с рабочей точностью равен />, остальные являютсясобственными значениями оставшейся матрицы n-1-го порядка. Тогда,если QL–алгоритм выполнен без ускорения сходимости, то все равно />, и поэтому автоматическиможно выделить величину сдвига ks.

Если матрица А эрмитова, то очевидно, что и всематрицы Аs эрмитовы; если Адействительная и симметричная, то все Qsортогональны и все Аsдействительны и симметричны.

/>2.2. Реализацияразложения

Таким образом, разложение /> производится в два этапа.Сначала матрица А посредством двух конечных последовательностейпреобразований Хаусхолдера где />, приводится к верхнейдвухдиагональной форме следующего вида:

/>

Далее реализуется итерационный процесс приведениядвухдиагональной матрицы J0кдиагональной форме, так что имеет место следующая последовательность: /> где /> а Siи Ti  – диагональные матрицы.

Матрицы Tiвыбираются так, чтобы последовательность матриц /> сходилась кдвухдиагональной матрице. Матрицы же Si выбирают так, чтобывсе Ji  сохраняли двухдиагональную форму. Переход /> осуществляетсяс помощью плоских вращений (10) – преобразований Гивенса. Отсюда, /> где

/>

а матрица /> вычисляется аналогично сзаменой /> на/>.

Пусть начальный угол /> произволен, однакоследующие значения угла необходимо выбирать так, чтобы матрица Ji+1 имела ту же форму,что и Ji. Таким образом /> неаннулирует  ни одного элемента матрицы, но добавляет элемент />; /> аннулирует/> нодобавляет />;/> аннулирует/> нодобавляет /> ит.д., наконец, /> аннулирует/> и ничегоне добавляет.

Этот процесс часто называют процессом преследования. Так как/>, то />, и Mi+1 – трехдиагональнаяматрица, точно так же, как и Mi.Начальный угол /> можновыбрать так, чтобы преобразование /> было QR–преобразованиемсо сдвигом, равным s.

Обычный QR–алгоритм со сдвигом можно записать вследующем виде:

/>

где /> –верхняя треугольная матрица. Следовательно, />. Параметр сдвига s определяется собственным значением нижнего минора (размерности 2´2) матрицы Mi.При таком выборе параметраs метод обладает глобальной и почти всегдакубичной сходимостью.

/>2.3. Пример сингулярногоразложения

Проведем преобразование Хаусхолдера на матрице />,

К первой компоненте первого столбца прибавляем норму первогостолбца, получим />. Пусть />

Преобразованная матрица A2вычисляется следующим образом. Для первого столбца имеем:

/>

так как />

Таким образом, в первый столбец были введены нули и егодлина не изменилась. Получим второй столбец:

  />

для третьего столбца:

/>

окончательно,

/>

Столбцы матрицы A2получаются вычитанием кратных вектора v1из столбцов A1. Эти кратныепорождаются скалярными произведениями, а не отдельными элементами, как вгауссовом исключении.

Прежде чем вводить дальнейшие нули под диагональю,преобразованием вида A3=A2Q1,получают нули в первой строке. Нули уже стоящие в первом столбце, не должныбыть испорчены, длина первого столбца должна быть сохранена; поэтому привнесении нулей в первую строку нельзя менять первый элемент строки, изменяемвторой элемент и зануляем третий. Для матрицы большего размера на этом шагебыло бы получено n–2 нуля.

Преобразование порождается первой строкой A2:

/>

Строка матрицы A3с номером i получается по формуле

/>.

Таким образом, из каждой строки A2вычитается надлежащее кратное />. Этодает матрицу

/>

Поскольку первая компонента />нулевая,то нули первого столбца A2сохраняются в A3, Так как Q1 ортогональная, то длина каждой строки в A3 равна длине соответствующей строки в A2.

Теперь можно добиться новых нулей под диагональю, неиспортив полученных ранее:

/>

Поскольку ранг этой матрицы равен лишь 2, то теперь третийстолбец имеет на диагонали и под диагональю  элементы порядка ошибкиокругления. Эти элементы обозначены в матрице через 0.000, чтобы отличить их отэлементов, в точности равных нулю. Если бы матрица имела полный ранг, то нужнобыло бы выполнить еще одно преобразование, чтобы получить нули в третьемстолбце:

/>

Если бы не ошибки округлений, то в данном примере третийдиагональный элемент был бы точным нулем. Элементы под диагональю во всехстолбцах указаны как точные нули, потому что преобразования так и строились,чтобы получить там нули. Последнее преобразование H3в этом примере могло бы быть тождественным, однако тогда оно не было быхаусхолдеровым отражением. Фактически использование H3попутно изменяет знак  элемента – 1.080 в матрице A4.

Получилась искомая двухдиагональная матрица, и первый этапзакончен. Прямое использование ортогональных преобразований не позволяетполучить какие–либо новые нули. Для общего порядка n нужно nпреобразований H и n–2  преобразований Q, чтобы достигнутьданного места. Число преобразований не зависит от строчной размерности m,но от m зависит работа, затрачиваемая на выполнение каждогопреобразования.

/>глава 3. Использованиесингулярного разложения  в методе наименьших квадратов

При использовании метода сингулярного разложения (SVDSingular Value Decomposition) мы проводим разложение дляматрицы плана. При этом основное уравнение y=Xb приобретает вид y=USVTb.Отсюда следует, что коэффициенты b можно получить решая уравнение UTy=SVTb. Т.е. если все sj,j=1,…,n, являющиесядиагональными элементами S отличны отнуля, то последнее уравнение разрешимо и

/>, где />.

Однако такое решение не всегда желательно, если некоторые sjмалы. Для правильного использования метода SVD мы должны ввести границу t отражающую точность входных данных иточность использованной плавающей арифметики. Всякое sj, большее, чем t, приемлемо, и соответствующее /> вычисляется по (1.20).Любое sj,меньшее, чем t, рассматриваетсякак пренебрежимо малое, и соответствующему /> можетбыть придано произвольное значение. С этой произвольностью значений связана не единственностьнабора коэффициентов, получаемых методом наименьших квадратов. Изменениявходных данных и ошибки округлений, меньшие, чем t, могут привести ксовершенно другому набору коэффициентов, определяемых этим методом. Посколькуобычно желательно, чтобы эти коэффициенты были по возможности малы, то полагаем/>=0, если sj £t.

Отбрасывание чисел sj, меньших, чем t,приводит к уменьшению числа обусловленности с /> до/>. Поскольку числообусловленности является множителем в увеличении ошибки, то следствием будетболее надежное определение коэффициентов />.

Продемонстрируем использование метода на следующем примере:

t Y 1900 75994575 1910 91972266 1920 105710620 1930 123203000 1940 131669275 1950 150697361 1960 179323175 1970 203211926

Следует определить значение Y при X =1980.

Если аппроксимировать эти данные квадратичным многочленом /> и использовать двойнуюточность, то в результате получим следующие коэффициенты />. При одинарной точностивычислений коэффициенты будут иметь значения />.У этих двух наборов коэффициентов не совпадают даже знаки. Если такую модельиспользовать для предсказания, то для коэффициентов, вычисленных с двойнойточностью, прогноз будет Y=227780000, а дляобычной точности Y=145210000. Ясно, что второйнабор коэффициентов бесполезен. Исследуем достоверность результатов. Матрицаплана для данного примера имеет размеры 8 ´3

/>

Рис. 2. Численныйпример

/>

Ее сингулярные числа: />.

Число обусловленности равно />,что говорит о близости базисных функций 1, t и t2к линейной зависимости. Для того, чтобы исправить ситуацию можно предпринятьодну из следующих мер.

Во–первых, можно выбрать границу для относительнойошибки, которая бы отражала точность данных и точность арифметики. Если взятьграницу в интервале />, то отбросимтретье сингулярное число. При этом получим следующие наборы коэффициентов длядвойной и обычной точности:

/>

Теперь коэффициенты находятся в гораздо лучшем согласии другс другом. Кроме того, коэффициенты стали существенно меньше, а это значит, чтоне будет столь большого, как прежде, взаимного уничтожения слагаемых привычислении квадратичного многочлена. Прогнозное значение Y(1980)будет соответственно 212910000 и 214960000. Эффект обычной точности ещезаметен, однако результаты уже не являются катастрофическими.

Можно также определить набор нулевых коэффициентов,соответствующих пренебрежимо малому сингулярному числу. Вот эти коэффициенты: />. Для значений t от1900 до 1970 величина функции /> непревосходит 0.0017, поэтомупри любом a коэффициенты можноизменить />, и при этом значения,выдаваемые моделью изменятся не более чем на 0.0017a. Любой из четырех перечисленных нами наборов коэффициентовможно получить из другого подобным изменением.

Во–вторых, можно улучшить ситуацию заменой базиса.Модели

/>

гораздо более удовлетворительны. Важно при этом то, чтонезависимая переменная преобразуется из интервала [1900, 1970] в какой–нибудьболее приемлемый интервал вроде [0, 70] или, еще лучше, [–3.5, 3.5]. Числаобусловленности при этом равны 5750 и 10.7 соответственно. последнее значениеболее чем приемлемо даже при счете с обычной точностью.

Удобнее всего воспользоваться стандартными способамистатистического анализа, т.е. матрицу плана преобразуем к стандартизованномуварианту Матрица стандартизованных данных есть матрица наблюдений с нулевым средними дисперсией 1. Это означает, что данные берутся в виде отклонений от среднего,которое мы считаем равным 0, вводим нормировку деля каждый член столбца матрицына корень квадратный из суммы квадратов отклонений.

Во втором случае, после преобразования матрицы плана ееобусловленность сильно уменьшается, и, соответственно, повышается точностьрасчетов.

Данную программу можно использовать и при решении системылинейных уравнений вместо методов Гаусса, Жордана, Холесского и пр. Вприложении 2 приведен пример расчета линейной системы, которая изначально неможет быть решена этими методами вследствие вырожденности матрицыкоэффициентов. Тем не менее, исследуемый метод дает нам правильное решение.

/>/>/>ЗАКЛЮЧЕНИЕ

В работе описаны компьютерные методы решения задачинаименьших квадратов. Для использования данных методов составленасоответствующая программа на алгоритмическом языке FORTRAN.Программа апробирована, результаты тестирования показывают работоспособностьпрограммы.

Результаты данной разработки могут быть использованы в самыхразнообразных расчетах, где необходимо провести аппроксимацию данных заданнымифункциями.

/>/>/>ЛИТЕРАТУРА

1.      Беллман Р. Введение в теорию матриц. -М.: Наука, 1969, 368с.

2.      Гантмахер Ф.Р. Теория матриц. -М.: Наука, 1988, 548с.

3.      Ланкастер П. Теория матриц. -М.: Наука, 1982, 387с.

4.      Лоусон Ч., Хенсон Р. Численное решение задач наименьших квадратов. М.:Статистика, 1979, 447с

5.      Марчук Г.И. Методы вычислительной математики. М.: Наука, 1980

6.      Мэйндоналд Дж. Вычислительные алгоритмы в прикладной статистике. М.:Финансы и статистика, 1988, 350с

7.      Стренг Г. Линейная алгебра и ее применения. М.: Мир, 1980, 454с

8.      Уилкинсон Дж., Райнш К. Справочник алгоритмов на языке АЛГОЛ. Линейнаяалгебра, М.: Машиностроение, 1976, 390с

9.      Фаддеев Д.К., Фаддеева В.Н. Вычислительные методы линейной алгебры. -М.:Физматгиз, 1963, 536с.

10.   Форсайт Дж., Малькольм М., Моулер К. Машинные методы математическихвычислений. М.: Мир, 1980, 279с

11.   Харебов К.С. Компьютерные методы решения задачи наименьших квадратов ипроблемы собственных значений. Владикавказ.: Изд-во СОГУ, 1995, 76 с.

/>/>/>ПРИЛОЖЕНИЕ 1. Исходные текстыпрограммы

                REAL A(3,3), U(3,3), V(3,3),SIGMA(3), WORK(3),Y(3),C(3),Y0(3)

                INTEGER I,IERR, J, M, N, NM

                OPEN(6,FILE=«SVD.OUT»,STATUS=«UNKNOWN»,FORM=«FORMATTED»)

        OPEN (5,FILE= «SVD.IN»,STATUS=«UNKNOWN»,FORM=«FORMATTED»)

140     FORMAT(3I5)

150     FORMAT(4E15.7)

        READ(5,140) NM,M,N

        DO 131 I=1,M

        READ(5,150) (A(I,J),J=1,N)

131     CONTINUE

        READ (5,150) (Y(I),I=1,M)

                CALLSVD(NM,M,N,A,SIGMA,.TRUE.,U,.TRUE.,V,IERR,WORK)

                IF(IERR.NE.0) WRITE (6,2) IERR

2             FORMAT(15H TROUBLE.IERR=,I4)

        WRITE(6,120)

120     FORMAT(/'МАТРИЦА А')

        DO 121 I=1,M

        WRITE(6,130) (A(I,J),J=1,N)

130     FORMAT(8E15.7)

121     CONTINUE

        WRITE (6,160) (Y(I),I=1,N)

160     FORMAT(/'ПРАВЫЕ ЧАСТИ'/8E15.7)

210     FORMAT(/'СИНГУЛЯРНЫЕ ЧИСЛА')

        WRITE(6,210)

                DO 3 J=1,N

                WRITE(6,6) SIGMA(J)

3             CONTINUE

        SMA=SIGMA(1)

        SMI=SIGMA(1)

        DO 211 J=2,N

        IF(SIGMA(J).GT.SMA) SMA=SIGMA(J)

        IF(SIGMA(J).LT.SMI.AND.SIGMA(J).GT.0.)SMI=SIGMA(J)

211     CONTINUE

        OBU=SMA/SMI

230     FORMAT(/'ЧИСЛО ОБУСЛОВЛЕННОСТИ=',E15.7)

        WRITE(6,230) OBU

        SIGMA1=0.

        DO 30 J=1,N

        IF(SIGMA(J) .GT. SIGMA1) SIGMA1=SIGMA(J)

        C(J)=0.

30      CONTINUE

        TAU=SIGMA1*0.1E-6

        DO 60 J=1,N

        IF(SIGMA(J).LE.TAU) GO TO 60

        S=0.

        DO 40 I=1,N

        S=S+U(I,J)*Y(I)

40      CONTINUE    

        S=S/SIGMA(J)

        DO 50 I=1,N

        C(I)=C(I) + S*V(I,J)

50      CONTINUE

60      CONTINUE

                write (6,560)

        WRITE (6,6) (C(I),I=1,3)

        DO 322 J=1,N

        SS=0.

        DO 321 I=1,M

321     SS=A(J,I)*C(I)+SS

322     Y0(J)=SS

                write (6,570)

        WRITE (6,6) (Y0(I),I=1,3)

C             WRITE(6,7)

C             DO 4 I=1,M

C             WRITE(6,6) (U(I,J),J=1,N)

C4          CONTINUE

C             WRITE(6,7)

C             DO 5 I=1,N

C             WRITE(6,6) (V(I,J),J=1,N)

C5          CONTINUE

6             FORMAT(3E15.7)

560         format(2x,'roots')

570         format(2x,'right')

7             FORMAT(1H )

                STOP

                E N D

SUBROUTINE SVD(NM,M,N,A,W,MATU,U,MATV,V,IERR,RV1)

        REAL A(NM,N),W(N),U(NM,N),V(NM,N),RV1(N)

        LOGICAL MATU,MATV

        IERR=0

        DO 100 I=1,M

        DO 100         J=1,N                    

        U(I,J)=A(I,J)

100     CONTINUE

                G=0.0

                SCALE=0.0

                ANORM=0.0

                DO 300 I=1,N

                L=I+1

                RV1(I)=SCALE*G

                G=0.0

                S=0.0

                SCALE=0.0

                IF(I.GT.M) GO TO 210

                DO 120 K=I,M

120         SCALE=SCALE+ABS(U(K,I))

                IF(SCALE.EQ.0.0) GO TO 210

                DO 130 K=I,M

                U(K,I)=U(K,I)/SCALE

                S=S+U(K,I)**2

130     CONTINUE

                F=U(I,I)

                G=-SIGN(SQRT(S),F)

                H=F*G-S

                U(I,I)=F-G

                IF(I.EQ.N) GO TO 190

                DO 150 J=L,N

                S=0.0

                DO 140 K=I,M

140         S=S+U(K,I)*U(K,J)

                F=S/H

                DO 150 K=I,M

                U(K,J)=U(K,J)+F*U(K,I)

150         CONTINUE

190     DO 200 K=I,M

200         U(K,I)=SCALE*U(K,I)

210         W(I)=SCALE*G

                G=0.0

                S=0.0

                SCALE=0.0

                IF(I.GT.M.OR.I.EQ.N) GO TO 290

                DO 220 K=L,N

220         SCALE=SCALE+ABS(U(I,K))

                IF(SCALE.EQ.0.0) GO TO 290

                DO 230 K=L,N

                U(I,K)=U(I,K)/SCALE

                S=S+U(I,K)**2

230         CONTINUE

                F=U(I,L)

                G=-SIGN(SQRT(S),F)

                H=F*G-S

                U(I,L)=F-G

                DO 240 K=L,N

240         RV1(K)=U(I,K)/H

                IF(I.EQ.M) GO TO 270

                DO 260 J=L,M

                S=0.0

                DO 250 K=L,N

250         S=S+U(J,K)*U(I,K)

                DO 260 K=L,N

                U(J,K)=U(J,K)+S*RV1(K)

260         CONTINUE

270         DO 280 K=L,N

280         U(I,K)=SCALE*U(I,K)

290         ANORM=AMAX1(ANORM,ABS(W(I))+ABS(RV1(I)))

300         CONTINUE

                IF(.NOT.MATV) GO TO 410

                DO 400 II=1,N

                I=N+1-II

                IF(I.EQ.N) GO TO 390

                IF(G.EQ.0.0) GO TO 360

                DO 320 J=L,N

320         V(J,I)=(U(I,J)/U(I,L))/G

                DO 350 J=L,N

                S=0.0

                DO 340 K=L,N

340         S=S+U(I,K)*V(K,J)

                DO 350 K=L,N

                V(K,J)=V(K,J)+S*V(K,I)

350         CONTINUE

360         DO 380 J=L,N

                V(I,J)=0.0

                V(J,I)=0.0

380         CONTINUE

390         V(I,I)=1.0

                G=RV1(I)

                L=I

400         CONTINUE

410         IF(.NOT.MATU) GO TO 510

                MN=N

                IF(M.LT.N) MN=M

                DO 500 II=1,MN

                I=MN+1-II

                L=I+1

                G=W(I)

                IF(I.EQ.N) GO TO 430

                DO 420 J=L,N

420         U(I,J)=0.0

430         IF(G.EQ.0.0) GO TO 475

                IF(I.EQ.MN) GO TO 460

                DO 450 J=L,N

                S=0.0

                DO 440 K=L,M

440         S=S+U(K,I)*U(K,J)

                F=(S/U(I,I))/G

                DO 450 K=I,M

                U(K,J)=U(K,J)+F*U(K,I)

450         CONTINUE

460         DO 470 J=I,M

470         U(J,I)=U(J,I)/G

                GO TO 490

475         DO 480 J=I,M

480         U(J,I)=0.0

490         U(I,I)=U(I,I)+1.0

500         CONTINUE

510         DO 700 KK=1,N

                K1=N-KK

                K=K1+1

                ITS=0

520         DO 530 LL=1,K

                L1=K-LL

                L=L1+1

                IF(ABS(RV1(L))+ANORM.EQ.ANORM) GOTO 565

                IF(ABS(W(L1))+ANORM.EQ.ANORM) GOTO 540

530         CONTINUE

540         C=0.0

        S=1.0

        DO 560 I=L,K

                F=S*RV1(I)

                RV1(I)=C*RV1(I)

                IF(ABS(F)+ANORM.EQ.ANORM) GO TO565

                G=W(I)

                H=SQRT(F*F+G*G)

                W(I)=H

                C=G/H

                S=-F/H

                IF(.NOT.MATU) GO TO 560

                DO 550 J=1,M

                Y=U(J,L1)

                Z=U(J,I)

                U(J,L1)=Y*C+Z*S

                U(J,I)=-Y*S+Z*C

550         CONTINUE

560         CONTINUE

565         Z=W(K)

                IF(L.EQ.K) GO TO 650

                IF(ITS.EQ.30) GO TO 1000

                ITS=ITS+1

                X=W(L)

                Y=W(K1)

                G=RV1(K1)

                H=RV1(K)

                F=((Y-Z)*(Y+Z)+(G-H)*(G+H))/(2.0*H*Y)

                G=SQRT(F*F+1.0)

                F=((X-Z)*(X+Z)+H*(Y/(F+SIGN(G,F))-H))/X

                C=1.0

                S=1.0

                DO 600 I1=L,K1

                I=I1+1

                G=RV1(I)

                Y=W(I)

                H=S*G

                G=C*G

                Z=SQRT(F*F+H*H)

                RV1(I1)=Z

                C=F/Z

                S=H/Z

                F=X*C+G*S

                G=-X*S+G*C

                H=Y*S

                Y=Y*C

                IF(.NOT.MATV) GO TO 575

                DO 570 J=1,N

                X=V(J,I1)

                Z=V(J,I)

                V(J,I1)=X*C+Z*S

                V(J,I)=-X*S+Z*C

570         CONTINUE

575         Z=SQRT(F*F+H*H)

                W(I1)=Z

                IF(Z.EQ.0.0) GO TO 580

                C=F/Z

                S=H/Z

580         F=C*G+S*Y

                X=-S*G+C*Y

        IF(.NOT.MATU) GO TO 600

                DO 590 J=1,M

                Y=U(J,I1)

                Z=U(J,I)

                U(J,I1)=Y*C+Z*S

                U(J,I)=-Y*S+Z*C

590         CONTINUE

600         CONTINUE

        RV1(L)=0.0

                RV1(K)=F

                W(K)=X

                GO TO 520

650         IF(Z.GE.0.0) GO TO 700

                W(K)=-Z

                IF(.NOT.MATV) GO TO 700

                DO 690 J=1,N

690         V(J,K)=-V(J,K)

700         CONTINUE

                GO TO 1001

1000      IERR=K

1001      RETURN

                E N D

/>ПРИЛОЖЕНИЕ 2. контрольныйпримерВходные данные

/>

(матрица изначально сингулярна – первая строка равна сумме второй итретьей с обратным знаком)

3    3    3

   .3200000E 02   .1400000E 02   .7400000E02

 -0.2400000E 02 -0.1000000E 02 -0.5700000E02

 -0.8000000E 01 -0.4000000E 01 -0.1700000E02

 -0.1400000E 02  0.1300000E 02  0.1000000E01

Полученный результат

МАТРИЦА А

   .3200000E+02   .1400000E+02   .7400000E+02

  -.2400000E+02  -.1000000E+02  -.5700000E+02

  -.8000000E+01  -.4000000E+01  -.1700000E+02

ПРАВЫЕ ЧАСТИ

  -.1400000E+02   .1300000E+02   .1000000E+01

СИНГУЛЯРНЫЕ ЧИСЛА

   .1048255E+03

   .7310871E-06

   .1271749E+01

ЧИСЛО ОБУСЛОВЛЕННОСТИ=   .1433830E+09

  Корни

   .1215394E+01   .1821742E+01  -.1059419E+01

  Правые корни после проверки

  -.1400000E+02   .1300000E+02   .1000001E+01

Видно, что правые частисоответствуют начальным данным. Решение верно.

еще рефераты
Еще работы по математике