Расчет матричной экспоненты на основе преобразования подобия с использованием функций MathCad

В составе встроенных функций векторно-матричной алгебры MathCad имеются функции eigenvals(A) и eigenvecs(A).

Функция eigenvals(A), принимая в качестве аргумента квадратную матрицу A, возвращает набор ее собственных чисел, сформированный в виде вектора.

Функция eigenvecs(A), принимая матрицу A, возвращает матрицу собственных векторов T.

Рассмотрим пример использования этих функций. Зададим квадратную матрицу A в виде матричной константы:

Вызовем функцию расчета собственных чисел этой матрицы и поместим результат в вектор Λ:

, (2.27)

и выведем значение вектора Λ в численной форме:

Вызовем функцию расчета собственных векторов матрицы A и поместим результат в матрицу T:

, (2.28)

и выведем значение матрицы T в численной форме:

Таким образом, встроенные функции MathCad Позволяют получить собственные функции и собственные вектора матрицы A.

Функции eigenvals(A) и eigenvecs(A) используют для расчета собственных чисел и собственных векторов численные методы, что неизбежно вызывает, хотя и не большую, ошибку. Оценим уровень ошибок, вычислив матрицу T^-1∙A∙T, которая теоретически должна быть диагональной:

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

Для вычисления матричной экспоненты e^At на MathCad нам потребуется сформировать диагональную матрицу вида:

, (2.29)

причем с учетом физической реализуемости, нужно доопределить экспоненты в диагонали так, чтобы они были равны нулю при отрицательном времени, т.е. DD(Λ, t) должна быть нулевой при t<0.

Сделать это можно, определив функцию пользователя:

EX(a,t) := if(t<0,0,exp(a∙t)) (2.30)

Введенная таким образом функция двух аргументов совпадает с e^at только для t ≥ 0, для t<0 она тождественно равна нулю.

В составе палитры матричных операций MathCad имеется оператор векторизации функций (см.п.2.7 в [1]). В нашем случае мы подвергнем векторизации введенную нами функцию EX(a,t) по параметру a.

Функция – это векторная функция с элементами:

(2.31)

П

(2.32)

рименив к векторной функции функцию diag из состава встроенных функций MathCad, получим диагональную матрицу вида:

Матрица (2.32) подобна физически реализуемой матричной экспоненте e^At. Остается только воспользоваться обратным преобразованием подобия, чтобы получить выражение для матричной экспоненты в исходной системе координат:

(2.33)

Выражение (2.33) – оператор MathCad, задающий функцию расчета матричной экспоненты (обозначена Ф(Λ,t)) , как функцию вектора собственных значений Λ и времени t.

Итак, порядок расчета матричной экспоненты можно представить в виде следующей последовательности операций:

Определяем собственные числа матрицы A с помощью функции eigenvals(A).

Определяем матрицу собственных векторов матрицы A с помощью функции eigenvecs(A).

Создаем скалярную функцию – «физически реализуемую» экспоненту EX(a,t) := if(t<0,0,exp(a∙t)).

С помощью оператора векторизации превращаем скалярную функцию EX(a,t) в векторную функцию и, используя функцию diag, формируем диагональную матрицу .

Формируем оператор вычисления матричной экспоненты,.: