Программная реализация метода дифференциальной прогонки для решения краевой задачи

Министерство образования и науки РФ

Федеральное государственное бюджетное образовательное учреждение

высшего профессионального образования

«БРЯНСКИЙ ГОСУДАРСТВЕННЫЙ УНИВЕРСИТЕТ

ИМЕНИ АКАДЕМИКА И.Г.ПЕТРОВСКОГО»

(БГУ)

 

 

 

Естественно-научный институт

Физико-математический факультет

 

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

 

 

Курсовая работа

«Программная реализация метода дифференциальной прогонки для решения краевой задачи»

 

 

 

 

 

 

 

 

 

 

 

 

Выполнил:

студент 5 курса 1 группы

специальности «Прикладная математика и информатика»

Денежкин Д. В.

Руководитель:

кандидат

физико-математических наук,

доцент

Трубников С.В.

 

 

 

 

 

 

 

Оглавление

Введение …………………………..………………………………3

Глава 1. Итерационный метод вращений …………………….....4

1.1. Определения и основные  факты ……..……………………...4

1.2. Постановка задачи……………….……..…………..................6

1.3. Суть Метода…………….. …………………………………....8

Глава 2. Программная реализация итерационного метода вращений для получения собственных значений эрмитовой матрицы………………..11

2.1. Методы оптимизации при написании программы…………11

2.2 Описание некоторых  процедур и функций………………….12

2.3 Тестирование…………………………………………………..13

2.4 Листинг………………………………………………………...14

Заключение………………………………………………………...19

Список литературы………………………………………………..20

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

Введение

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

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

Матрицы возникающих систем могут иметь различные структуры и свойства. Уже сейчас имеется потребность в решении систем линейных алгебраических уравнений с матрицами полного заполнения порядка нескольких тысяч. При решении ряда прикладных задач методом конечных элементов в ряде случаев появляются системы, обладающие симметричными положительно определёнными ленточными матрицами порядка несколько десятков тысяч с половиной ширины ленты до тысячи. И, наконец, при использовании в ряде задач метода конечных разностей необходимо решить системы разностных уравнений с разрежёнными матрицами порядка миллион.

Одним из самых распространенных методов решения систем линейных уравнений является метод вращений Якоби. Этот метод (который также называют методом простых итераций) известен в различных вариантах уже более 200 лет.

 

 

Глава 1. Итерационный метод вращений.

1. Определения и основные факты

 

Собственный вектор — понятие в линейной алгебре, определяемое для квадратной матрицы или произвольного линейного преобразования как вектор, умножение матрицы на который или применение к которому преобразования даёт коллинеарный вектор — тот же вектор, умноженный на некоторое скалярное значение, называемое собственным числом матрицы или линейного преобразования.

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

Матрица называется симметричной, если она совпадает с транспонированной (со своей транспозицией): A = AT, или ai,j = aj,i.Она называется эрмитовой, или самосопряженной, если она равна комплексно-сопряженной матрице от своей транспозиции (эрмитово-сопряженной от нее): A = AH, или ai,j = aj,i*.Она называется ортогональной, если ее транспозиция является обратной к ней: AAT = ATA = 1.Она называетсяунитарной, если эрмитово-сопряженная от нее равна обратной. Наконец, матрица называется нормальной, если она является перестановочной с эрмитово-сопряженной: AAH = AHA.Для действительных матриц эрмитовость означает симметрию, унитарность совпадает с ортогональностью, и оба этих класса являются нормальными.

 При поиске собственных значений матрицы, "эрмитовость" является весьма важной концепцией: все собственные значения эрмитовых матриц действительны. С другой стороны, собственные значения действительных несимметричных матриц могут быть либо действительными, либо парами комплексно - сопряженных чисел. Собственные значения комплексной неэрмитовой матрицы в общем случае комплексные.

Концепция "нормальности" важна при поиске собственных векторов. Система собственных векторов нормальной матрицы с невырожденными (несовпадающими) собственными значениями является полным и ортогональным базисом N-мерного векторного пространства. Для нормальной матрицы с вырожденными собственными значениями имеется свобода в определении собственных векторов, соответствующих вырожденным собственным значениям, связанная с заменой их любой линейной комбинацией. Это означает, что мы всегла можем провести процесс ортогонализации Грама - Шмидта и найти полный ортогональный набор собственных векторов, как и в невырожденном случае. Очевидно, что матрица с колонками из ортонормированного множества собственных векторов является унитарной. Для матрицы собственных векторов, полученных из действительной симметричной матрицы, выполняется свойство ортогональности.

Если матрица не является нормальной, как, например, любая действительная несимметричная матрица, то в общем случае нельзя отыскать ортонормированный набор собственных векторов, нельзя даже гарантировать ортогональности любой пары из них (кроме редких случаев). В общем случае эти N собственных векторов будут образовывать неортогональный базис в N-мерном пространстве (но не всегда). Если собственные вектора не образуют N-мерный базис, то матрицу будем называть дефектной.

1.2 Постановка задачи           Эти задачи решены разными математиками в разное время и много раз. Но прежде, чем рассказывать о методах решения (точнее, здесь пойдет речь об одном методе – методе вращений, предложенным Якоби) нужно рассказать о самой задаче.

Итак, число   называется собственным значением матрицы  , если существует не нулевой вектор  , такой что

(1.1)

При этом вектор   называется собственным вектором матрицы  .

Выражение (1.1) можно переписать так:

(1.2)

Здесь   - единичная матрица.

Выражение в скобках определяет матрицу:

(1.3)     То есть, мы получили матричное уравнение:

(1.4)

В выражениях (1.2) и (1.4) под нулем в правой части понимается нуль-вектор.

Матричное уравнение (1.4) представляет собой систему линейных однородных уравнений. Оно имеет решение, если главный определитель равен нулю:

(1.5)

Когда мы решим это уравнение, получим n значений  . Это и будут собственные значения. Для каждого   составим матрицу (1.3) и решим уравнение (1.4). В конце концов, получим n собственных векторов  .

Вот здесь и возникли проблемы, о которых идет речь в названии статьи. Ведь для решения нужно составить уравнение (1.5) и решить его (что для больших порядков уже не просто). А потом еще и для каждого   решить (1.4). То есть нужно решить n матричных уравнений. Задача получается достаточно объемной.

ПРИМЕЧАНИЕ: Коэффициенты   характеристического многочлена (1.5) оператора   представляют собой инварианты относительно выбора базиса. В частности:

(1.6)

Это выражение называется следом и обозначается Sp. В литературе можно еще встретить обозначение tr. А почему бы нам просто не сменить базис, таким образом, чтобы в новой матрице все элементы, кроме стоящих на главной диагонали, обнулились? Если нам это удастся сделать, то вектора задающие новый базис будут собственными векторами, а элементы стоящие на главных диагоналях – собственными значениями. 1.3 Суть Медота

Якоби (математик XIX века) рассматривал вещественную симметричную матрицу  . Почему симметричную? Потому, что с несимметричными матрицами дело обстоит гораздо сложнее. К тому же, симметричные матрицы встречаются в приложениях гораздо чаще. Для перехода к новому базису, Якоби предложил следующее преобразование.

Построим матрицу  :

(2.1)

Возьмем теперь нашу матрицу  и выполним преобразование:

(2.2)

Такое преобразование называется вращением (отсюда и название метода – метод вращений), а матрица   – матрицей вращения.

Преобразования (2.2) обладают следующими свойствами:

1. Сохраняют свойства  симметрии:  ;

2. Изменяются только элементы, стоящие в i- ой и j- ой строке  и в j- м и i- м столбце;

3. Сохраняются собственные  значения и собственные векторы.

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

Чтобы построить эту матрицу   нужно ответить на два вопроса:

Во-первых, каково значение угла  ?

И, во-вторых, куда мы будем ставить   и  ?

Ответ на второй вопрос дает свойство 2. Разумеется,   и  нужно ставить так, чтобы изменился максимальный недиагональный элемент  .

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

(2.3)

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

Итак, после одного преобразования, мы получили новую матрицу  , у которой на месте максимального элемента матрицы   стоит 0. Самое главное, что при таком преобразовании сохранились собственные значения и вектора (см. свойство 3). Совершенно верно, Производим те же операции с новой матрицей  . Ищем максимальный элемент в матрице   и уничтожаем его описанным выше способом.

Элемент, убитый на предыдущем шаге может возникнуть, но уже меньший по модулю! В книге “Ильин В.А., Позняк Э.Г. Линейная алгебра М.: 1974, «Наука»” авторами доказывается сходимость данного метода. Доказательство основано, на том, что сумма квадратов недиагональных элементов матрицы   (  - номер шага) стремится к нулю.

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

 

Глава 2. Программная реализация итерационного метода вращений для получения собственных значений матрицы на языке Visual Basic.

2.1 Методы оптимизации  при написании программы

Переходим к реализации алгоритма. Так как матрица симметричная, то нет необходимости хранить ее полностью. Достаточно хранить только верхний или нижний треугольник (часть, матрицы лежащая выше или ниже главной диагонали, соответственно). Реализуем это следующим образом: создаем одномерный массив размеров K = ((1 + N) / 2) * N, где N – это размерность матрицы. Во всех местах, где мы обращаемся к элементам матрицы, будем вызывать специальную функцию GetIndex, которая по двум индексам квадратной матрицы будет давать нам индекс в одномерном массиве.

Такая оптимизация даст нам, например, для матрицы размер 20 на 20, следующие преимущества:

Массив размером 20 на 20 типа Double, будет занимать в памяти   байтов, т.к. число типа Double занимает 8 байт.

Массив размером   типа Double, будет отнимать 1680 байт. И так, мы сэкономили   памяти.

Нужно ли выполнять полное умножение матриц при вращении? Так, как матрица (2.1) почти полностью состоит из 0, то нет смысла выполнять лишние операции. Все равно получится 0. Для вращения воспользуемся следующими преобразованиями:

(5.1)

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

2.2 Описание некоторых процедур и функций

В модуле mdlMain находится несколько процедур и функция.   

Clear() - Процедура отчистки рабочей области. Вызывается из меню и из главной процедуры. К алгоритму никакого отношения не имеет;

   Show(Row As Long) - Процедура вывода матрицы на рабочий лист. Row параметр, номер строки откуда выводить. К алгоритму никакого отношения не имеет;

   MultMatrix(FirstMatr() As Double, SecondMatr() As Double, ResMatrix() As Double) - Процедура умножения матриц. FirstMatr() – первый массив, SecondMatr() – второй массив, ResMatrix() – результирующий массив. Все массивы размерности N на N;

   Transp(InputMatrix() As Double) - Процедура транспонирования матрицы. Входной параметр InputMatrix() – массив опять же размерности N на N;

   Swap(A As Double, B As Double) - Процедура , меняющая значения две переменных A и B

   MyGenerate() - Процедура формирующая диагональную матрицу и вращающая ее произвольным образом. Результатом работы этой процедуры будет симметричная матрица;

   Main() - Из названия видно, что это главная процедура. Здесь реализован алгоритм!

   Sp() As Double - Функция вычисления следа матрицы. Напомню, след должен сохраняться (см. ПРИМЕЧАНИЕ)

   Private Declare Sub Sleep Lib "kernel32" (ByVal dwMilliseconds As Long) - Это API-шная функция, используется для задержки выполнения.

2.3 Тестирование

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

2.4 Листинг

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

Заключение

Рассмотренный метод вращения решает полную проблему собственных

значений и собственных векторов матриц (симметрических) в том смысле, что

определяются все собственные значения и собственные векторы.

При решении полной проблемы собственных значений для

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

Кроме того, вычислительная техника может иметь далеко не рекордное быстродействие и объём оперативной памяти, и заведомо конечную разрядность двоичного представления чисел и связанные с этим ненулевые вычислительные погрешности. Поэтому итерационные методы, и метод вращений в частности, получили большое применение в решении систем линейных алгебраических уравнений. Современная вычислительная техника позволяет проводить исследование устойчивости и сходимости итерационного метода в зависимости от параметров задачи.

 

 

 

 

 

 

 

 

 

 

 

 

Список литературы

1 Вержбицкий В.М. Численные методы (математический анализ и обыкновенные дифференциальные уравнения). – М.: Высшая школа, 2001.

2.Демидович Б.П., Марон  И.А. Численные методы анализа. –  М.: Наука, 1970.

3.Ильин В.А., Позняк Э.Г. Линейная алгебра М.: 1974, «Наука»