voeto.ru страница 1
скачать файл



4. Алгоритм решения системы дифференциальных уравнений, имеющей действительные кратные корни характеристического уравнения.
Пусть имеется характеристическое уравнение некоторого дифференциального уравнения

с заранее известными корнями (три действительных корня, имеющих кратности 1, 2 и 3

соответственно). Например:

expand((p+2)^2*(p+1)^3*(p+3));



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

(СДУ) в нормальной форме. Зададим порядок СДУ

n:=6:

и матрицу А правых частей СДУ (вариантом задания матрицы А может быть также

процедура: seq(seq(a[i,j],j=1..n),i=1..n) )

A:=matrix(n,n,[0,1,0,0,0,0,0,0,1,0,0,0,0,0,0,1,0,0,0,0,0,0, 1,0,0,0,0,0,0,1,-12,-52,-91,-82,-40,-10]):

Запишем вектор решений:

X:=vector(n,[seq(x[j](t),j=1..n)]): (22)

Укажем вектор начальных условий Х0:

X0:=vector(n,[1,0,1,0,1,0]):

По матрице А и вектору решений Х приведем развернутую запись СДУ в нормальной форме

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

(23) seq(x[j](0)=X0[j],j=1..n);





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

k матрицы А:

ke:=eigenvals(A);



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

информация)

kev:=eigenvects(A);


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

векторами k и r в строгом соответствии с номерами компонент. Определение числа

J компонент векторов k и r.


Введем единичную матрицу

E:=matrix(n,n,[]):

for i from 1 to n do

for j from 1 to n do

if i=j then E[i,j]:=1

else E[i,j]:=0 fi:

od:od:

Введем промежуточный вектор, в который запишем корни

kk:=vector(n,[seq(ke[j],j=1..n)]):

Введем вспомогательную матрицу,

K:=matrix(n,n,[]):

в которую будем записывать нули и единицы в соответствии с ниже приводимым

алгоритмом вычисления кратности корня.

Приведем алгоритм вычисления кратности корня в предположении, что сами корни

известны (например, корни определяются процедурой eigenvals(A)) и представлены в

виде списка или вектора.

for i from 1 to n-1 do

for j from i+1 to n do

if kk[i]=kk[j] and kk[i]<>999 then

K[i,j]:=1: kk[j]:=999:

else K[i,j]:=0 fi:

od:od:

evalm(K): визуализация полученной матрицы К.

Кратность корня запишем в вектор

rv:=vector(n,[]):

По матрице К определяем кратность корня kk[i]

for i from 1 to n do

if kk[i]<>999 then

rv[i]:=add(K[i,j], j=i+1..n)+1:

else rv[i]:=0 fi:

od:

В результате получили вектор, компоненты которого указывают кратность корня: номер

компоненты это номер корня, а её значение - кратность корня (берутся только те компо-

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

eval(rv); rv[1]:



Вычисление числа J различных компонент вектора kk по вектору rv

J:=0:

for i from 1 to n do

if rv[i]<>0 then J:=J+1 fi:

od:

Получили число J разных корней. Если J<n, то имеются кратные корни.

J;



Причем сумма всех кратностей rv[i] должна быть равна числу корней n, то есть





Операция логического сравнения

evalb(sr=n);



Упростим задачу представления корней и их кратностей.

Различные (без повторений) компоненты вектора kk перепишем в вектор k, а заодно

и значимые позиции rv[i]>0 также перепишем в вектор rr

k:=vector(J,[]):

по следующему алгоритму

ir:=1:

for j from 1 to n do

if rv[j]=0 then next else k[ir]:=kk[j]: rr[ir]:=rv[j]

fi: ir:=ir+1:

od:

r:=convert(rr,list):

Проведем визуализацию полученных таким образом корней k[i] и кратностей их

повторения r[i], перечисленных ниже слева направо, i=1,...,J

eval(k); eval(r);





Перейдем к вычислению решения СДУ, соответствующего данному корню. Запишем харак-теристическую матрицу M=A-kE для каждого корня k[i]:

for j from 1 to J do

M[j]:=evalm(A-k[j]*E)

od:

Визуализация характеристической матрицы

eval(M[1]):

Примечание. Покажем, что матрица, возведенная в степень, равную кратности корня

k[i], совпадающей с количеством корней (J=1, r[1]=n), равна нулю - только в этом случае!

if r[2]-n<>0 then print(`не тот случай`) fi:



evalm(M[2]^r[2]):

Для кратного корня решение СДУ связано с вычислением матрицы b, описание которой

приводится ниже. Введем матрицу b

for i from 1 to J do

b[i]:=matrix(n,r[i],[])

od:

evalm(b[1]): evalm(b[2]): col(b[2],2)[1]:

Для вычисления собственных векторов корней с единичной кратностью зададим

ноль-вектор правых частей СЛАУ

p0:=vector(n,[]):

for i from 1 to n do

p0[i]:=0

od:

Итак, в случае кратного корня решение надо искать, как указывалось ранее, в следующем

виде

(24)

Если, например, r[1]=2 , то из системы (24) надо взять только одно верхнее уравнение или

воспользоваться формулой

(25)




при том, что . (26)
Из равенства M*col(b,r)=0, учитывая ранг матрицы М, находим собственный вектор

col(b[i],r[i]), решая систему линейных уравнений. Однако, имеет смысл вычислить все

собственные векторы матрицы А (как решение СЛАУ таким образом linsolve(M[i],p0 )),

представляющие собой крайние правые столбцы матриц b[i] (в частном случае матрица

b[i] может иметь один столбец).

for i from 1 to J do

col(b[i],r[i]):=linsolve(M[i],p0) linsolve

od:

Визуализация полученного собственного вектора р=col(b,r) и других.

evalm(col(b[1],r[1])); evalm(col(b[2],r[2]));evalm(col(b[2],1)):





Определение собственного вектора проводится с точностью до константы, поэтому здесь

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

решение. Практика показала, что лучше всё-таки эту константу оставить.

А сейчас заменим системную переменную _t на переменную с другим именем,

например, g. Введем матрицу g

g:=matrix(J,n,[]):

Запишем алгоритм замены системной переменной _t

for i from 1 to J do

if r[i]=1 then

col(b[i],r[i]):=subs(seq(_t[j]=g[i,j],j=1..r[i]),col(b[i],r[i]))

else col(b[i],r[i]):=subs(seq(_t[j]=g[i,r[i]-j+1],j=1..r[i]-1),col(b[i],r[i])) fi:

od:

Получили собственные векторы, стоящие в одиночных и крайних правых столбцах

матриц b[i], которые можно сравнить с приведенными ранее kev.

evalm(col(b[1],r[1])): evalm(col(b[2],r[2])): evalm(col(b[3],r[3])):

Теперь на основании формулы (25) можно провести вычисление остальных вектор-

столбцов матриц b[i], у которых r[i]>1, решая соответствующие СЛАУ. Начать следует

с крайних левых вектор-столбцов р=col(b[i],1).

for i from 1 to J do

if r[i]>1 then linsolve

col(b[i],1):=linsolve(1/(r[i]-1)!*M[i]^(r[i]-1),col(b[i],r[i])) fi:

od:

Все крайние левые столбцы матриц b[i], r[i]>1 получены, в том числе собственные векторы.

col(b[1],1): col(b[2],1);



При решении СЛАУ все переменные делятся в зависимости от ранга матрицы М на

свободные и базисные. Свободные переменные будем записывать в матрицу g с компо-

нентами g[i,j], у которых номер строки i совпадает с номером матрицы b[i]. Причем

будем полагать, что число компонент g[i,j] в строке i не может превышать r[i], то есть

j=1,...,r[i]. По-прежнему в полученном выражении перейдем от системного параметра _ t

к параметру с именем g.

for i from 1 to J do

if r[i]>1 then

col(b[i],1):=subs(seq(_t[j]=g[i,j],j=1..n),col(b[i],1)) fi:

od:

Визуализация результата подстановки

col(b[3],1): col(b[3],2):

А дальше, если r[i]>2 , вычисляем все остальные столбцы col(b[i],j), j=2,..,(r[i]-1) как

функции 1-го столбца col(b[i],1) по формуле (25) , то есть - ещё (r[i]-2)- раз. Найдем вектор

р=col(b[i],j) из уравнения (25) для j=r-1, решая СЛАУ

for i from 1 to J do

for j from 2 to r[i]-1 do

if r[i]>2 then

col(b[i],j):=1/(j-1)!*evalm(M[i]^(j-1)&*col(b[i],1)):

col(b[i],j):=convert(col(b[i],j),vector) fi:

od:od:

evalm(col(b[2],2)[1]): J:

Таким образом матрица b получена. Остается только записать её компоненты, необходи-

мые для дальнейших вычислений.

for s from 1 to J do

for i from 1 to n do

for j from 1 to r[s] do

b[s][i,j]:=evalm(col(b[s],j)[i])

od:od:od:

Выборочная визуализация вычисленных компонент матрицы b

evalm(b[1]): evalm(b[2]): evalm(b[3]): b[2][1,1]:

Столбцы col(b,u), u=1...r матрицы b можно рассматривать как векторы, компоненты

которых зависят от параметров, определяемых из начальных условий . Введем вектор Т

for j from 1 to J do

if r[j]>1 then

T[j]:=vector(r[j],[seq(t^(u-1),u=1..r[j])]) fi:

od:

eval(T[2]); eval(T[1]):



Рассмотрим корень k[j], j=1,..,J - кратности . Соответствующее этому корню

решение СДУ ищется в виде вектора-столбца Q, умноженного на exp(k[j]*t) и некоторую

константу. Введем матрицу Q, столбцы которой привязаны к корню k[j]

Q:=matrix(n,J,[]):

Тогда

for j from 1 to J do

if r[j]>1 then col(Q,j):=evalm(b[j]&*T[j])

else col(Q,j):=evalm(b[j]):

col(Q,j):=convert(col(Q,j),vector) fi:

od:

type(Q,matrix); col(Q,2): col(Q,2)[2]: col(Q,3):



Таким образом, решение исходной СДУ получено в общем виде. Его можно записать как

взвешенную сумму вектор-столбцов col(Q,j), если каждый из них домножить на

C[j]*exp(k[j]*t). Константы в этом решении определяются из начальных условий. Для этого

подставим в матрицу Q начальное значение t=0, необходимое для вычисления её

параметров. Результат подстановки отразим в матрице Q0

Q0:=matrix(n,J,[]):

Причем столбцы матрицы Q0 также привязаны к корню k[j].

for j from 1 to J do

if r[j]>1 then col(Q0,j):=evalm(subs(t=0,col(Q,j)))

else col(Q0,j):=col(Q,j) fi:

od:

Визуализация матрицы Q0

evalm(col(Q0,1)): evalm(col(Q0,2)): evalm(col(Q0,3)):

Запишем матрицу AQ0[i] коэффициентов матрицы Q0[i], соответствующую корню k[i] и

представляющую i-ю матрицу СЛАУ

for i from 1 to J do

AQ0[i]:=matrix(n,r[i],[])

od:

evalm(AQ0[2]):

При её вычислении будем использовать то обстоятельство, что пара r[s] и col(Q0,s) позво-

ляет выделять (для r[s]>1) в col(Q0,s) коэффициенты при g[s,j]. Вначале найдем коэф-

фициенты в одиночных столбцах матриц b[s], r[s]=1 или col(Q0,s)

for s from 1 to J do

if r[s]=1 then col(AQ0[s],s):=evalm(col(Q0,s)/g[s,s]):

AQ0[s]:=convert(col(AQ0[s],s),matrix) fi:

od:

col(AQ0[1],1);



Ниже изменение индекса j до значения j=r[s], (r[s]>1) определяется кратностью корня k[j].

Крайний правый столбец с номером j=r[s] в матрице b[s] является собственным векто-

ром, для которого общая константа равна g[s,r[s]]. Следовательно, этот столбец можно

вычислять отдельно, например, так col(AQ0[s],r[s]):=col(b[s],r[s]). Однако, ниже в

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

for s from 1 to J do

for j from 1 to r[s] do

for i from 1 to n do

if r[s]>1 then AQ0[s][i,j]:=coeff(col(Q0,s)[i],g[s,j]) fi:

od:od:od:

rank(AQ0[2]): evalm(AQ0[2]): матрица СЛАУ и её ранг.

Полученные матрицы AQ0[s] теперь надо объединить в одну AS, которая и будет итоговой

матрицей СЛАУ: . Именно из этого уравнения мы найдем вектор констант Rh,

численные значения компонент которого подставим в матрицу g и получим окончательное

решение Х(t). Итак

AS:=concat(seq(AQ0[j],j=1..J)):

Rh:=vector(n,[]): rank(AS);



Тогда, исходя из начальных условий Х0, компоненты g[s,i], s=1,..,J, i=1,..,n матрицы Q (Q0)

определим как решение СЛАУ

Rh:=linsolve(AS,X0); Rh[4]: linsolve



Осуществим подстановку вычисленных компонент Rh[i] в матрицу Q. Напомним, что

в матрице Q столбцы привязаны к корню k[s], а сам корень к номеру s матрицы b[s],

s=1,..,J. Напомним, что матрица g имеет компоненты g[s,j], у которых номер строки s

совпадает с номером матрицы b[s]. Причем будем полагать, что число компонент g[s,j] в

строке s не может превышать r[s], то есть j=1,...,r[s].

for s from 1 to J do

for i from 1 to n do

if s=1 then

col(Q,s)[i]:=subs(seq(g[s,j]=Rh[j],j=1..r[s]),col(Q,s)[i])

else col(Q,s)[i]:=subs(seq(g[s,j]=Rh[(sum(r[m], m=1..s-1))+j],j=1..r[s]),col(Q,s)[i]) fi:

od:od:

Визуализация полученного результата. Выпишем все вектор-столбцы матрицы Q

evalm(col(Q,1)): evalm(col(Q,2)): evalm(col(Q,3)):

Так как вектор констант С уже содержится в g[s,j] и мы их определили в Rh[i], то дополни-

тельно умножать решение на эти константы нет надобости. Поэтому сразу запишем

решение (1) системы (2) в виде:



Где константы C[i] также определяются из начальных условий, а здесь они равны единице.

for s from 1 to J do

C[s]:=1

od:

Более наглядной является следующая покоординатная запись решения:

for i from 1 to n do

X[i]:=add(C[j]*exp(k[j]*t)*col(Q,j)[i], j=1..J)

od:

Выборочное представление полученного решения:

seq(x[j](t)=X[j],j=2..3);



Решения представим графически на рис. 1.

plot({seq(X[j],j=1..3)},t=0..8,color=[black,blue,red,green],title=`Решение СДУ`); `Рис.1`;





Проверка. Решение задачи штатным ПО Maple.
n:=6: X0:=vector(n,[1,0,1,0,1,0]):

A:=matrix(n,n,[0,1,0,0,0,0,0,0,1,0,0,0,0,0,0,1,0,0,0,0,0,0, 1,0,0,0,0,0,0,1,-12,-52,-91,-82,-40,-10]):

Вектор решений :

X:=vector(n,[seq(x[j](t),j=1..n)]):

По матрице А и вектору решений Х запишем саму систему дифференциальных уравнений

в векторно-матричной нормальной форме с указанием начальных условий:

seq(diff(x[j](t),t)=evalm(row(A,j)&*X),j=1..n): seq(x[j](0)=X0[j],j=1..n):

Для решения СДУ штатным ПО неоходимо записать её формульную постановку в

следующем виде

sys:=seq(diff(x[j](t),t)=evalm(row(A,j)&*X),j=1..n); fcns:={seq(x[j](t),j=1..n)}: simplify(dsolve({sys,seq(x[j](0)=X0[j],j=1..n)},fcns)):



Вычисление собственных чисел и соответствующих им собственных векторов.

p := [eigenvectors(A)]:

xm:=charmat(A,p[1][1]):

Свойства собственных векторов: умножение характеристической матрицы на собствен-

ный вектор должно дать нулевой вектор-столбец.

evalm(xm&*p[1][3][1]);



Ч.Т.Д.


Константин Титов стр. 5/10/2015
скачать файл



Смотрите также:
4. Алгоритм решения системы дифференциальных уравнений, имеющей действительные кратные корни характеристического уравнения
242.11kb.
Метод численного решения, дифференциальные уравнения, трубопроводные системы
899.52kb.
Специальность «Дифференциальные уравнения, динамические системы и оптимальное управление» область математики, посвященная изучению дифференциальных уравнений
2944.3kb.
При необходимости более детального просмотра увеличьте масштаб документа
209.85kb.
Академик Власов В. В. (Бибу, Балаково)
30.83kb.
Исследование машинного агрегата
89.93kb.
1. Делимость целых чисел. Деление с остатком. Алгоритм Евклида нахождения наибольшего общего делителя. Решение линейных уравнений в целых числах
12.21kb.
Тема «Решение однородных тригонометрических уравнений»
23.46kb.
Задачи управления для гидродинамических моделей механики сплошной среды
16.22kb.
Решение целых уравнений
49.27kb.
Программа вступительных испытаний в магистратуру по направлению 010100. 68 Математика Программа обсуждена на заседании кафедры ит
150.53kb.
Лабораторная работа №8 " Критерии устойчивости " студентка гр. Ас27 Позняк О. В. Проверил: Иванюк Д. С
15.93kb.