Работа опубликована в интернет-журнале "Математика в Вузе", №3

Вернуться на страницу <Методические разработки>

Задана матрица A (7 на 7). Требуется найти ее Жорданову форму J, построив базис "сверху вниз" (как башенную диаграмму). По виду J определить поведение Aⁿ и AⁿX при неограниченном росте n. С помощью экспоненты матрицы найти решение системы обыкновенных дифференциальных уравнений dZ/dt = AZ (Z=Z(t), столбец).

Для решения был использован Maple 6 (^tm).

> restart:with(linalg):

Warning, the protected names norm and trace have been redefined and unprotected

> (ввод матрицы A) A:=matrix([[-17,25,-7,-12,-2,-5,-22],[-9,12,4 и т.д. по строкам]]);

> charpoly(A,lambda); lambda:=eigenvals(A);

> E:=diag(1,1,1,1,1,1,1): B1:=evalm(A-lambda[1]*E);

> C:=vector(7,0):

Решаем систему B1*X=0, решение - собственный вектор e1.

> e1:=linsolve(B1,C);

Используем найденный базисный вектор e1 как первый столбец матрицы перехода S.

> S1:=matrix(7,1,0):

> for i from 1 to 7 do S1[i,1]:=subs(_t[1]=1,e1[i]); end do: evalm(S1);

Теперь обрабатываем корневое подпространство второго собственного числа, - строим его Жорданову "башню".

> B2:=evalm(A-lambda[2]*E);

> rank(B2); rank(B2^2); rank(B2^3); rank(B2^4); rank(B2^5); rank(B2^6);

Выводы:

Вторая башня имеет 3 этажа, в ней всего 6 векторов (6= алг. кратность lambda[2]). Почему? Да ведь мы искали сначала, "сколько" векторов обнуляются онократным применением B2. Оказалось, они составляют 7 - 4 = 3 мерное пространство (столько векторов на первом этаже башни). Затем искали, "сколько" векторов обнуляютcя двукратным применением B2. Их, разумеется, больше: 7-2 = 5-мерное пространство. Столько векторов у нас будет на ДВУХ этажах. Затем искали "сколько" векторов обнуляются трехкратным применением B2. Таких - еще больше: 7-1 = 6 - мерное пространство. Столько у нас векторов должно быть на ТРЕХ этажах. А потом мы ничего нового уже не получили: ранг матрицы с ростом степени прекратил уменьшаться и мы не можем обнаружить векторов, для обнуления которых понадобилось бы более. чем 3 воздействия матрицей B2.

Итог:

В трех этажах векторы число л.н. векторов = (7 неизв. - 1 ур-ие) = 6

В первых двух этажах (B2^2*X=0) число л.н. векторов = 7-2 = 5

В первом этаже (B2*X=0) число л.н. векторов (собственных) = 7-4 = 3

Следовательно, на третьем этаже 6-5=1 вектор, на втором 5-3=2 вектора, на первом - 3 вектора.

Теперь можем нарисовать СТРУКТУРУ БАШНИ.

*

**

***

000

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

2000000

0110000

0011000

0001000

0000110

0000010

0000001

Теперь матрицу перехода S можно искать, решая систему линейных уравнений SJ=AS

(7*7 = 49 неизвестных Sij и столько же уравнений).

Но нам задано построить Жорданову башню (базис, обладающий Жордановой структурой) "сверху вниз".

Для связи с дальнейшим присвоим векторам второй башни (см. ее структуру выше) имена:

e4

e3 e6

e2 e5 e7

0 0 0

Напоминаю, что все они под воздействием матрицы B2 отображаются вниз, на векторы

предыдущего этажа, а первый этаж - в 0.

Ищем векторы на этажах башни. Решаем методом Гаусса системы (B2^k)*X=0, k=1..3.

6-мерная линейная оболочка векторов всех трех этажей второй башни:

> kernB23:=linsolve(B2^3,C);

5 - мерная линейная оболочка векторов первых двух этажей второй башни:

> kernB22:=linsolve(B2^2,C);

3 - мерная линейная оболочка векторов первого этажа второй башни (подпр-во собств. векторов lambda [2]):

> kernB21:=linsolve(B2,C);

Это сделано так же, как делается вручную. Получены решения систем в параметрической форме. Но имеется и готовая команда kernel(A) (или nullspace(A)), - она дает базисные комбинации значений неизвестных для системы AX=0. Для наглядности я здесь буду действовать, как при решении методом Гаусса вручную. Поставим базисные векторы полученных общих решений в матрицы

> z3:=matrix(7,6): for j from 1 to 6 do h:=[seq(E[j,k],k=1..6)]: _t:=h:

for i from 1 to 7 do   z3[i,j]:=subs(_t=h,kernB23[i]) end do: end do:

> z2:=matrix(7,5): for j from 1 to 5 do h:=[seq(E[j,k],k=1..5)]: _t:=h:

for i from 1 to 7 do z2[i,j]:=subs(_t=h,kernB22[i]) end do: end do:

> z1:=matrix(7,3): for j from 1 to 3 do h:=[seq(E[j,k],k=1..3)]: _t:=h:

for i from 1 to 7 do z1[i,j]:=subs(_t=h,kernB21[i]) end do: end do:

Теперь надо выделить векторы 3-го этажа, они лежат в kernB23, но не лежат в KernB22. Составляем из матриц z2,z3 (базисных векторов-столбцов пр-в KernB22 и KernB23) матрицу M23

> M23:=concat(z2,z3);

Здесь первые 5 столбцов - базис линейной оболочки первых двух этажей (пр-ва KernB22),

последние 6 столбцов - базис линейной оболочки трех этажей (всей башни, пространства KernB23) Подпр-во kernB22 содержится в kernB23 (т.к. из равенства B2^2*X=0 следует равенство B2^3*X=0).

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

Команда gaussjord(A) действует по строкам (проверено опытом), - столбцы не смешиваются друг с другом, и мы можем по единичкам в итоговом виде матрицы выяснить, какие столбцы в M23 надо брать. Поскольку нет полной гарантии, что Maple не обнулит именно ту группу векторов-столбцов, которую мы хотим оставить в живых, то, применяя к матрице M23 команду gaussjord, следим за смыслом результата.

> M23a:=gaussjord(M23);

Судя по "ступенькам" ранг =6; первые 6 столбцов матрицы M23 линейно независимы (напоминаю, команда gaussjord не смешивала столбцы при преобразовании M23 в M23a), причем 6-ой столбец не зависел от первых пяти. Значит он ДОПОЛНЯЕТ базис пр-ва kernB22 до базиса пр-ва kernB23 и поэтому является верхним этажом (или, если угодно, крышей) башни. Теперь получим все векторы башни, расположенные на вертикали под ним. Заодно убедимся, что подвал вертикали - нулевой, как и полагается (проверка).

> e4:=submatrix(M23,1..7,[6]): e3:=evalm(B2&*e4): e2:=evalm(B2&*e3): evalm(transpose(B2&*e2));

Наращиваем матрицу перехода S векторами e2, e3, e4 (снизу вверх по вертикали в башне ).

> S1_4:=concat(S1,e2,e3,e4);

Для проверки правильности S1_4 можно вычислить (знак # - признак комментария, - не команды, - уберите его, если надо выполнить).

> # evalm(B2&*S1_4); evalm(B2^2&*S1_4); evalm(B2^3&*S1_4):

Теперь надо найти вершину второй вертикали (см. схему башни вверху). Векторы первых ДВУХ этажей - базис подпространства kern (B2^2). Но векторы ВТОРОГО этажа удовлетворяют еще одному требованию: при однократном воздействии матрицы B2 они ДОЛЖНЫ НЕ обнуляться. Множество всех векторов, удовлетворяющих таким ДВУМ требованиям, не является подпространством, т.к. их линейная комбинация может нарушить второе требование. Это надо помнить. Поступим так: уже известный нам вектор e2 дополним до найденного базиса трехмерного пр-ва kernB21 (первого этажа) (этот базис нас не устраивает как раз тем, что в нем нет вектора e2).

> M12:=concat(e2,z1);

> M12a:=gaussjord(M12);

Т.е. первые три вектора из M12 нам годятся как базис пр-ва kernB21, имеющий в своем составе вектор e2. Добавим к ним вектор e3 из второго этажа первой вертикали и дополним этот набор векторов до базиса 5 - мерного пр-ва kernB22 (двух этажей)

> M22:=concat(e3,submatrix(M12,1..7,1..3),z2); gaussjord(M22);

Очевидно, 5-ый столбец MM не лежит в kernB21, лежит в kernB22, отличен от e3. Т.е. это - еще один вектор из второго этажа. Используем его как e6, - вершину второй вертикали. Воздействуя на него матрицей B2, найдем вектор на 1 этаже под ним (вектор e5)

> e6:=submatrix(M22,1..7,[5]): e5:=evalm(B2&*e6):

Теперь осталось дополнить векторы e2, e5 первого этажа до базиса kernB21 (всего первого этажа).

> M11:=concat(e2,e5,z1): gaussjord(M11);

Годится третий из M11. Ставим их в матрицу перехода S.

> e7:=submatrix(M11,1..7,[3]):S:=concat(S1_4,e5,e6,e7);

Теперь вычисляем S^(-1)*A*S и получаем J

> J:=evalm(inverse(S)&*A&*S);

Теперь ПРИЛОЖЕНИЯ.

Maple, кажется, не знает формулу для J^(n). Поэтому просто вычислим ее по биному Ньютона для n>1.

> M:=diag(2,1,1,1,1,1,1); N:=matrix(7,7,0): for i from 1 to 6 do N[i,i+1]:=J[i,i+1] end do: evalm(N);

Теперь раскладываем (M+N)^(n) по биному Ньютона и помним, что N^3=0.

> n:='n': Jn:=evalm(diag(2^n,1,1,1,1,1,1)+ binomial(n,1)*diag(2^(n-1),1,1,1,1,1,1)&*N+ expand(binomial(n,2))*diag(2^(n-2),1,1,1,1,1,1)&*N^2);

Очевидно, что при неограниченном росте n первая ячейка неограниченно экспоненциально растет, во второй ячейке имеется линейный (Jn[2,3] и Jn[3,4]) и квадратичный (Jn[2,4]) рост. В третьей ячейке видим линейный рост (Jn[5,6]). Не меняется только последняя ячейка = Jn[7,7]. В матрице J^n есть неограниченно растущие элементы. A^(n) = S*J^(n)*S^(-1) и S не вырождена. Поэтому c ростом n в матрице A^(n) будут неограниченно расти некоторые элементы.

Теперь о поведении A^(n)*X с ростом n. Поскольку A^(n)*X = S*J^(n)*S^(-1)*X, то, обозначив Y=S^(-1)*X, потребуем, чтобы вектор Y имел нулевыми некоторые элементы так, чтобы в произведении J^(n)*Y был уничтожен вклад от "плохих" ячеек.

Итог: векторы X=S*Y имеют образом A^(n)*X, которое вообще не меняется с ростом n.

> Y:=matrix([[0],[t1],[0],[0],[t2],[0],[t3]]): X:=evalm(S&*Y); AnX:=evalm(inverse(S)&*Jn&*Y);

То же самое можно было получить в терминах собственных векторов - только собственные векторы с собственным числом=1 не меняются под воздействием A.

Теперь решаем систему дифф. уравнений dZ/dt = AZ, Z(0)=Z0. Решение Z=exp(At)*Z0 и Maple умеет вычислять экспоненту. (Уберите знак # если хотите увидеть ответ). Нам важно понимать, почему она такая.

> # expAt:=exponential(A*t): Z0:=vector(7): evalm(expAt&*Z);

Если подставить A*t = S*J*S^(-1) в ряд Тэйлора - Маклорена для экспоненты, то получится exp(A*t) = =S*exp(J*t)*S^(-1). Експоненту Жордановой клетки можно вычислить в формульном виде, если разложить J*t = lambda*E*t + N (где N имеет единички над диагональю), воспользоваться тем, что E и N коммутируют, поэтому exp(lambda*E*t+N*t) = exp(lambda*E*t)*exp(N*t). Вычислить последние две экспоненты можно с помощью формулы Тэйлора - Маклорена, если заметить, что E^(n) =E, а при n = m (порядок ячейки) получаем N^(m)=0. Так получится формула для J^n, которую можно запомнить. Мы напишем ее для ячейки 4 порядка. Предварительно освободим переменную lambda от ее предыдущего смысла:

> lambda:='lambda': t:='t': J:='J': J:=evalm(lambda*diag(1,1,1,1)):  for i from 1 to 3 do J[i,i+1]:=1 end do: J:=evalm(J); expJt:=exponential(J*t);

Зная эту формулу, можно вычислить exp (A*t) = S *exp(J*t)*S^(-1) и затем получить решение системы dZ/dt=AZ в виде Z=exp(A*t)*Z0, где Z0 - вектор начального состояния системы (начальных условий).

Наверх

Вернуться на страницу <Методические разработки>