Размещено на http://www.allbest.ru/

Размещено на http://www.allbest.ru/

Содержание

Введение

1. Постановка задачи

2. Программная реализация методов на языке программирования Си

3. Код программы

4. Проверка работы программы

Заключение

Список использованных источников



Введение

Большое количество задач с механики, физики и техники требует нахождение собственных значений и собственных векторов матриц, т.е. таких значений л, для которых существует нетривиальное решение однородной системы линейных алгебраических уравнений. Тут А-действительная квадратичная матрица порядка n с элементами ajk, а --вектор с компонентами x1, x2,…, xn Каждому собственному значению лi соответствует хотя бы одно нетривиальное решение. Если даже матрица А действительная, ей собственные числа (все или некоторые) и собственные векторы могут быть недействительными. Собственные числа являются корнями уравнения, где Е - единичная матрица порядка n.

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

ЗАДАНИЕ НА КУРСОВУЮ РАБОТУ

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

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



1. Постановка задачи

Большое количество задач науки и техники требует отыскания собственных чисел (СЧ) и собственных векторов (СВ). Широкий выход на практику проблемы СЧ и СВ породил много синонимов для понятия этой проблемы.

Введем основные обозначения и названия:

E - единичная матрица размерности [n Ч n],

A - искомая квадратная матрица размерности [n Ч n],

x - n-мерный столбец,

л - скаляр, используется как переменная характеристического уравнения и как собственное число,

A ? лE - характеристическая матрица,

|A ? лE| - характеристический определитель или вековой определитель. Рассматриваемый как многочлен D(л) характеристический определитель называется характеристическим многочленом.

D(л) = 0 - характеристическое уравнение или вековое уравнение.

Ненулевой столбец x, удовлетворяющий уравнению Ax = лx называется собственным вектором (столбцом), а число л - собственным числом (значением), соответствующим СВ x. Иногда соответствие СЧ и СВ выражают словами "СВ принадлежит СЧ" или "СВ для СЧ". Методы вычисления собственных значений делятся на две принципиально различные группы. В одной находятся методы развертывания вековых определителей:

методы Данилевского, Крылова, Леверрье и др., с последующим применением какого-либо итерационного метода решения алгебраического уравнения. В другой группе - итеративные методы без развертывания вековых определителей: методы Якоби, скалярных произведений и др. Задача нахождения СВ при найденных СЧ значительно проще, чем нахождение СЧ. Ядро оператора A ? лE есть множество СВ для СЧ. Ядро может быть найдено, например, идеальным методом Гаусса. Однако почти в каждом итеративном методе вычисления СЧ есть своего рода льготная добавка для вычисления СВ без решения системы (A ? лE) x = 0.

Метод Данилевского

Этот метод самый эффективный по числе арифметических операций, используемых для развертывания векового определителя. Так, для определителей 5-го порядка по методу Данилевского требуется более чем в три раза меньшее число арифметических действий, по сравнению с непосредственным развертыванием, а для определителей 9-го порядка -- в тысячу раз меньшее.

Метод заключается в приведении подобным преобразованием матрицы к нормальному виду Фробениуса

Вековой определитель для матрицы F имеет вид

После разложения D(л) по первой строке получаем

Таким образом, знание матрицы Фробениуса является знанием коэффициентов характеристического уравнения |F ? лE| = 0. Оказывается, коэффициенты характеристического многочлена |A ? лE| совпадают с коэффициентами многочлена D(л). Действительно, в силу подобия A = S?1P S, и далее

Нахождение формы Фробениуса

Начнем с последней строки.

Основное преобразование выполняется, когда . Тогда поделим в A предпоследний столбец на , назовем его ведущим и, умножая на , вычтем его из i-го столбца, i = 1,..., n?2, n. Результатом будет матрица A, последняя строка которой будет совпадать с последней строкой формы Фробениуса.

Покажем, что этому вычислительному алгоритму соответствует некоторое преобразование подобия. Положим матрицу такой, чтобы выполнялось равенство . Очевидно строки , кроме (n ? 1)-й, совпадают со строками единичной матрицы, предпоследняя такова:

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

.

Далее, если , то над матрицей A1 произведем аналогичные операции, взяв за основу (n?1)-ю строку. В результате получим, в которой две последние строки будут совпадать с формой Фробениуса. Продолжая процесс, получим форму Фробениуса матрицы A:

(1)

если, конечно, для всех было верно .

Остановимся теперь на исключительных случаях: равенство нулю элемента для какого-либо i. Возможны два типа исключительных случаев.

A. Пусть для всех . Этот случай облегчает дальнейшие вычисления. Действительно, положим k:= n ? i ? 1, и

Тогда

Матрица F1 - матрица Фробениуса. Остается применить метод Данилевского к матрице B.

B. Пусть и . Тогда, переставляя в матрице местами j-ю и k-ю строки, а затем j-й и k-й столбцы, получим матрицу , годную для основного преобразования метода Данилевского, и как нетрудно показать, подобную матрице . Действительно, = . Переставим строки j, k в обеих частях равенства. На месте левой матрицы E появится , получаемая перестановкой j-й и k-й строки. Аналогично, переставляя столбцы, получим матрицу , в которой, в отличие от E переставлены j-й и k-й столбцы. Легко заметить, что = . Следовательно, = . Непосредственно умножая, легко убедиться, что = E, откуда = и подобна .

Вычисление собственных столбцов

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

Найдем сначала y - собственный столбец матрицы P:

(P ? лE)y = 0 (2)

с точностью до постоянного множителя . Из формул (1) и (2) получим расчетную формулу для x - собственного столбца матрицы A:

(3)

где , либо единичная матрица, либо получена из единичной подстановкой i-го и k-го столбцов i < k. При вычислениях на ЭВМ или вручную целесообразно учесть особую структуру матриц преобразований и вычисления производить последовательно со столбцами

,

не пользуясь общими правилами умножения матрицы на столбец.

Квазиединичное преобразование меняет местами две координаты столбца . Каждое преобразование меняет в столбце только i-ю координату: .

Этот алгоритм естественно наводит на мысль об экономии памяти. Вместо запоминания матриц можно запомнить массив размерности (n ? 1) Ч n из строк . Вместо запоминания матриц можно запомнить целочисленный массив длины n ? 1.

Если перед нами стоит задача вычисления собственных столбцов, и мы желаем это сделать, используя информацию, полученную при вычислении собственных чисел, то в первом исключительном случае надо продолжить работу со всей матрицей A, а не только с B. Алгоритм вычисления собственных столбцов немного усложнится. Одному СЧ может принадлежать несколько СВ. Каков алгоритм в общем случае вполне ясно из следующего примера.

Рассмотрим вариант с одним исключительным случаем n = 6, i = 2. Тогда получим на последнем шаге (четвертом, так как один пропущен):



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

Первое строка-уравнение относительно c таково:

Будучи совместным, оно либо определяет постоянную c (множитель при c отличен от нуля) и, соответственно, один вектор y, либо не определяет. Во втором случае получаем два собственных вектора матрицы :

Несовместность (4) (множитель при c равен нулю, свободный член в (4) не нуль) дает один собственный вектор v1.

Формула (3) для рассматриваемого примера имеет вид программа матрица данилевский

Погрешность метода

При точном выполнении арифметических операций (умножении, делении, сложении,вычитании) все погрешности возникают из приближенного вычисления корней векового уравнения. Анализ этих ошибок совпадает с анализом ошибок метода, использованного для решения векового уравнения. При неточном выполнении арифметических операций сам метод может давать очень грубые ошибки в коэффициентах векового уравнения. Таков случай паталогической близости к нулю элемента (” ? ” - знак приближенности). Когда нельзя быть уверенным в знаке и в том, что он отличен от нуля, столбец n ? i ? 1 нельзя брать ведущим. Ведущим желательно сделать такой j-й столбец, что, используя методику второго исключительного случая.

2. Программная реализация методов на языке программирования Си

Головная программа

Сначала задается количество переменных и количество ограничений.

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

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

3. Код программы

#include <stdio.h>

#include <stdbool.h>

#include <stdlib.h>

// вспомогательные типы

// тип-обёртка для массива, чтобы не таскать размер отдельным параметром

Typedef

struct array

{

int size;

float* data;

} array;

// массив массивов

// capacity - максимальный размер (указывается при создании),

// size - текущий

typedef

struct vector

{

int size;

int capacity;

array* data;

} vector;

// все объявления функций

// создаём и уничтожаем массив

array newArray(int n);

void deleteArray(array* a);

// создаём вектор вместимости cap

vector newVector(int cap);

//.. и удаляем

void deleteVector(vector* v);

// добавляем в вектор новый массив

void push_back(vector* v, array a);

// прототип функции для метода Данилевского

void danilevsky(int n, float a[n][n], float q[n + 1]);

// функции для работы с массивами и векторами

array newArray(int n)

{

array a;

a.size = n;

a.data = malloc(n * sizeof(float));

return a;

}

void deleteArray(array* a)

{

if (a) {

free(a->data);

}

}

vector newVector(int cap)

{

vector v;

v.size = 0;

v.capacity = cap;

v.data = malloc(cap * sizeof(array));

return v;

}

void deleteVector(vector* v)

{

int i;

if (v && v->data) {

for (i = 0; i < v->size; ++i) {

deleteArray(&v->data[i]);

}

free(v->data);

}

}

void push_back(vector* v, array a)

{

if (v && v->size < v->capacity) {

v->data[v->size++] = a;

}

}

// реализация метода Данилевского

// min и max есть в stdlib.h

int max(int a, int b)

{

return (a > b)? a: b;

}

int min(int a, int b)

{

return (a < b)? a: b;

}

// переставляем r и c строки и столбцы в матрице a

void swapRowsAndColumns(int n, float a[n][n], int r, int c)

{

int k;

for (k = 0; k < n; ++k) {

float tmp = a[r][k];

a[r][k] = a[c][k];

a[c][k] = tmp;

}

for (k = 0; k < n; ++k) {

float tmp = a[k][r];

a[k][r] = a[k][c];

a[k][c] = tmp;

}

}

// собственно реализация метода Данилевского

// если будет вырожденный случай, мы изменим n,

// чтобы можно было запустить функцию ещё раз

// если вырожденного случая нет, n = 0 на выходе

// N - настоящий размер матрицы, n - текущий

array DanilevskyImpl(int* pn, int N, float a[N][N])

{

int z,j,i,k, n = *pn;

if (n == 1) {

array p = newArray(1);

p.data[0] = -a[0][0];

*pn = 0;

return p;

}

array row = newArray(n);

for (z = 1; z < n; ++z) {

for (j = 0; j < n; ++j) {

row.data[j] = a[n - z][j];

}

float t = a[n - z][n - z - 1];

if (t == 0) {

int j0 = 0;

bool all_is_zero = true;

for (; j0 < n - z - 1; ++j0) {

if (a[n - z][j0]!= 0) {

all_is_zero = false;

break;

}

}

if (all_is_zero) { // первый вариант вырожденного случая

// к оставшейся подматрице мы потом применим метод снова

n -= z;

*pn -= z;

array p = newArray(z);

for (j = 0; j < z; ++j) {

p.data[j] = -a[n][j + n];

}

deleteArray(&row);

return p;

} else { // второй вариант

// переставляем j0 и n - z - 1 строки...

swapRowsAndColumns(n, a, j0, n - z - 1);

// и нам нужно сделать итерацию ещё раз

--z;

}

} else { // вырожденного случая нет, следуем основному алгоритму

// делим весь столбец

for (i = 0; i < n; ++i) {

a[i][n - z - 1] /= t;

}

// вычитаем из всех столбцов с номером j ведущий, умноженный на a[n - z][j]

for (j = 0; j < n - z - 1; ++j) {

t = a[n - z][j];

for (i = 0; i < n; ++i) {

a[i][j] -= a[i][n - z - 1] * t;

}

}

// кроме самого столбца с номером n - z - 1

for (j = n - z; j < n; ++j) {

t = a[n - z][j];

for (i = 0; i < n; ++i) {

a[i][j] -= a[i][n - z - 1] * t;

}

}

// домножаем слева на M^(-1), как в книжке

for (j = 0; j < n; ++j) {

t = 0;

for (k = 0; k < n; ++k) {

t += row.data[k] * a[k][j];

}

a[n - z - 1][j] = t;

}

}

}

// вектор коэффициентов

array p = newArray(n);

for (j = 0; j < n; ++j) {

// мы возвращаем сами коэффициенты, а в матрице Фробениуса они с минусом

p.data[j] = -a[0][j];

}

// если случай невырожденный, мы обработали всю матрицу

*pn = 0;

deleteArray(&row);

return p;

}

// перемножаем многочлены с коэффициентами p1 и p2

array mult(array p1, array p2)

{

int j,i;

int n1 = p1.size, n2 = p2.size;

int n = n1 + n2; // степень итогового многочлена

array res = newArray(n);

for (i = 0; i < n; ++i) {

// коэффициент при x^i

float t = 0;

for (j = max(0, i - n2 + 1); j <= min(i, n1 - 1); ++j) {

t += p1.data[n1 - 1 - j] * p2.data[n2 - 1 - (i - j)];

}

if (n - 1 - i < n2) {

t += p2.data[n - 1 - i];

}

if (n - 1 - i < n1) {

t += p1.data[n - 1 - i];

}

res.data[n - 1 - i] = t;

} return res;

}

// из векторов коэффициентов многочлена составим коэффициенты характеристического

array collect(vector v)

{

int i,k;

int n = v.size;

array coeff = newArray(v.data[0].size);

array res;

// копируем v.data[0] в coeff

for (i = 0; i < v.data[0].size; ++i) {

coeff.data[i] = v.data[0].data[i];

}

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

for (k = 1; k < n; ++k) {

res = mult(coeff, v.data[k]);

deleteArray(&coeff);

coeff = res;

}

return coeff;

}

// метод Данилевского

void danilevsky(int n, float a[n][n], float q[n + 1])

{

int i;

vector v = newVector(n);

int m = n;

while (m) {

// на каждом шаге получаем новый вектор коэффициентов

push_back(&v, DanilevskyImpl(&m, n, a));

}

// перемножаем все многочлены

array coeff = collect(v);

// у нас коэффициенты идут в обратном порядке

for (i = 0; i < n; ++i) {

q[i] = coeff.data[n - 1 - i];

}

q[n] = 1; // старший коэффициент всегда равен 1

deleteArray(&coeff);

deleteVector(&v);

}

int main(int argc, char **argv)

{

int n; // размер исходного массива

int i,j, rezult; // рабочие переменные

char c;

printf("\n n=");

scanf("%d", &n);

float a[n][n]; // исходная матрица, после вычислений - матрица Фробениуса

float q[n + 1]; // вектор коэффициентов ее характеристического полинома по возрастающим степеням

// ввод исходной матрицы

printf("\n Vvesti matricu\n");

for (i = 0; i < n; i++) {

printf("\n stroka %d (%d chisel):\n", i + 1, n);

for (j = 0; j < n; j++)

scanf("%f", &a[i][j]);

}

danilevsky(n, a, q);

printf ("\nMatrica Frobeniusa:\n");

for (i = 0; i < n; i++) {

for (j = 0; j < n; j++)

printf(" %10.5f", a[i][j]);

printf("\n");

}

printf("\n Koeff. harakterist. polinoma:\n");

for (i = 0; i < n + 1; i++)

printf(" %12.6f", q[i]);

printf("\n\n");

return 0;

}

4. Проверка работы программы

Пример: вычислить матрицу Фробениуса и коэффициенты характеристического полинома для следующей матрицы (n=4):

-4.0 -3.0 1.0 1.0

a[4][4]= 2.0 0.0 4.0 -1.0

1.0 1.0 2.0 -2.0

1.0 1.0 -1.0 -1.0

Результат программы:



Заключение

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

Список использованных источников

1. Мехеев С.Е.- Численные методы / Мехеев С.Е.- Новосибирск, 2013. С. 73-77.

Размещено на Allbest.ru