Системы линейных алгебраических уравнений

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

Высшего образования

«ФИНАН

«ФИНАНСОВЫЙ УНИВЕРСИТЕТ ПРИ ПРАВИТЕЛЬСТВЕ

РОССИЙСКОЙ ФЕДЕРАЦИИ»

Департамент анализа данных, принятия решений и финансовых технологий

 

С.А.Зададаев

 

 

Матричные уравнения (RStudio)

Учебно-методические рекомендации для проведения

семинара №29 по компьютерному практикуму

С индивидуальными заданиями

 

Для бакалавров направления 38.03.01 «Экономика»

 

 

Электронное издание

 

 

Москва 2018

Матричные уравнения (RStudio)

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

 

Системы линейных алгебраических уравнений

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

                                                              (1)

эквивалента матричному уравнению

,

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

,   и .

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

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

Задание 1. Решить приближенно систему (1)

                          

Решение.  Введем соответствующие матрицы  и  на листе Excel и скопируем их через буфер обмена в R (Рис.1):

Рис.1

Здесь мы объединили на листе Excel две матрицы  и  в одну, чтобы их было проще скопировать в буфер обмена и считать в промежуточную матрицу W.

Далее считываем буфер обмена в R и образуем матрицы  и :

# В excel сформировать таблицу чисел и скопировать в буфер обмена 

Data <- read.table("clipboard", h=FALSE, dec=",", sep = "\t") # Чтение из буфера обмена excel-формата

W <- data.matrix(Data)      # Объявить таблицу чисел Data матрицей W в R

#-------------------------------------------------------------------------------------

A <- W[, 1:3]                # Задать матрицу A как первые 3 столбца таблицы W

B <- W[, 4]                   # Задать матрицу B как 4-ой столбец таблицы W

A                                   # Посмотреть матрицу A

B                                   # Посмотреть матрицу B

 

Рис.2

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

X <- solve(A, B)       # Найти приближенное решение системы уравнений AX=B

X                                 # Посмотреть решение X

с результатом:

> X <- solve(A,B) # Найти приближенное решение системы уравнений AX=B> X              # Посмотреть решение XV1 V2 V3 -3 2 5

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

Матричные уравнения

Задание 2. Для матрицы A из предыдущего примера

решить приближенно матричное уравнение .

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

E <- diag(3)          # Вводим единичную матрицу E (3x3)

X <- solve(A, E)   # Решаем уравнение AX= E

X                           # Выводим результат в консоль   

с результатом:

> E <- diag(3)          # Вводим единичную матрицу (3x3)> X <- solve(A, E)      # Решаем уравнение AX=E> X                     # Выводим результат в консоль       [,1] [,2]   [,3]V1 0.2727273 0.1818182 -0.09090909V2 0.1818182 0.4545455 0.27272727V3 1.0000000 1.0000000 1.00000000

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

X <- solve(A)      # Находим обратную матрицу к A

X                          # Смотрим результат

действительно, приводит к тому же результату.

Задание 3. Для матрицы Q

Где А*- номер по классному журналу.

найти обратную .

Решение.  Точное значение определителя матрицы , следовательно, матрица вырождена и её обратная матрица не определена (не существует). Кстати, «машинный» ноль определителя здесь вовсе не ноль, хотя и очень маленький, порядка :

det(Q)           # Определитель матрицы Q      

 > det(Q) # Определитель матрицы Q

[1] 6.661338e-16

Посмотрите, как отреагирует процедура solve на нашу задачу:

Q <- matrix(1:9, nrow = 3, ncol = 3, byrow = TRUE); Q # Задаем и отображаем матрицу Q

solve(Q)                                                                                        # Находим обратную матрицу к Q                                                                      

В результате компилятор выдаёт вполне ожидаемую ошибку:

> Q <- matrix(1:9, nrow = 3, ncol = 3, byrow = TRUE); Q [,1] [,2] [,3] [1,] 1 2 3 [2,] 4 5 6 [3,] 7 8 9 > solve(Q) Error in solve.default(Q): system is computationally singular: reciprocal condition number = 2.59052e-18

Замечание. В версии 3.4.2 и выше в данном сообщении будет указано на точную вырожденность матрицы Q.

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

MySolve <- function(A){
tryCatch(solve(A), error = function(x) print("Матрица необратима!"))
}

 

# Пример использования функции MySolve:

MySolve(Q)     # Находим обратную матрицу или получаем сообщение о вырожденности

MySolve(A)      # Находим обратную матрицу или получаем сообщение о вырожденности

Задание 4. Рассмотрим в качестве следующего примера задачу о нахождении матрицы X, удовлетворяющую уравнению

,

где , а   E – единичная (3x3).

Где А*- номер по классному журналу.

Решение. Для нахождения матрицы X необходимо обе части уравнения слева домножить на матрицу A:

и полученное уравнение справа домножить на квадрат обратной матрицы к B:

В итоге нам необходимо реализовать вычисление для матрицы

= .

Введем данные в R из Excel аналогично через буфер обмена (Рис.3):

 

Рис.3.

Здесь мы также временно объединили две матрицы в одну W и собираемся в R левую половину объявить матрицей A, а правую – матрицей B.

# В excel сформировать таблицу чисел и скопировать в буфер обмена

Data <- read.table("clipboard", h=FALSE, dec=",", sep = "\t") # Чтение из буфера обмена

W <- as.matrix.data.frame(Data) # Объявить таблицу чисел матрицей в R

A <- W[, 1:3]                              # Объявить матрицей A первые три столбца матрицы W

B <- W[, 4:6]                             # Объявить матрицей B последние три столбца матрицы W

Далее нам остается вызвать соответствующее матричное произведение

=  

X <- A%*%(solve(B)%*%solve(B)); X # Вычислить ответ A*((B^-1)*(B^-1))

library(expm)                                       # Подключаем библиотеку для степеней матриц %^%

X <- A%*%solve(B%^%2); X          # То же самое, но короче... A*((B^-1)^2)

Eсли требуется выгрузить ответ X обратно в Excel, то необходимо добавить:

# Скопировать в буфер обмена матрицу X в excel-формате:

write.table(X,"clipboard", quote=FALSE, col.names = FALSE, row.names = FALSE,

                                 sep = "\t",dec=",")

(и не забудьте выгрузить ответ X обратно в Excel по сочетанию Ctrl + v)

Результат вычислений представлен на рис.4

Рис.4

    Обратите внимание на то, что оператор solve(B) возвращает обратную матрицу, а оператор solve(B%^%2) – квадрат обратной матрицы. Вообще, используя определение обратной матрицы, несложно доказать, что квадрат обратной матрицы совпадает с обратной для квадрата матрицы:

,

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

 

Системы нелинейных алгебраических уравнений*

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

                                                                                     (1)

методом Ньютона-Рафсона.

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

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

Успехом применения данного метода служит удачный выбор начального приближения решения системы.

 

Задание 5. Найти точки пересечения прямой  и окружности .

Где А*- номер по классному журналу

Решение. В первую очередь для уверенного выбора начального приближения построим графики прямой и окружности:

x.line <- seq(-1, 5, length.out = 50) # Для прямой разбиваем отрезок [1, 5] на 50 точек x.line

y.line <- 2-x.line                            # Вычисляем соответствующие орлинаты прямой y.line

t <- seq(0, 2*pi, length.out = 50)        # Задаем параметр-угол t, описывающий окружность

x.circle <- 2 + 2*cos(t)                   # Вычисляем абсциссы окружности x.circle

y.circle <- 2*sin(t)                          # Вычисляем ординаты окружности y.circle

plot(x.line, y.line, type = "l", lwd = 2, col = "blue")     # Строим прямую

abline(h = 0, v = 0, lty = "84")                                      # Отображаем оси координат

lines(x.circle, y.circle, type = "l", lwd = 2, col = "red")  # Достраиваем окружность

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

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

Из графиков на полученном рисунке 5

Рис.5.

хорошо видны точки пересечения кривых. Левый корень системы имеет координаты примерно , а правый – :

text(0.5, 1.5, "A")

text(3.5, -1.5, "B")

    Найдем их уточненные значения. Объявим нашу систему уравнений

в виде специальной функции f, где вместо переменной  будем использовать , а вместо переменной  соответственно :

f <- function(x) { c(Line = x[1] + x[2] - 2, Circle = (x[1] - 2)^2 + x[2]^2 - 4) }

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

    Далее активируем библиотеку rootSolve:

library(rootSolve)    # Активация пакета "rootSolve"

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

В окрестности точки A:

multiroot(f, start = c(0.5, 1.5), ctol = 1e-8) # Метод Ньютона-Рафсона

> multiroot(f, start = c(0.5, 1.5), ctol = 1e-8) # Метод Ньютона-Рафсона с начальным приближением (0.5, 1.5)$root[1] 0.5857864 1.4142136  # Это найденное решение системы (левая точка пересечения) $f.root    Line   Circle -9.392487e-14 9.046097e-12 # Насколько отличаются от нулей правые части уравнений $iter                  # Количество проделанных итераций[1] 4 $estim.precis[1] 4.570011e-12       # Достигнутая точность

 

В окрестности точки B:

multiroot(f, start = c(3.5, -1.5), ctol = 1e-8) # Метод Ньютона-Рафсона

> multiroot(f, start = c(3.5, -1.5), ctol = 1e-8) # Метод Ньютона-Рафсона с начальным приближением (3.5, -1.5)$root[1] 3.414214 -1.414214  # Это найденное решение системы (правая точка пересечения) $f.root    Line   Circle -1.976197e-14 8.952838e-12 # Насколько отличаются от нулей правые части уравнений $iter[1] 4                 # Количество проделанных итераций $estim.precis[1] 4.4863e-12        # Достигнутая точность