\documentclass{report}
\usepackage{amsmath,amsthm,amssymb}
\usepackage{mathtext}
\usepackage[T1,T2A]{fontenc}
\usepackage[utf8]{inputenc}
\usepackage[english,russian]{babel}
\usepackage[left=25mm, top=20mm, right=10mm, bottom=20mm, nohead, nofoot]{geometry}
\usepackage{listings}


\begin{document}

\begin{center}
\normalsize{\bfМОСКОВСКИЙ ГОСУДАРСТВЕННЫЙ УНИВЕРСИТЕТ \\
им М.В. Ломоносова \\}
\
\\
\normalsize{ Факультет вычислительной математики и кибернетики}\\
\noindent\rule{\textwidth}{1pt} %линия

\topskip=0pt
\vspace*{\fill}
    \normalsize{\bf
        Компьютерный практикум по учебному курсу \\
        <<ВВЕДЕНИЕ В ЧИСЛЕННЫЕ МЕТОДЫ>> \\
        ЗАДАНИЕ №1 \\
        \
        \\
        ОТЧЕТ \\
        о выполненном задании \\
    }
    \normalsize{
        студента \underline{206} учебной факультета ВМК МГУ \\
        \underline{Арутюнова Эдгара Арамовича }
        \\
    }
    {\tiny
        (фамилия имя отчество студента)
    }
\vspace*{\fill}


гор.Москва \\ 2018 \\

\thispagestyle{empty}

\end{center}

\newpage %титульный

\tableofcontents

\newpage %начало

\parindent=1cm

\chapter{Постановка задачи и ее целей}
\section{Цель работы}
\

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


Дана система уравнений $Ax=f$ порядка $n*n$ с невырожденной матрицей $A$.
Написать программу,
решающую систему линейных алгебраический уравнений
заданного пользователем
размером ($n$ - параметр программы) методом Гаусса и
методом Гаусса с выбором главного
элемента).
 Предусмотреть возможность задания элементов матрицы системы
 и её правой части как
во входном файле данных, так и путём задания специальных формул.

\section{Цели и задачи практической работы}
\
\begin{enumerate}
    \item Решить заданную СЛАУ методом Гаусса.
    и методом Гаусса с выбором главного элемента.
    \item Вычислить опеределитель матрицы $det(A)$.
    \item Вычислить обратную матрицу $A^{-1}$.
    \item Опеределить число обусловленности $M_A=||A|| * ||A_{-1}||$
    \item Исследовать вопрос вычислительной устойчивости метода
    Гаусса (при больших значениях параметра $n$);
    \item Правильность решения СЛАУ подтвердить системой
    тестов используя ресурсы on-line системы \\
    http://www.wolframalpha.com.
\end{enumerate}

\chapter{Описание методов}
\section{Решение СЛАУ методом Гаусса}
\

Метод Гаусса — классический метод решения системы линейных алгебраических уравнений
(СЛАУ).Рассмотрим систему линейных уравнений с действительными постоянными коэффициентами:

\begin{equation*}
 \begin{cases}
   a_{11} \cdot x_{1} +
   a_{12} \cdot x{2} +
   \ldots +
   a_{1n} \cdot x_n = y_{1}
   \\
   a_{11} \cdot x_{1} +
   a_{12} \cdot x{2} +
   \ldots +
   a_{1n} \cdot x_n = y_{1}
   \\
    \vdots
   \\
   a_{n1} \cdot x_{1} +
   a_{n2} \cdot x{2} +
   \ldots +
   a_{nn} \cdot x_n = y_{1}
 \end{cases}
\end{equation*}

Или в матричной форме:

\[
\begin{bmatrix}
    a_{11} \ldots a_{1,n}
    \\ \vdots \ddots \vdots
    \\ a_{n1} ... a_{n,n}
\end{bmatrix}
    \cdot
\begin{bmatrix}
    x_{1}
    \\ \vdots
    \\ x_{n}
\end{bmatrix}
    =
\begin{bmatrix}
    y_{1}
    \\ \vdots
    \\ y_{n}
\end{bmatrix}
\]
 Метод Гаусса заключается в применении следующих элементарных преобразований:
\begin{enumerate}
    \item Перестановка двух уравнений в системе.
    \item Умножение уравнения на число, отличное от нуля.
    \item Прибавлением к одному уравнению другого,
    домноженного на коэффициент, отличныйот нуля.
\end{enumerate}
\
\\
Метод Гаусса решения системы линейных уравнений включает в себя 2 стадии:
\begin{enumerate}
    \item Прямой ход.
    \item Обратной ход.
\end{enumerate}

\section{Прямой ход}
\
Метод Гаусса условно можно разделить на 2 части: прямой ход и обратный ход.
Прямой ход:
На первом этапе матрица коэффициентов с помощью элементарных преобразований
приводится к верхнетреугольному или ступенчатому виду, и устанавливается ее
совместность. Среди элементов первого столбца выбирается ненулевой элемент и строка
с этим элементом меняется с первой строкой матрицы, далее все элементы первой строки
делятся на первый элемент этой строки. Затем первая строка помноженная на первые
элементы строк вычитается из всех последующих строк, тем самым мы обнуляем столбец
под первым элементом. После того, как указанные преобразования были завершены мы
переходим к нижнему правому минору размера на единицу меньше текущей матрицы и
повторяем с ним вышеизложенный алгоритм, пока не достигнем минора размера 0.
Обратный ход:
На втором этапе из каждого уравнения(i-го), начиная с последнего, выражается i
переменная через i элемент столбца b и последующие переменные, которые были
вычислены на предыдущих шагах.
Заметим, что при решении СЛАУ методом Гаусса, если ведущие элементы окажутся очень
маленькими по абсолютной величине, то при делении на такие ведущие элементы
получается большая погрешность округления (вычислительная погрешность). Чтобы
избежать сильного влияния вычислительной погрешности на решение, применяется метод
Гаусса с выбором главного элемента.
Метод Гаусса с выбором главного элемента отличается от метода Гаусса лишь тем, что
при прямом ходе мы выбираем не произвольные ненулевой элемент в первом столбце, а
наибольший.

\section{Число операций(сложение, умножение, деление) }

Общее число операций, необходимых для решения задачи по методу
Гаусса, имеет порядок
$O(n^{3})$.

\chapter{Программа}

\section{Структура}
\
В программе во всех функциях следующие имена переменных обозначают:
$n$ - размер матрицы коэффицентов
$A$ - массив размера $n \cdot (n+1)$, где последний столбец -
столбец свободных
членов, а оставшаяся матрица - матрица коэффициентов СЛАУ
$rev\_matrix$ - обратная матрица

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

void creat\_matrix(double ***A, int n)

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

\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_

void read\_matrix(double ***A, int n, const char *filename)

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

\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_

double *gauss (double **A, double **rev\_matrix, int n)

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

\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_

double *gauss\_with\_main (double **A, double **rev\_matrix, int n)

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

\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_

void init\_rev\_matrix(double ***rev\_matrix, int n)

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

\chapter{Тесты}
\

\section{Ручные тесты}
\
\begin{equation*}
 \begin{cases}
    3 \cdot x_{1} -
    2 \cdot x_{2} +
    2 \cdot x_{3} +
    2 \cdot x_{4}
    =
    8
    \\
    2 \cdot x_{1} -
    1 \cdot x_{2} +
    2 \cdot x_{3}
    =
    4

    \\
    2 \cdot x_{1} +
    1 \cdot x_{2} +
    4 \cdot x_{3} +
    8 \cdot x_{4}
    =
    -1
    \\
    1 \cdot x_{1} +
    3 \cdot x_{2} -
    6 \cdot x_{3} +
    2 \cdot x_{4}
    =
    8
    \\
\end{cases}
\end{equation*}

\begin{equation*}
 \begin{cases}
    2 \cdot x_{1} -
    3 \cdot x_{2} +
    1 \cdot x_{3} +
    2 \cdot x_{4}
    =
    4
    \\
    4 \cdot x_{1} +
    3 \cdot x_{2} +
    1 \cdot x_{3} +
    1 \cdot x_{4}
    =
    5

    \\
    1 \cdot x_{1} -
    7 \cdot x_{2} -
    1 \cdot x_{3} -
    2 \cdot x_{4}
    =
    7
    \\
    2 \cdot x_{1} +
    5 \cdot x_{2} +
    1 \cdot x_{3} +
    1 \cdot x_{4}
    =
    1
    \\
\end{cases}
\end{equation*}


\begin{equation*}
 \begin{cases}
    1 \cdot x_{1} -
    1 \cdot x_{2} +
    1 \cdot x_{3} -
    1 \cdot x_{4}
    =
    0
    \\
    4 \cdot x_{1} -
    1 \cdot x_{2} -
    1 \cdot x_{4}
    =
    0

    \\
    2 \cdot x_{1} +
    1 \cdot x_{2} -
    2 \cdot x_{3} +
    1 \cdot x_{4}
    =
    0
    \\
    5 \cdot x_{1} +
    1 \cdot x_{2} -
    4 \cdot x_{4}
    =
    0
    \\
\end{cases}
\end{equation*}

\section{Генерируемый тест}
\

\begin{equation*}
 A_{ij} =
 \begin{cases}
    q_{M}^{i+j} + 0.1 * (j - i), \ \ i \neq j
    \\
    (q_{M} - 1)^{i+j}, \ \ \ \ \ \ \ \ \ \ \ i = j
\end{cases}
где\ q_{M}=1.001 - 2 \cdot M \cdot 10^{-3}, i, j = \overline{1, ..., n}
\end{equation*}
\newpage

\chapter{Описание программы и её исходный код}
\section{main.c}
\



\begin{lstlisting}

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
double m = 1; //начальное значение определителя
// Заполнение матрицы коэффициентов по формуле
void
creat_matrix(double ***A, int n)
{
    *A = calloc(n, sizeof(**A));
    for(int i = 0; i < n; i++)
    (*A)[i] = calloc(n + 1, sizeof(***A));
    int M = 6;
    double q = 1.001 - 2 * M * 0.001;
    for (int i = 1; i <= n; i++) {
        for (int j = 1; j <= n; j++) {
            if(i == j) {
                (*A)[i - 1][j - 1] = pow(q - 1, i + j);
            }
            else {
                (*A)[i - 1][j - 1] = pow(q, i + j) - 0.1 * (j - i);
            }
        }
        (*A)[i - 1][n] = exp(1 / i) * cos(1 / i);
    }
}
// Считывание матрицы коэффициентов из файла
void
read_matrix(double ***A, int n, const char *filename)
{
    FILE *file = fopen(filename, "r");
    *A = calloc(n, sizeof(**A));
    for(int i = 0; i < n; i++)
    (*A)[i] = calloc(n + 1, sizeof(***A));
    for (int i = 0; i < n; i++) {
        for (int j = 0; j < n + 1; j++) {
            fscanf(file, "%lf", &((*A)[i][j]));
        }
    }
}
// Метод Гаусса
double
*gauss (double **A, double **rev_matrix, int n)
{
    double *x = calloc(n, sizeof(*x));
    double tmp;

    // Прямой ход

    for (int i = 0; i < n; i++)
    {
        int p;
        for (p = 0; p < n && A[p][i] == 0; p++); // Поиск не нулевого
        элемента в столбце
        if (p == n) {
            printf("Система несовместна или имеет более одного решения
            \ndet(A)=0\nОбратная матрица не существует\n");
            free(x);
            return NULL;
        }
        else { // Перестановка строки с ненулевым элементом в начало
            double *c = A[i];
            A[i] = A[p];
            A[p] = c;
            c = rev_matrix[i];
            rev_matrix[i] = rev_matrix[p];
            rev_matrix[p] = c;
        }
        tmp = A[i][i];
        m *= tmp; // Домножаем определитель
        A[i][n] /= tmp;
        for (int j = n - 1; j >= 0; j—) {
            // Деление строки на первый элемент
            A[i][j] /= tmp;
            rev_matrix[i][j] /= tmp;
        }
        for (int j = i + 1; j < n; j++) {
        // Вычитание строки из последующих
            tmp = A[j][i];
            A[j][n] -= tmp * A[i][n];
            for (int k = n - 1; k >= 0; k--) {
                A[j][k] -= tmp * A[i][k];
                rev_matrix[j][k] -= tmp * rev_matrix[i][k];
            }
        }
    }

    // Обратный ход

    x [n - 1] = A[n - 1][n];
    for (int i = n - 2; i >= 0; i--)
    {
        x[i] = A[i][n];
        for (int j = i + 1;j < n; j++) {
            x[i] -= A[i][j] * x[j];
            for (int k = 0; k < n; k++) {
                rev_matrix[i][k] -= A[i][j] * rev_matrix[j][k];
            }
        }
    }
    return x;
}

// Метод Гаусса с выбором главного элемента
double
*gauss_with_main (double **A, double **rev_matrix, int n)
{
    double *x = calloc(n, sizeof(*x));
    double tmp;
    // Прямой ход

    for (int i = 0; i < n; i++)
    {
        double max = fabs(A[i][i]);
        int max_i = i;
        for (int j = i; j < n; j++) {
            // Поиск главного элемента
            if(max < fabs(A[j][i])) {
                max = fabs(A[j][i]);
                max_i = j;
            }
        }
        if(A[max_i][i] == 0)
        {
            printf(
            "Система несовместна или имеет более одного решения
            \ndet(A)=0\nОбратная матрица не существует\n");
            free(x);
            return NULL;
        }
        if(i != max_i) {
            // Перестановка строки с главным элементом в начало
            double *c;
            c = A[i];
            A[i] = A[max_i];
            A[max_i] = c;
            c = rev_matrix[i];
            rev_matrix[i] = rev_matrix[max_i];
            rev_matrix[max_i] = c;
            m *= -1; // Меняем знак определителя
        }
        tmp = A[i][i];
        m *= tmp; // Домножаем определитель
        A[i][n] /= tmp;
        for (int j = n - 1; j >= 0; j--) {
            // Деление строки на первый элемент
            A[i][j] /= tmp;
            rev_matrix[i][j] /= tmp;
        }
        for (int j = i + 1; j < n; j++) {
            // Вычитание строки из последующих
            tmp = A[j][i];
            A[j][n] -= tmp * A[i][n];
            for (int k = n - 1; k >= 0; k--) {
                A[j][k] -= tmp * A[i][k];
                rev_matrix[j][k] -= tmp * rev_matrix[i][k];
            }
        }
    }

    // Обратный ход

    x [n - 1] = A[n - 1][n];
    for (int i = n - 2; i >= 0; i--)
    {
        x[i] = A[i][n];
        for (int j = i + 1;j < n; j++) {
            x[i] -= A[i][j] * x[j];
            for (int k = 0; k < n; k++) {
                rev_matrix[i][k] -= A[i][j] * rev_matrix[j][k];
            }
        }
    }
    return x;
}

// Создание единичной матрицы
void
init_rev_matrix(double ***rev_matrix, int n)
{
    *rev_matrix = calloc(n, sizeof(**rev_matrix));
    for (int i = 0; i < n; i++) {
    (*rev_matrix)[i] = calloc(n, sizeof(***rev_matrix));
    (*rev_matrix)[i][i] = 1;
    }
}

int
main(int argc, const char * argv[]) {
    int n;
    double **A;
    double **rev_matrix;
    double *x;

    /************************************
    Выбор способа задания матрицы коэффициентов:
    1) Через формулу - 1 аргумент программы (размер матрицы)
    2) Матрица коэффициентов записана в файле - 2 аргумента (имя файла,
    размер матрицы)
    ************************************/

    if (argc == 1) {
        printf("Некорректный запуск программы, укажите аргументы");
        return 1;
    }
    if(argc == 2) {
        n = atoi(argv[1]);
        creat_matrix(&A, n);
    }
    else {
        n = atoi(argv[2]);
        read_matrix(&A, n, argv[1]);
    }
    init_rev_matrix(&rev_matrix, n);

    // Выбор метода решения СЛАУ

    printf("Выберите метод:\n1)Метод Гаусса\n2)Метод Гаусса с выбором
    главного элемента\n");
    int id;
    scanf("%d", &id);
    if(id == 1) {
        x = gauss(A, rev_matrix, n);
    }
    else {
    x = gauss_with_main(A, rev_matrix, n);
    }
    if (x == NULL) {
        return 0;
    }
    // Вывод полученной информации

    printf("x = (");
    for (int i = 0; i < n - 1; i++) {
        printf("%.g, ", x[i]);
     }
    printf("%.g)\n", x[n - 1]);
    printf("det(A) = %.10f\n", m);
    printf("Обратная матрица:\n");
    for (int i = 0; i < n; i++) {
        for (int j = 0; j < n; j++) {
            printf("%.g ", rev_matrix[i][j]);
        }
        printf("\n");
    }
    free(x);
    for (int i = 0; i < n; i++) {
        free(A[i]);
        free(rev_matrix[i]);
    }
    free(A);
    free(rev_matrix);
    return 0;
}
\end{lstlisting}


\newpage
\chapter{Вывод программы}
\

>> cat m.txt

<< 3 -2 2 -2 8

<< 2 -1 2 4

<< 2 1 4 8 -1

 << 1 3 -6 2 3

 >> ./m m.txt 4

 << Выберите метод:

 << 1)Метод Гаусса

 << 2)Метод Гаусса с выбором главного элемента

 >> 1

 << x = (3.8571428571, 1.2857142857, -1.0714285714, -0.5714285714)

 << det(A) = 461.99999999999988631316227838397026062011718750000000

 << Обратная матрица:

 << 0.4 0.4 -0.1 -0.3

 << 0.1 0.2 -0.02 -0.3

 << -0.1 -0.1 0.1 0.2

 << -0.1 0.2 -0.02 -0.004

\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_

 >> cat m.txt

 << 2  3  1  2  4

 << 4  2  1  1  5

 << 1 -7 -1 -2  7

 << 2  5  1  1  1

 >> ./m m.txt 4

 << Выберите метод:

 << 1)Метод Гаусса

 << 2)Метод Гаусса с выбором главного элемента

 >> 1

 x = (17.0000000000, 10.0000000000, -106.0000000000, 23.0000000000)

 << det(A) = -1.00000000000000088817841970012523233890533447265625

 << Обратная матрица:

3 -4 3 4

2 -3 2 3

-2e+01 3e+01 -2e+01 -3e+01

5 -6 4 5

\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_

>> cat m.txt

<< 1 -1  1 -1 0

<< 4 -1  0 -1 0

<< 2  1 -2  1 0

<< 5  1  0 -4 0

<< Выберите метод:

<< 1)Метод Гаусса

<< 2)Метод Гаусса с выбором главного элемента

<< 1

<< Система несовместна или имеет более одного решения

<< det(A)=0

<< Обратная матрица не существует



\newpage
\chapter{Исследование вычислительной устойчивости метода
Гаусса}
\

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

Это можно легко проверить, для этого была написана специальная функция, которая
генерирует матрицу Гильберта заданного размера и свободный вектор состоящий из 1:

\begin{lstlisting}
void
g_matrix(double ***A, int n)
{
    *A = calloc(n, sizeof(**A));
    for(int i = 0; i < n; i++)
        (*A)[i] = calloc(n + 1, sizeof(***A));
    for (int i = 1; i <= n; i++) {
        for (int j = 1; j <= n; j++) {
            (*A)[i - 1][j - 1] = 1. / (i + j - 1);
        }
        (*A)[i - 1][n] = 1;
    }
}
\end{lstlisting}

И уже на матрице размера 10 мы получаем существенные отклонения от правильного
результата даже для метода Гаусса с выбором главного элемента.

Результат работы метода Гаусса:
$ x = (-9.9973706711, 989.7722374956, -23755.1401375090,
240195.7620502072, -1261048.7910691723, 3783198.9648144487,
-6725766.1627083644, 7000357.8211527597, -3937735.6951384256,
923673.4643101940)$

Правильный ответ:

$х = (-10, 990, -23760, 240240, -1261260, 3783780, -6726720, 7001280,
-3938220, 923780) $

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


\newpage

\chapter{Вывод}
\

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

\newpage

\chapter{Литература}
\
<<Вводные лекции по численным методам>> Д. П. Костомаров, А. П. Фаворский.

\end{document}