Решение систем линейных уравнений с невырожденной матрицей

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

Тамбовский государственный университет им. Г.Р. Державина

Институт математики, физики и информатики

Кафедра математического анализа

КУРСОВАЯ РАБОТА НА ТЕМУ:

«РЕШЕНИЕ СИСТЕМ ЛИНЕЙНЫХ  УРАВНЕНИЙ С НЕВЫРОЖДЕННОЙ МАТРИЦЕЙ»

                   Выполнили:

                   студенты 25 гр.  Мешкова Татьяна

                                                                                        Сторожев Андрей

                   Научный руководитель:

                   доктор ф-м. наук, проф. Г.И. Малашонок

Тамбов 2012

Содержание

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

1.Цели и задачи…………………………………………………………………4

2. Решение систем линейных уравнений………………………………..…….5

3.Число обусловленности..………………………………...............................6-9

4.Алгоритм решения систем линейных уравнений.………………….….10-15

5.Примеры.....……………………………….................................................16-19

6.Заключение………………………………………………………………...…20

Список использованной литературы…………………………………….……21

Ведение

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

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

              К основным задачам линейной алгебры можно отнести задачи:

              1.Решения систем линейных алгебраических уравнений.

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

              3.Нахождение собственных значений и собственных функций матриц.

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

Цели и задачи

Цель: Разработать  программы реализующие  решение систем линейных  уравнений с невырожденной матрицей.

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

Задачи:

      Разработать  программу для   реализации  решений систем линейных уравнений с невырожденной матрицей

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

Решение систем линейных уравнений

Приступаем к рассмотрению квадратной матрицы

Пусть также нам задан вектор-столбец

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

где неизвестным является вектор-столбец

В развернутой форме задача может быть записана следующим образом

Решение системы - совокупность n чисел c1, c2, …, cn, таких что подстановка каждого ci вместо xi в систему обращает все её уравнения в тождества.

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

Число обусловленности матрицы - числовая характеристика вычислительной сложности обращения матрицы.

21

Число обусловленности.

              Число обусловленности квадратной матрицы A определяется, как

k(A) = ||A||·||A^(-1)||

Число обусловленности имеет следующее значение: если машинная точность, с которой совершаются все операции с вещественными числами, равна ε, то при решении системы линейных уравнений Ax = b результат будет получен с относительной погрешностью порядка ε·k(A). Хотя число обусловленности матрицы зависит от выбора нормы, если матрица хорошо обусловлена, то её число обусловленности будет мало при любом выборе нормы, а если она плохо обусловлена, то её число обусловленности будет велико при любом выборе нормы. Таким образом, обычно норму выбирают исходя из соображений удобства. На практике наиболее широко используют 1-норму, 2-норму и ∞-норму, задающиеся формулами:

Описание алгоритма

              Если матрица A задана непосредственно, алгоритм оценки числа обусловленности имеет вид:

   вычисление 1-нормы (или inf-нормы) матрицы A прямым алгоритмом за время O(N2)

   построение треугольной факторизации матрицы A за время O(N3)

   получение нижней границы нормы матрицы A-1, используя итеративный алгоритм с трудоемкостью O(N2)

   получение оценки - произведения двух норм. Общая трудоемкость: O(N3)

              Если матрица A задана своей треугольной факторизацией (LU-разложением), то нам нет необходимости осуществлять шаг (2). Однако на шаге (1) мы уже не можем использовать прямой алгоритм - для него требуется знание матрицы A, а не её факторизации. Поэтому используется следующий алгоритм:

   вычисление 1-нормы (или inf-нормы) матрицы A итеративным алгоритмом за время O(N2)

   получение нижней границы нормы матрицы A-1 итеративным алгоритмом за время O(N2)

   получение оценки - произведения двух норм. Общая трудоемкость: O(N2)

              Недостатком итеративного алгоритма является то, что он позволяет получить не само число обусловленности, а лишь его оценку снизу, причем довольно-таки грубую. Обычно оценка меньше самого числа обусловленности лишь на 5-10 процентов, но бывают случаи, когда они отличаются очень значительно (в численных экспериментах оценка в среднем была меньше числа обусловленности на 8%, максимум - на 87%). Обычно когда оценка какой-то величины в десять раз меньше самой величины, эту оценку нельзя назвать точной. Тем не менее, даже такая оценка бывает достаточной.

Для вычисления  1-нормы используется метод rmatrixrcond1 function

Оценка числа обусловленности матрицы (1-норме)

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

Входные параметры:
- Матрица А. Массив, индексы изменяются в пределах [0 .. N-1, 0 .. N-1].
N - размер матрицы A.

Результат: 1/LowerBound (условие (А)) //обратное к нижней границе

ПРИМЕЧАНИЕ:
Если k(А) очень велико, то матрица считается вырожденной, k(А) = INF,
0.0 возвращается в таких случаях.

public static double rmatrixlurcond1(Double[][] lua,

            int n)

        {

            double result = 0;

            v = 0.;

            result = rmatrixrcondluinternal(lua, n, true, false, 0);

            return result;

        }

Аналогично для вычисления inf-нормы используется метод  rmatrixlurcondinf function

              Оценка числа обусловленности матрицы (бесконечность нормы).
Алгоритм расчета нижней границы числа обусловленности. В этом случае,
алгоритм не возвращает нижнюю границу числа обусловленности, но возвращает
обратное число (чтобы избежать переполнения в случае особой матрицы).
Входные параметры:
- Матрица А. Массив,индексы изменяются в пределах [0 .. N-1, 0 .. N-1].
N - размер матрицы A.

Результат: 1/LowerBound (условие (А))  //обратное к нижней границе
ПРИМЕЧАНИЕ:
Если k(А) очень велико, то матрица считается вырожденным, k(А) = INF,
0.0 возвращается в таких случаях.

   public static double rmatrixlurcondinf(Double[][] lua,

            int n)

        {

            double result = 0;

            v = 0.;

            result = rmatrixrcondluinternal(lua, n, false, false, 0);

            return result;

        }

              В данных методах ( rmatrixlurcondinf() и rmatrixlurcond1() ) для вычислинений используется метод  rmatrixrcondluinternal();

Далее в функции rmatrixrcondluinternal(...) используя LU разложение мы получаемся det( A )

и det ( A-1 ) .таким образом мы получили некоторое значение k( A ) :

k(A) =

              Умножив полученное значение k( A ) на ε, получим погрешность решения системы линейных уравнений, заданную в виде A*x=B.

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

Алгоритм функции:

* Автоматическое определение вырожденных случаях

* Число обусловленности оценки

* O (N ^ 3) сложности

ВХОДНЫЕ ПАРАМЕТРЫ

    - Массив [0 .. N-1, 0 .. N-1], матрицы системы

    N - размер

    B - массив [0 .. N-1], правая часть

ВЫХОДНЫЕ ПАРАМЕТРЫ

    Информация - то же, что и в RMatrixSolve

    Вес - то же, что и в RMatrixSolve

    X - то же, что и в RMatrixSolve

public static void rmatrixsolve(double[][] a,

            int n,

            double[] b,

            Integer info,

            densesolverreport rep,

            Double x[])

        {

            xa = new double[n+1];

            xb = new double[n+1];

            tx = new Double[n+1];

            bm = new double[n][n];

            xm=new double[n][n];

           // int i_ = 0;

            info = 0;

            if( n<=0 )

            {

                info = -1;

                return;

            }

            for(int i_=0; i_<=n-1;i_++)

            {

                bm[i_][0] = b[i_];

            }

            rmatrixsolvem(a, n, 1, true, info, rep,x);

            if(!main.fl)

            {

                for(int i_=0; i_<=n-1;i_++)

                {

                    x[i_] = xm[i_][0];

                }

            }

        }

Реализованный метод основан, на представление матрицы А в форме LU.

Сперва, наверно, надо пояснить, что такое LU разложение.

LU-разложение — представление матрицы A в виде LU, где L — нижнетреугольная матрица с единичной диагональю, а U — верхнетреугольная. LU-разложение еще называют LU-факторизацией.

              Матрица L является нижнетреугольной с единичной диагональю, поэтому ее определитель равен 1. Матрица U — верхнетреугольная матрица, значит ее определитель равен произведению элементов, расположенных на главной диагонали (I,j=1…n).

      Будем использовать следующие обозначения для элементов матриц L = (lij), U = (uij), ; причем диагональные элементы матрицы L: lii = 1,                        .   Тогда, если известно LU-разложение матрицы, её определитель можно вычислить по формуле det(A) = det(LU) = det(L)det(U) = det(U)

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

Для

1.                                                                            (1)

2.                                                                                           (2)

              В итоге мы получим матрицы — L и U. В программной реализации данного метода (компактная схема Гаусса) для представления матриц L и U можно обойтись всего одним массивом, в котором совмещаются матрицы L и U. Например вот так(для матрицы 3х3 размером ):

              Фрагмент программы на Java для нахождения матриц L и U.

//переменная n(размерность исходной ''квадратной'' матрицы) должна получить значение до этого момента

            double[][] A = new double[n][n];

            double[][] L = new double[n][n];

            double[][] U = new double[n][n];

//до этого момента массив A должен быть полностью определен

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

            {

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

                {

                    U[0][i] = A[0][i]

                    L[i][0] = A[i][0] / U[0][0];

                    double sum = 0;

                    for (int k = 0; k < i; k++)

                    {

                        sum += L[i][k] * U[k][j];

                    }

                    U[i][j] = A[i][j] - sum;

                    if (i > j)

                    {

                        L[j][i] = 0;

                    }

                    else

                    {

                        sum = 0;

                        for (int k = 0; k < i; k++)

                        {

                            sum += L[j][k] * U[k][i];

                        }

                        L[j][i] = (A[j][i] - sum) / U[i][i];

                    }

                }

            }

              Если матрица A исходной системы разложена в произведение треугольных L и U, то, значит, вместо мы можем записать эквивалентное  уравнение

Lux=b.

              Введя вектор вспомогательных переменных y=(y1,y2,….,yn)T, последнее можно переписать в виде системы

Ly=b

                     (3)

Ux=y

Решение системы, реализованной на JAVA:

y[0]=b[0];

    for(int i=1;i<n;i++){

        y[i]=b[i];

        for(int j=0;j<i;j++){

           y[i]-= lua[i][j]*y[j];

        }

    }

    x[n-1]=y[n-1];

    for(int i=n-2;i>=0;i--){

        x[i]=y[i];

        for(int j=n-1;j>=i;j++){

           x[i]-= lua[i][j]*x[j];

        }

    }

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

Получим сначала формулы для вычисления элементов yi.

Очевидно, все yi могут быть последовательно найдены.

Развернём теперь Ux=y.

Итак, решение СЛАУ посредством LU-факторизации сводится к:

1.      представление матрицы А в виде  LU-факторизации,

2.      решение системы (3)