close

Вход

Забыли?

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

?

FreeBASIC6

код для вставкиСкачать
Пособие для начинающих программировать на языке высокого уровня FreeBASIC. Работа состоит из нескольких глав. В главе 6 "Линейная алгебра" – продолжение пособия, которое будет интересно для учащихся школ, студентов институтов, а также преподавателей
 1
Евгений Рыжов, инженер
Программирование на языке FreeBASIC
пособие для начинающих
Пособие для начинающих программировать на языке высокого уровня FreeBASIC.
Работа состоит из нескольких глав. В главе 6 "Линейная алгебра" –
продолжение
пособия, которое
будет интересно для учащихся школ, студентов институтов,
а также преподавателей.
Фрагмент 6. Линейная алгебра
"Простите, мне очень тяжело читать по своим лекциям...
Иногда мне кажется, что это записки сумасшедшего..."
Профессор математики
1.
Из записок сумасшедшего
Вы думаете, как устроена Вселенная?
Вот откровение: Вселенная состоит из определенного вида голограмм, называемых "Эгрегоры". Они являются сложными энергоинформационными образованиями, рожденными и поддерживаемыми мыслями 2
и действ
иями людей находящимися под покровительством различных планет. Эгрегоры состоят из одного или нескольких кристаллических модулей взаимодействующих между собой по иерархическому принципу. Кристаллические модули эгрегоров разных стран питаются Биоплазмой кот
орую получают из разных фрактальных планетарных каналов. Есть различные методики подключения к эгрегорам и правила работы с ними. Подключение к эгрегору страны и работа с ним облегчает достижение целей, которые ставит перед собой человек
. Одной из методик подключения к эгрегору страны является использование ключей доступа к эгрегору в виде числовых космоматриц
.
Видите, все очень просто –
достаточно нащупать нужную числовую космоматрицу!
Но кроме шуток –
в этом что
-
то есть. Вот посмотрите на досуге лекции г
енерал
-
майора Петрова Константина Павловича (1945 –
2009) –
очень познавательно:
http://kob.su/petrov
Константин Павлович Петров
На данный момент не существует какой
-
либо даже слабой теории, объясняющей работу моз
га и структуру мысли в частности. Мысль -
это наиболее заметный элемент мышления, т.е. процесса происходящего в мозгу; элемент, с которым каждый день имеет дело любой человек. Мысль имеет разные проявления и формы:
воспоминания, фантазии, ассоциации, чувст
ва, убеждения и гипотезы. Какая мысль приходит Вам в голову глядя на работу Эдуарда Путерброда
, представленную в проекте "Неожиданный Путерброт"?
Фотография работы приведена на следующем листе.
Выставка работ Путерброт
а -
счастливый случай выйти не только
на новое понимание искусства, но и приобщиться к совершенно неожиданному направлению в искусстве третьего тысячелетия
, к искусству Миллениума или искусству информационники
! Экспозиция "Неожиданный Путерброт" родилась не как экспромт, а как итог продуманно
го стремления утвердить идеи творчества позднего периода этого замечательного, трагически погибшего дагестанского театрального художника. 3
Эдуард Путерброт, начиная с 90
-
х годов, пытался воплотить в своих работах неведомое прежде людям Земли, знания, получе
нные им необычайным образом. Путём трансакции настоящего и будущего, реальности и виртуальности... Вот
–
вот –
теплее!
Одна из работ Эдуарда Путерброда
Помните диалог из фильма "Осенний марафон" (1979):
Билл:
Дом, где я спал, как называется? Тресвевател
?
Бузыкин:
Вытрезвитель.
Билл:
Андрей, там было много новых слов... Андрей, я алкач?
Бузыкин:
Алкач, алкач...
Билл:
Андрей, а ви —
ходок?
Бузыкин:
Ходок, я ходок!..
Билл:
Андрей.
..
Хорошо сидим...
Итак, мы гурьбой уже почти написали трактат
и хотим опубли
ковать его. Можно поступить следующим образом:
Если мы понимаем, что написали, и можем доказать это -
посылаем в какой
-
нибудь математический журнал
!
Если мы понимаем, что написали, но не можем это доказать -
посылаем в журнал по физике
!
Если мы не понимаем
, что написали, но можем доказать это -
посылаем в журнал по экономике
!
Если мы ни хрена не понимаем и не можем доказать это -
посылаем в журнал по философии
!
В противном случае (
вдруг найдется еще вариант
!) посылаем работу на чудесный сайт публикации доку
ментов:
http://www.docme.ru/
4
Хотел вот написать очередной "Фрагмент" про алгоритмы "аппроксимации" и "оптимизации"
, но посмотрел доступные "сетевые ресурсы" и ужаснулся –
даже примитивные программки решения систем линейных алгебраических уравнений куда
-
то
исчезли, а деловые молодые люди пытаются заработать на этом несчастье! Они не понимают, что научив абсолютно даром десяток людей ремеслу, научивший, в конечном счете, получает существенно больше, чем "сдирая трояки" за чужие программки... Таким посредника
м не место в цивилизованном сообществе! Платить следует только за произведенный продукт или предоставленную услугу.
Вот почему этот выпуск "Фрагментов" будет посвящен линейным уравнениям и матрицам.
И вот первый вопрос к читателям:
одно уравнение с одним неизвестным определяет одну точку,
два уравнения с двумя неизвестными определяют две прямые,
три уравнения с тремя неизвестными определяют три плоскости,
а как просто и внятно ответить, что определяют четыре уравнения с четырьмя неизвестными?
Может быть, э
то частный случай четырех уравнений с пятью неизвестными как в замечательном рассказе Никиты Владимировича Разговорова
(1920
—
1982) "Четыре четырки" (1963), опубликованном в 14 томе "Библиотеки современной фантастики" (1967 -
Антология советской фантастики)
.
Цитирую:
Наконец, на последней из вставленных страниц профессор Ир сразу же разобрал еще одну константу: "Бусука —
лучший друг тида". Несколько секунд он сидел совершенно неподвижно. Потом, преодолев оцепенение, снова весь погрузился в лежавшую перед ним
страницу. Одна за другой он выписал еще три константы, подчеркивая незнакомые слова.
Бусука воет, ветер носит.
Любить, как Бусука палку.
Четыре четырки, две растопырки, седьмой вертун
-
бусука.
На этом вставные константы заканчивались. Профессор Ир с раздра
жением посмотрел на следующую страницу. Увы! Здесь уже снова шли знакомые физические постоянные. Удельный заряд электрона!.. Все, что он смог узнать о таинственной Бусуке, свелось пока к этим выписанным строчкам, четырем уравнениям с пятью неизвестными...
5
Не волнуйтесь –
в этом выпуске "Фрагментов" речь будет идти только о "тривиальном" случае N линейных алгебраических уравнений с N неизвестными
, хотя имеют место быть и случаи, когда количество неизвестных меньше количества уравнений и когда количество н
еизвестных больше количества уравнений. Великая вещь Математика! Она
-
то точно знает, что Бусука не есть слон –
Бусука это собака
. Только про собаку говорят "Собака воет, ветер носит!"
2.
Традиционный экскурс в историю
Алгебра (от араб. "аль
-
джабр" —
во
сполнение) —
раздел математики, который можно грубо охарактеризовать как обобщение и расширение арифметики. Слово "алгебра" также употребляется в названиях различных алгебраических систем. В более широком смысле под алгеброй понимают раздел математики, пос
вящённый изучению операций над элементами множества произвольной природы, обобщающий обычные операции сложения и умножения чисел. Истоки алгебры уходят к временам глубокой древности. Ещё 4000 лет назад вавилонские учёные могли решать квадратные уравнения. Тогда никаких обозначений не было, и уравнения записывались в словесной форме. Первые обозначения появились в Древней Греции благодаря ученому Диофанту. Диофант Александрийский
уже был упомянут в 5
-
ом выпуске "Фрагментов".
Алгебру можно грубо разделить на
следующие категории:
Элементарная алгебра (школьная алгебра),
Абстрактная алгебра (современная алгебра),
Линейная алгебра (алгебра векторных пространств),
Универсальная алгебра (свойства алгебраических структур),
Алгебраическая теория чисел (расширение ал
гебры),
Алгебраическая геометрия (алгебра решения проблем геометрии),
Алгебраическая комбинаторика (абстрактной алгебры в комбинаторике).
Не будем уходить в глубь веков, но нельзя не упомянуть французского математика Франсуа Виета
(Francois Viete; 1540 —
1603) —
основоположника символической алгебры. Виет четко представлял себе конечную цель —
разработку нового языка, своего рода обобщенной арифметики, которая дала бы возможность проводить математические исследования с недостижимыми ранее глубиной и общно
стью: "Все математики знали, что под их алгеброй... были скрыты несравненные сокровища, но не умели их найти; задачи, которые они считали наиболее трудными, совершенно легко решаются десятками с помощью нашего искусства, представляющего поэтому самый верны
й путь для математических изысканий".
Огромное влияние на развитие математики оказала "Геометрия" Декарта, вышедшая в 1637 году. Декарт
включает в геометрию широкий класс кривых, в том числе "механические" (трансцендентные, вроде спирали), и провозглашает
, что у каждой кривой есть определяющее уравнение. Он строит такие уравнения для 6
алгебраических кривых, проводит их классификацию (позже основательно переделанную Ньютоном). Декарт подчёркивает, хотя и не доказывает, что основные характеристики кривой не з
ависят от выбора системы координат. Система координат у Декарта была перевёрнута по сравнению с современной (ось ординат горизонтальна), и отрицательные координаты не рассматривались. Термины "абсцисса" и "ордината" изредка встречаются у разных авторов, хо
тя в широкое употребление их ввёл только Лейбниц
в конце XVII века, вместе с термином "координаты".
На рисунке показано, что с
истема линейных уравнений трех переменных определяет набор из трех плоскостей
и т
очка пересечения является решением.
Рене Де
карт
(Renatus Cartesius —
Картезий; 1596 —
1650) —
французский математик, философ, физик и физиолог, создатель аналитической геометрии и современной алгебраической символики, автор метода радикального сомнения в философии, механицизма в физике, предтеча ре
флексологии.
7
Готфрид
Вильгельм
Лейбниц
(Gottfried Wilhelm von Leibniz; 1646 —
1716) —
немецкий
философ
, логик
, математик
, физик
, историк
, дипломат
, изобретатель
и
языковед
. Основатель и первый президент Берлинской Академии наук, иностранный член Французской Академии наук.
Так начиналась линейная алгебра —
важная в приложениях часть алгебры, изучающая векторы, векторные, или линейные пространства, линейные отображения и системы линейных уравнений. Векторные пространства встречаются в математи
ке и её приложениях повсеместно. Линейная алгебра широко используется в абстрактной алгебре и функциональном анализе и находит многочисленные приложения в естественных науках. Исторически первым вопросом линейной алгебры был вопрос о линейных уравнениях.
П
остроение теории систем таких уравнений потребовало таких инструментов, как теория матриц и определителей, и естественно привело к появлению теории векторных пространств.
В линейной алгебре рассматриваются четыре класса основных задач:
решение систем лине
йных алгебраических уравнений (СЛАУ),
вычисление определителей,
нахождение обратных матриц,
определение собственных значений и собственных векторов матриц.
Все эти задачи имеют важное прикладное значение при решении различных проблем науки и техники. Кром
е того, задачи линейной алгебры являются вспомогательными при реализации многих алгоритмов
вычислительной математики, математической физики, обработки результатов экспериментальных исследований.
3.
Некоторые предварительные сведения
Система M линейных алгебраических уравнений с N неизвестными
(или, линейная система, также употребляется аббревиатура СЛАУ) это система уравнений вида:
A(1,1)*X(1)+A(1,2)*X(2)+...+A(1,N)*X(1)=B(1)
A(2,1)*X(1)+A(2,2)*X(2)+...+A(2,N)*X(2)=B(2)
(1)
...
A(M,1)*X(1)+A(M,
2)*X(2)+...+A(M,N)*X(M)=B(M)
Здесь
:
M
—
количество уравнений, а N
—
количество неизвестных.
X(1),X(2),...,X(N)
—
неизвестные, которые надо определить.
A(1,1),A(1,2),...,A(M,N)
—
коэффициенты системы
B(1),B(2),...,B(M)
—
свободные члены системы.
8
Индексы коэ
ффициентов A(I,J)
обозначают номера уравнения (
I
) и неизвестного (
J
), при котором стоит этот коэффициент, соответственно.
Все элементы A(I,J)
и B(I)
предполагаются известными.
Система называется однородной, если все её свободные члены равны нулю:
B(1)=B(2)
=...=B(M)=0
, иначе —
неоднородной.
Система называется квадратной, если число M
уравнений равно числу N
неизвестных.
В настоящее время с помощью компьютеров численно решаются СЛАУ высокого порядка (более чем N = 10^6). Такие решения осуществляются при помо
щи прямых или итерационных численных методов
. Прямые методы позволяют в предположении отсутствия ошибок округления (при проведении расчетов на идеальном, т.е. бесконечноразрядном компьютере) получить точное решение задачи за конечное число арифметических д
ействий. Итерационные методы, или методы последовательных приближений, позволяют вычислить последовательность, сходящуюся к решению задач при бесконечном выполнении итераций (на практике, разумеется, ограничиваются конечным значением итераций в зависимости
от требуемой точности).
К решению систем линейных алгебраических уравнений сводятся многочисленные практические задачи (по некоторым оценкам более 75% всех задач
). Можно с полным основанием утверждать, что решение линейных систем является одной из самых распространенных и важных задач вычислительной математики. Конечно, существует много методов и современных пакетов прикладных программ для решения СЛАУ, но для того, чтобы их успешно использовать, необходимо разбираться в основах построения методов и алгор
итмов, иметь представления о недостатках и преимуществах используемых методов.
3.1.
Метод последовательных приближений
Для решения СЛАУ итерационным методом необходимо преобразовать квадратную систему уравнений (1) к виду:
X(1)=(B(1)
-
A(1,1)*X(1)
-
A(1,2)
*X(2)
-
...
-
A(1,N)*X(1))/A(1,1)+X(1)
X(2)=(B(2)
-
A(2,1)*X(1)
-
A(2,2)*X(2)
-
...
-
A(2,N)*X(2))/A(2,2)+X(2)
(2)
X(N)=(B(N)
-
A(N,1)*X(1)
-
A(N,2)*X(2)
-
...
-
A(N,N)*X(N))/A(N,N)+X(N)
Задавая столбец начальных приближений X0(1),X0(2),...,X
0(N)
и подставляя их в правую часть системы (2) можно вычислить новые приближения неизвестных:
X1(1),X1(2),...,X1(N)
, которые опять подставляются в правую часть системы (2). Таким образом осуществляется итерационный процесс (метод простых итераций). Разраб
отчиком алгоритма принято считать Зейделя.
Филипп Людвиг Зейдель
(Philipp Ludwig von Seidel; 1821 —
1896) —
немецкий математик и астроном. Для сходимости итерационных методов необходимо, чтобы значения диагональных элементов матрицы СЛАУ были преоблада
ющими по абсолютной 9
величине по сравнению с другими элементами. Условие сходимости можно обеспечить преобразованием исходной матрицы путем перестановки уравнений и неизвестных. Методы простых итераций Зейделя имеют разные области сходимости. Эти методы мож
но применять и к решению систем нелинейных уравнений. Ниже приведен текст программы, соответствующей указанному алгоритму.
' P R O G R A M "MatrSeid"
' 08.10.1991
'
-------------------------------------
------------------------
' Решение системы линейных алгебраических уравнений
'
-------------------------------------------------------------
' Используется итерационный метод Зейделя
' N -
количество неизвестных
' A(I,J) -
элемент расши
ренной матрицы ' I -
номер строки расширенной матрицы
' J -
номер столбца расширенной матрицы
' E -
требуемая погрешность вычисления
' S -
суммарная текущая погрешность
' IC -
счетчик итераций
' IM -
мак
симально допустимое число итераций
#Lang "FB" ' режим FreeBASIC
-
совместимости
Dim As Integer N, I, J, K, IM, IC
Read N
Dim As Single A(N,N+1), X(N), E, DS, S Data 3, 100, 0.000001 ' исходные данные
Data 1, 4, 2, 3 ' Вариант 1
Data 2,
7, 3, 4 ' итерации сходятся!
Data 0, 8, 1, -
5
'Data 10, -
7, 0, 7 ' Вариант 2
'Data -
3, 2, 6, 4 ' итерации расходятся!
'Data 5, -
1, 5, 6
Read IM, E
Open "CON" For Output As #2
'Open "MatrSeid.res" For Output As #2
' В
вод расширенной матрицы по строкам и обнуление неизвестных
For I = 1 To N For J = 1 To N+1: Read A(I,J): X(I) = 0: Next J
Next I
' Выполнение итерационного цикла методом Зейделя
For K = 1 To IM: IC = K ' повторять итерации
S = 0 ' сумм
арная
погрешность
For I = 1 To N: DS = A(I,N+1) For J = 1 To N: DS = DS
-
A(I,J)*X(J) : Next J
DS = DS/A(I,I): X(I) = X(I)+DS: S = S +Abs(DS)
Next I If S < E Then Exit For ' Точность
достигнута
Next K If K > IM Then
Print #2, "Error!" ' Итерации
исчерпаны
Else
10
Print #2, " Vector X(I):" ' Вывод
вектора
решения
For I = 1 To N
Print #2, Using " ###.###### "; X(I);
Next I
EndIf
Print #2,
Print #2, "Iteration Counter = "; IC
Close #2
Sleep
В качестве исходных данных рассматриваются два варианта, пе
рвый из которых позволяет найти решения методом итераций:
Vector X(I):
1.000000 -
1.000000 3.000000 Iteration Counter = 6
Второй вариант решить методом итераций не получается:
Error!
Iteration Counter = 100
Условие сходимости можно попробовать обеспечить преобразованием исходной матрицы путем перестановки уравнений и неизвестных
. Кроме того, количество итераций зависит от выбранных начальных значений неизвестных.
3.2.
Вычисление определителя матрицы
Понятие определителя вводится как одно из основных в линейной алгебре. Это понятие распространяется только на квадратные матрицы, у которых порядки равны: N = M), т.е. квадратные таблицы из чисел (как, например, 3). Говорят об определителях матриц, элементами которых являются действительные (или к
омплексные) числа и тогда определитель -
действительное (или комплексное) число. Определитель матрицы A(I,J) обозначается как det(A). Здесь, как и раньше, первый индекс I означает номер строки, а второй индекс J —
номер столбца. В вычислительной практике и
спользуются два способа вычисления определителя матрицы:
1. по определению
-
через разложение по строке или столбцу (для N <=3),
2. по методу Гаусса
-
приведение матрицы к треугольному виду (для N > 3).
A(1,1), A(1,2), ... A(1,N)
A(2,1), A(2,2), ... A(2,
N)
(3)
...
A(N,1), A(N,2), ... A(N,N)
Для матриц вида (3) вводятся понятия главной и побочной диагоналей. Главной диагональю матрицы A(I,J)
называется диагональ, идущая из левого верхнего угла этой матрицы в правый нижний ее уго
л. Побочной диагональю той же матрицы называется диагональ, идущая из левого нижнего угла в правый верхний угол. Среди квадратных матриц выделяют класс так называемых диагональных матриц, у каждой из которых элементы, расположенные вне главной диагонали, р
авны нулю.
11
Если порядок N матрицы A равен единице, то эта матрица состоит из одного элемента A(1,1) и определителем первого порядка, соответствующим такой матрице, называют величину самого этого элемента. Далее, если порядок матрицы A равен двум, т.е. есл
и эта матрица имеет вид:
A(1,1), A(1,2)
A(2,1), A(2,2)
то определителем второго порядка, соответствующим такой матрице, называют число, равное A(1,1)*A(2,2) -
A(1,2), A(2,1) (4)
Формула (4) представляет собой правило составления опред
елителя второго порядка по элементам соответствующей ему матрицы. Словесная формулировка этого правила такова: определитель второго порядка, соответствующий матрице A, равен разности произведения элементов, стоящих на главной диагонали этой матрицы, и прои
зведения элементов, стоящих на побочной ее диагонали.
Перейдем теперь к выяснению понятия определителя любого порядка N
, где N>=2
. Понятие такого определителя вводится индуктивно, в предположении, что уже введено понятие определителя порядка N
-
1
, соответс
твующего произвольной квадратной матрице порядка N
-
1
. При этом договариваются называть минором любого элемента A(I,J)
матрицы N
-
го порядка определитель порядка N
-
1
, соответствующий той матрице, которая получается из исходной матрицы A
в результате вычеркив
ания I
-
ой строки и J
-
го столбца (той строки и того столбца, на пересечении которых стоит элемент A(I,J)
). Минор элемента A(I,J)
обычно обозначают символом М(I,J)
.
Отсюда рождается идея алгоритма первого способа вычисления определителя матрицы –
"по опреде
лению" или иначе через разложение по строке или столбцу. И этот алгоритм, конечно, должен быть "индуктивно рекурсивным"! Из свойств определителя следует, что определитель матрицы порядка N
может быть представлен в виде суммы N
определителей N
-
1
порядка (разложение по строке или столбцу)... Рекурсия —
процесс повторения элементов самоподобным образом
. В программировании рекурсия —
вызов подпрограммы (процедуры или функции) из неё же самой, непосредственно (простая рекурсия) или через другие функц
ии (сложная или косвенная рекурсия), например, функция FA
вызывает функцию FB
, а функция FB
—
функцию FA
. Количество вложенных вызовов функции или процедуры называется глубиной рекурсии. Использование рекурсии в некоторой программе позволяет описать "беско
нечное вычисление", причем без явных повторений частей программы. К сожалению не все трансляторы допускают рекурсивные вызовы...
Ниже представлен текст программы для компилятора FreeBASIC.
' P R O G R A M "MatrDet_Rec"
' 08.12.2012
'
-------------------------------------------------------------
' Рекурсивная процедура вычисления определителя матрицы
'
-------------------------------------------------------------
' Для вычисления определителя квадра
тной матрицы используется
' разложение по столбцу -
рекурсивная процедура.
' A(I,J) -
элементы матрицы,
' N -
размерность матрицы.
'
#Lang "FB" ' режим FreeBASIC
-
совместимости
12
Sub GetMatr(A() As Single, B() As Single, M As integer,_
I As Integer, J As Integer)
' Вычеркивание из матрицы строки и столбца
Dim As Integer ki, kj, di, dj
di = 0
For ki = 1 To M -
1
If (ki = I) Then di = 1
dj = 0
For kj = 1 To M -
1
If (kj = J) Then dj = 1
B(ki,kj) = A(ki+di,kj+dj)
Next kj
Next ki
E
nd Sub
'
Function Determinant(A() As Single, N As Integer) As Single
' Вычисление определителя матрицы
Dim As Integer i, j, k
Dim As Single b(N, N), d
d = 0: k = 1
If (N < 1) Then
Print "Determinant: Cann't run. N ="; N
GoTo Ex
EndIf
If (N = 1) Then d = A(1,1): GoTo Ex
If (N = 2) Then d = A(1,1)*A(2,2)
-
A(2,1)*A(1,2): GoTo Ex
For i = 1 To N
GetMatr(A(), b(), N, i, 1)
d = d + k*A(i,1)*Determinant(b(), N
-
1)
k = -
k
Next i
Ex: Determinant = d
End Function
'
' Основной модуль тестовой програ
ммы
Dim As Integer I, J, N
Read N ' чтение заданной размерности
' Элементы квадратной матрицы
Dim As Single A(N, N), DET Data 3 ' текущая размерность массива
'Data 1, 4, 2 ' Вариант
1
'Data 2, 7, 3 ' DET = 7.00
0
'Data 0, 8, 1
Data 10, -
7, 0 ' Вариант
2
Data -
3, 2, 6 ' DET = -
155.000
Data 5, -
1, 5
For I = 1 To N
For J = 1 To N: Read A(I, J): Next J
Next I
Open "CON" For Output As #2
'Open "MatrDet_Rec.res" For Output As #2
Print #2, " Matr
ix A(I,J):"
13
For I = 1 To N
For J = 1 To N
Print #2, Using " ###.###### "; A(I, J);
Next J
Print #2,
Next I
' Вычисление и вывод определителя
DET = Determinant(A(), N)
Print #2, Using " DET = #####.##### "; DET
Close #2
Sleep
Ниже приводится текст программ
ы, соответствующей второму способу вычисления определителя матрицы –
"по методу Гаусса", т.е. методом приведения матрицы к треугольному виду. Операции по выполнению прямого хода метода Гаусса в соответствии с теоремами линейной алгебры не изменяют величины
определителя. Очевидно, что определитель треугольной матрицы равен произведению ее диагональных элементов. Таким образом, для вычисления определителя квадратной матрицы необходимо привести ее к треугольному виду и найти произведение ее диагональных элемен
тов. При этом следует учесть, что выбор главного элемента осуществляется перестановкой строк матрицы, каждая из перестановок изменяет знак определителя на противоположный. При малых диагональных элементах треугольной матрицы перемножение их может привести к исчезновению порядка, при больших величинах диагональных элементов -
к переполнению.
' P R O G R A M "MatrDet_Gauss"
' 08.10.1991
'
-------------------------------------------------------------
' Вычисление определителя матрицы методом Гаусса
'
-------------------------------------------------------------
' Для вычисления определителя квадратной матрицы она
' приводится к треугольному виду и находится произведение
' ее диагональных элементов.
' N -
размерность квадратной матрицы;
' A() -
элементы исходной матрицы;
' S -
величина определителя матрицы.
'
#Lang "FB" ' режим FreeBASIC
-
совместимости
Sub MatDet (N As Integer, A() As Single, ByRef S As Single)
Dim As Single R, P
Dim As Inte
ger I, J, K, K1
P = 1
For K = 1 To N
-
1
K1 = K + 1: S = A(K,K): J = K
For I = K1 To N
R = A(I, K)
If Abs(R) > Abs(S) Then
S = R: J = I
EndIf
Next I
'
14
If S = 0 Then Exit Sub ' нулевой
элемент
If J <> K Then
For I = K To N
R = A(K,I)
A(K, I) = A(J,I)
A(J, I) = R
Next I
P = -
P
EndIf
For J = K1 To N: A(K,J) = A(K,J)/S: Next J
For I = K1 To N: R = A(I,K)
For J = K1 To N: A(I,J) = A(I,J) -
A(K,J)*R: Next J
Next I
P = P * S
Next K
S = P * A(N,N)
End Sub
'
' Основной модуль тестовой програ
ммы
Dim As Integer N, I, J
Read N
Dim As Single A(N,N), DET
' Исходные данные:
Data 3 ' текущая размерность массива
'Data 1, 4, 2 ' Вариант
1
'Data 2, 7, 3 ' DET = 7.000
'Data 0, 8, 1
Data 10, -
7, 0 ' Вариант
2
Da
ta -
3, 2, 6 ' DET = -
155.000
Data 5, -
1, 5
For I = 1 To N
For J = 1 To N: Read A(I,J): Next J
Next I
Open "CON" For Output As #2
'Open "MatrDet_Gauss.res" For Output As #2
Print #2, " Matrix A(I,J):"
For I = 1 To N
For J = 1 To N
Print #2, Usin
g " ###.###### "; A(I,J);
Next J
Print #2,
Next I
'
MatDet(N, A(), DET)
'
Print #2, Using " DET = #####.##### "; DET
Close #2
Sleep
15
3.3.
Решение системы методом Крамера
Сегодня трудно указать один или несколько способов наиболее эффективных в смысле скорости решения системы уравнений с нужной точностью и минимальным использованием имеющихся ресурсов, например, памяти компьютера. Для это требуется большая и тщательная теоретическая и экспериментальная работа по сравнительной оценке многочисленных извес
тных способов. Теперь, мой дорогой читатель, ты к этому готов?
Стоит еще напомнить известное из курса высшей алгебры правило Крамера для решения систем линейных алгебраических уравнений, которое сегодня признается практически невыгодным, так как требует слишком большого количества арифметических операций и операций чтения
-
записи. Но и не привести программу, реализующую этот метод, было бы неправильно! Ведь проблема численного решения линейных уравнений интересует математиков уже несколько столетий, а перв
ые реальные математические результаты появились в XVIII веке во многом благодаря работам Крамера. В 1750 году Габриэль Крамер
(1704
-
1752) опубликовал свои труды по детерминантам квадратных матриц и предложил алгоритм нахождения обратной матрицы, известный как правило Крамера.
Метод Крамера требует вычисления N+1
определителей. Каждое значение X(I)
вычисляется как отношение определителей:
X(I) = DET(I)/DET(0)
, где
DET(0)
–
основной определитель системы,
DET(I)
–
определители, вычисляемые для матриц, получае
мых заменой столбца левой части уравнений на столбец свободных членов.
В приведенном ниже тексте программы используются уже рассмотренные выше подпрограммы -
процедура Sub GetMatr()
и функция Function Determinant()
.
' P R O G R A M "MatrKramer"
' 08.12.2012
'
-------------------------------------------------------------
' Решение системы линейных уравнений Методом Крамера
'
-------------------------------------------------------------
' A(I,J) -
исхо
дная расширенная матрица
' B(I,J) -
модифицированная квадратная матрица' Det(0) -
главный определитель матрицы системы;
' DET(I) -
определитель матрицы c заменой столбцов
' N -
количество неизвестных и уравнений в системе
' X(I) -
вектор решения = D
ET(I)/DET(0)
#Lang "FB" ' режим FreeBASIC
-
совместимости
Sub GetMatr(A() As Single, B() As Single, M As integer,_
I As Integer, J As Integer)
'
16
' Вычеркивание из матрицы строки и столбца
Dim As Integer ki, kj, di, dj
di = 0
For ki = 1 To M -
1
If
(ki = I) Then di = 1
dj = 0
For kj = 1 To M -
1
If (kj = J) Then dj = 1
B(ki,kj) = A(ki+di,kj+dj)
Next kj
Next ki
End Sub
'
Function Determinant(A() As Single, N As Integer) As Single
' Вычисление определителя матрицы
Dim As Integer i, j, k
Dim As Single b(N, N), d
d = 0: k = 1
If (N < 1) Then
Print "Determinant: Cann't run. N ="; N
GoTo Ex
EndIf
If (N = 1) Then d = A(1,1): GoTo Ex
If (N = 2) Then d = A(1,1)*A(2,2)
-
A(2,1)*A(1,2): GoTo Ex
For i = 1 To N
GetMatr(A(), b(), N, i, 1)
d = d + k*A(i,1)*Determinant(b(), N
-
1)
k = -
k
Next i
Ex: Determinant = d
End Function
'
' Основной модуль тестовой программы
Dim As Integer I, J, N
N = 3 ' текущая размерность массива
' Элементы квадратной матрицы
Dim As Single A(N,N+1)
, B(N,N), JD, DET(N) Data 1, 4, 2, 3
Data 2, 7, 3, 4
Data 0, 8, 1, -
5
'DET(0) -
определитель матрицы = 7.000
'Data 10, -
7, 0, 7
'Data -
3, 2, 6, 4
'Data 5, -
1, 5, 6
'DET(0) -
определитель матрицы = -
155.000
For I = 1 To N
For J = 1 To N+1: Read A(I,J): Next J
Next I
Open "CON" For Output As #2
'Open "MatrKramer.res" For Output As #2
Print #2, " Matrix A(I,J):"
'
17
For I = 1 To N
For J = 1 To N+1
Print #2, Using " ###.###### "; A(I,J);
Next J
Print #2,
Next I
' Вычисл
ение и вывод всех определителей
Print #2, " Matrix B(I,J):"
For JD = 0 To N ' DET(0) ' определитель матрицы = -
155.000
For J = 1 To N
For I = 1 To N: B(I,J) = A(I,J): Next I ' весь
столбец
If J = JD Then For I = 1 To N: B(I,J) = A(I,N+1): Next I
Next J
For I = 1 To N
For J = 1 To N+1
Print #2, Using " ###.###### "; B(I,J);
Next J
Print #2,
Next I
DET(JD) = Determinant(B(), N)
Print #2, Using " DET = #####.##### "; DET(JD)
Next JD
Print #2, " Vector X(I):"
For I = 1 To N
Print #2, Using " ###.###
### "; DET(I)/DET(0)
Next I
Close #2
Sleep
3.4.
Решение системы методом Гаусса
Гаусс в 1809 году опубликовал работу, посвященную движению небесных тел, в которой был изложен метод для решения линейных систем, известный как метод исключения. Иоганн Кар
л Фридрих Гаусс
(Johann Carl Friedrich Gauss; 1777 —
1855) —
немецкий математик, астроном и физик, считается одним из величайших математиков всех времён, "королём математиков". Лауреат медали Копли (1838), иностранный член шведской (1821) и российской (182
4) Академий наук, английского Королевского общества.
Текст программы, пожалуй, не требует дополнительных пояснений –
все предельно ясно! Никакие подпрограммы не используются. С нее можно начинать повествование о линейных уравнениях...
18
' P R O G R A M "Matr_Gauss"
' 08.10.1991
'
-------------------------------------------------------------
' Решение уравнений методом последовательного исключения
'
-----------------------------------------------------
--------
' Метод Гаусса основан на приведении исходной матрицы
' к треугольному виду. В программе обозначено:
' N -
порядок системы,
' A(I,J) -
коэффициенты левой части системы
' B(I) -
свободные члены
#Lang "FB" ' режим FreeBASIC
-
совмест
имости
Data 3
Data 10, -
7, 0, 7
Data -
3, 2, 6, 4
Data 5, -
1, 5, 6
Dim As Integer N, I, J, K
Read N
Dim As Single A(N, N), B(N), X(N), H
' Ввод исходных коэффициентов A(I,J) и B(I) For I = 1 To N
For J = 1 To N: Read A(I,J): Next J
Read B(I)
Next I
' Прямой ход исключения переменных
For I = 1 To N
-
1
For J = I+1 To N
A(J,I) = -
A(J,I)/A(I,I)
For K = I+1 To N
A(J,K) = A(J,K)+A(J,I)*A(I,K)
Next K
B(J) = B(J)+A(J,I)*B(I)
Next J
Next I
X(N) = B(N)/A(N,N)
' Обратный х
од вычисления неизвестных
For I = N
-
1 To 1 Step -
1
H = B(I)
For J = I+1 TO N
H = H
-
X(J)*A(I,J)
Next J
X(I) = H/A(I,I)
Next I
Print " VECTOR X(I)"
For I = 1 TO N: PRINT Using " ###.###### "; X(I)
Next I
Sleep
'
' Vector X(I):
' -
0.000004
' -
1.000005
' 1.0
00000
19
При желании получить большую точность вычисления следует заменить строку:
Dim As Single A(N, N), B(N), X(N), H
на
строку
:
Dim As Double A(N, N), B(N), X(N), H
3.5.
Решение с использованием подпрограмм DECOMP и SOLVE
В очередной раз призываю доро
гих читателей скачать из сети и полистать замечательную книжку с примерами на Фортране:
Форсайт Дж., Малькольм М., Моулер К.,
Машинные методы математических вычислений. М., Мир, 1980, 280 с., илл
.
Далее цитирую эту книгу со страницы 43:
Глава 3. Линейные системы уравнений.
Одна из задач, наиболее часто встречающихся в научных вычислениях, решение системы линейных уравнений; при этом обычно число уравнений равно числу неизвестных. Такую систему можно записать в виде:
Ax = b,
где
А
-
заданная квадратная матр
ица порядка n,
b
-
заданный вектор столбец с n компонентами,
x
-
неизвестный вектор
-
столбец с n компонентами.
Источники линейных уравнений включают аппроксимацию непрерывных дифференциальных или интегральных уравнений конечными, дискретными алгебраическим
и системами, локальная линеаризация систем нелинейных уравнений, построение полиномов или кривых какого либо иного специального вида по заданной информации. Некоторые из этих приложений будут обсуждены в последующих главах.
И далее со страницы 61:
Парагра
ф 3.3. Подпрограммы DECOMP и SOLVE
Почти в любой машинной библиотеке имеются подпрограммы, основанные на вариантах гауссова исключения с частичным выбором ведущего элемента для решения систем линейных уравнений. Детали реализации разных подпрограмм очень р
азличны. Эти детали могут серьезно отразиться на времени выполнения данной подпрограммы, но они не должны сильно влиять на ее точность, если подпрограмма правильно составлена. В этом параграфе мы опишем две такие подпрограммы DECOMP и SOLVE.
DECOMP выполня
ет ту часть гауссова исключения, которая зависит лишь от матрицы. Она сохраняет множители и информацию о ведущих элементах. Подпрограмма SOLVE использует эти результаты, чтобы получить решение для произвольной правой части.
' P R O G R A M "Matr_Fors
"
' 04.12.2012
'
-------------------------------------------------------------
' Иллюстрирующая программа для подпрограмм DECOMP и SOLVE
'
-------------------------------------------------------------
' Рассматривается решение системы линейных алгебраических
' уравнений вида: A*X = B.
'
#Lang "FB" ' режим FreeBASIC
-
совместимости
20
Sub SOLVE(NDIM As Integer, N As Integer, A() As Single,_
B() As Single, IPVT() As Integer)
' Подпрограмму не с
ледует использовать, если
' подпрограмма DECOMP обнаружит вырожденность!
'
' Входные данные:
' NDIM -
заявленная строчная размерность массива
' N -
порядок матрицы
' A -
факторизованная матрица, полученная DECOMP
' B -
вектор правых частей у
равнений
' IPVT -
вектор ведущих элементов, полученных DECOMP
'
' Выходные данные:
' B -
вектор решения X
'
Dim As Integer KB, KM1, NM1, KP1, I, K, M
Dim As Single T
'
' Прямой ход вычислений:
If(N <> 1) Then
NM1 = N
-
1
For K = 1 To NM1
KP1 = K+1
M = IPVT(K)
T = B(M)
B(M) = B(K)
B(K) = T
For I = KP1 To N
B(I) = B(I)+A(I,K)*T
Next I
Next K
'
' Oбратная подстановка:
For KB = 1 To NM1
KM1 = N
-
KB
K = KM1+1
B(K) = B(K)/A(K
,K)
T = -
B(K)
For I = 1 To KM1
B(I) = B(I)+A(I,K)*T
Next I
Next KB
EndIf
B(1) = B(1)/A(1,1)
End Sub
'
'
Sub DECOMP(NDIM As Integer, N As Integer, A() As Single,_
ByRef COND As Single, IPVT() As Integer, WORK() As Single)
' п
одпрограмма вычисляет разложение действительной матрицы
' методом гауссова исключения и оценивает ее обусловленность
'
21
' Входные данные:
' NDIM -
заявленная строчная размерность массива
' N -
порядок матрицы
' A -
матрица подлежащая разложению
'
' Выходные данные:
' A -
верхняя треугольная матрица
' COND -
оценка "обусловленности" матрицы A
' если COND+1=COND, то A в пределах машинной точности
' является вырожденной. COND принимается 10^32
' IPVT -
вектор ведущих элементов:
' IPVT(K) -
индекс K
-
ой ведущей строки
' IPVT(N) -
(
-
1)^<число перестановок>
'
Dim As Single EK, T, ANORM, YNORM, ZNORM
Dim As Integer NM1, I, J, K, KP1, KB, KM1, M
IPVT(N) = 1
If(N > 1) Then
NM1 = N -
1
'
' вычисление 1
-
нормы матрицы A
ANORM = 0.0
For J =1 To N
T = 0.0
For I = 1 To N
T = T+Abs(A(I,J))
Next I If(T > ANORM) Then ANORM = T
Next J
'
' Гауссово исключение с частичным выбором ведущего элемента
For K = 1 To NM1
KP1=K+1
'
' Поиск ведущего элемента
M = K
For I = KP1 To N
If(Abs(A(I,K)) > Abs(A(M,K))) Then M = I
Next I
IPVT(K) = M
If(M <> K) Then IPVT(N) = -
IPVT(N)
T = A(M,K)
A(M,K) = A(K,K)
A(K,K) = T
'
' Пропустить этот шаг, если ведущий элемент равен нулю
If(T <> 0) Then
' Вычи
слить множители
For I = KP1 To N
A(I,K) = -
A(I,K)/T
Next I
22
' Переставлять и исключать по столбцам
For J = KP1 To N
T = A(M,J)
A(M,J) = A(K,J)
A(K,J) = T
If(T <> 0) Then
For I = KP1 To N
A(I,J) = A(I,J)+A(I,K)*T
Next I
E
ndIf
Next J
EndIf
Next K
' COND = <1
-
норма A>*<оценка для 1
-
нормы обратной к A>
' ESTIMATE = <1
-
норма Z>/<1
-
норма Y>
'
' Решить систему <транспонированная для A>*Y = E
For K = 1 To N
T = 0.0
If(K <> 1) Then
KM1 = K
-
1
For I = 1 To KM1
T = T+A(I,K)*WORK(I)
Next I
EndIf
EK = 1.0
If(T < 0.0) Then EK = -
1.0
If(A(K,K) = 0.0) Then GoTo Ex
WORK(K) = -
(EK+T)/A(K,K)
Next K
For KB = 1 To NM1
K = N
-
KB
T = WORK(K)
KP1 = K+1
For I = KP1 To N
T = T+A(I,K)*WORK(I)
Next I WORK(K) = T
M = IPVT(K)
If(M <> K) Then
T = WORK(M)
WORK(M) = WORK(K)
WORK(K) = T
EndIf
Next KB
'
YNORM = 0.0
For I = 1 To N
YNORM = YNORM+Abs(WORK(I))
Next I
'
'
23
' Р
ешить
систему
A
*
Z
= Y
SOLVE(NDIM, N, A(), WO
RK(), IPVT())
'
ZNORM = 0.0
For I = 1 To N
ZNORM = ZNORM + Abs(WORK(I))
Next I
'
' Оценить обусловленность
COND = ANORM*ZNORM/YNORM
If(COND < 1.0) Then COND=1.0
Return
'
' Случай матрицы размера 1 х 1
EndIf
COND = 1.0
If(A(1,1
) <> 0.0) Then Return
'
' Точная вырожденность
Ex: COND = 1.0E+32
Return
End Sub
'
' Основной модуль тестовой программы
Dim As Integer NDIM
NDIM = 10 ' резервируемая размерность массива
Dim As Single A(NDIM, NDIM), B(NDIM), WORK(NDIM), COND, CONDPI
Dim As Integer IPVT(NDIM), I, J, N
Dim As Single DET N = 3 ' текущий порядок системы
' Элементы квадратной матрицы
Data 10, -
7, 0
Data -
3, 2, 6
Data 5, -
1, 5
For I = 1 To N
For J = 1 To N: Read A(I, J): Next J
Next I
Open "
CON" For Output As #2
'Open "Matr_Fors.res" For Output As #2
Print #2, " Matrix A(I,J):"
For I = 1 To N
For J = 1 To N
Print #2, Using " ###.###### "; A(I, J);
Next J
Print #2,
Next I
' Вызов
процедуры
DECOMP
DECOMP (NDIM, N, A(), COND, IPVT(), WORK())
Pri
nt #2, Using " COND = ###.###### ";COND
CONDPI = COND + 1
'
'
24
If (CONDPI = COND) Then
Print #2, " Matrix is Singular!"
GoTo St
EndIf
' Вектор правой части уравнений
B(1) = 7: B(2) = 4: B(3) = 6
Print #2, " Vector B(I):"
For I = 1 To N
Print #2, Using
" ###.###### "; B(I);
Next I
Print #2,
' Вызов
процедуры
SOLVE
SOLVE (NDIM, N, A(), B(), IPVT())
Print #2, " Vector X(I):"
For I = 1 To N
Print #2, Using " ###.###### "; B(I);
Next I
Print #2,
' Вычисление определителя матрицы A по формуле:
' DET(A) = IPV
T(N)*A(1,1)*A(2,2)*...*A(N,N)
DET = IPVT(N)
For I = 1 To N
DET = DET*A(I,I)
Next I
Print #2, Using " DET = #####.##### "; DET
St:
Close #2
Sleep
Думаю, что указанная программа будет актуальной еще, по крайней мере, лет десять –
двадцать. Сохраните ее и об
язательно поэкспериментируйте!
К сожалению временные рамки и объем работы не позволили даже обозначить проблему решения системы линейных уравнений методом Монте
-
Карло
. Будем живы
–
здоровы –
обсудим и этот интереснейший метод решения –
скорострельность совре
менных компьютеров делает его весьма привлекательным...
На этом позвольте закруглиться.
До новых встреч!
Пишите: eugene_r@mail.ru
Документ
Категория
Информатика
Просмотров
292
Размер файла
609 Кб
Теги
freebasic алгоритмы программирование
1/--страниц
Пожаловаться на содержимое документа