Семинар 2 Уравнение Сильвестра. II

A, B - матрицы размера n x n, m x m соответственно, правая часть Y - матрица размера m x n. Неизвестной является матрица H размера m x n. В теоретическом курсе был установлен следующий результат (см. § 3). Теорема. Если спектры матриц A и B не пересекаются, то для любой Y существует единственное решение уравнения Сильвестра (1). Уравнение Сильвестра (1)

Вопрос: Как найти решение уравнения (1) ? На прошлом семинаре предлагалось сводить решение уравнения (1) к нахождения решения системы алгебраических уравнений Ch = D, (2) где C – матрица, составленная из элементов заданных матриц A и B, вектор D составлен из элементов заданной матрицы Y, вектор h – искомый, его компонентами являются элементы искомой матрицы H. Для нахождения решения системы (2) предлагалось использовать пакет программ, разработанный в Институте математики им. С.Л. Соболева. Уравнение Сильвестра (1)

В настоящее время существуют различные математические пакеты (Maple, Mathematica и др.), которые позволяют находить решение системы (2). Цель: Используя современные пакеты, найти решение системы (2), по нему восстановить решение уравнения Сильвестра (1) и провести сравнительный анализ решений, полученных на прошлом семинаре с использованием пакета программ, разработанного в Институте математики им. С.Л. Соболева. (1) (2)(2)

В зависимости от уровня подготовки студентов предлагаются разные варианты организации вычислительного процесса. Уровень 1. Предлагается готовая программа (в качестве примера ниже дается такая программа для Maple), в которой осуществлены переход от уравнения (1) к системе (2), вычисление решения системы (2), восстановление решения уравнения (1) по решению системы (2), вычисление невязки. Уровень 2. Студентам предлагается самостоятельно провести все расчеты, с использованием одного из пакетов. (1) (2)(2)

Пример программы (Maple) Задаются матрицы A, B, Y так, чтобы можно было следить за спектрами матриц A, B: with(LinearAlgebra): DA:=, >;TA:=, >; detTA:=Determinant(TA); condTA:=ConditionNumber(TA); A:=TA.DA.MatrixInverse(TA); DB:=,,, >; TB:=,,, >; detTB:=Determinant(TB); condTB:=ConditionNumber(TB); B:=TB.DB.MatrixInverse(TB); Y:=,,, >;

Переход от уравнения (1) к системе (2), нахождение решения системы (2), восстановление решения уравнения (1), вычисление невязки. H:=Matrix(4,2,symbol=h): Y1:=H.A-B.H: sys:=[Y1[1,1]=Y[1,1],Y1[1,2]=Y[1,2],Y1[2,1]=Y[2,1],Y1[2,2]=Y[2,2],Y 1[3,1]=Y[3,1],Y1[3,2]=Y[3,2],Y1[4,1]=Y[4,1],Y1[4,2]=Y[4,2]]: var:=[h[1, 1],h[1, 2],h[2, 1],h[2, 2],h[3, 1],h[3, 2],h[4, 1],h[4, 2]]: (HH, f) := GenerateMatrix( sys, var ); condHH:=ConditionNumber(HH); h_sol:=Vector(LinearSolve(HH,f)): H_sol:=,,, >; MatrixNorm(Y-H_sol. A+B.H_sol,2);