«Р. Г. Козин АЛГОРИТМЫ ЧИСЛЕННЫХ МЕТОДОВ ЛИНЕЙНОЙ АЛГЕБРЫ И ИХ ПРОГРАММНАЯ РЕАЛИЗАЦИЯ Учебно-методическое пособие Москва 2012 УДК [512.64:004.4]:519.6(07) ББК 22.143я7+22.19я7+32.973.202-018.я7 К59 Козин Р.Г. Алгоритмы ...»

МИНИСТЕРСТВО ОБРАЗОВАНИЯ И НАУКИ РФ

НАЦИОНАЛЬНЫЙ ИССЛЕДОВАТЕЛЬСКИЙ ЯДЕРНЫЙ

УНИВЕРСИТЕТ «МИФИ»

Р. Г. Козин

АЛГОРИТМЫ ЧИСЛЕННЫХ МЕТОДОВ

ЛИНЕЙНОЙ АЛГЕБРЫ

И ИХ ПРОГРАММНАЯ РЕАЛИЗАЦИЯ

Учебно-методическое пособие

Москва 2012

УДК [512.64:004.4]:519.6(07) ББК 22.143я7+22.19я7+32.973.202-018.я7 К59 Козин Р.Г. Алгоритмы численных методов линейной алгебры и их программная реализация: Учебно-методическое пособие. – М.: НИЯУ МИФИ, 2012. – 192 с.

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

Пособие предназначено для студентов группы К4-331, обучающихся по специальности «Прикладная математика», для методической поддержки лекций, практических занятий и лабораторного практикума по курсу «Специальные главы математики».

Книга будет полезно инженерам и аспирантам при решении конкретных практических задач линейной алгебры.

Рецензент канд. физ.-мат. наук, доцент Н. В. Овсянникова © Национальный исследовательский ISBN 978-5-7262-1711- ядерный университет «МИФИ», Оглавление Предисловие

Введение

Глава 1. Нормы вектора и матриц, мера обусловленности матрицы........ 1.1. Норма вектора

1.2. Норма матрицы

1.3. Мера обусловленности и ранг матрицы

Глава 2. Методы решения систем линейных уравнений

2.1. Итерационные методы решения систем линейных уравнений... 2.2. Прямые методы решения системы линейных уравнений............ 2.3. Метод регуляризации для решения плохо обусловленных систем.

2.4. Решение систем с прямоугольными матрицами

Глава 3. Методы решения задачи на собственные значения................. 3.1. Собственные значения и собственные вектора матрицы........... 3.2. Решение частичной проблемы собственных чисел для симметричной матрицы.

3.3. Решение задачи на собственные значения методом Данилевского

3.4. Метод Крылова

3.5. Определение собственных векторов в общем случае................. 3.6. Метод вращения

Глава 4. Алгоритмы для работы с полиномами

4.1. Подсчет значения полинома и деление на множители............... 4.2. Локализация корней полинома

4.3. Метод парабол для нахождения всех корней полинома............ Лабораторный практикум

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

Приложение. Некоторые сведения о математическом пакете wxMaxima

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

В книге приведены методы решения следующих задач линейной алгебры:

• определение нормы вектора и матрицы, меры обусловленности матрицы;

• нахождение решения системы линейных уравнений;

• вычисление определителя матрицы;

• определение обратной матрицы;

• нахождение максимального и минимального собственных чисел матрицы;

• нахождение характеристического многочлена матрицы;

• определение всех собственных чисел и векторов матрицы;

• выделение линейного и квадратичного полиномов из заданного полинома;

• локализация и последующее уточнение корней полинома.

Каждая тематический раздел завершается несколькими вопросами и заданиями для самоконтроля. В заключении приведены темы заданий для лабораторного практикума и приложение с описанием некоторых встроенных команд пакета wxMaxima.

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

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

Наиболее востребованными задачами линейной алгебры являются следующие:

1. Найти вектор х с координатами хi, i = 1 n, являющейся решение системы линейных уравнений где А ( Aij, i, j = 1 n ) – исходная матрица; b ( bi, i = 1 n ) – заданный вектор правых частей системы.

2. Для заданной матрицы А вычислить определитель det(A).

3. Найти матрицу А–1, обратную заданной А, для которой выполняется матричное равенство A1 A = E, где Е – единичная матрица ( Eij = ij ).

4. Для заданной матрицы А определить собственные числа i и соответствующие им собственные векторы хi, которые удовлетворяют матричному уравнению Ах = х.

5. Нахождение корней полинома:

выделение из полинома линейных и квадратичных множителей.

МЕРА ОБУСЛОВЛЕННОСТИ МАТРИЦЫ

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

Норма вектора ||x|| – положительная скалярная величина, вычисляемая через его компоненты xi, i = 1n, и явно или косвенно характеризующая его длину. Как и длина вектора, она должна удовлетворять соотношениям:

Приведем три способа определения нормы вектора.

1. Кубическая норма 2. Октаэдрическая норма 3. Сферическая норма Концы множества векторов, нормы которых удовлетворяют равенствам x куб = 1, x окт = 1 и x сф = 1, нарисуют в п-мерном пространстве соответственно поверхности п-мерного куба, п-мерного октаэдра и п-мерной сферы. Это обстоятельство и объясняет их название.

Легко проверить, что все указанные нормы удовлетворяют трем соотношениям (1.1).

Из рис. 1.1, на котором для двумерного случая построены куб, октаэдр и сфера, соответствующие уравнениям: x куб = 1, x окт = = 1, видно, что для любого вектора х между этими нормами выполняются соотношения Справедливость соотношения легко подтверждается на примере конкретного векторов не следует равенство самих векторов, в то время как равные вектора всегда имеют равные нормы.

Приведем алгоритмы вычисления двух норм x куб и x окт (поскольку они одновременно демонстрируют, как можно найти максимальное значение и сумму любой последовательности чисел):

Норма матрицы A – положительная скалярная величина, которая устанавливает соответствие между нормами исходного вектора x и вектора, являющегося результатом воздействия на исходный вектор матрицы Ax. При этом в зависимости от критерия, используемого для нахождения нормы матрицы:

определяются два вида норм: согласованная и подчиненная (наименьшая из всех согласованных).

Выведем нормы матрицы подчиненные, соответственно, кубической, октаэдрической и сферической нормам вектора.

Матрица AT A симметричная и положительно определенная (так как ( AT A)T = AT ( AT )T = AT A и ( AT Ax, x) = ( Ax, Ax) > 0 ). Поэтому она имеет полную линейно-независимую систему собственных векторов x i, i = 1 n, а ее собственные числа действительные и положительные: 1 2... n > 0.

Представим произвольный вектор х в виде разложения по этим векторам x = ci x i и подставим его в верхнее скалярное произвеi дение Отсюда Замечание. Поскольку для симметричной матрицы AT A = A2, то В качестве примера нормы матрицы согласованной со сферической нормой вектора приведем евклидову норму Она связана со сферической нормой следующим образом:

A сф A евк (рис. 1.2).

Замечание. Приведенные ниже неравенства устанавливают нижнюю и верхнюю грани для всех указанных норм где 1 – максимальное по модулю собственное число матрицы, На практике, если значение нормы матрицы выходит за пределы этого интервала, то это указывает о том, что норма матрицы определена неверно. Кроме того, эти границы позволяют грубо оценить само значение нормы матрицы.

В заключение приведем алгоритмы вычисления кубической нормы матрицы (для октаэдрической нормы в этом алгоритме необходимо поменять местами индексы в операторах цикла) и произведения B = AT A, которое используется при определении сферической нормы матрицы:

На рис. 1.2 показана экранная форма программы, которая для введенной матрицы вычисляет нормы и меру обусловленности ( ( A) = A1 * A ). Результаты расчета подтверждают, что A сф A евк и для любой из указанных норм матрицы справедливо неравенство A < n max Aij = 5 0,92685 = 4,63425.

Для расчета сферической нормы необходимо предварительно умножить матрицу на транспонированную, а затем методом вращения найти все собственные числа полученной симметричной положительно определенной матрицы. Собственные числа этой матрицы формируются на месте ее диагональных компонент. (Метод вращения описывается в п 3.6.) Дополнительно приведем листинг функций и подпрограмм для математического пакета wxMaxima, которые дублируют вычисления, отображенные на рис. 1.1. Результаты их исполнения представлены на рис. 1.3 и 1.4.

/* функции вычисляют нормы матрицы: кубическая, октаэдрическая и евклидова */ nrm1(A):=lmax(create_list(sum(abs(A[i,j]),j,1,length(A)),i,1,length(A)));

nrm2(A):=lmax(create_list(sum(abs(A[i,j]),i,1,length(A)),j,1,length(A)));

nrm4(A):=sqrt(sum(sum(A[i,j]^2,j,1,length(A)),i,1,length(A)));

/* подпрограмма умножает транспонированную матрицу на исходную /* результат возвращает в матрице C2 */ Ct_C(C):=block [i,j,k,n,Ct], /* локальные переменные подпрограммы */ n:length(C), /* определям размерность матрицы */ C2:zeromatrix(n,n), /* формируем нулевую матрицу, размерности n */ Ct:transpose(C), /* трансформируем матрицу */ for j thru n do for k thru n do C2[i,j]:C2[i,j]+Ct[i,k]*C[k,j] /* подпрограмма находит положение максимального по модулю внедиагонального */ /* элемента матрицы; результат возвращается в переменных: lmax, kmax, Cmax,fi */ maxElem(C):=block [i,j,n,buf], /* локальные переменные */ Cmax:-1,lmax:0, kmax:0, n:length(C), buf:abs(C[i,j]), if buf>Cmax then (Cmax:buf,lmax:i,kmax:j) if C[lmax,lmax]=C[kmax,kmax] then fi:signum(C[lmax,kmax])*%pi/ else fi:0.5*atan(2*C[lmax,kmax]/(C[lmax,lmax]-C[kmax,kmax])) /* подпрограмма реализует метод вращения для симметричной матрицы /* iter, eps – число итераций и точность */ lmd(A,iter,eps):=block [i,j,n,it,buf], /* локальные переменные */ n:length(A),maxIt:0, for it thru iter do maxElem(A), if Cmax 0, > 0, из схемы (2.1) мы получим схемы, соответствующие различным, известным итерационным методам. Сейчас же используем ее для решения следующих вопросов, касающихся всех итерационных схем: условие сходимости итерационного процесса, количественная оценка скорости сходимости, выбор критерия окончания итераций и признак расходимости итерационного процесса.

1. Первую проблему решает следующая теорема.

Теорема. Для сходимости процесса (2.1) необходимо и достаточно, чтобы все собственные значения матрицы В были меньше единицы: i ( B ) < 1.

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

Следствие. Для сходимости итерационного процесса (2.1) достаточно выполнения неравенства B < 1.

Доказательство. Пусть Bx = x. Тогда B x Bx = x. Отсюда следует неравенство 1 > B, которое согласно предыдущей теореме гарантирует сходимость итерационного процесса.

2. Получим оценку скорости сходимости итерационного процесса при условии B < 1.

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

4. Итерационный процесс x ( k +1) = Bx ( k ) + b можно считать расходящимся, если неравенство выполняется подряд для трех–пяти последовательных значений k.

Справедливость этого утверждения для одномерного случая (п = 1) демонстрирует рис. 2.1.

Рассмотрим итерационные схемы, которые вытекают из общей схемы (2.1).

1. Метод Якоби получается при выборе параметров: = 1, C = D, где D – диагональная составляющая исходной матрицы A = AL + D + AH, где AL, AH – соответственно ее нижняя и верхняя составляющие.

Отсюда следует С учетом очевидного вида обратной диагональной матрицы легко выводится покоординатная запись схемы метода Якоби Достаточное условие сходимости B < 1 для метода имеет вид этом имеем покомпонентная запись схемы и достаточное условие сходимости метода Если матрица А симметрична и положительно определена (ее собственные числа 1 2... n > 0 ), то для нее можно указать оптимальный параметр при котором метод простой итерации всегда сходится, причем с максимальной скоростью, так как в этом случае сферическая норма матрицы В принимает минимально возможное значение Покажем справедливость этих утверждений.

Если матрица А – симметричная, то B = E A также симметричная, и поэтому Поскольку А еще и положительно определенная матрица, то ее собственные числа i ( A) > 0, и поэтому все функции fi () = = 1 i ( A) будут располагаться внутри зоны, очерченной двумя ломаными линиями 1 1 ( A) и 1 n ( A) (рис. 2.2). Сферической норме B () = max 1 i ( A) как функции параметра на этом рисунке будут соответствовать верхние участки ломаных линий. Видно, что норма принимает минимальное значение при опт, которое задается координатой точки пересечения ломаных линий жения (2.14) и (2.15).

Замечание. Если матрица плохо убусловлена, то 1 >> n. Тогда ции для системы будет сходиться очень медленно.

3. В методе верхней релаксации принимается > 0, C = D +AL.

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

Однако доказано, что если матрица А симметричная и положительно определенная, то метод верхней релаксации всегда сходится Частный вариант метода верхней релаксации при = 1 известен как метод Зейделя. В том случае выражения (2.16) и (2.18) преобразуются к виду 2.1.3. Алгоритмическая реализация методов Рассмотрим особенности реализации алгоритмов.

1. В качестве нулевого вектора берется вектор правой части матричного уравнения x (0) = b.

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

3. В алгоритме автоматически отлавливается «несходимость»

итерационного процесса для заданной системы с помощью критерия – увеличение нормы norm = x ( k +1) x ( k ) в течение пяти покуб следовательных итераций. Если последнее имеет место, то переменной code присваивается значение двойки и осуществляется выход из алгоритма.

4. Для окончания процесса используются два критерия: число итераций ограничено сверху заданным максимальным числом K и norm <, где K, задаются пользователем. Окончание процесса согласно одному из этих критериев соответствует code = 0. При этом в переменных K, возвращаются значения использованного числа итераций и достигнутой точности в оценке нормы разности соседних приближений, а в векторе x – полученное решение.

5. Используются: буферный вектор x0i, i = 1,n; переменные norm_old и fl для хранения соответственно нормы, вычисленной на предыдущей итерации, и числа итераций, при которых последовательно возрастает значение этой нормы.

6. exit sub – обозначает выход из подпрограммы.

Алгоритм метода Якоби: xi( k +1) = Входные параметры процедуры: A, b, n, x, K,..

for i = 1, n нормировка уравнений, задание x if Aii < 1030 ( 0) then { code = 1, exit sub} for k = 1, K цикл по итерациям Алгоротм простой итерации:

Входные параметры процедуры: A, b, n, x,, K,.

for i = 1, n нормировка уравнений, задание x Алгоритм метода верхней релаксации:

Входные параметры процедуры: A, b, n, x,, K,. Используются буферные переменные: 1, buf, nt.

for i = 1, n нормировка уравнений, задание x if Aii < 1030 ( 0) then { code = 1, exit sub} for k = 1, K цикл по итерациям Была написана программа, которая позволяет применить любой итерационный метод для решения заданной пользователем системы.

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

Чтобы матрицу системы сделать симметричной и положительно определенной, систему следует умножить на исходную транспонированную матрицу. Если система не является плохо обусловленной, то такое преобразование обеспечит сходимость итерационного процесса. Для запуска выбранного итерационного метода необходимо задать параметр, максимальное число итераций и требуемую точность > x ( k ) x ( k 1). По окончании итерационного процесса в таблице выдается решение и невязка ( = b Ax ( k ) ).

Замечание. В принципе, в качестве критерия окончания итерационного процесса можно было бы использовать условие = b Ax ( k ) < eps. Этот критерий является более качественным, но на каждой итерации он потребует значительного объема дополнительных операций. К тому же, как видно из таблицы на рис. 2.4., применяемый критерий обеспечивает выполнение и этого критерия.

С помощью программы было проведено исследование сходимости методов простой итерации и верхней релаксации для модельной системы (см. таблицу на рис. 2.3). Результаты исследования представлены в табл. 2.1 и 2.2.

Таблица 2.1. Зависимость предельного числа итераций от параметра (при = 0, 0000001 ) для метода верхней релаксации Таблица 2.2. Зависимость предельного числа итераций от параметра (при = 0, 0000001 ) для метода простой итерации при итерации 5905 3310 1733 1419 1204 1105 1359 расходится Характер зависимости, зафиксированной в табл. 2.2, хорошо объясняет рис. 2.5, где показано местоположение функций f i () = 1 i ( A) и опт для метода простой итерации.

Для преобразованной системы метод Якоби также сходится.

При этом ее решение определяется с заданной точностью за полученной с учетом (2.8), 1,2 = ±0,99 < 1 < B куб = 1, 4, т.е. «необходимое и достаточное условие» выполнено, а «только достаточное условие» – нет.

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

0,3 1 x = 2. Результаты, полученные с помощью программы для данной системы, приведены в табл. 2.3.

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

Словесная (вербальная) формулировка задачи: при заданной плотности источников q(x,y) найти распределение температуры T ( x, y ) в плоской прямоугольной области, две границы которой теплоизолированы, а через другие происходит теплообмен.

Математически эта задача записывается следующим образом:

Здесь k, T0, x, y – соответственно коэффициент теплопроводности, температура окружающей среды и коэффициенты поверхностной теплопроводности.

Для уменьшения зависимости задачи от исходных параметров перейдем к новым величинам:

С учетом этих изменений исходная задача примет вид Приведенную задачу будем решать разностным методом. Для этого область покроем прямоугольной сеткой (рис. 2.6), а уравнение теплопроводности и граничные условия запишем через значения Tij искомой функции T(x,y) в узлах сетки Разностная задача аппроксимирует исходную с точностью порядка O( hx2, hy ). Последнее легко проверить, если в разностные аппроксимации всех производных подставить соответствующее разложение Фурье, подобное данному Исключив с помощью граничных условий из разностных уравнений значения Tij, связанные с внешними узлами, и разрешив после этого все линейные уравнения относительно диагональных элементов разреженной матрицы, получим систему уравнений в форме, соответствующей методу Зейделя.

Процедура реализация метода Зейделя для данной системы выглядит следующим образом:

- все узлы сетки обрабатываются слева направо слой за слоем, начиная с нижнего слоя (i = 0);

- выход из итерационного процесса происходит по критерию max Tijk Tijk 1 < eps, где k – номер итерации, или по исчерпание заданного числа итераций.

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

kill(all);

/* подпрограмма рассчитывает коэффициенты разностной разреженной системы */ /* линейных уравнений */ coef(hx,hy,alfx,alfy):=block M:zeromatrix(9,3), M[1,3]:hx^2*hy^2/(2*(hx^2+hy^2)), M[1,1]:M[1,3]*2/hx^2, M[1,2]:M[1,3]*2/hy^2, M[2,3]:M[1,3], M[2,1]:M[2,3]/hx^2, M[2,2]:M[1,2], M[3,3]:hx^2*hy^2/(2*(hx^2+hy^2+alfx*hx*hy^2)), M[3,1]:M[3,3]*2/hx^2, M[3,2]:M[3,3]*2/hy^2, M[4,3]:M[1,3], M[4,1]:M[1,1], M[4,2]:M[4,3]/hy^2, M[5,3]:M[1,3], M[5,1]:M[5,3]/hx^2, M[5,2]:M[4,2], M[6,3]:M[3,3], M[6,1]:M[3,1], M[6,2]:M[6,3]/hy^2, M[7,3]:hx^2*hy^2/(2*(hx^2+hy^2+alfy*hy*hx^2)), M[7,1]:M[7,3]*2/hx^2, M[7,2]:M[7,3]*2/hy^2, M[8,3]:M[7,3], M[8,1]:M[8,3]/hx^2, M[8,2]:M[7,2], M[9,3]:hx^2*hy^2/(2*(hx^2+hy^2+alfx*hx*hy^2+alfy*hy*hx^2)), M[9,1]:M[9,3]*2/hx^2, M[9,2]:M[9,3]*2/hy^ define(f(x,y),(2-x)*(1-y)); /* заданная функция тепловыделения */ /* подпрограмма рассчитывает значения плотности тепловыделения в узлах сетки */ create_q(Nx,Ny,hx,hy):=block [i,j,x,y], q:zeromatrix(Nx+1,Ny+1), x:hx*(i-1), for j thru Ny+1 do (y:hy*(j-1), q[i,j]:f(x,y)) /* подпрограмма реализует итерационный метод Зейделя для разностной системы */ /* линейных уравнений */ zadel(Nx,Ny,iter,eps):=block [i,j,k,buf,nt], it:0, T:zeromatrix(Nx+1,Ny+1), /* обработка нижнего слоя сетки */ buf:M[1,1]*T[2,1]+M[1,2]*T[1,2]+M[1,3]*q[1,1], nt:abs(T[1,1]-buf), T[1,1]:buf, if nt>norm then norm:nt, buf:M[2,1]*(T[i-1,1]+T[i+1,1])+M[2,2]*T[i,2]+M[2,3]*q[i,1], nt:abs(T[i,1]-buf), T[i,1]:buf, if nt>norm then norm:nt buf:M[3,1]*T[Nx,1]+M[3,2]*T[Nx+1,2]+M[3,3]*q[Nx+1,1], nt:abs(T[Nx+1,1]-buf), T[Nx+1,1]:buf, if nt>norm then norm:nt, /* обработка промежуточных слоев сетки */ buf:M[4,1]*T[2,j]+M[4,2]*(T[1,j-1]+T[1,j+1])+M[4,3]*q[1,j], nt:abs(T[1,j]-buf), T[1,j]:buf, if nt>norm then norm:nt, buf:M[5,1]*(T[i-1,j]+T[i+1,j])+M[5,2]*(T[i,j-1]+T[i,j+1])+M[5,3]*q[i,j], nt:abs(T[i,j]-buf), T[i,j]:buf, if nt>norm then norm:nt buf:M[6,1]*T[Nx,j]+M[6,2]*(T[Nx+1,jT[Nx+1,j+1])+M[6,3]*q[Nx+1,j], nt:abs(T[Nx+1,j]-buf), T[Nx+1,j]:buf, if nt>norm then norm:nt /* обработка последнего слоя сетки */ buf:M[7,1]*T[2,Ny+1]+M[7,2]*T[1,Ny]+M[7,3]*q[1,Ny+1], nt:abs(T[1,Ny+1]-buf), T[1,Ny+1]:buf, if nt>norm then norm:nt, buf:M[8,1]*(T[iNy+1]+T[i+1,Ny+1])+M[8,2]*T[i,Ny]+M[8,3]*q[i,Ny+1], nt:abs(T[i,Ny+1]-buf), T[i,Ny+1]:buf, if nt>norm then norm:nt buf:M[9,1]*T[Nx,Ny+1]+M[9,2]*T[Nx+1,Ny]+M[9,3]*q[Nx+1,Ny+1], nt:abs(T[Nx+1,Ny+1]-buf), T[Nx+1,Ny+1]:buf, if nt>norm then norm:nt, if A max < nul ( 0) then det = 0 ( матрица вырожденная) выход из процедуры if k max k then for j = k,(n + 1) (переставляем строки ) det = det* Akk (вычисление определителя) Матрица перестановок обладает следующими свойствами.

1. Матричная операция T (i, j ) A переставляет строки, например:

2. Матричная операция AT (i, j ) переставляет столбцы, например:

3. T (i, j )T (i, j ) = E, например, После приведения системы уравнений к системе с треугольной матрицей легко записать алгоритм решения последней (учтено, что правые части уравнений хранятся в дополнительном столбце расширенной матрицы – Ai ( n +1) ) if Ann < nul then {det = 0, exit sub} else det = det* Ann Здесь первый условный оператор завершает проверку матрицы и подсчет значения определителя.

2.2.2. Метод оптимального исключения Гаусса-Жордана В отличие от метода Гаусса в данном методе исходная система Ax = b с помощью п последовательных преобразований сводится к системе с единичной матрицей, решение которой располагается на месте правого крайнего столбца расширенной матрицы (добавлен столбец, состоящий из компонентов вектора b) преобразованной системы Ex LAx = Lb, где L = L( n )...L(1).

Приведем формулы, по которым выполняется k-е преобразование. К данному моменту текущая расширенная матрица исходной системы имеет следующей вид Первоначально компонент матрицы Аkk приводится к единице путем деления k-й правой подстроки на компонент Аkk:

Затем обнуляются все остальные компоненты k-го столбца расширенной матрицы, путем вычитания из каждой i-й строки k-й строки, умноженной на Aik:

Замечание. Легко проверить, что преобразования матрицы (2.27) и (2.28) эквивалентно матричной операции L( k ) A( k 1), где A( k 1) – исходная матрица после (k – 1)-го преобразования, а матрица L( k ) имеет вид (отличается от единичной матрицы, размерности пп, только k-м столбцом):

Приведем алгоритм, реализующий данный метод:

det = 1 (переменная для хранения значения определителя ) вызов процедуры выбора ведущего элемента Здесь:

1) в пределах каждого k-го преобразования k-й столбец явно не обрабатывается, поскольку он в последующих преобразованиях не используется;

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

(Программа, использующяя метод оптимального исключения, приведена в п. 2.2.4.) det = det* Akk (вычисление определителя) 2.2.3. Нахождение обратной матрицы методом оптимального Полное преобразование матрицы методом оптимального исключения может быть записано в матричном виде из которого следует матричная форма нахождения обратной матрицы (так как по определению A1 A = E ) где T(i) – матрица перестановок при i-ом преобразовании, а L( i ) задано (2.29).

Согласно выражению (2.31), если исходную матрицу А дополнить справа единичной матрицей Е размерности пп, то в качестве алгоритма нахождения обратной матрицы можно использовать алгоритм (2.30), заменив в нем правые пределы циклов с п + 1 на п + п.

При этом ту же замену необходимо выполнить и в процедуре выбора ведущего элемента (2.24) при перестановке строк.

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

Обсудим детали оптимального алгоритма расчета обратной матрицы. Обозначим соответствующие матрицы после (k – 1)-го преобразования Они имеют вид Обратим внимание на следующие обстоятельства.

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

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

Все это позволяет сделать вывод, что обе матрицы можно хранить на месте исходной и во время k-го преобразования обрабатывать по единым формулам (2.32), которые получены в результате обобщения матричных операций A(k ) = L(k ) ( A(k 1)T (k ) ), B(k ) = L(k ) B(k 1) :

выбор ведущего элемента в к ой правой подстроке ( 0) Ниже приводится оптимальный алгоритм в развернутом виде (если возвращается нулевое значение определителя матрицы, то матрица – вырожденная или плохо обусловлена):

det = 1, m = 0 (инициализация переменных) вызов процедуру выбора ведущего элемента в подстроке вызов процедуры, переставляющей строки в обработанной матреце (см. ниже) Алгоритм процедуры выбора ведущего элемента. Здесь для хранения номеров переставляемых столбцов используется массив M ij, i = 1, 2; j = 1, n, а в переменной т хранится число выполненных перестановок:

Процедура, переставляющая строки в преобразованной матрице:

if m = 0 then выход из процедуры for l = m,1 step = 1 ( переставляем строки в обратном порядке) В заключение в качестве примера, приведем экранную форму программы (рис. 2.11), написанной на Qt, на которой сформирована обратная матрица для исходной матрицы 4 5 0.

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

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

Далее полученную матрицу можно протестировать, проверев выполнение соотношения A(обр)*А=E. Результаты тестирования полученной обратной матрицы показаны на рис. 2.12.

Для сравнения приведен листинг аналогичной консольной программы, состаленной на языке Python:

# -*- coding: cp1251 -*import random # --- процедура выбора ведущего элемента --def select_Amax(k):

global m,det,n for i in xrange((n-1),-1,-1):

print "Определитель исходной матрицы = ","%.5f" % det print "Число использованных перестановок столбцов = ",n print "Получена следующая обратная матрица А(обр):" print "%.5f" % A[i][j], #формат – вещественное, 5 знаков print "Матрица ВЫРОЖДЕННАЯ" # --- Главная программа --print "Нахождение обратной матрицы методом Гаусса-Жордано с выбором ведущего элемента" fl_end= while fl_end:

A=[] #исходная матрица – !!! это надо делать при каждом повторении B=[] #для хранения копии матрицы E=[] #используется при тестировании точности обратной матрицы (буфер для единичной матрицы) m=raw_input("Введите размерность матрицы А: ") #завершается нажатием – ENTER ans=raw_input("Компоненты матрицы задать случайным образом (1-ДА/2-НЕТ)? ") for i in xrange(0,m): #ввод произвольной матрицы пользователем str1="Введите компанент матрицы – А("+str(i)+","+str(j)+"):

ans=raw_input("Показать сформированную матрицу (1-ДА/2-НЕТ)? ") if ans=="1":

print "Сформирована следующая матрица А:" for i in xrange(0,m):

print "%.5f" % A[i][j], #формат – вещественное, 5 знаков после точки det= L=[] #задание вспомогательного массива для хранения номеров переставляемых столбцов for i in xrange(0,2):

for j in xrange(0,m):

ar1.append(-1) L.append(ar1) ans=raw_input("Протестировать точность обратной матрицы [т.е. ? – А(обр)*А=Е] (1-ДА/2-НЕТ)? ") if ans=="1":

print "Произведение матриц А(обр)*А (с точностью до 7-го знака):" for i in xrange(0,m): # расчитываем произведение А(обр)*А print "%.7f" % E[i][j], #формат – вещественное, 7 знаков после точки ans=raw_input("Выйти из программы (1-ДА/2-НЕТ)? ") if ans=="1": fl_end= Ниже показан результат работы консольной программы. Обратите внимание, что в языке Python для объединения команд в блоки используется одинаковое число (обычно четыре) левых пробелов.

Кстати, это соответствует правилу «хорошего тона» в программировании.

Нахождение обратной матрицы методом Гаусса–Жордано с выбором ведущего элемента Введите размерность матрицы А: Компоненты матрицы задать случайным образом (1-ДА/2-НЕТ)? Введите компанент матрицы – А(0,0): Введите компанент матрицы – А(0,1): Введите компанент матрицы – А(0,2): Введите компанент матрицы – А(1,0): Введите компанент матрицы – А(1,1): Введите компанент матрицы – А(1,2): Введите компанент матрицы – А(2,0): Введите компанент матрицы – А(2,1): Введите компанент матрицы – А(2,2): Показать сформированную матрицу (1-ДА/2-НЕТ)? Сформирована следующая матрица А:

0.00000 2.00000 3. 4.00000 5.00000 0. 0.00000 6.00000 8. Определитель исходной матрицы = 8. Число использованных перестановок столбцов = Получена следующая обратная матрица А(обр):

5.00000 0.25000 -1. -4.00000 0.00000 1. 3.00000 0.00000 -1. Протестировать точность обратной матрицы [т.е. ? – А(обр)*А=Е] (1ДА/2-НЕТ)? Произведение матриц А(обр)*А (с точностью до 7-го знака):

1.0000000 -0.0000000 -0. 0.0000000 1.0000000 0. 0.0000000 0.0000000 1. Выйти из программы (1-ДА/2-НЕТ)? Наконец, аналогичная работа выполнена в пакете wxMaxima.

Листинг команд пакета wxMaxima, выполняющих расчет обратной матрицы и ее тестирование:

numer:true; /* выдавать числа в десятичном виде */ fpprintprec:5; /* в выводимом числе 5 цифр, не считая точки */ M:matrix([0,2,3],[4,5,0],[0,6,8]); /* задаем исходную матрицу */ mulMatr(A,B):=block /* подпрограмма перемножает матрицы A*B */ [i,j,k,n], /* локальные переменные подпрограммы */ E:zeromatrix(n,n), /* генерируем нулевую «глобальную» матрицу */ for k thru n do E[i,j]:E[i,j]+A[i,k]*B[k,j] mulMatr(OM,M); /* вызов подпрограммы */ Результаты выполнения команд, представленных выше, показаны на рис. 2.13. Сравните их с рис. 2.11 и 2.12 и резултатами работы программы на языке Python.

1. В чем суть метода исключения Гаусса?

2. Для каких целей в прямых методах используется процедура «выбора ведущего элемента»?

3. Запишите алгоритм для k-го преобразования исходной матрицы в методе Гаусса.

4. Запишите алгоритм нахождения решения системы с нижней треугольной матрицей.

5. Чем метод Гаусса отличается от оптимального метода исключения Гаусса-Жордана?

6. Что такое матрица перестановок и каковы ее свойства?

7. Напишите алгоритм, который меняет местами i-ю и k-ю строки матрицы А.

8. Каким образом выбирается ведущий элемент при нахождении обратной матрицы методом Гаусса–Жордана?

9. Почему при нахождении обратной матрицы необходимо запоминать номера переставляемых столбцов при выборе ведущего элемента?

Представим систему Ах = b в виде где Ain +1 = bi.

Если ввести (п + 1)-мерные векторы Ai = ( Ai1,.., Ain, Ain +1 ), i = 1, n ;

y = ( x1,..., xn,1), то эти уравнения можно заменить системой скалярных произведений Соотношения (2.34) позволяют дать другую формулировку для исходной задачи: в (п + 1)-мерном пространстве R n +1 найти (п + 1)мерный вектор у с уп+1 = 1, перпендикулярный п-мерной гиперплоскости, «натянутой» на векторы Ai, i = 1, n. Если удастся построить такой вектор, то первые п компонент этого вектора будут являться компонентами решения исходной системы уравнений.

Предлагается следующая схема решения этой задачи.

1. К исходной расширенной матрице добавляем снизу строку с компонентами (0,0,…,1). Если матрица исходной системы невырожденная, то строки этой новой матрицы в пространстве R n +1 образуют полную линейно-независимую систему векторов. Это следует из того обстоятельства, что ее определитель равен определителю исходной матрицы, который для невырожденной системы отличен от нуля.

2. Применяя к указанной полной линейно-независимой системе векторов процедуру ортогонализации Грама–Шмидта (2.35), строим ортонормированный базис для пространства R n +1 (рис. 2.14).

Последний орт этой системы (последняя строка расширенной матрицы) будет перпендикулярен п-мерной гиперплоскости, «натянутой» на векторы Ai, i = 1, n.

Детализируем отдельные операции алгоритма процедуры ортогонализации:

- расчет скалярного произведения sk (i, k ) = ( Ai, Ak ) = Aij Akj :

- нормировка вектора A k = k Нулевое значение подкоренного выражения означает линейную зависимость системы исходных векторов, т.е. вырожденность матрицы системы Ax = b.

- корректировка вектора A i = Ai ( Ai, Ak ) Ak 3. Поделив последнюю строку матрицы на ее последнюю компоненту, получаем решение исходной системы:

Замечание. Алгоритм (2.36) необходимо поставить вместо последнего оператора в схеме (2.35), которая является общей записью алгоритма ортогонализации.

На рис. 2.15 приведено основное окно экранной формы программы, которая с помощью рассмотренных выше алгоритмов (пп. 2.2. и 2.2.4) находит решение и компоненты вектора невязки ( i = bi Aij x j ) для заданной системы уравнений Ax = b.

Здесь предусмотрены возможность корректировки значения машинного нуля (переменной nul ), используемого для сранения действительных чисел с нулем, а также включение и отключение дополнительного критерия плохо обусловленной матрицы куб < nul (наряду с критериями Akk < nul или norm < nul соответственно для методов Гаусса–Жордана или метода ортогонализации).

На рис. 2.16 приведена экранная форма справки для указанной программы.

Рис. 2. Рис. 2. На рис. 2.17 демонстрируется применение метода оптимального исключения к системе с вырожденной матрицой (ее определитель равен нулю) при nul = 1E 30. Благодаря «пропуску» ошибок округления, удалось получить для данной системы одно из ее частных решений x = (0,3,0) (например, другое – x = (1,1,1) и т.д.).

Приведенные на рис. 2.15–2.17 результаты использования дополнительных возможностей для системы уравненей с вырожденной матрицей. Сранение этих результатов наглядно показывает, к каким существенным измениниям в решении х приводит незначительная входная погрешность в правой части системы b.

Результат применения метода исключения к системе с вырожденной матрицей и отключенным вторым критерием по невязке показан на рис. 2.18. Небольшое изменение правой части системы b привело к полной перестройке решения системы х (сравните с рис.

2.17).

В заключение для сравнения приведено решение предыдущих задач в рамках пакета wxMaxima.

Листинг программы для wxMaxima:

fpprintprec:5;

A:matrix([0,1,2],[4,5,0],[0,7,8]); /* "хорошая" матрица системы */ /* находим решение "хорошей" системы */ solve([x2+2*x3=3,4*x1+5*x2=6,7*x2+8*x3=9],[x1,x2,x3]);

B:matrix([1,2,3],[4,5,6],[7,8,9]); /* "плохая" матрица системы */ determinant(B);

/* находим решение "плохой" системы */ solve([x1+2*x2+3*x3=6,4*x1+5*x2+6*x3=15,7*x1+8*x2+9*x3=24],[x1,x2,x 3]);

Полученные результаты показаны на рис. 2.19. Из общего решения вырожденной системы легко получить любое частное решение, напрмер: при %r4=0 x=(0,3,0), при %r4=1 x=(1,1,1) и т.д. (сравните с рис. 2.17).

вырожденной системы 2.2.5. Второй метод ортогонализации Этот метод применим к системам с симметричной положительно определенной матрицей. Поскольку только с помощью такой матрицы можно определить следующее скалярное произведение ( x, Ay ), удовлетворяющее всем свойствам, присущим скалярному произведению:

Например, проверим свойство 1, так как выполнение остальных свойств очевидно:

Видно, что для симметричной матрицы правые, а следовательно, и левые части этих выражений равны.

Теперь рассмотрим, как строится решение системы Ax = b данным методом.

1. Берем в пространстве Rn полную линейно-независимую систему векторов Ej, j = 1, n, в виде столбцов единичной матрицы (ее определитель равен 1) 2. Применяя процедуру ортогонализации Грама–Шмидта к данной системе векторов, строим в пространстве Rn ортонормированный базис в смысле скалярного произведения (х, Ау):

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

а перед вычислением корня необходимо проверять, чтобы подкоренное выражение положительное. Невыполнение этого условия означает, что матрица системы Ax = b не является положительно определенной.

Приведем алгоритм расчета скалярного произведения sk ( j, k ) ( E j, AE k ) = Eij Aim Emk с учетом хранения векторов в виде столбцов единичной матрицы:

3. Ищем решение системы в виде разложения по векторам построенного базиса, сформированного на месте единичной матрицы или через компоненты вектора х Для определения неизвестных коэффициентов разложения подставим разложение (2.37) в исходную систему Ax = b, а затем умножим ее скалярно последовательно на вектора Ei, i = 1, n:

В виду ортонормированности нового базиса в левой части уравнения (2.39) останется только одно слагаемое с j = i, равное ci. Отсюда следует, что Видно, что метод потребует значительного количества операций, но в методическом плане он представляет определенный интерес. По этой причине с ним полезно познакомиться.

2.2.6. Метод квадратного корня (метод Холецкого) Метод применим только для систем с симметричной положительно определенной матрицей.

Согласно этому методу исходная матрица представляется в виде произведения двух матриц где B – нижняя треугольная матрица.

Такое разложение возможно только для симметричной положительно определенной матрицы:

Чтобы найти рекуррентные формулы для определения компонент матрицы В, запишем матричное соотношение (2.41) для четырехмерного случая B11 B31 B31 B21 + B32 B22 B31 + B32 + B B11 B41 B41 B21 + B42 B22 B41 B31 + B42 B32 + B43 B33 B41 + B42 + B43 + B Отсюда, сравнивая компоненты левой и правой матриц, получим Здесь в двух верхних строчках под корнем явно присутствуют определители первых миноров исходной матрицы, а согласно критерию Сильвестра для положительно определенной матрицы они положительны. Последнее утверждение справедливо и для остальных подкоренных выражений. Это еще раз подтверждает необходимость того, чтобы матрица была положительно определенной.

Эти формулы по индукции позволяют записать следующий алгоритм расчета всех компонент матрицы В для п-мерного случая:

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

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

После того как разложение матрицы А выполнено система Ax = b распадается на две системы с треугольными матрицами или в покомпонентной записи Их решение рассчитывается элементарно Последовательное вычисление компонент промежуточного вектора у можно выполнять в алгоритме (2.42) после определения компонент соответствующей строки матрицы В (см. приведенный ранее алгоритм). При этом для экономии памяти матрицу В можно хранить на месте исходной матрицы А, а вектор у формировать на месте вектора х. Необходимо помнить, что все суммы, встречающиеся в алгоритмах, вычисляется путем введения соответствующих внутренних циклов.

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

Function kvroot () As Integer Dim i%, j%, k% kvroot = 0 'матрица положительно-определенная If Not (a(0, 0) > 0) Then kvroot = 1 'матрица не положительно-определенная Exit Function a(0, 0) = Sqr(a(0, 0)) x(0) = a(0, n) / a(0, 0) a(1, 0) = a(1, 0) / a(0, 0) x(1) = a(1, n) a(1, 1) = a(1, 1) – a(1, 0) * a(1, 0) If Not (a(1, 1) > 0) Then kvroot = 1 'матрица не положительно-определенная Exit Function a(1, 1) = Sqr(a(1, 1)) x(1) = (x(1) – a(1, 0) * x(0)) / a(1, 1) If n = 2 Then GoTo cont x(k) = a(k, n) – a(k, 0) * x(0) If Not (a(k, k) > 0) Then kvroot = 1 'матрица не положительно-определенная Exit Function a(k, k) = Sqr(a(k, k)) x(k) = x(k) / a(k, k) cont:

x(i) = x(i) / a(i, i) End Function Здесь систама хранится в расширенной матрице с компонентами: a(0,0),…a(n-1,n).

На рис. 2.20 представлен скриншот программы, которая позволяет найти решение произвольной системы линейных уравнений как вторым методом ортогонализации, так и методом квадратного корня. Для преобразования произвольной системы в равносильную систему с симметричной и положительно определенной матрицей следует использовать кнопку «Умножить систему на А(т)» – умножить систему на транспонированную исходную матрицу. Имеется возможность задать систему уравнений случайным образом. Это позволяет провести сравнительный анализ точности (по норме невязки полученного решения) используемых методов.

Метод прогонки является самым экономичным способом решения систем с трех диагональной матрицей Здесь для хранения системы использовано четыре вектора a, b, c, d, причем компоненты a1 = cn = 0 не используются.

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

Введем рекуррентные зависимости Найдем формулы для расчета входящих в (2.44) неизвестных параметров k, k.

Сравнивая первое уравнение в (2.43) с соотношением (2.44), сразу получаем Далее, заменим в k-ом уравнении системы (2.43) неизвестное xk 1 согласно формуле (2.44):

Отсюда легко выводим рекуррентные формулы для последовательного определения остальных неизвестных параметров k, k :

Вычисления по формулам (2.45), (2.46) называется прямым ходом прогонки.

Теперь, подставляя в последнее уравнение (2.43) соотношение (2.44) с k = n 1, получим Формулы (2.47) и (2.44) (с k = n 1, n 2,...,1 ) позволяют рассчитать все компоненты хk неизвестного решения системы (2.43). Эти расчеты называются обратным ходом прогонки.

Покажем методом индукции, что для устойчивости метода прогонки ( bk + ak k 1 > 0 ) достаточно выполнения следующих соотношений между параметрами системы При выполнении условия (2.48) имеем (см. (2.45)) b1 > и 1 < 1. Предположим, что k 1 < 1. Тогда bk + ak k 1 > 0. Кроме того, bk + ak k 1 > ck и, следовательно, k < 1 и т.д.

Заметим, что метод прогонки легко обобщается на системы с пятью (и более) диагональными матрицами. В этом случае исходное рекуррентное соотношение следует брать в виде Замечание. При программной реализации метода в целях экономии памяти параметры k, k можно сохранять на месте компонент векторов bk, d k, а решение вернуть в векторе с. В этом случае алгоритм метода будет иметь вид В заключение рассмотрим пример использования метода прогонки для нахождения решения краевой задачи для диффренциального уравнения второго порядка.

Пример. Распределение стационарной температуры в тепловыделяющем стержне с теплоизолированными боковыми стенками, один торец которого теплоизолирован, а через другой происходит теплообмен, в безрамерном виде задается уравнением и граничными условиями Разностный аналог этой задачи, аппроксимирующий ее с квадратичной точность относительно шага разностной сетки h =, заn писывается в виде Разностная сетка, используемая в разностной задаче, показана на рис. 2.21.

Отсюда для определения неизвестных значений yi для i = 0, n, получаем следующую систему линейных уравнений с трех диагональной матрицей Листинг программы в wxMaxima, которая решает систему методом прогоки для q( x) = x и = 1 :

kill(all);

/* подпрограмма, реализующая метод прогонки; a,b,c,d – списки с коэффициентами */ /* системы, n – порядок системы; решение возвращается в списке с */ progonka(a,b,c,d,n):=block /* прямой ход прогонки */ d[1]:d[1]/b[1], b[1]:-c[1]/b[1], for k:2 thru n-1 do (t:b[k]+a[k]*b[k-1], d[k]:(d[k]-a[k]*d[k-1])/t, b[k]:c[k]/t), /* обратный ход прогонки */ c[n]:(d[n]-a[n]*d[n-1])/(b[n]+a[n]*b[n-1]), for k:(n-1) thru 1 step -1 do c[k]:b[k]*c[k+1]+d[k] /* главная программа */ numer:true;

fpprintprec:6;

define(q(x),x); /* правая часть исходного уравнения */ /* формирование коэффициентов системы */ n:11; h:1/(n-1); alf:1;

a:[]; b:[]; c:[]; d:[]; t:[]; /* инициализация списков */ t:append(t,[0]); a:append(a,[0]); b:append(b,[-1]); c:append(c,[1]);

d:append(d,[-h^2*q(0)/2]);

for i:2 thru n-1 do t:append(t,[h*(i -1)]), a:append(a,[1]), b:append(b,[-2]), c:append(c,[1]), d:append(d,[-h^2*q(t[i])]) t:append(t,[1]); a:append(a,[1]); b:append(b,[-(1+h*alf)]); c:append(c,[0]);

d:append(d,[-h^2*q(1)/2]);

progonka(a,b,c,d,n);

c; /* выводим список значений сеточного решения, полученного методом прогонки */ r(x):=(4-x^3)/6; /* точное решение исходной задачи */ c1:[]; for i thru n do c1:append(c1,[r(t[i])]);

c1; /* выводим список значений точного решения */ /* строим график точного решений и значения сеточного уравнения */ plot2d([[discrete, t,c], r(x)], [x,0,1], [style, points, lines], [color, red, blue], [point_type, asterisk], [legend, "разностное решение", "точное"], [xlabel, "длина – x"], [ylabel, "температура"]);

Результаты выполнения программы даны на рис. 2.21 и 2.22.

Список значений решения разностного уранения:

Список значений точного решения в точках сетки:

Из рис. 2.22 видно, что значения разностной задачи хорошо согласуются с точным решением (графически и численно).

2.3. Метод регуляризации для решения Решение системы уравнений Ax = b с плохо обусловленной матрицей существенно зависит от ошибок в исходных данных и ошибок округления. На рис. 2.23 для таких систем схематично показано соотношение между областью ошибок в задании исходной системы уравнений и областью, в которой распологаются ее решения.

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

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

В заключение укажем способ выбора параметра. Обычно на практике проводится ряд вычислений с различными значениями.

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

На рис. 2.26 дано нормальное решение системы с плохо обусловленной матрицей. Заметим, что точное решение исходной системы равно (1; 1; 1). Видно, что решить непосредственно исходную систему методом Гаусса–Жордано не удалось.

2.4. Решение систем с прямоугольными матрицами Системы уравнений Ax = b, в которых матрицы A m*n имеют прямоугольную форму, делятся на недоопределенные (m < n) и переопределенные (m > n). Решение таких систем имеет свою специфику.

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

Итак, решений много, но рассмотрим, какое из них следует считать оптиальным.

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

Легко проверить, что оно является решением исходного уравнения В качестве доказательства оптимальности этого решения (в указанном выше смысле) приведем результаты вычислительного эксперимента, в котором расчитывались нормы нормального псевдорешения и совокупности частных решений. Последние строились путем обнуления всевозможных комбинаций «лишних» компонент вектора х. Напомним, что число таких комбинаций равно числу сочетаний из п элементов по (п – т):

В эксперименте использовалась недоопределенная система со случайной матрицей рамерности m = 3, n = 5, для которой можно построить 10 частных решений.

Листинг тестирующей программы на wxMaxima:

kill(all);

numer:true;

fpprintprec:6;

ratprint:false;

/* недоопределенная система */ b:matrix([1],[2],[3]);

g[i,j]:=random(1.0); A:genmatrix(g,3,5); /* генерация случайной матрицы A(3,5) */ rank(A);

At:transpose(A);

xo:At.invert(A.At).b; /* находим нормальное псевдорешение и его норму */ L_norm:[sqrt(sum(xo[i,1]*xo[i,1],i,1,5))];

for i thru 4 do /* определение всех возможных частных решений */ B:submatrix(A,i,j), x1:zeromatrix(5,1),x1[i,1]:1, x1[j,1]:1, x:invert(B).b, norm:sqrt(sum(x[i,1]*x[i,1],i,1,3)), L_norm:append(L_norm,[norm]), for k thru 5 do if x1[k,1]=0 then (x1[k,1]:x[k1,1],k1:k1+1) else x1[k,1]:0, xo:addcol(xo,x1) print("Список норм нормального псевдорешения и10 частных решений:"); L_norm;

print("Минимальная норма в списке = ",min_norm:lmin(L_norm));

print("Матрица, содержащая нормальное псевдорешение и 10 частные решения:");xo;

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

2.4.2. Решение переопределенных систем Переопределенная система может быть несовместной, т.е. не иметь точного решения. В этом случае в качестве ее решения используется обобщенное псевдорешение, которое строится таким образом, чтобы обеспечить минимальное значение нормы невязки для данной системы уравнений, т.е. минимизирует следующий функционал Согласно формуле (2.5) обобщенное псевдорешение является решением уранения Чтобы показать оптимальность обобщенного псевдорешения (в указанном выше смысле) был поставлен вычислительный эксперимент. В нем норма невязки обобщенного псевдорешения сравнивается с соответствующими нормами частных решений всевозможных систем, которые получаются из исходной системы путем отбрасывания «лишних» уравнений.

В эксперименте использовалась переопределенная система со случайной матрицей рамерности m = 5, n = 3. Она допускает построение 10 частных решений.

Листинг тестирующей программы на wxMaxima:

kill(all);

/* подпрограмма формирует случайную матрицу размерности m*n */ rndA(m,n):=block A:zeromatrix(m,n), for j thru n do A[i,j]:random(1.0) numer:true;

fpprintprec:6;

ratprint:false;

/* переопределенная система */ b:matrix([1],[2],[3],[4],[5]);

rndA(5,3); A;

rank(A); At:transpose(A);

xo:invert(At.A).At.b;

nvo:b-A.xo;

print("Норма невязки для обобщенного псевдорешения = ",sqrt(sum(nvo[i,1]*nvo[i,1],i,1,5)));

L_norm:[]; /* задаем нулевой список для сбора норм невязки частных решений */ for i thru 4 do /* определение норм невязки для всех возможных частных решений */ for j:i+1 thru 5 do B:submatrix(i,j,A), bt:submatrix(i,j,b), x:invert(B).bt, nv:b-A.x, norm:sqrt(sum(nv[i,1]*nv[i,1],i,1,5)), L_norm:append(L_norm,[norm]), xo:addcol(xo,x) print("Список норм невязки для частных решений:"); L_norm;

print("Минимальная норма невязки для частных решений = ",min_norm:lmin(L_norm));

print("Матрица, содержащая обобщенное псевдорешение и частные решения:"); xo;

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

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

1. При каком условии вектор y = ( x1,..., xn,1), используемый в первом методе ортогонализации, будет линейно-независимым от векторов-строк расширенной матрицы А?

2. Какие манипуляции с векторами выполняет процедура Грама– Шмидта?

3. Как с помощью процедуры Грама–Шмидта можно проверить: является ли система векторов линейно независимой?

4. Для каких систем линейных уравнений можно использовать второй метод ортогонализации и почему?

5. В чем сущность метода квадратного корня и почему он применяется только для систем с симметричной положительно определенной матрицей?

6. Запишите алгоритм модифицированного метода ортогонализации Грама–Шмидта.

7. Как подбирается параметр в методе регуляризации, который используется при решении плохо обусловленных систем линейных уравнений?

8. Запишите условие сходимости метода прогонки.

9. Составьте алгоритм метода прогонки для системы с пятидиагональной матрицей.

10. Распишите через компоненты скалярное произведение двух столбцов единичной матрицы sk ( j, k ) ( E j, AE k ) и составьте алгоритм для вычисления его значения.

11. Сравните норму нормального псевдорешения и нормы всех частных решений, полученных из исходной системы после приравнивания единице «лишних» компонент вектора х. Выполните задание в пакете wxMaxima.

НА СОБСТВЕННЫЕ ЗНАЧЕНИЯ

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

Для того чтобы система (3.1) имела нетривиальное решение, должно выполнятся равенство Уравнение (3.2), называемое характеристическим полиномом матрицы А, представляет собой полином степени п относительно и имеет п корней.

Приведем некоторые свойства собственных значений и векторов.

1. Матрица А порядка пп имеет п собственных значений.

2. Для каждого уникального существует, по крайней мере, один собственный вектор (так как не всегда k-кратное собственное число имеет k собственных векторов). Каждый собственный вектор определяет для данной матрицы только собственное направление в пространстве Rn.

3. Собственные векторы, соответствующие различным, линейно-независимы.

Доказательство этого утверждения проведем «по принципу индукции».

но-независимы.

Докажем, что векторы x1,..., xk, xk +1 также линейно-независимые.

Доказательство проведем «от противного». Предположим, что это не так, т.е. что и не все ai = 0. Умножим равенство (3.3) на матричный оператор Но по предположению вектора x1,..., xk – линейно-независимые и поэтому так как ( i k +1 ) 0 для i = 1, k. Отсюда согласно равенству (3.3) ak +1 = 0 и, следовательно, в (3.3) все ai = 0, i = 1, k + 1. А это означает, что все векторы x1,..., xk, xk +1 – линейно-независимые.

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

5. Собственные числа действительной симметричной матрицы являются действительными числами. Кроме того, если матрица положительно определенная, то ее собственные числа положительные.

Пусть x – собственный вектор матрицы А. Тогда Отсюда, так как (Ах, х) для симметричной матрицы и (х, х) – действительные числа, следует, что также действительное. Если матрица положительно определенная, то в (3.5) ( Ax, x) > 0, ( x, x) > 0, и поэтому > 0.

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

7. Матрицы А, А–1 и B = ci Ai имеют одни и те же собственные векторы, а собственные числа двух последних матриц ( A1 ), ( B) связаны с собственными числами матрицы А ((А)) следующими соотношениями:

8. Преобразование матрицы А вида C = B 1 AB, где матрица В невырожденная ( det( B ) 0 ), называется подобным, а сами матрицы А и С подобными. Подобные матрицы имеют одинаковые характеристические полиномы (так как учтено, что det( B 1 ) = ) и, следовательно, одинаковые собстdet( B) венные числа, а их собственные вектора (соответственно х, у) связаны соотношением Если B = B, то преобразование, выполняемое над А с помоT щью такой матрицы, называется ортогональным подобным преобразованием.

9. Для любой нормы матрицы справедливо неравенство i A. Оно следует из цепочки соотношений: Ax = x В заключение приведем два полезных соотношения, которые на практике могут быть использованы для нахождения двух последних неизвестных собственных значений или для проверки надежности рассчитанной совокупности собственных значений:

Подумайте: будут ли эти соотношения выполняться для действительной матрицы при наличии у нее комплексных собственных значений?

3.2. Решение частичной проблемы собственных чисел Рассмотрим итерационный метод нахождения максимального по модулю собственного числа симметричной матрицы. Все собственные числа такой матрицы действительные, а собственные вектора ei, i = 1, n, образуют в пространстве R n ортонормированный базис.

Замечание. Этот метод применим и для обычной матрицы, если у нее существует полный набор собственных векторов и собственные числа действительные. В этом случае может пригодиться теорема Пиррона: если все компоненты матрицы положительные ( Aij > 0 ), то максимальное собственное число матрицы простое и положительное 1 > 0.

произвольный вектор x (0) = ci ei и будем последовательно умноi = жать его на матрицу А:

Для полученных таким образом векторов можно построить следующие скалярные произведения (учтено, что вектора ei, i = 1, n, ортонормированны):

порядка» (когда 1 < 1 ):

Такая запись справедлива, так как Случай 2 не возможен для положительно определенной матрицы, так как у нее все i > 0.

где a = const.

Видно, что в этом случае при k ( s1 ) 2 s2 и вектор e( k ) поочередно меняет свое «значение». Оба собственных вектора е1 и е можно определить из соотношений (здесь k – четное) следующим образом (так как собственные вектора определяются с точность до константы):

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

для k = 1, K max - максимальное число итераций если (1k ) ) 2 1k ) = 0,0001, то 1 = 1k ) имеет место случай иначе 1 = 1k ) имеет место случай Посмотрим, какую информацию относительно собственных значений матрицы А можно получить, если применить процедуру (3.14) к симметричным матрицам вида собственных значений матрицам будем обозначать следующим образом: i ( A), i ( B ), i (C ).

1. Определение минимального по модулю собственного числа матрицы А. Его можно рассчитать как или, если матрица А положительно определенная ( A > 0 ), то матрица C 0, и поэтому 2. Определение границ отрезка, содержащего все собственные значения матрицы А.

- Если при реализации процедуры (3.14) имеет место случай - Если матрица А положительно определенная, то На рис. 3.1 и 3.2 приведены скриншоты программы, в которой использован алгоритм (3.14). Программа позволяет задать «в ручную» или случайным образом (см. рис. 3.1) произвольную исходную матрицу А, найти для нее максимальное по модулю собственное число, с помощью этих данных сформировать два варианта вспомогательной матрицы В (см. рис. 3.1) и определить для нее максимальное по модулю собственное число. Используя эту информацию, пользователь может рассчитать минимальное по модулю собственное число матрицы А и границы отрезка, содержащего все ее собственные значения. Например, согласно результатам, представленным на рисунках, следует, что все собственные числа матрицы лежат на отрезке В табл. 3.1 приведены результаты небольшого исследования, выполненного с помощью указанной программы. Строгая регулярность в представленных результатах отсутствует.

Таблицы 3.1. Зависимость числа итераций от размерности случайной матрицы для eps = 0, 000001 при величине выборки 1. Для чего производят нормировку текущего вектора e1( k ) = при реализации итерационного метода нахождения максимального собственного значения?

2. Что такое характеристический полином матрицы?

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

4. Как с помощью рассмотренного итерационного метода найти минимальное по модулю собственное значение матрицы?

5. Как связаны между собой собственные значения и собственные векk торы матриц А и B = Ai ?

6. Составьте алгоритмы для вычисления компонент матриц 7. Какие матрицы называются подобными?

8. Чем можно объяснить отсутствие строгой регулярности в результатах, представленных в табл. 3.1?

3.3. Решение задачи на собственные значения Сущность метода Данилевского состоит в приведении матрицы А с помощью (п – 1) подобных преобразований к подобной ей матрице Фробениуса Р:

Для последней характеристический полином записывается просто путем разложения по элементам первой строки Подобное преобразование матрицы А представляет собой последовательный перевод ее строк, начиная с п-й, в строки матрицы Фробениуса.

Детально рассмотрим k-е ( k = 1, 2,...,(n 1) ) преобразование, формирующее ( n (k 1) )-ю строку матрицы Р.

После (k 1) -го преобразования матрица имеет вид ( A(0) = A ) Над ней выполняются следующие операции.

1. «Сажаем» единицу на место компонента Ank k 1) n k путем деления на него всего (п – k )-го столбца^ 2. Обнуляем все остальные компоненты приводимой строки. Для этого вычитаем из каждого j-го столбца (п – k)-й столбец, умноженный на соответствующий компонент An k k 1) j (п – (k –1)-й строки Эти операции равносильны матричной операции где матрица M ( n k ) отличается от единичной только (п – k)-й строкой, компоненты которой равны Чтобы завершить k-е подобное преобразование, необходимо умножить промежуточную матрицу A( k ) на матрицу ( M ( n k ) ) 1 :

Последняя получается из единичной матрицы подстановкой в ее (п – k)-ю строку компонент (п – (k – 1))-й строки матрицы A( k 1).

Легко проверить, что в этом случае будет выполняться необходимое соотношение Матричная операция (3.20) приводит к изменению компонент только одной (п – k)-й строки A( k ) или с учетом имеющихся нулевых компонент в нижних подстолбцах матрицы A( k ) После построения характеристического полинома и нахождения его корней – собственных чисел матрицы А – легко определяются собственные векторы матрицы Р(y), а затем и матрицы А, так как Для определения векторов у распишем систему ( P E ) y = 0 :

Положив уп = 1 (это можно сделать, так как собственные векторы определяются с точностью до постоянной), из последнего и последующих (до первого) уравнений легко находим При таких значениях уi первое уравнение системы (3.25) удовлетворяется автоматически, поскольку оно совпадает с характеристическим уравнением (3.2).

Теперь с учетом вида матриц легко рассчитываются компоненты собственных векторов х (первоначально устанавливается х = у):

Остановимся на особенностях программной реализации данного метода.

1. Обрабатываемая матрица хранится на месте исходной.

2. Для повышения устойчивости алгоритма перед выполнением очередного k-го подобного преобразования необходимо выполнять процедуру выбора «ведущего» элемента: найти среди компонент An k k 1) j, j = 1,(n k ), максимальный по модулю элемент и, если jmax (n k ), то переставить местами верхние части столбцов и строки с номерами (п – k) и jmax. Такая перестановка ( TA( k 1)T ) будет подобной, поскольку для матрицы перестановок Т справедливо соотношение T 1 = T. Если окажется, что «ведущий» элемент равен нулю, то имеет место сингулярный случай, когда характеристический полином распадается на произведение двух характеристических полиномов:

причем для одного из них матрица Фробениуса P ( k ) уже получена, а для подматрицы A( n k ) ее предстоит еще найти с помощью того же метода. В этом особом случае собственные векторы матрицы следует искать вне данного метода (см. п. 3.5).

3. Строку M n n kj), j = 1, n, (см. формулы (3.20)) матрицы M ( n k ) можно первоначально формировать в дополнительном векторе т, а в конце k-го преобразования сохранять на месте (п – (k – 1))-й обработанной строки исходной матрицы. При этом формулы (3.18), (3.19) перепишутся следующим образом Здесь исключена явная обработка (п – (k – 1))-й приводимой строки, поскольку ее компоненты используются при завершении kго подобного преобразования (см. формулы (3.23)).

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

Скриншоты этой программы представлены на рис. 3.3–3.7. Так, на рис. 3.3 дана форма с описанием порядка работы с программой, на рис. 3.4 – главная форма программы, на которой приведен пример матрицы, для которой матрица Фробениуса распадается на две подматрицы (приведен этап получения первой подматрицы).

На рис. 3.5 показаны корни характеристического полинома, полученные методом парабол для первой подматрицы Фробениуса.

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

Обратите внимание, что, как и положено, сумма следов двух матриц Фробениуса равна следу исходной матрицы (21+7=28), а суммы корней каждой из подматриц совпадает с ее следом.

Ниже приведены листинги на Visual Basic алгоритма, вычисляющего матрицу Фробениуса, и прицедуры выбора ведущего элемента, используемой на каждом шаге алгоритма.

If v_element(k) = 1 Then ‘выделена подматрица Фробениуса ma(j) = a(n – k, j) 'сохраняем обрабатываемую строку 'ставим 1 в обрабатываемой строке 'ставим 0 вместо остальных элементов обрабатываемой строки 'умножаем на обратную матрицу (М) Next k Function v_element (k As Integer) As Integer Dim amax#, jmax%, i%, j%, buf#, cd% amax = Abs(a(n – k, n – 1 – k)) buf = Abs(a(n – k, j)) If Not (k = jmax) Then If amax < 1E-30 Then cd = v_element = cd End Function Здесь ma(n) – вспомогательный массив для временного хранения обрабатываемой строки; n – размерность обрабатываемо матрицы; v_element – функция, реализующая процедуру выбора ведущего элемента.

В данном методе нахождение коэффициентов характеристического полинома базируется на использовании матричного тождества Гамильтона – Кели где рi – коэффициенты характеристического полинома Возьмем произвольный вектор y (0) 0 и умножим на него (3.30). В результате получим векторное уравнение относительно неизвестных рi:

или в новых обозначениях y ( k ) = Ay ( k 1) = Ak y (0) Эту систему линейных уравнений можно переписать в привычной форме где Приведем алгоритм для расчета матрицы Y и вектора y(п) (не забудте, что правый столбец матрицы Yi n 1 yi(0) задается пользователем):

Система (3.32) решается любым известным методом, например, методом исключения Гаусса. Если матрицы Y окажется особенной (это имеет место, если векторы y (0), y (1),..., y ( n 1) линейно-зависимые), то следует сменить начальный вектор y (0).

Чтобы не произошло переполнения в процессе вычисления векторов y ( k ), исходную матрицу целесообразно предварительно отAi j нормировать: Ai j =. После этого собственные значения отнормированной и исходной матриц будут связаны соотношением После определения собственных значений матрицы А ее собственные векторы можно искать в виде разложения по полной линейно-независимой системе векторов y (0), y (1),..., y ( n 1) :

Для нахождения неизвестных коэффициентов i подставим (3.33) в уравнение Ax = x. Тогда с учетом процедуры построения векторов y ( k ), получим Заменим в (3.34) y ( n ) выражением (3.31) и соберем подобные слагаемые Поскольку векторы y (0), y (1),..., y ( n 1) – линейно-независимые, то коэффициенты в этом равенстве при всех y ( k ) должны быть равны нулю. Отсюда получаем Здесь в последнем соотношении в скобках стоит характеристический полином, который равен нулю. Поэтому, не нарушая это равенство, можно положить 1 = 1. Тогда остальные неизвестные коэффициентов можно рассчитать с помощью рекуррентной формулы Для кратных недостающие собственные векторы (если они существуют) следует искать, изменяя вектор y (0). В этом случае необходимо проверить, что полученные собственные векторы x ( i ) являются линейно-независимыми. Для этой цели можно использовать процедуру ортогонализации, которая проходит до конца только для линейно-независимых векторов.

На рис. 3.8 и 3.9 приведены скриншоты программы, написанной на Visual Basic, которая ипользует метод Крылова для построения характеристического полинома для производьной матрицы, введенной пользователем или случайным образом сформированной самой программой. На рис. 3.8 показана главная экранная форма программы для нахождения коэффициентов характеристического полинома произвольной матрицы методом Крылова, а на рис. 3.9 – вспомогательная, которая вызывается пунктом меню «метод парабол» на главной форме (см. рис. 3.8).

После ввода матрицы можно построить «систему Крылова» (см.

формулу (3.32)), а затем найти коэффициенты характеристического полинома исходной матрицы. Система Крылова решается методом оптимального исключения Гаусса–Жордана (см. п. 2.2.2).

Далее методом парабол (см. п. 4.3) можно определить все корни построенного полинома (собственные значения исходной матрицы). Вместе с корнями программа выдает значения следа матрицы ( Aii ) и сумму этих корней. Сравненительный анализ этих велиi = чин позволяет оценить корректность выполнения всех промежуточных этапов решения задачи нахождения собственных значений исходной матрицы.

Дополнительно на рис 3.10 представлены листинг и результаты аналогичного поиска собственных значений матрицы методом Крылова в рамках пакета wxMaxima. Видно, что на всех этапах решения задачи данные представленные на рис. 3.8–3.10 совпадают.

При этом трудоемкость получения этих результатов для пакета существенно меньше.

Листинг программы для пакета wxMaxima:

/* формируем систему Крылова для матрицы A и начального вектора y0 */ formKrilov(A,y0):=block [k,i,m,n], /* локальные переменные */ n:length(A), Y:zeromatrix(n,n-1), /* генерация нулевой матрицы размерности n*(n-1) */ Y:addcol(Y,y0), /* добавляем к матрице Y справа столбец y0 */ for m thru n do Y[i,n+1-k]:Y[i,n+1-k]+A[i,m]*Y[m,n+2-k], yn:zeromatrix(n,1), for m thru n do yn[i]:yn[i]-A[i,m]*Y[m,1] /* главная программа */ numer:true;

fpprintprec:6;

A:matrix([1,2,3],[4,5,6],[7,8,9]);

y0:matrix([1],[0],[1]);

formKrilov(A,y0);

print("Матрица и правая часть системы Крылова");

Y_obr:invert(Y); /* находим обратную матрицу системы Крылова */ n:length(A);

p:zeromatrix(n,1);

/* решаем систему Крылова через обратную матрицу p=Y_obr*yn */ for i thru n do for m thru n do p[i]:p[i]+Y_obr[i,m]*yn[m];

c:[1]; for i thru n do c:append(c,p[i]);

print("Список коэффициентовхарактеристического полинома = ",c);;

solve(x^3+p[1]*x^2+p[2]*x+p[3],[x]); /* определение корней полинома */;

Результаты выполнения программы даны на рис. 3.10.

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

для случайных матриц, компонентами которых являлись псевдослучайные равномерно распределенные числа из интервала [-0,5; 0,5].

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

4. После нахождения уточненного значения корня, выделяем его из текущего полинома и тем самым понижаем его порядок. Коэффициенты нового полинома, полученного после выделения действительного корня из предыдущего полинома, приведены в нижней таблице рис. 4.5. Используя кнопку «обновить исходный полином», коэффициенты нового полинома ставим в верхнюю таблицу на место их предыдущих значений (рис. 4.6).

5. Новый полином имеет второй порядок, поэтому его корни определяются по известным формулам – кнопка «найти корни квадратичного полинома». Экранная форма со значениями этих корней приведена на рис. 4.7.

6. С помощью кнопки «найти значение полинома» – можно рассчитать значение полинома для найденного корня и убедится в правильности полученных его действительных и мнимых значений.

7. Пункт меню «графическая локализация», расположенный на главной форме, загружает форму программы для локализации корней на комплексной плоскости и последующего уточнения локализованных корней. На рис. 4.8 показан результат графической локализации методом перебора всех корней исходного полинома. Затем методом наискорейшего спуска произведено уточнение действительного корня (его точное значение равно 3). При этом использован алгоритм одномерного спуска с отступом назад.

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

Это позволяет произвести графическую локализацию местоположения корней полинома на комплексной плоскости. Затем значение каждого локализованного корня можно уточнить, минимизируя выпуклый квадратичный функционал методами наискорейшего или покоординатного спуска. Согласно первому методу спуск к нулевому значению функционала осуществляется против направления, заданного ортами вектора градиента функционала в текущей начальной точке комплексной плоскости.

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

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

Приведем соответствующие алгоритмы.

1. Алгоритмы метода наискорейшего спуска для двух стратегий одномерного спуска.

Алгоритм it = 0 число использованных итераций Алгоритм it = 0 число использованных итераций do while h > 0, 2. Алгоритм метода покоординатного спуска:

it = 0 число использованных итераций do while it < it max ' определение направления спуска для данной координаты ' одномерный спуск вдоль данной координаты Здесь it max - максимальное число допустимых итераций (спусков), h – параметр спуска, esp – требуемая точность.

На рис. 4.9 приведена экранная форма с описанием порядка работы с программой.

Далее показано, как в рамках пакета wxMaxima решаются некоторые из ренее рассмотренных задач.

/* строит отображение, которое дает полином для заданной окружности */ /* здесь: с-список коэффициентов полинома,x0,y0,r-центр окружности и ее радиус */ /* n-число исследуемых точек на окружности */ otobragenieP(c,x0,y0,r,n):=block dfi:2*%pi/n, /* шаг по углу */ lx:[],ly:[], /* инициализируем два списка */ xt:x0+r*cos(dfi*i),yt:y0+r*sin(dfi*i), /* функция введена ранее – вычисляет значение полинома в заданной точке */ p(c,xt,yt), lx:append(lx,[p1]),ly:append(ly,[p2]) /* добавляем элементы в списки /* использованием 2-х списков строим в отдельном окне график отображения*/ plot2d([discrete,lx,ly]) numer:true;

fpprintprec:5; /* в выводимом числе 5 цифр, не считая точки */ c:[1,-5,8,6];

/* локализация корней полинома на плоскости */ n:length(c);

c1:makelist(abs(c[i]),i,2,n); /* формирование 2-х вспомогательных списков */ c2:makelist(abs(c[i]),i,1,n-1);

/* оценка местоположения корней полинома на комплексной плоскости */ print(" Mаксимальный R = ",1+lmax(c1)/abs(c[1]));

print(" Mинимальный r = ",1/(1+lmax(c2)/abs(c[n])));

r:9; x0:0; y0:0;

n:960; otobragenieP(c,x0,y0,r,n); /* строим отображение */;

pol(x):=x^3-5*x^2+8*x-6;

/* проверяем – есть у полинома действительные корни на отрезке [строим график функции в основном окне */ wxplot2d(pol(t),[t,-9,9],[color,red]);

wxplot2d(pol(t),[t,2.5,3.5],[color,red]); /* уточняем положение действ.корня */ /* внутренняя функция находит корни полинома */ to_poly_solve([x^3-5*x^2+8*x-5=0], [x]);

Результаты выполнения команд и подпрограмм листинга показаны на рис. 4.10.

На рис. 4.11 дано отображение полиномом (x^3-5*x^2+8*x-6) на комплексной плоскости окружности радиуса 9 (построено пакетом wxMaxima). Проверка наличия у полинома дейсивительного корня на отрезке [-9,9] (построено пакетом wxMaxima) показана на рис.

4.12, а уточнение положения дейсивительного корня на отрезке [2.5,3.5] – на рис. 4.13.

Рис. 4. Рис. 4. Метод парабол является наиболее надежным алгоритмом поиска всех корней полинома. Метод не требует задания приближенных значений корней и определяет их последовательно один за другим.

Сходимость метода теоретически не доказана, но не известно ни одного случая, когда бы метод не сходился.

Для запуска метода произвольно выбираются три начальные точки на комплексной плоскости z = x + iy z ( x, y ), обычно z0 = (1,0), z1 = (1,0), z2 = (0,0). По значениям полинома в этих точках Pn ( z0 ), Pn ( z1 ), Pn ( z2 ) (для их расчета используется схема (4.3)) строится интерполяционный многочлен Лагранжа второго порядка:

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

Процесс повторяется до тех пор, пока значение полинома для новой точки z2 (или расстояние между старой и новой точками z2) не будет меньше заданной малой величины P ( z2 ) < eps.

После того, как корень найден, порядок полинома понижается с использованием схем (4.6) или (4.9) соответственно для действительного корня (мнимая часть корня равна нулю) и для случая комплексно сопряженных корней (мнимая часть корня отлична от нуля). Затем поиск корня продолжается уже для полинома меньшего порядка.

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

Замечание. Все переменные (кроме коэффициентов полинома) в алгоритме должны быть описаны как комплексные. Соответственно, и все операции над ними должны учитывать их природу.

На рис. 4.14 приведены корни полинома, вычисленные методом парабол в рамках программы, представленной в п. 4.2. При реализации алгоритма метода парабол для комплексных чисел был определен новый тип переменной – комплексное число, состоящее из действительной и комплексной частей. Все операции с комплексными числами, используемые в алгоритме, такие, как были оформлены в виде соответствующих подпрограмм. Листинги на Visual Basic этих подпрограмм и прицедуры, реализующей метод парабол, следующие:

‘складывает два комплексных числа Sub addc (a As complex, b As complex, c As complex) End Sub ‘вычитание для двух комплексных чисел Sub subc (a As complex, b As complex, c As complex) a.r = b.r – c.r a.c = b.c – c.c End Sub ‘умножение двух комплексных числа Sub mulc (aa As complex, b As complex, c As complex) aa.r = b.r * c.r – b.c * c.c aa.c = b.r * c.c + b.c * c.r End Sub ‘деление для двух комплексных чисел Sub divc (a As complex, b As complex, c As complex) Dim d# d = c.r * c.r + c.c * c.c If d = 0# Then Exit Sub End If a.r = (b.r * c.r + b.c * c.c) / d a.c = (c.r * b.c – b.r * c.c) / d End Sub ‘корень квадратный из комплексного числа Sub sqrtc (a As complex, b As complex) Dim r# r = Sqr(b.r * b.r + b.c * b.c) a.r = Sqr((r + b.r) / 2) a.c = Sqr((r – b.r) / 2) If b.c < 0 Then a.c = -a.c End If End Sub Sub muller () ReDim zn(3) ReDim l(3) Dim dz#, dz1#, dz2#, normp# Dim it%, txt$ Dim z0 As complex Dim z1 As complex Dim z2 As complex Dim z3 As complex Dim z4 As complex Dim z5 As complex Dim aa As complex Dim bb As complex Dim cc As complex Dim dd As complex Dim l0 As complex Dim l1 As complex Dim l2 As complex nr = 0 'число найденных корней Do While n > ‘наличие нулевого корня If Abs(a(n)) <.0001 Then root(nr).r = root(nr).c = zn(0).r = - zn(0).c = zn(1).r = zn(1).c = For i = 0 To 2 'значения полинома в трех точках Call polinomc(zn(i)) l(i).r = p.r l(i).c = p.c Next i normp = Sqr(l(2).r * l(2).r + l(2).c * l(2).c) Do While normp > eps Call subc(z0, zn(0), zn(2)) 'zn(0)-zn(2) -> z Call subc(z1, zn(0), zn(1)) 'zn(0)-zn(1) -> z Call mulc(z3, z0, z1) 'z0*z1 -> z Call divc(l0, l(0), z3) 'l(0)/z3 -> l0 L Call subc(z0, zn(1), zn(0)) 'zn(1)-zn(0) -> z Call subc(z1, zn(1), zn(2)) 'zn(1)-zn(2) -> z Call mulc(z3, z0, z1) 'z0*z1 -> z Call divc(l1, l(1), z3) 'l(1)/z3 -> l1 L Call subc(z0, zn(2), zn(1)) 'zn(2)-zn(1) -> z Call subc(z1, zn(2), zn(0)) 'zn(2)-zn(0) -> z Call mulc(z3, z0, z1) 'z0*z1 -> z Call divc(l2, l(2), z3) 'l(2)/z3 -> l2 L Call addc(z0, l0, l1) Call addc(aa, z0, l2) 'A aa=l0+l1+l Call addc(z3, zn(2), zn(1)) 'zn(2)+zn(1) -> z Call mulc(z0, l0, z3) 'l0*z3 -> z0 L0*(z(2)+z(1)) Call addc(z3, zn(0), zn(2)) 'zn(0)+zn(2) -> z Call mulc(z1, l1, z3) 'l1*z3 -> z1 L1*(z(0)+z(2)) Call addc(z3, zn(1), zn(0)) 'zn(1)+zn(0) -> z Call mulc(z2, l2, z3) 'l2*z3 -> z2 L2*(z(1)+z(0)) Call addc(z3, z0, z1) Call addc(bb, z3, z2) 'B bb=z0+z1+z Call mulc(z3, zn(2), zn(1)) 'zn(2)*zn(1) -> z Call mulc(z0, l0, z3) 'L0*z3 -> z Call mulc(z3, zn(2), zn(0)) 'zn(2)*zn(0) -> z Call mulc(z1, l1, z3) 'L1*z3 -> z Call mulc(z3, zn(0), zn(1)) 'zn(0)*zn(1) -> z Call mulc(z2, l2, z3) 'L2*z3 -> z Call addc(z3, z0, z1) Call addc(cc, z3, z2) 'C cc=z0+z1+z Call mulc(z4, aa, cc) 'A*C -> z z4.r = 4 * z4.r Call subc(dd, z3, z4) 'D dd=z3-z Call sqrtc(z3, dd) 'sqr(dd) -> z Call addc(z4, bb, z3) Call subc(z5, bb, z3) Call divc(z0, z4, aa) 'z0 два корня Call divc(z1, z5, aa) 'z 'находим ближайший корень Call subc(z3, zn(2), z0) 'zn(2)-z0 -> z Call subc(z4, zn(2), z1) 'zn(2)-z1 -> z dz1 = Sqr(z3.r * z3.r + z3.c * z3.c) dz2 = Sqr(z4.r * z4.r + z4.c * z4.c) MsgBox "При поиске корня исчерпаны все итерации !", 48, "Сообщение" Call polinomc(zn(2)) l(2).r = p.r 'значение полинома в третей точке normp = Sqr(l(2).r * l(2).r + l(2).c * l(2).c) root(nr) = zn(2) If Abs(zn(2).c) <.0001 Then Call extr_root1(zn(2).r) Call extr_root2(zn(2)) root(nr) = zn(2) 'сохраняем второй комплексно-сопряженный корень root(nr).c = -root(nr).c Lab:

Loop 'корни линейного или квадратичного полинома If n = 1 Then root(nr).r = -a(1) / a(0) root(nr).c = 0# Else Call fine_root root(nr).r = rr root(nr).c = rc root(nr).r = rr root(nr).c = rc End If End Sub ‘вычисление значения полинома в точке х Sub polinomc (x As complex) Dim i%, r1 As Double p.r = a(0) p.c = p.r = x.r * p.r – x.c * p.c + a(i) p.c = x.r * p.c + x.c * r Next i End Sub Здесь используются переменными типа – comlex, введенного следующим оператором:

Type complex End Type На рис. 4.14 даны корни полинома x 3 5 x 2 + 8 x 6, определенные методом парабол. Форма вызывается пунктом меню «метод парабол» на главной форме (см. рис. 4.2).

Заметим, что в рамках пакета wxMixima результаты, отображенные на рис. 4.15, получются с помощью одной внутренней команды (здесь %i – обозначение мнимой единицы).

1. Почему при вычислении значения полинома схема Горнера является предпочтительной по сравнению с другими способами?

2. В каких случаях используются схемы деления полинома на линейный и квадратичный множитель?

3. Каким образом можно грубо оценить местоположение корней полинома на комплексной плоскости?

4. Опишите, каким образом можно с помощью приведенной выше программы последовательно найти все корни полинома.

5. В чем отличие методов наискорейшего и покоординатного спусков?

6. Составьте алгоритм, реализующий метод парабол.

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

Привем список этих программ с кратким описанием их возможностей.

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

2. ITER_MET – в приложении реализованы три итерационных метода, с помощью которых можно найти решения произвольной системы уравнений. Исходную систему можно автоматически преобразовать в равносильную систему с положително определенной симметричной матрицей. В приложении исходную систему можно ввести вручную или сформировать случайным образом с помощью датчика равномерно распределенных псевдослучайных чисел.

3. CONVERT – приложение позволяет находить для произвольной матрицы обратную и значение определителя.

4. ORTGAUSS – приложение позволяет найти решение произвольной системы уравнений методом оптимального исключения или первым методом ортогонализации.

5. KV_ROOT – приложение позволяет найти решение произвольной системы уравнений с симметричной положительно определенной матрицей методом квадратного корня или вторым методом ортогонализации. В приложении можно преобразовать обычную систему в равносильную с симметричной положительно определенной матрицей.

6. REGUL – приложение строит решение для произвольной системы и ее регуляризованного прототипа. Может использоваться для поиска нормального решения для плохо обусловленных систем.