close

Вход

Забыли?

вход по аккаунту

?

FreeBASIC17

код для вставкиСкачать
Пособие для начинающих программировать на языке высокого уровня FreeBASIC. Брошюра 17 "Вычислительные задачи алгебры" – продолжение пособия, которое будет интересно для учащихся школ, студентов институтов, а также преподавателей.
 1
Евгений Рыжов, инженер
Программирование на языке FreeBASIC
пособие для начинающих
Пособие для начинающих программировать на языке высокого уровня FreeBASIC.
Работа состоит из нескольких глав. В главе 17 "Вычислительные задачи алгебры" –
продолжени
е пособия, которое будет интересно для учащихся школ,
студентов
институтов, а также преподавателей.
Фрагмент 17. Вычислительные задачи алгебры
Он стал поэтом —
для математика
у него не хватало фантазии.
Давид Гильберт об одном из своих учеников...
На экзаменах глупцы предлагают вопросы,
На которые мудрецы не могут ответить.
Оскар Уайльд
2
Но вот мехмат. Евгений мнется,
В замочну скважину глядит,
Вздохнув, рукой за дверь берется
И в залу входит. Страшный вид!
Доцент, профессор —
колет, режет,
Вой, крик
и, плач, зубовный скрежет,
Пытаемых студентов стон —
И смерть, и ад со всех сторон.
Нет, не могу. Кладу перо я,
Чтоб описать кровавый пир,
Потребен минимум Омир...
Роман в стихах "Евгений Неглинкин"
Песнь четвертая. Глава XII.
Шуточная поэма "Евгений Не
глинкин" занимала такое же место в студенческом фольклоре математиков, как породившая ее поэма Пушкина "Евгений Онегин" в русской поэзии, выгодно отличаясь от других перлов этого жанра, как мастерством исполнения, так и грандиозностью замысла. Написанная е
ще в предвоенные годы двумя студентами мехмата Леонидом Трудлером и Александром Штерном под собирательным псевдонимом Аллеон Труште, она (по воспоминаниям профессора Б.В.Гнеденко) сразу очаровала не только студентов мехмата, н
о и преподавателей.
1.
Почт
овый ящик. Введение.
Ждал сообщений от читателей и дождался! Наш давний читатель "NiV"
знакомый по письму в брошюре "Фрагмент 8. Решение нелинейных уравнений" пишет:
"...понравилась брошюра "Фрагмент 14. Классические ортогональные полиномы" и хотелось бы продолжить разговор про полиномы Лагерра..."
Что еще о них можно написать не слишком утомляя аудиторию? Конечно, писать о них можно бесконечно, но будет ли это интересно всем? Однако попробую... ведь "Благосостояние народов зависит от прогресса математики"
! Это утверждение принадлежит Наполеону. Наполеон I Бонапарт (Napoleon Bonaparte; 1769 —
1821) —
император Франции в 1804
—
1815 годах, великий полководец и государственный деятель, заложивший основы современного французского государства. С ним нельзя не сог
ласиться -
всегда требуются вычисления (неважно чего) и именно математическое исчисление лежит в основе всех наук. Люди все стараются посчитать –
сколько ступенек нам следует пройти, чтобы забраться наверх, сколько облаков плывут по небу, какое расстояние преодолеть, чтобы не попасть под машину, сколько соли положить в суп, чтобы он был не пересолен... короче –
поверить алгеброй гармонию (А.С.Пушкин "Моцарт и Сальери", 1830)... Но есть и серьезные научные исследования в этой области. Феномену математики, на
пример, посвящен очерк Владимира Андреевича Успенского "Апология математики, или о математике как части духовной культуры".
3
Владимир Андреевич Успенский
, родился 27 ноября 1930 в Москве, российский математик, лингвист, публицист, доктор физико
-
математиче
ских наук (1964), профессор. Автор работ по математической логике, лингвистике, а также художественных произведений в жанре мемуарной прозы. За книгу "Апология математики" В.А. Успенский получил главный приз премии "Просветитель
-
2010" в области естественны
х и точных наук.
Его страница в сети:
http://lpcs.math.msu.su/~uspensky/
Книгу "Апология математики" можно скачать по ссылке:
http
://
knigosite
.
org
/
library
/
books
/29069
Главные тезисы книги:
1.
Математика, вне зависимости от её практического использования, принадлежит духовной культуре.
2.
Отдельные фрагменты математики входят в общеобязательную часть этой культуры.
Наша сегодняшняя задача гора
здо прозаичнее –
удовлетворить, насколько это возможно, просьбу читателя "NiV". Иными словами, речь пойдет о вычислительных задачах алгебры
, как части вычислительной математики
—
области математики, посвященной приближённому решению математических и физиче
ских задач, аналитическое решение которых представляется невозможным или затруднительным.
Вычислительная математика обладает широким кругом прикладных применений для проведения научных и инженерных расчётов. На её основе в последнее десятилетие образовали
сь такие новые области естественных наук, как вычислительная химия, вычислительная биология и так далее. Детали можно узнать на страницах:
http://ru.wikipedia.org/wiki/Вычислительная_м
атематика
http://en.wikipedia.org/wiki/Computational_mathematics
В вычислительной математике выделяют следующие направления:
анализ математических моделей, разработка методов и алго
ритмов решения стандартных математических задач, автоматизация программирования.
Вычислительные задачи алгебры [01]
обладают обманчивой простотой, поскольку исследователь по самой постановке вопроса имеет дело с конечномерными объектами —
многочленами, ве
кторами, матрицами
. Здесь не приходится заниматься тонкими и глубокими вопросами дискретизации и, казалось бы, имеется алгоритм —
бери и решай численно задачу. Но эта простота чисто внешняя, и при численном решении задач алгебры приходится сталкиваться с т
рудными ситуациями, 4
так что во многих случаях вычислителю приходится проявлять подлинное искусство, а не холодную рутинную выучку. В брошюре используется стандартная терминология
и факты, известные читателю из курсов общей алгебры и линейной алгебры [02]
.
2.
Традиционная "художественная галерея".
Когда б вы знали, из какого сора
Растут стихи, не ведая стыда,
Как желтый одуванчик у забора,
Как лопухи и лебеда.
"Мне ни к чему одические рати...", 1940
Анна Андреевна Ахматова (1889 —
1966)
Буквально на дн
ях по телевизору показывали документальный фильм "Я пришел дать вам сказку. Художник Ефим Честняков"
-
автор сценария Татьяна Басова, автор идеи Савва Ямщиков, режиссер Алексей Илюхин, произведено: Россия, телекомпания "Голд Медиум" 2008, продолжительность
00:44:37. Если не смотрели –
посмотрите, например, на сайте ежегодного Международного фестиваля кинофильмов и телепрограмм "Радонеж" по ссылке:
http://fest.radonezh.ru/video/52212
Фильм посвящен русск
ому художнику Ефиму Честнякову, ученику Ильи Репина, которому прочили большое будущее. Но он уехал в свою родную деревню Шаблово Костромской области и всю жизнь прожил там в безвестности и бедности, рисуя соседских ребятишек, сочиняя им стихи и сказки. Уже
после его смерти Саввой Ямщиковым была начата работа по реставрации картин Честнякова и возрождению забытого имени замечательного художника и светлого человека, образца бескорыстного служения своему народу.
Ефим Васильевич Честняков
(Евфимий Самойлов) (1874 —
1961) —
русский художник (портреты и сказочные сюжеты в русле псевдонаивного искусства), писатель (сказки, сказания, роман в стихах, стихи, размышления), скульптор (мелкая глиняная пластика), создатель детского театра в Шаблово. Пик творчества приш
ёлся на первую четверть XX века. Умер художник 27 июня 1961 года в своем доме
-
мастерской прозванным Честняковым "шалашка" (очень современно -
шарашка —
жаргонное название НИИ и КБ закрытого типа)...
5
Картина "Город Всеобщего Благоденствия"
После сме
рти Честнякова многое было разобрано односельчанами "на память", а главная картина "Город Всеобщего Благоденствия" была даже поделена на куски. Из глиняного "Кордона", состоявшего более чем из 800 фигур, до настоящего времени дошло ничтожное количество —
о
коло сорока произведений; что
-
то было растащено детьми из опустевшего дома Ефима. Картины, которые с 1970
-
х годов начала восстанавливать группа реставраторов, "дошли до нас в тяжёлом состоянии". Но сейчас не об этом. В фильме после фрагмента о "шалашке" (0
0:26:45)
зав.отделом Костромского художественного музея Татьяна Сухарева интересуется: куда уходит крестьянство
на картине "Вход в город всеобщего благоденствия"?
Картина "Вход в город всеобщего благоденствия"
6
Очевидно, Ефим Васильевич Честняков знал отв
ет на этот вопрос!
Он получил его в свой первый "петербургский" период (1899
—
1905), после общения с "интересными" людьми, помнившими члена петербургской академии наук, великого математика Леонарда Эйлера
(1707 —
1783), одного из "теоретиков" "Полой Земли".
У Эйлера были знаменитейшие предшественники: древнегреческий философ Платон
(428
\
427 —
348
\
347 до н.э., ученик Сократа, учитель Аристотеля), английский астроном и математик Эдмунд Галлей
(1656 —
1742)... Другой серьезный ученый —
шотландский физик и матем
атик Джон Лесли
(1766
—
1832) —
развил теорию Эйлера, предположив, что внутри планеты находится сразу два солнца —
Плутон и Прозерпина. В 19 веке подобные гипотезы вышли за пределы научных кругов и превратились в идефикс для всякого рода мечтателей, жуликов и шизофреников. Трудно сказать, к какой из этих категорий принадлежал американец Джон Клив Саймс
(1779
—
1829), но этот бывший вояка впервые сделал теорию "Полой Земли" достоянием самой широкой мировой общественности. Саймс не публиковал книг, так как был сл
ишком занят чтением лекций. Схема устройства планеты выглядела примерно так:
Предоставляя богатое поле для воображения, теории полой Земли нашли своё отражение в художественной литературе, кино, видеоиграх; часто эта тема встречается в эзотерической ли
тературе...
7
3.
Сопутствующая тема. Два слова о "математической физике".
Пьяный физик отдыхает на скамейке,
прикрывшись учебником "Теория поля".
К нему подходит прохожий и говорит:
-
Агроном, а нажрался как физик...
Так вот плавненько мы перешли к устр
ойству Мира. Есть такое научное направление "Математическая физика"
, в его русле развивается теория математических моделей физических явлений. Направление это относится к математическим наукам: критерий истины в ней —
точное математическое доказательство. Однако, в отличие от чисто математических наук, в математической физике исследуются физические задачи на математическом уровне, а результаты представляются в виде теорем, графиков, таблиц и т.д. и получают доступную физическую интерпретацию. Редакционная к
оллегия журнала "Journal of Mathematical Physics" (англоязычный рецензируемый научный журнал, издаваемый Американским институтом физики с 1960; официальный сайт журнала: http://jmp.aip.org/) определяет математическую физику как "применение математики к физ
ическим задачам и разработка математических методов, подходящих для таких приложений и для формулировок физических теорий". Первоначально математическая физика сводилась к краевым задачам для дифференциальных уравнений. Это направление составляет предмет "
классической" математической физики...
Никакими интеллектуальными доводами вы не сможете передать глухому ощущение музыки. Точно так же никакими интеллектуальными доводами нельзя передать понимание природы человеку "другой культуры". Философы пытаются ра
ссказать о природе без математики. Я пытаюсь описать природу математически. Но если меня не понимают, то не потому, что это невозможно. Может быть, моя неудача объясняется тем, что кругозор этих людей чересчур ограничен и они считают человека центром Вселе
нной. Так писал в своей лекции "Связь математики с физикой" Ричард Филлипс Фейнман
(Richard Phillips Feynman; 1918 —
1988) —
выдающийся американский учёный, основные достижения относятся к области теоретической физики.
Для желающих поближе познакомиться с "физической" стороной рассматриваемой сегодня темы, рекомендую книгу:
Фейнман Р. Характер физических законов, М., "Наука", 1987. 164 с.
Библиотечка Квант, выпуск 62, ссылка на книгу:
http://vivovoco.ibmh.msk.su/VV/Q_PROJECT/FEYNMAN/CONT.HTM
Лекция 1. Пример физического закона –
Закон тяготения.
Лекция 2. Связь математики с физикой.
Лекция 3. Великие законы сохранения.
Лекция 4. Симметрия физических законов.
8
Лекция 5. Различие про
шлого и будущего.
Лекция 6. Вероятность и неопределенность.
Лекция 7. В поисках новых законов.
Чтобы немного укротить математический оптимизм Ричарда Фейнмана ниже приведена ссылка на мультимедийную статью Николая Викторовича Левашова
(1961 –
2012; русско
го писателя, ученого
-
исследователя) "Теория Вселенной и объективная реальность" (2006). Статья состоит из 6 небольших частей примерно по 10 минут каждая. Ссылка на статью:
http://levasho
v
-
media.com/articles_content.php?id=25
"...Человек последние несколько тысяч лет постоянно пытался осмыслить окружающий Космос. Создавались разные модели Вселенной и представления о месте человека в ней. Постепенно эти представления сформировались в, та
к называемую, научную теорию Вселенной. Эта теория была окончательно сформирована в середине двадцатого века...
На основе созданных человеком представлений об окружающей природе, создаются технологии, приборы и машины. И от того, какими они создаются –
за
висит и то, будет ли существовать земная цивилизация или нет. Если эти представления не правильны или не точны, подобное может обернуться катастрофой и гибелью не только цивилизации, но и самой жизни на прекрасной планете, которую, мы называем Землёй. И та
ким образом, из понятий чисто теоретических, представления о природе Вселенной переходят в категорию понятий, от которых зависит будущее цивилизации и будущее жизни на нашей планете.
То, какими будут эти представления, должно волновать не только философов и учёных естественных наук, но и каждого живущего человека. Таким образом, представления о природе Вселенной, если они правильные, могут стать ключом к невиданному прогрессу цивилизации, и, если они неправильные –
привести к гибели и цивилизации..."
Все ф
илософы мечтали описать мир одним уравнением.
Оно должно быть красивое, не очень сложное и небольшое, такое, чтобы помещалось на футболке. Когда его найдут, физики будут довольны. Они назовут это "уравнением Господа" и выйдут на пенсию. До этого, к счастью
, еще далеко.
Призвал Господь Бог архангела Гавриила и говорит: "Гавриил, узнай, как там люди поживают, чем занимаются, не нужно ли им чего". Возвращается Гавриил через некоторое время. Бог его спрашивает: "Что так невесел, Гавриил?" Гавриил отвечает: "Го
споди, лучше бы я не видел, что они там творят. То, что с моралью у них проблемы, это еще ничего. Но у них столько оружия! Ядерного, термоядерного... и им все мало. Господь ему отвечает: "
Ничего, до морали дело еще дойдет, а ты добавь
-
ка им пока еще один н
елинейный член в уравнение Шредингера
".
9
Вы помните фильм "Зеркало для героя" (1987)?
Фантастическая притча Владимира Ивановича Хотиненко по мотивам повести Святослава Рыбаса. Там Сергей уезжает на поезде, несколько месяцев путешествует по стране и возвращ
ается с убеждением, что разорвать петлю можно только одним способом —
ни во что не вмешиваться.
Андрей его позицию не принимает:
...Я раньше как думал –
надо показать, убедить... Так нет: "товарищ Тюкин, вам телеграмма!" От гиппопотама... Ну, я и запил..
. Решил взорвать "Пьяную". В церкви даже был... Мертвые оживают, живые как мертвые. Всё можно, понимаешь? А назавтра –
всё по новой. Устаю очень... Но ты не видел главного: они меняются!
Если изо дня в день повторять одно и то же –
они запоминают. Нужно то
лько ежедневно и с полной отдачей...
Наверное, неслучайно в начале фильма звучит серенада Смита из оперы Бизе "Пертская красавица" (1866) в исполнении Геннадия Михайловича Пищаева:
http://www.you
tube.com/watch?v=32LGJQOsH2w
10
4.
Основная часть. Вычислительные задачи алгебры.
Казалось бы, ничего интересного не может скрываться за утверждениями:
всякая числовая матрица имеет характеристический многочлен...
всякий характеристический многочлен им
еет сопровождающую матрицу...
Но уже при беглом ознакомлении с литературой по обсуждаемой теме бросается в глаза, что авторов книг можно разделить на две группы, исповедующие различные подходы и старательно рассматривающие:
методы перехода от многочленов к сопровождающим матрицам и нахождение их корней;
методы перехода от матриц к характеристическим многочленам и нахождение их корней.
Потому и родилась идея не просто в очередной раз рассмотреть ортогональные полиномы Лагерра, а связать рассмотрение с испо
льзованием матриц
.
4.1.
Ортогональные полиномы Лагерра.
Основные сведения об ортогональных полиномах приведены в брошюре "
Фрагмент 14
. Классические ортогональные полиномы". Здесь сделаем лишь несколько добавлений. Понятие "ортогональные многочлены" был
о введено в конце XIX века в работах П.Л. Чебышёва по непрерывным дробям (о нем написано в п.п. 2.2.2. брошюры "Фрагмент 3. Некое техническое устройство"), позднее развито Александром Александровичем Марковым (1937
—
1994) и Томасом Иоаннесом Стилтьесом (185
6 —
1894). Ортогональные многочлены нашли множество применений. В понятии "ортогональность" нет ничего таинственного. Слово "ортогональность" происходит от греческих слов "прямой, правильный" + "угол", математическое понятие, является обобщением перпендику
лярности
–
"если скалярное произведение двух элементов пространства равно нулю, то они называются ортогональными друг другу". Согласно Беллу, термин "ортогональные" был введён Р. Мёрфи только в 1833
–
1835.
Теория специальных функций, зародившаяся в работах
Эйлера, Гаусса, Лапласа, Якоби, Римана и Чебышева, давно превратилась в классический раздел математики, глубоко проникающий в анализ, теорию функций комплексной переменной, теоретическую и математическую физику, теорию представлений групп, вычислительную математику [3, 4]
. Классические ортогональные полиномы (непрерывного аргумента) —
полиномы Якоби, Лагерра и Эрмита —
образуют самый простой класс специальных функций. В то же время их теория допускает естественное обобщение и на другие специальные функции математической физики. Кроме того, используя ту же логическую схему рассуждений, что и при рассмотрении полиномов Якоби, Лагерра и Эрмита, можно естественным образом ввести определение классических ортогональных полиномов дискретной переменной и построить их теорию как на равномерных, так и неравномерных сетках. Основные свойства полиномов Якоби, Лагерра и Эрмита проще всего вывести исходя из линейного дифференциального уравнения второго порядка
. Так полиномы Лагерра удовлетворяют дифференциальному уравнени
ю
(
являются решением уравнения на полуоси 0
x
):
11
0
ny
'
y
)
1
x
(
"
xy
.
(1)
В общем случае независимая переменная x
может принимать любые комплексные значения. Мы будем рассматривать только действительные пер
еменные. Переменная y
-
решение дифференциального уравнения (1), n
-
порядок рассматриваемого полинома и -
свободный параметр. Полиномы определяются уравнением (1) с точ
ностью до постоянны
х множителей. Поэтому в
практике встречаются три нормировки
задающие
:
L
s
–
стандартизованн
ые полином
ы;
L
e
-
полинома с единицей при старшей степени
;
L
n
–
ортонормированн
ые полином
ы.
Обобщенные полинома Лагерра можно вычислять по формуле:
))
x
exp(
x
(
dx
d
)
x
exp(
x
!
n
1
)
;
x
(
L
n
n
n
n
,
(2)
но обычно используются всевозможные рекуррентные формулы.
4.2.
Две полезные подпрограммы.
Рассмотренные ниже две подпрограммы только на первый взгляд лежат в стороне от обсуждаемой проблемы. На самом же деле они очень даже будут полезны в дальнейшем, как впрочем и все другие подпрограммы из уникального справочника [05]
. В книге [06]
рассматриваются этапы конструирования практических вычислительных алгоритмов на примере решения систем линейных уравнений, приводится описание нескольких паке
тов и библиотек программ, созданных в США и Великобритании и нашедших широкое практическое применение. В приложении [06, стр. 193]
рассмотрено типовое математическое обеспечение задач линейной алгебры:
подпрограммы линейной алгебры
BLAS
;
подпрограммы паке
та LINPACK;
подпрограммы решения проблемы собственных значений
EISPACK
;
п
одпрограммы решения линейных алгебраических уравнен
ий библиотеки IMSL;
п
одпрограммы решения совместных ли
нейных уравнений библиотеки NAG;
п
одпрограммы для метода наименьших квадратов из книги Лоусона и Хэнсона.
Позже с целью замены существующих пакетов программ LINPACK и EISPACK была создана библиотека LAPACK, предназначенная для решения следующих классов задач:
анализ и решение систем линейных алгебраических уравнений;
решение станда
ртной и обобщенной задачи на собственные значения;
декомпозиция особенных матриц.
В приложении E [06, стр. 248]
рассмотрены подпрограммы пакета EISPACK
для решения алгебраической проблемы собственных значений. Пакет EISPACK представляет собой первый приме
р того, как следует воплощать знания экспертов в вычислительных программах для последующего их широкого применения. Подпрограммы пакета созданы на основе набора алгольных процедур, разработанных Уилкинсоном и Райншем в конце 1960
-
х годов [05]
. Использовани
е языка Алгол несколько ограничивало применение процедур. Поэтому был разработан проект их 12
реализации на языке Фортран с целью разработки высокоэффективных версий для различных машин и их детальной проверки. Описание подпрограмм пакет приведено в лекциях п
о информатике [07]
.
4.2.1.
Подпрограмма "MBLN.bas".
Процедура "Balancing a Matrix for Calculation of Eigenvalues and Eigenvectors" описана в справочнике [05, стр. 277]
как алгоритм II.11 "Масштабирование матриц при вычислении собственных значений и век
торов". В
библиотеке
EigenSystem Subroutine Package [07, стр
. 200]
EISPACK: F269
-
2 –
BALANC (July 1975).
A Fortran IV Subroutine to Balance a Real General Matrix.
Поскольку время, затрачиваемое на эту процедуру, мало по сравнению со временем решения всей задачи на собственные значения, то рекомендуется всегда прибегать к процедуре balance.
Кроме того, balance выделяет "изолированные" собственные значения, т.е. такие, которые определяются без вычислений. Причем они определяются абсолютно точно, хотя и могут
быть плохо обусловленными. Ниже приведен текст программы на языке FreeBASIC, иллюстрирующей работу процедуры.
' P R O G R A M "MBLN"
' 22.09.2013
'
----------------------------------------------------
---------
' Подпрограмма уравновешивания действительной матрицы.
'
-------------------------------------------------------------
' ординарная переработка программы из книги:
' Уилкинсон Дж., Райнш Ч. Справочник алгоритмов на языке Алгол
' Стр. 277: Алго
ритм II.11. Масштабирование матриц.
' A() -
исходная и результирующая числовые матрицы;
' N -
порядок квадратной матрицы;
' FPBase -
основание системы счисления с плавающей запятой;
' LOW,IGH -
индикаторы границ индексов;
' SCALE -
массив с ин
формацией о преобразованиях.
#Lang "FB" ' режим совместимости с FreeBASIC
#Include "windows.bi" ' логические переменные
' Объявление подпрограммы уравновешивания матрицы
Declare Sub BALANC(N As Integer, LOW As Integer, IGH As Integer) ' Объяв
ление общих областей данных
Common Shared A() As Double
Common Shared SCALE() As Double
' Основной модуль программы
' Тестовая матрица A() для подпрограммы MBLN при N = 7
Data 7
' 1 2 3 4 5 6 7
Data 6E0, 0, 0, 0
, 0, 1E0, 0
Data 0, 4E0, 0, 3E
-
4, 1E
-
2, 2E
-
2, 1E
-
1
Data 1E0, 1E2, 7E0, 0, 0, -
2E0, 2E1
Data 0, 2E4, 0, 1E0, -
4E2, 3E2, -
4E3
Data -
2E0, -
3E2, 0, 1E
-
2, 2E0, 2E0, 4E1
Data 0, 0, 0, 0, 0, 0, 0
Data 0, 1E1, 0, 4E
-
3, 1E
-
1, -
2E
-
1, 3E0
13
'
Dim As Integer N, LOW, IGH, I, J
Read N
Dim As Integer Cnt(N), ER
Dim As Double A(N,N),SCALE(N)
Dim As String F
For I = 1 To N
For J = 1 To N: Read A(I,J): Next J
Next I
'
Open Cons For Output A
s #2
'Open "MBLN.res" For Output As #2
Print #2,: Print #2, " Subroutine to Balance a Real General Matrix"
'
BALANC(N, LOW, IGH)
'
Print #2, " Transformed the square Matrix A(I,J):": Print #2,
For I = 1 To N
For J = 1 To N: Print #2, Using " ######.# "
; A(I,J);: Next J
Print #2,
Next I
'
Close #2
Sleep
'
Sub BALANC(N As Integer, LOW As Integer, IGH As Integer)
Dim As Integer I, J, K, L, M, JJ, IEXC, FPBase
Dim As Double C, F, G, R, S, B2
Dim As Bool NoConv
FPBase = 10 ' машиннозависимая
константа
B2 = FPBase*FPBase
K = 1: L = N
GoTo L100
' Внутренняя процедура exc() перестановки строк и столбцов
L20: SCALE(M) = J ' коэффициенты масштабирования
If J <> M Then
For I = 1 To L
F = A(I,J): A(I,J
) = A(I,M): A(I,M) = F
Next I
For I = K To N
F = A(J,I): A(J,I) = A(M,I): A(M,I) = F
Next I
EndIf
'
If IEXC = 2 Then GoTo L130
' Поиск строк с изолированными собственными значениями и сдвиг вниз
L1:
L80: If L = 1 Then GoTo L28
0
L = L -
1
14
L100:For J = L To 1 Step -
1
For I = 1 To L
If I = J Then GoTo L110
If A(J,I) <> 0.0 Then GoTo L120
L110:Next I
M = L
IEXC = 1
GoTo L20
L120:Next J
GoTo L140
' Поиск столбцов с изолированными собственными значениями и сдвиг влево
L130:K = K + 1
'
L2:
L140:For J = K To L
'
For I = K To L
If I = J Then GoTo L150
If A(I,J) <> 0 Then GoTo L170
L150:Next I
'
M = K
IEXC = 2
GoTo L20
L170:Next J
' Масштабирование подматрицы в строках от K до L
For I = K To L
SCALE(I) = 1.0
L180:Next I
' Начало итерацион
н
ого цикла
Iteration:
L190:NoConv = FALSE
'
For I = K To L
C = 0.0: R = 0.0
'
For J = K To L
If J = I Then GoTo L200
C = C + Abs(A(J,I))
R = R + A
bs(A(I,J))
L200:Next J
' Защита от исчезновения C или R
If C = 0 Or R = 0 Then GoTo L270
G = R/FPBase: F = 1: S = C + R
L3: If C < G Then F = F*FPBase: C = C*B2: GoTo L3
G = R*FPBase
L4: If C >= G Then F = F/FPBase: C = C/B2: GoTo L4
' Пр
оведение масштабирования
L240:If (C + R)/F < 0.95*S Then
G = 1/F: SCALE(I) = SCALE(I)*F: NoConv = TRUE
For J = K To N: A(I,J) = A(I,J)*G: Next J
15
For J = 1 To L: A(J,I) = A(J,I)*F: Next J
EndIf
'
L270:Next I
'
If NoConv Then GoTo L1
90
'
L280:LOW = K
IGH = L
End Sub
Результат выполнения программы:
' Subroutine to Balance a Real General Matrix
' Transformed the square Matrix A(I,J):
' 7.0 0.1 0.2 0.0 0.0 1.0 -
2.0 ' 0.0 4.0 1.0 3.0 1.
0 0.0 20.0 ' 0.0 1.0 3.0 4.0 1.0 0.0 -
20.0 ' 0.0 2.0 -
4.0 1.0 -
4.0 0.0 30.0 ' 0.0 -
3.0 4.0 1.0 2.0 -
20.0 20.0 ' 0.0 0.0 0.0 0.0 0.0 6.0 1.0 ' 0.0 0.0 0.0 0.0 0.0 0.0 0.0 А вот у автора появилось ощущение, что выполнение этой процедуры не слишком благотворно!
4.2.2.
Подпрограмма "MHQR.bas".
Процедура "The QR Algorithm for Real Hessenberg Matrices" описана в справочнике [05, стр. 316]
как алгоритм II.14 "QR
-
алгоритм для вычисления собственных значений действительной матрицы в форме Хессенберга". В
библиотеке
EigenSystem Subroutine Package [07, стр
. 330]
EISPACK: F286
-
2 -
HQR (July, 1975).
A Fortran IV Subroutine to Determine the Eigenvalues of a Real Upper Hessenberg Matr
ix. При решении проблемы собственных значений, как правило, в классе подобных матриц выбирают матрицу, имеющую по возможности более простой вид. Самым простым видом является -
диагональный, однако не каждая матрица подобна диагональной. Для ряда алгоритмов
решения полной проблемы собственных значений, объем вычислительных операций можно существенно уменьшить, если исходную матрицу А первоначально преобразовать к верхней форме Хессенберга Н, элементы которой hij с номерами i > j+1 равны нулю. Процесс приведе
ния можно выполнить с помощью элементарных устойчивых или унитарных преобразований. Ниже приведен текст программы на языке FreeBASIC, иллюстрирующей работу процедуры.
' P R O G R A M "MHQR"
' 12.08.19
94
'
-------------------------------------------------------------
' Демонстрационная программа для подпрограммы "MHQR"
'
-------------------------------------------------------------
'
' Исследуются два примера из указанной ниже книги.
' A() -
исход
ная квадратная матрица
' N -
размер квадратной матрицы
16
' WRe -
действительные части собственных значений
' WIm -
мнимые части собственных значений
' Cnt -
счетчик итераций
' ER -
код завершения подпрограммы
' EPS -
машинная точность
#L
ang "FB" ' режим совместимости с QBasic
' Объявление подпрограммы вычисления собственных значений
Declare Sub MHQR(N As Integer, A() As Double, WRe() As Double, WIm() As Double, Cnt() As Integer, ER As Integer, EPS As Double) ' Основной модуль про
граммы
'
Dim As Integer NM, I, J
Read NM
Dim As Integer Cnt(NM), ER
Dim As Double A(NM,NM), WRe(NM), WIm(NM), EPS
Dim As String F
'
EPS = 2.22045E
-
016 ' значение машинного "Эпсилон" Open Cons For Output As #2
'Open "MHQR.res" For Output As #2
'
' Пример
1.
Матрица для изу
чения возможности зацикливания.
' Исходные данные со стр. 324:
Data 8
Data 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0
Data 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0
Data 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0
Data 0.0,
0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0
Data 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0
Data 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0
Data 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0
Data 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0
'
For I = 1 To NM
For J = 1 To NM: Read A(I,J): Next J
Next I
'
MHQR(NM, A(), WRe(), WIm(), Cnt(), ER, EPS)
'
If ER <> 0 Then Print #2, " Error executing MHQR!"
Print #2,: Print #2, " Computed Eigenvalues of the Matrix1:": Print #2,
Print #2, " Re(W) Im(W) IN"
Print #2,
F = " ##.######## ##.######## ###"
For I = 1 To NM
Print #2, Using F; WRe(I); WIm(I); Cnt(I)
Next I
'
' Пример
2.
Проверка пр
именимости расщепления матрицы.
' Исходные данные со стр. 324:
17
Data 8
Data 3.0, 2.0, 1.0, 2.0, 1.0, 4.0, 1.0, 2.0 Data 2.0, 1.0, 3.0, 1.0, 2.0, 2.0, 1.0, 4.0 Data 0.0, 3.0, 1.0, 2.0, 1.0, 2.0, 1.0, 3.0 Data 0.0, 0.0, 1.0, 1.0, 2.0, 1.0, 3.0, 1.0 Data 0.0, 0.0, 0.0, 1.E
-
7, 3.0, 1.0, 4.0, 2.0
Data 0.0, 0.0, 0.0, 0.0, 1.E
-
6, 2.0
, 1.0, 4.0
Data 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 2.0, 3.0 Data 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 3.0, 2.0 '
Read NM
For I = 1 To NM
For J = 1 To NM: Read A(I,J): Next J
Next I
'
MHQR(NM, A(), WRe(), WIm(), Cnt(), ER, EPS)
'
If ER <> 0 Then Print #2, " Error e
xecuting MHQR!"
Print #2,: Print #2, " Computed Eigenvalues of the Matrix2:": Print #2,
Print #2, " Re(W) Im(W) IN"
Print #2,
F = " ##.######## ##.######## ###"
For I = 1 To NM
Print #2, Using F; WRe(I); WIm(I); Cnt(
I)
Next I
'
Close #2
Sleep
'
Sub MHQR
(NK As Integer, H() As Double, WRe() As Double, WIm() As Double, NT() As Integer, ER As Integer, EPS As Double)
' 11.09.1987
'
-------------------------------------------
------------------
' Вычисление собственных чисел матрицы Хессенберга
'
-------------------------------------------------------------
' Ординарная переработка программы из книги:
' Уилкинсон Дж., Райнш Ч., Справочник алгоритмов на языке
' Алгол, лине
йная алгебра, М., Машиностроение, 1976.
' 390 с., ил. (QR
-
алгоритм для матрицы, заданной в верхней
' форме Хессенберга, с. 316).
' H(N,N) -
исходная матрица в форме Хессенберга;
' WRe(N) -
действительные части собственных значений;
' WIm(N) -
мнимы
е части собственных значений.
'
Dim As Integer N, NA, ITS, L, K, I, J, M
Dim As Double T, X, S, Y, W, Z, P, Q, Q1, R, R1
ER = 0: T = 0: N = NK
'
NextW:If N <= 0 Then GoTo Fin
ITS = 0: NA = N -
1
'
18
NextIt:
If N >= 2 Then
For L = N To 2 Step -
1
X = Abs(H(L,L
-
1))
S = EPS*(Abs(H(L
-
1,L
-
1)) + Abs(H(L,L)))
If X <= S Then GoTo Cont1
Next L
EndIf
L = 1
Cont1:
X = H(N,N)
If L = N Then GoTo L400
Y = H(NA,NA): W = H(N,NA)*H(NA,N)
If L = NA Then GoTo L500
If IT
S >= 30 Then GoTo L600 ' Ошибка
> 30 итераций
If ITS = 10 Or ITS = 20 Then
' Формирование дополнительного смещения
T = T + X
For I = 1 To N
H(I,I) = H(I,I) -
X
Next I
S = Abs(H(N,NA)) + Abs(H(NA,N
-
2))
X = 0.75*S: Y
= X: W = -
0.4375*S*S
EndIf
ITS = ITS + 1
If N -
2 < L Then GoTo Cont2
' Поиск двух соседних малых элементов
For M = N
-
2 To L Step -
1
Z = H(M,M): R = X -
Z: S = Y -
Z
P = (R*S -
W)/H(M+1,M) + H(M,M+1)
Q = H(M+1,M+1) -
Z -
R -
S
R =
H(M+2,M+1)
S = Abs(P) + Abs(Q) + Abs(R)
P = P/S: Q = Q/S: R = R/S
If M = L Then GoTo Cont2
R1 = Abs(H(M,M
-
1))*(Abs(Q) + Abs(R))
Q1 = EPS*ABS(P)*(ABS(H(M
-
1,M
-
1)) + ABS(Z) + ABS(H(M+1,M+1)))
If R1 <= Q1 Then GoTo Cont2
Next M
Cont2:
If M+2 <= N Then
For I = M+2 To N: H(I,I
-
2) = 0: Next I
EndIf
If M+3 <= N Then
For I = M+3 To N: H(I,I
-
3) = 0: Next I
EndIf
' Двойной шаг QR
-
алгоритма над строками и столбцами
If M > NA Then GoTo NextIt
For K = M To NA
If K <> M Th
en
P = H(K,K
-
1): Q = H(K+1,K
-
1)
19
R = 0
If K <> NA Then R = H(K+2,K
-
1)
X = Abs(P) + Abs(Q) + Abs(R)
If X = 0 Then GoTo Cont3
P = P/X: Q = Q/X: R = R/X
EndIf
S = SQR(P*P + Q*Q + R*R)
If P < 0 Then S = -
S
If K = M And L = M Then G
oTo L8
H(K,K
-
1) = -
S*X
If K = M Then H(K,K
-
1) = -
H(K,K
-
1)
L8:P = P + S: X = P/S: Y = Q/S: Z = R/S
Q = Q/P: R = R/P
' Изменение строки
If K > N Then GoTo L24
For J = K To N
P = H(K,J) + Q*H(K+1,J)
If K = NA Then GoTo L11
P = P + R*H
(K+2,J)
H(K+2,J) = H(K+2,J) -
P*Z
L11:H(K+1,J) = H(K+1,J) -
P*Y
H(K,J) = H(K,J) -
P*X
Next J
L24:J = N
If K+3 < N Then J = K + 3
If L > J Then GoTo Cont3
' Изменение
столбца
For I = L To J
P = X*H(I,K) + Y*H(I,K+1)
If K = NA Then G
oTo L13
P = P + Z*H(I,K+2)
H(I,K+2) = H(I,K+2) -
P*R
L13:H(I,K+1) = H(I,K+1) -
P*Q
H(I,K) = H(I,K) -
P
Next I
Cont3:
Next K
GoTo NextIt
onew:
L400:WRe(N) = X + T: WIm(N) = 0
NT(N) = ITS
N = NA
GoTo NextW
iwow:
L500: P = 0.5*(Y -
X)
Q = P*P + W
Y = Sqr(Abs(Q))
NT(N) = -
ITS: NT(NA) = ITS: X = X + T
'
If Q > 0 Then
' Пара действительных значений
If P < 0 Then Y = -
Y
20
Y = P + Y
WRe(NA) = X + Y: WIm(NA) = 0
WRe(N) = X -
W/Y: WIm(N) = 0
Else
' Пара комплек
сно
-
сопряженных собственных значений
WRe(NA) = X + P: WIm(NA) = Y
WRe(N) = X + P: WIm(N) = -
Y
EndIf
N = N -
2
GoTo NextW
'
Fin:
L600:ER = 1
End Sub
Результат
выполнения
программы
:
' Computed Eigenvalues of the Matrix1:
' Re(W) Im(W) IN
' -
1.00000000 0.00000000 0
' -
0.70710678 0.70710678 3
' -
0.70710678 -
0.70710678 -
3
' -
0.00000000 1.00000000 3
' -
0.00000000 -
1.00000000 -
3
' 1.00000000 0.00000000 2
' 0.70710678 0.70710678 17
' 0.70710678 -
0.70710678 -
17
' Computed Eigenvalues of the Matrix2:
' Re(W) Im(W) IN
' -
2.42951247 0.00000000 0
' 1.94495965 0.00000000 2
' 0.86636719 0.00000000 -
2
' 2.99999999 0.00000000 2
' 5.61818555 0.00000000 -
2
' 0.35424730 0.00000000 0
' 0.00000117 0.00000000 1
' 5.64575163
0.00000000 7
Результат выполнения совпадает с приведенным в книге!
Процедура MHQR может быть использована для отыскания всех собственных значений действительной матрицы, заданной в верхней форме Хессенберга. В случае произвольной действитель
ной матрицы ее следует преобразовать к такому виду с помощью одной из процедур elmhes, dirhes или orthes [5, стр. 298]
.
21
4.3.
Проблема собственных значений.
Когда алгебра была еще скудна теоремами, следующее утверждение получило название основной теоре
мы алгебры
: "Многочлен степени n с комплексными коэффициентами имеет ровно n корней (с учетом их кратностей)". Впервые это утверждение сформулировал Альбер де Жирар в 1629, но он даже не пытался его доказывать. Первым осознал необходимость доказательства о
сновной теоремы алгебры Даламбер, но его доказательство (1746) не было признано убедительным. Свои доказательства предложили Эйлер (1749), Фонсене (1759) и Лагранж (1771), но и эти доказательства были небезупречны. Первым удовлетворительное доказательство основной теоремы алгебры получил Гаусс, который привел три разных доказательства (1799, 1815 и 1816), а в 1845 опубликовал еще и уточненную версию своего первого доказательства.
Для поиска корней полиномов разработаны особые методы. В этом случае могут бы
ть найдены все корни. После того как один из корней полинома найден, степень полинома может быть понижена, после чего поиск корня повторяется. Способ нахождения корней многочленов первой и второй степени, то есть способ решения линейных и квадратных уравне
ний, был известен ещё в древнем мире. То, что корни общего уравнения пятой степени и выше не выражаются при помощи рациональных функций и радикалов от коэффициентов, было доказано математиком Абелем в 1824 году [08]
. Норвежский математик Нильс Генрик Абель
(1802
–
1829) доказал следующую теорему: Общее алгебраическое уравнение с одним неизвестным степени выше 4
-
й неразрешимо в радикалах, т.е. не существует формулы, выражающей корни общего уравнения степени выше 4
-
й через коэффициенты с помощью операций сложен
ия, вычитания, умножения, деления, возведения в натуральную степень и извлечения корней натуральной степени. Из книги [08]
читатели узнают, как решать алгебраические уравнения 3
-
й и 4
-
й степени с одним неизвестным и почему для решения уравнений более высок
ой степени не существует общих формул (в радикалах). При этом будет возможность познакомиться с двумя очень важными разделами современной математики —
теорией групп и теорией функций комплексного переменного.
Но эта теорема совсем не означает, что корни т
аких уравнений не могут быть найдены с необходимой точностью. Довольно часто применимы на практике метод Лагерра (Laguerre's Method), как самый "быстрый" и метод Собственных Значений Матриц (Eigenvalues of Companion Matrix), как самый "точный". Чем
-
то сред
ним между указанными методами является метод Дженкинса
-
Трауба (Jenkins
-
Traub Algorithm, 1970). Но вернемся к ортогональным полиномам Лагерра. В брошюре "Фрагмент 14. Классические ортогональные полиномы" для вычисления комплексных корней полинома Лагерра ис
пользовалась подпрограмма PLRT()
от 12.08.1994, реализующая метод последовательных приближений Берстоу
-
Ньютона. Попробуем что
-
нибудь новенькое...
4.3.1
.
От многочлена к сопровождающей матрице.
Метод перехода от многочленов к сопровождающим матрицам и н
ахождение их корней проиллюстрируем расчетами для ортогональных полиномов Лагерра 7
-
ой степени. Метод собственных значений сопровождающих матриц (Eigenvalues of 22
Companion Matrix) -
очень точный метод вычисления нулей многочленов. Напомним, что многочлен (п
олином) одной переменной называется унитарным или приведенным, если его коэффициент при старшей степени равен единице. Везде, где это возможно, будем использовать "приведенные" полиномы. Ниже приведена программа на языке FreeBASIC и результаты ее выполнени
я с необходимыми комментариями.
' P R O G R A M "MPOL"
' 02.09.2013
'
-------------------------------------------------------------
' Всякое про полиномы Лагерра степени 7 и квадратной матрице
'
--
-----------------------------------------------------------
'
' NP -
степень полинома;
' AP() -
коэффициенты полинома;
' AM() -
коэффициенты сопровождающей матрицы;
' AT() -
коэффициенты транспонированной матрицы.
'
#Lang "FB" ' реж
им совместимости с FreeBASIC
' Объявление подпрограммы вычисления коэффициентов полинома:
Declare Sub CPLG(N As Integer, Alph As Single, CP() As Double, ER As Integer)
' Объявление подпрограммы вычисления собственных значений:
Declare Sub MHQR(N As Integ
er, A() As Double, WRe() As Double, WIm() As Double, Cnt() As Integer, ER As Integer, EPS As Double)
' Объявление подпрограммы вычисления методом Горнера Declare Sub PGorn(N As Integer, A() As Double, ByRef X As Double, ByRef Y As Double, ByRef P As Doub
le)
'
' Основной модуль программы
'
' Массивы переменных (нижний индекс по умолчанию равен 0):
Dim As Integer NP
NP = 7 ' заданная степень полинома Лагерра
Dim As Integer N, I, J, K, ER, Cnt(1 To 10)
Dim As Single Alph
Dim As Double EPS, Sum,
X, X1, Y, P
Dim As Double AP(0 To 10), AP0, A1(0 To 10), WRe(1 To 10), WIm(1 To 10)
Dim As Double AM(1 To 10,1 To 10), AT(1 To 10,1 To 10)
Dim As String F
Open Cons For Output As #2
'Open "MPOL.res" For Output As #2
' Исходные данные
EPS = 2.22045E
-
016 ' значение машинного "Эпсилон"
Alph = 0.0 ' параметр полинома
Print #2,
Print #2, Using " Parameter polynomial Alph = ##.##"; Alph
Print #2, Using " Coefficients of the Polynomial N = ##"; NP
CPLG(NP, Alph, AP(), ER)
23
Print #2,: Print #2, " K
CPLG"
For K = 0 To NP
Print #2, Using " ## ##.######^^^^"; K; AP(K)
Next K
Print #2,
' Приведение к единице старшего коэффициента
AP0 = AP(0)
For J = 0 To NP: A1(J) = AP(J)/AP0: Next J
' Формирование сопровождающей матрицы
For I = 1 To NP
For J
= 1 To NP
AM(J,I) = 0
If I <> NP And J = I+1 Then AM(J,I) = 1
If I = NP Then AM(J,I) = -
A1(NP
-
J+1)
Next J
Next I
' Транспонирование сопровождающей матрицы
For I = 1 To NP
For J = 1 To NP: AT(I,J) = AM(J,I): Next J
Next I
' Вывод полученных результатов
P
rint #2,
Print #2, " Accompanying the Matrix Polynomial AM(I,J):": Print #2,
For I = 1 To NP
For J = 1 To NP: Print #2, Using " ######.# "; AM(I,J);: Next J
Print #2,
Next I
Print #2,
Print #2, " Transposed accompanying Matrix AT(I,J):": Print #2,
F
or I = 1 To NP
For J = 1 To NP: Print #2, Using " ######.# "; AT(I,J);: Next J
Print #2,
Next I
Print #2,
' Вычисление собственных чисел матрицы (корней полинома)
Print #2,: Print #2, " Computed Eigenvalues of the Matrix": Print #2,
MHQR(NP, AM(), WRe()
, WIm(), Cnt(), ER, EPS)
Print #2, " Re(W) Im(W) Cnt"
Print #2,
F = " ##.######## ##.######## ###"
For I = 1 To NP
Print #2, Using F; WRe(I); WIm(I); Cnt(I)
Next I
Print #2,
' Прямая проверка корней
Print #2, " Root
s calculated according to the Program:": Print #2,
24
For I = 1 To NP ' перебор всех корней
Sum = 0
For K = 0 To NP ' перебор всех коэффициентов
Sum = Sum + AP(NP
-
K)*WRe(I)^K
Next K
Print #2, Using " I = ### Sum = ##.####^^^^"; I; Sum
Next I
' Из книги Крылов В.И. Приближенное вычисление интегралов:
Data 0.193043676560, 1.026664895339, 2.567876744951, 4.900353084526
Data 8.182153444563, 12.734180291798, 19.395727862263
For I = 1 To NP ' перебор
всех
коэффициентов
Read WRe(I)
Next I
'
Print #2,: Print #2
, " Roots in a book published Krylov:": Print #2,
For I = 1 To NP ' перебор всех корней
Sum = 0
For K = 0 To NP ' перебор всех коэффициентов
Sum = Sum + AP(NP
-
K)*WRe(I)^K
Next K
Print #2, Using " I = ### Sum = ##.####^^^^"; I; Sum
Next I
'
Print #2,
Print #2,: Print #2, " Value and Derivative of a Polynomial"
Print #2, " X YG PG"
For X = 0 To 5.0 Step 0.5
PGorn(NP, AP(), X, Y, P)
Print #2, Using " #.## ###.####### ###.#######"; X; Y; P
Next X
' Попытка уточнения корней (нулей)
полинома:
Print #2,: Print #2, " Roots Sophisticated in this Program:": Print #2,
For I = 1 To NP ' перебор всех корней
X1 = WRe(I) ' значение первого приближения
'Print #2, Using " I = ### ##.########"; I; X1
J = 0
Do ' выполнять
цикл
X = X1
PGorn(NP, AP(), X, Y, P)
X1 = X -
0.1*Y/P ' итерационный метод Ньютона
J = J + 1
If J > 100 Then Exit Do
Loop While Abs(0.1*Y/P) > EPS
'Print #2, Using " I = ### ##.####^^^^"; I; WRe(I)
-
X
WRe(I) = X1
Next I
'
For I = 1 To NP ' перебор в
сех корней
Sum = 0
25
For K = 0 To NP ' перебор всех коэффициентов
Sum = Sum + AP(NP
-
K)*WRe(I)^K
Next K
Print #2, Using " I = ### Sum = ##.####^^^^"; I; Sum
Next I
Close #2
Sleep
'
Sub CPLG
(N As Integer, Alph As Single, CP() As Double, ER As Integer)
' Вычисление коэффициентов полинома Лагерра
'
-------------------------------------------------------------
'
' Подпрограмма соответствует явному выражению для полинома
' Лагерра, из книги: Суетин П.К. Классические ортогональные
' многочлены. М., Наука, 1979, 416 с. (стр. 220)
' N -
степень стандартизованного полинома Лагерра;
' Alph -
значение параметра полинома (Alph > -
1);
' CP() -
вектор коэффициентов полинома Лагерра;
' ER -
код завершения подпрограммы.
'
D
im As Integer I, M, Sign
Dim As Double X, S
Dim As Double FCN, FCL, FCM
ER = 0: CP(0) = 1 ' нет ошибок; свободный член
If Alph > -
1 And N > 0 And N <= 12 Then ' проверка
ограничений
'
FCN = 1
For I = 1 To N ' факториал
N
FCN = CDbl(I)*FCN
Next
I
'
Sign = 1
For M = 0 To N
'
FCM = 1
For I = 1 To M ' факториал
M
FCM = CDbl(I)*FCM ' преобразовать к двойной точности
Next I
'
FCL = 1
For I = 1 To N
-
M ' факториал
N
-
M
FCL = CDbl(I)*FCL ' преобразовать к двойной точности
Next I
'
X = N + Alph: S = 1
For I = 0 To N
-
M
-
1 ' члены
числителя
S = S*(X
-
CDbl(I))
Next I
CP(N
-
M) = Sign*S/(FCL*FCM)
Sign = -
Sign ' изменить знак
26
Next M
EndIf
End Sub
'
Sub MHQR
(NK As Integer, H() As Double, WRe() As Double, WIm() As Double, Cnt() As Integer, ER As Integer, EPS As Double)
' Вычисление собственных чисел матрицы Хессенберга
Текст подпрограммы опущен!
См. ранее пункт 4.2.2. Подпрограмма
"MHQR.bas".
End Sub
'
Sub PGorn
(NN As Integer, AA() As Double, ByRef XX A
s Double, ByRef YY As Double, ByRef PP As Double)
' Модуль вычисления значения и производной по схеме Горнера
'
-------------------------------------------------------------
' Полином по убывающим степеням аргумента:
' Y = A(0)*X^N+A(1)*X^(N
-
1)+.
..+A(N
-
1)*X+A(N)
' NN -
степень полинома;
' AA -
массив коэффициентов полинома;
' XX -
значение аргумента;
' YY -
значение полинома (в точке "ХX");
' PP -
значение производной (в точке "ХX").
'
Dim KK As Integer
PP = 0
YY = AA(0)
For KK = 1 To NN
PP = PP*XX + YY
YY = YY*XX + AA(KK)
Next KK
End Sub
Результаты
вычислений
:
' Parameter polynomial Alph = 0.00
' Coefficients of the Polynomial N = 7
' K CPLG
' 0 -
1.984127E
-
04
' 1 9.722222E
-
03
' 2 -
1.750000E
-
01
' 3
1.458333E+00
' 4 -
5.833333E+00
' 5 1.050000E+01
' 6 -
7.000000E+00
' 7 1.000000E+00
Если необходимо производить вывод коэффициентов полинома с единицей при старшем члене, то при выводе следует значения разделить на значение AP(0), как
это сделано в программе далее. Ниже распечатана сопровождающая матрица полинома и транспонированная сопровождающая матрица.
27
' Accompanying the Matrix Polynomial AM(I,J):
' 0.0 0.0 0.0 0.0 0.0 0.0
5040.0
' 1.0 0.0 0.0 0.0 0.0 0.0
-
35280.0
' 0.0 1.0 0.0 0.0 0.0 0.0 52920.0
' 0.0 0.0 1.0 0.0 0.0 0.0 -
29400.0
' 0.0 0.0 0.0 1.0 0.0 0.0 7350.0
' 0.0 0.0 0.0 0.0 1.0 0.0 -
882.0
' 0.0 0.0 0.0 0.0 0.0 1.0 49.0
' Transposed accompanying Matrix AT(I,J):
' 0.0 1.0 0.0 0.0 0.0 0.0 0.0
' 0.0 0.0 1.0 0.0 0.0 0.0 0.0
' 0.0 0.0
0.0 1.0 0.0 0.0 0.0
' 0.0 0.0 0.0 0.0 1.0 0.0 0.0
' 0.0 0.0 0.0 0.0 0.0 1.0 0.0
' 0.0 0.0 0.0 0.0 0.0 0.0 1.0
' 5040.0 -
35280.0 52920.0 -
29400.0 7
350.0 -
882.0 49.0
После выполнения подпрограммы MHQR на печать выводятся комплексные собственные числа сопровождающей матрицы.
' Computed Eigenvalues of the Matrix
' Re(W) Im(W) Cnt
' 0.19304368 0.00000000 0
' 1.02666490 0.00000000 0
' 2.56787674 0.00000000 3
' 4.90035308 0.00000000 1
' 8.18215344 0.00000000 3
' 12.73418029 0.00000000 1
' 19.39572786 0.00000000 7
Они одновременно являются корнями исследуемого полинома Лагерра степени 7 с параметром Alpha = 0. Для оценки результатов выводятся значения полинома при корнях, рассчитанных по программе.
' Roots calculated according to the Program:
' I = 1 Su
m = 2.5422E
-
15
' I = 2 Sum = 2.4241E
-
15
' I = 3 Sum = 3.9094E
-
14
' I = 4 Sum = 4.4025E
-
13
' I = 5 Sum = 2.8502E
-
12
' I = 6 Sum = 1.3763E
-
11
' I = 7 Sum = 9.1845E
-
11
Для сравнения "качества" выводятся значения полинома при корнях, взятых из книги Крылов В.И. Приближенное вычисление интегралов:
' Roots in a book published Krylov:
' I = 1 Sum = 1.2893E
-
12
' I = 2 Sum = -
2.9199E
-
13
' I = 3 Sum = -
4.1075E
-
13
' I = 4 Sum = -
1.5164E
-
12
28
' I = 5 Sum =
-
1.1171E
-
12
' I = 6 Sum = 1.6585E
-
11
' I = 7 Sum = -
5.6409E
-
10
Для наглядности распечатываются значения и производные полинома в точках из интервала 0 <= X <= 5 с шагом 0.5.
' Value and Derivative of a Polynomial
' X YG PG
' 0.00 1.0000000 -
7.0000000
' 0.50 -
0.5183392 -
0.1987196
' 1.00 -
0.0404762 1.5152778
' 1.50 0.5987584 0.8099609
' 2.00 0.6634921 -
0.5555556
' 2.50 0.1079567 -
1.5513238
' 3.00 -
0.7464286 -
1.7125000
' 3.50 -
1.45
58485 -
0.9909071
' 4.00 -
1.6285714 0.3777778
' 4.50 -
1.0369280 1.9912109
' 5.00 0.3253968 3.3819444
Для уточнения корней производится попытка уточнения корней Крылова методом Ньютона. Очевидно, это улучшает корни, делая их близкими к расчетным.
' Roots Sophisticated in this Program:
' I = 1 Sum = 6.6896E
-
15
' I = 2 Sum = -
3.1922E
-
15
' I = 3 Sum = -
5.2780E
-
15
' I = 4 Sum = -
4.2158E
-
14
' I = 5 Sum = -
6.2206E
-
13
' I = 6 Sum = 1.8217E
-
12
' I = 7 Sum =
-
3.6991E
-
11
Пожалуй, можно считать опыт удавшимся!
4.3.2
.
От матрицы к характеристическому многочлену.
Здесь будет рассмотрен методы перехода от матриц к характеристическим многочленам и нахождение их корней. Удивительно, но в доступных библиотеках
программ не смог обнаружить подпрограмму перехода от матрицы общего вида к характеристическому многочлену. Найти ее удалось только в учебном пособии по численным методам решения задач линейной алгебры [10]
, где детально представлено решение задачи на собс
твенные значения методом Данилевского, сущность которого состоит в приведении матрицы с помощью подобных преобразований к подобной ей матрице Фробениуса.
29
' P R O G R A M "MFRB"
' 21.10.2013
'
------
-------------------------------------------------------
' Приведение матрицы общего вида к форме Фробениуса
'
-------------------------------------------------------------
' ординарная переработка программы из книги:
' Козин Р.Г. Алгоритмы численных мет
одов линейной алгебры. М.:
' НИЯУ МИФИ, 2012. 141 с.
' Стр. 87: приведены листинги программы, вычисляющей матрицу
' Фробениуса, и процедуры выбора ведущего элемента.
'
#Lang "FB" ' режим совместимости с FreeBASIC
' Объявление подпрограммы вычисления
матрицы Фробениуса
Declare Sub Froben(N As Integer)
' Объявление процедуры выбора ведущего элемента
Declare Function LdElem(N As Integer, K As Integer) As Integer
' Объявление общей области данных
Common Shared A() As Double
'
Dim As Integer N, I, J, Nm
at
' Основной модуль программы
'
Data 4
' Тестовая матрица из работы Козина:
Data 1, 2, 3, 4
Data 5, 6, 7, 8
Data 0, 0, 9, 10
Data 0, 0, 11, 12
' Тестовая матрица из работы Леверье:
Data -
5.509882, 1.870086, 0.422908, 0.008814
Data 0.287865, -
11.811654, 5.711900, 0.058717
Data 0.049099, 4.308033, -
12.970687, 0.229326
Data 0.006235, 0.269851, 1.397369, -
17.596207
'
Read N
Dim As Double A(N,N)
Open Cons For Output As #2
'Open "MFRB.res" For Output As #2
Print #2,: Print #
2, " Subroutine for computing the matrix Frobenius"
' Ввод элементов исходной матрицы (по строкам)
For Nmat = 1 To 2
For I = 1 To N
For J = 1 To N: Read A(I,J): Next J
Next I
'
Froben(N)
' Вывод
полученной
матрицы
Print #2,: Print #2, Using " Resulting matrix ## Frobenius A(i,j):"; Nmat: Print #2,
For I = 1 To N
For J = 1 To N
30
Print #2, Using " ######.######"; A(I,J);
Next J
Print #2,
Next I
Print #2,
Next Nmat
Close #2
Sleep
'
Sub Froben(N As Integer
)
Dim As Integer K, I, J, M
Dim As Double MA(N), Amax, Tmp
For K = 1 To N
If LdElem(N, K) = 1 Then Exit For
' Сохранение обрабатываемой строки N
-
(K
-
1)
For J = 1 To N
MA(J) = A(N
-
(K
-
1),J)
Next J
' Деление обрабатываемого столбца (N
-
K) на A
(N
-
(K
-
1),N
-
K) For I = 1 To N
-
(K
-
1)
A(I,N
-
K) = A(I,N
-
K)/A(N
-
(K
-
1),N
-
K)
Next I
' Обнуление остальных элементов обрабатываемой строки
For I = 1 To N
-
(K
-
1)
For J = 1 To N
If J <> N
-
K Then
A(I,J) = A(I,J) -
A(I,N
-
K)*A(N
-
(K
-
1),
J)
EndIf
Next J
Next I
' Умножение промежуточной матрицы A на обратную матрицу М
For J = 1 To N
Tmp = 0
For M = 1 To N
Tmp = Tmp + MA(M)*A(M,J)
Next M
A(N
-
K,J) = Tmp
Next J
Next K
End Sub
'
Function LdEle
m(N As Integer, K As Integer) As Integer
Dim As Integer Jmax, I, J, cd, code
Dim As Double Amax, Tmp
cd = 0
' LdElem = leading element
Amax = Abs(A(N
-
(K
-
1),N
-
K))
Jmax = N -
K
' Поиск ведущего элемента в строке
For J = 1 To N
-
1
-
K
Tmp = Abs(A(N
-
(K
-
1),J))
If Tmp > Amax Then
31
Amax = Tmp
Jmax = J
EndIf
Next J
'
If K <> Jmax Then
For I = 1 To N
Tmp = A(I,N
-
K)
A(I,N
-
K) = A(I,Jmax)
A(I,Jmax) = Tmp
Next I
'
For J = 1 To N
Tmp = A(N
-
K,J)
A(N
-
K,J) = A(Jmax,J)
A(Jmax,J) = Tmp
Next J
End
If
' Установка признака малости максимального элемента
If Amax < 1E
-
30 Then cd = 1
LdElem = cd
End Function
Результат выполнения программы:
' Subroutine for computing the matrix Frobenius
' Resulting matrix 1 Frobenius A(i,j):
' 1.000000 2.000000 0.272727 0.727273
' 5.000000 6.000000 0.636364 0.363636
' 0.000000 0.000000 21.000000 2.000000
' 0.000000 0.000000 1.000000 0.000000
' Resulting matrix 2 Frobenius A(i,j):
' -
47.888430 -
797.278765 -
5349.455515 -
12296.550566
' 1.000000 0.000000 0.000000 0.000000
' 0.000000 1.000000 0.000000 0.000000
' 0.000000 0.000000 1.000000 0.000000
Александр Михайлович Дани
левский
(1906
-
1941), сотрудник Харьковского Электротехнического Института предложил элегантный и весьма эффективный метод вычисления коэффициентов характеристического полинома в 1937 году. Метод основан на приведении определителя уравнения частот к так наз
ываемой форме Фробениуса. Это приведение может быть выполнено самыми простыми средствами с достаточной для практики точностью и при большом числе степеней свободы требует в полтора раза меньше операций, чем знаменитый метод А.Н.Крылова. Данилевский погиб п
ри оккупации Харькова в годы Великой Отечественной войны.
Фердинанд Георг Фробениус
(нем. Ferdinand Georg Frobenius; 1849 —
1917) —
немецкий математик, получивший известность за вклад в теорию эллиптических функций, дифференциальных уравнений и теории груп
п. Он также был первым, кто ввёл понятие рациональной аппроксимации функций (ныне известный как аппроксимации Паде), и дал первое полное доказательство теоремы Гамильтона —
Кэли. Также он внёс свой вклад в определение дифференциально
-
геометрических 32
объекто
в современной математической физики, известных ныне как многообразия Фробениуса. К нашей сегодняшней тематике относится понятие "фробениусова нормальная форма линейного оператора"... Вопрос о нормальных формах матриц рассматривается в книге [11]
, достаточн
о детально представлены жорданова и фробениусовы нормальные формы.
Урбен Жан Жозеф Леверье
(фр. Urbain Jean Joseph Le Verrier, 1811
-
1877) -
французский математик, занимавшийся небесной механикой, большую часть своей жизни проработавший в Парижской обсерва
тории. В то время астрономов всего мира весьма сильно занимал вопрос о возмущениях планетной орбиты Урана. В 1845 и 1846 годах Леверье представил в Парижскую академию наук вычисления и установил предполагаемые элементы орбиты возмущающего тела. Иоганн Галл
е, который был тогда адъюнктом и наблюдателем в Берлинской обсерватории и имел в своем распоряжении хорошие звездные карты, получив письмо от Леверье 23 сентября 1846 года, немедленно начал наблюдения и в ту же ночь нашел неизвестную планету, возмущающую д
вижение Урана, весьма близко от места, указанного Леверье. Новая планета получила имя Нептун.
Используемая в программе тестовая матрица происходит из работы Леверье:
Leverrier U.J.J. Sur les variations seculafres des elements des orbites. J. Math., 1840.
и стала весьма популярной, особенно в книгах посвященных приведению векового уравнения к полиномиальному виду. В замечательной книге [12, 263]
при рассмотрении метода Алексея Николаевича Крылова (1863
–
1945) для нахождения всех собственных значений матрицы эта матрица приведена в качестве первого примера при определении коэффициентов характеристического полинома для матрицы. Напомним, что собственными значениями некоторой матрицы называются корни ее характеристического полинома. Далее рассматриваются метод Х
ессенберга, метод Самуэльсона, метод Данилевского, метод Леверье, эскалаторный метод... Согласитесь, приличная коллекция методов приведения к полиномиальному виду векового уравнения. Так в таблице IV.8 на [12, 288]
в правом столбце приведены коэффициенты х
арактеристического полинома с единицей при старшей степени:
t^4 + 47.888430t^3 + 797.27
875t^2 + 5349.4555t + 12296.550,
которые хорошо совпадает с коэффициентами, рассчитанными по программе.
Большинство методов, дающих решение полной проблемы собственны
х значений, включает предварительное вычисление коэффициентов характеристического полинома, которое осуществляется теми или иными средствами, минуя вычисление многочисленных миноров. Собственные значения вычисляются затем по какому
-
либо методу для приближе
нного вычисления корней полинома. Одним из лучших способов приближенного вычисления корней полинома признается способ Ньютона...
33
Использованная литература
01.
Бабенко К.И.
Основы численного анализа. Москва
-
Ижевск: РХД, 2002. 848 с.
02.
Крылов В.И
., Бобков В.В., Монастырный П.И.
Начала теории вычислительных методов. Линейная алгебра и нелинейные уравнения. Мн.: Наука и техника, 1985. 280 с.
03.
Суетин П.К.
Классические ортогональные многочлены. 3
-
е изд., М.: ФизМатЛит, 2005. 480 с.
04.
Никифоро
в А.Ф., Уваров В.Б.
Основы теории специальных функций. М.: Наука, 1974. 304 с.
05.
Уилкинсон Дж.Х., Райнш С.
Справочник алгоритмов на языке Алгол. Линейная алгебра. М.: Машиностроение, 1976. 390 с.
06.
Райс Дж.
Матричные вычисления и математическое обе
спечение. М
.: Мир
, 1984. 264 с
.
07.
Lecture Notes in Computer Science
(LNCS #6). 1976. 560 p.
Matrix Eigensystem Routines -
EISPACK Guide.
B.Smith, J.Boyle, J.Dongarra, B.Garbow, Y.Ikebe, V.Klema, C.Moler
08.
Алексеев В.Б.
Теорема Абеля в задачах и реш
ениях. М.: МЦНМО, 2001. 192 с.
09.
Крылов В.И.
Приближенное вычисление интегралов. М.: Наука, 1967. 500 с.
10.
Козин Р.Г.
Алгоритмы численных методов линейной алгебры. М.: НИЯУ МИФИ, 2012. 141 с.
11.
Милованов М.В., Толкачев, Р.И., Тышкевич Р.И., Фед
енко А.С.
Алгебра и аналитическая геометрия. Часть 2. Мн.: Выш. шк., 1987. 269 с.
12.
Фаддеев Д.К., Фаддеева В.Н.
Вычислительные методы линейной алгебры. М.: Наука, 1960. 655 с.
34
Как всегда расстаюсь не без сожаления
-
материал этой брошюры просто не
объятный, а тема весьма актуальная. Тут не брошюры, а книги писать нужно. И ведь что интересно –
пишут много, но, видимо, не те и не для тех! По пальцам можно пересчитать доступные издания для подготовки инженеров...
Но сегодня лимит времени и объема бр
ошюры исчерпан!
На этом позвольте закруглиться.
До новых встреч!
Чтобы не тратить попусту время, очень прошу, напишите, какую тему (программу) желательно рассмотреть в следующей брошюре. Обещаю, что тема будет назначена по двум полученным сообщениям,
содержащим близкие по смыслу запросы.
Вот адрес:
eugene_r@mail.ru
Документ
Категория
Информатика
Просмотров
181
Размер файла
988 Кб
Теги
многочлен, алгоритм, матрица, freebasic, программирование
1/--страниц
Пожаловаться на содержимое документа