Аналитически-численный расчет динамики выделенного класса моделей нелинейных систем с распределенными параметрами Текст научной статьи по специальности «Математика»

Теория сигналов

УДК 681.511.4

А. В. Прикота

Санкт-Петербургский государственный электротехнический университет "ЛЭТИ"

II

Аналитически-численный расчет динамики выделенного класса моделей нелинейных систем с распределенными параметрами

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

Системы квазилинейных уравнений, аналитически-численный расчет, распределенные параметры

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

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

В данной статье предлагается метод получения приближенных аналитических решений систем квазилинейных дифференциальных уравнений в частных производных гиперболического типа в случае двух независимых переменных, основанный на аналитически-численном методе решения систем обыкновенных нелинейных дифференциальных уравнений, изложенном в [1], [2]. Суть аналитически-численного метода заключается в нахождении коэффициентов ряда Тейлора для неизвестных точных решений преобразованием исходной системы уравнений, определением допустимого шага расчета по координате и перенесением текущей точки расчета на величину шага.

© А . В. Прикота, 2003

17

Известия вузов России. Радиоэлектроника. 2003. Вып. 1======================================

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

.9 u тад u ПЛ

A— + В— = c, (1)

д t д x

где A = A(t, x,u) и В = B(t, x,u) - матрицы с размером n х n; u = [«i, щ,..., un]T - вектор неизвестных; t - временная переменная; x - пространственная переменная; c=c (t,x,u) -вектор с размером n х 1.

Если матрица A неособая, то система (1) преобразуется к нормальному виду:

8 u .8 u ,

— + Ai — = b, (2)

8t 8x

где A = Ai(t,x,u) и b = b(t,x,u) - некоторые новые матрица и вектор соответственно.

Если матрица А особая, то система (1) приводится к нормальному виду после преобразования координат [3].

Систему (2) можно записать в характеристической форме

1

f8 u „ д u ^

= fk, k =1, ••■, n (3)

д t д x

и в виде инвариантов Римана

дкдх=м, к=1, ..., п,

дх

п

где 1 к = 1 к (},х,и) - левый собственный вектор матрицы А1; /к = ^ ¡к ¡Ь - компоненты

I =0

вектора - свободного члена Г; к = к , х, и) - соответствующее 1 к собственное значение матрицы А1; Гк = Гк (¿,х,и) - компоненты вектора инвариантов Римана; gк (¿,х,г) - компоненты вектора - свободного члена g.

Системы с п > 2 не всегда можно привести к инвариантному виду, но такое приведение всегда возможно для систем с п < 2 и для продолженных систем [4], [5].

В узком смысле система (2) называется гиперболической в односвязной области переменных (¿,х,и), если в каждой точке этой области собственные значения 2, . ., ^п матрицы А1 вещественны и различны [5].

Для выделения из бесконечного множества решений системы (2) единственного необходимо задать начальные условия и, если это необходимо, граничные. Доказано, что если начальные и/или граничные условия являются аналитическими функциями в окрестности точки (¿о, хо) и элементы матрицы А1 определены и аналитичны в окрестности точки

===================================== Известия вузов России. Радиоэлектроника. 2003. Вып. 1

(to, xo,u(t,x)), то система (2) имеет аналитическое решение в некоторой окрестности точки (¿o, xo), притом единственное в классе аналитических функций [6].

Таким образом, задача может быть поставлена следующим образом: для системы (2) и заданных аналитических начальных и/или граничных условий uk (t,o), uk (o, x),

(k = 1, ..., n) в области D (t >o, x>o) необходимо найти приближенные решения Uk(t,x), (k = 1, ..., n) в классе аналитических функций.

Нахождение аналитического решения. Запишем систему (2) в следующем виде:

A(D¿ ,D x )u(t, x)= G(D¿, D x )F(t, x) + H(t, x, u, F), (4)

где D¿,Dx - операторы частного дифференцирования по t и x соответственно; A(d¿ ,Dx) и G(d¿,Dx) - матрицы с размерами n x n и n x m соответственно, элементы которых представляют собой операторы частного дифференцирования или константы; F(t, x) - вектор внешних воздействий с размером 1 х m. В вектор H(t, x, u, F) попадают все оставшиеся члены уравнений, т. е. те, которые не попали в матрицы A(d¿,Dx) и G(Dt ,D x ).

Решение системы (2) ищется в области t >o, x >o; начальные условия Uk (o, x) = fk (x), k = 1, ..., n ; граничные условия щ (t,o) = gk (t), k = 1, ..., n [3]-[5].

На рис. 1 показана область существования решений. В подобласти 1 основное влияние на результат оказывают граничные условия, поэтому решение u(t,x) определяется в виде ряда Тейлора по переменной x:

/ \ ^^ Rk i (t) x ^ Uk\t,x)-У —--, k = 1, ..., n . В такой

zto i!

же форме представляются и все элементы H(t, x, u, F), зависящие от x. В области 2, напротив, основное влияние на результат оказывают начальные условия и u(t,x) определяется в виде ряда Тейлора по переменной t:

Í Ч Rk i (x) t Uk\t,x) => —!-, k = 1, ..., n . В такой

i=o i! ___

же форме представляются и все элементы O ^x

H(t, x, u, F), зависящие от t. Рис. 1

Подставив все ряды в матрицу H(t,x,u,F) и произведя все требуемые алгебраические действия над ними, получим следующие выражения: для области 1

strkH(t,x,u,F) = Y Tk,i(t)x , k = 1, ..., n ; (5а)

a i!

i=o

для области 2

strk H(t,x, u, F) = ¿

k = 1, ..., n .

(5б)

i=0

Преобразуем систему (4) с учетом (5а) и (5б) с использованием двойного преобразования Лапласа (^ p, Бх ^ ц) :

A(p, ц )u(p, ц)=G(p, ц) ^ ^ ц) + Н^ ц) + Q(p, ц),

где Q(p,q) - вектор, содержащий начальные и граничные условия. Запишем решение по правилу Крамера:

(6)

для области 1

для области 2

щ (p,q)=AApi, А(p,q) * о, k=1,

^-1

X В (РP)q-

uk(p, q)=-N-

X Ai (p )q-

i=0

N-1

X Вi (q ) p

uk(p, q)=-N-

X A- (q) pl

i=0

n

(7а)

(7б)

Предел суммы N в выражениях для щ (p, ц) ограничен конечным числом в связи с заменой рядов Тейлора, входящих в (6), полиномами и определяется заданной точностью расчета.

После выполнения деления в выражениях (7а), (7б) получим: для области 1

( ч

uk(p,ч)= Х )+1 ;

i=о q

для области 2

( ) (ц)

uк\P,ц) = Х, ■

i =0 Р

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

в области 1

uk (t,x) = , k = 1, -,

i=o i!

n;

(8а)

в области 2

uk(t,x)=E

Rk,i (x) ti

k = 1,

n.

i=o

i!

(8б)

Коэффициенты Я^ (х), Я^ (V) определяются рекуррентно. Если они определяются единственным образом, то полученное решение щ(V,х), (к = 1, ...,п) является единственным решением системы (2).

Выбор величины шага расчета. Критерием при выборе величины шагов расчета по временной переменной и по пространственной переменной Их является допустимая

погрешность расчета как локальная (погрешность на данном шаге), так и накопленная. Обе эти погрешности, очевидно, зависят от порядка полинома Тейлора N и от величины шага. Чем меньше шаг расчета, тем быстрее убывают члены рядов (8а), (8б) и, следовательно, тем меньше локальная погрешность, связанная с отбрасыванием всех членов ряда со степенью выше N.

Фундаментальным понятием в теории уравнений гиперболического типа является понятие характеристик. Уравнения характеристик вытекают из характеристической формы (3):

= £,кх,и), к = 1, ..., п . (9)

Система из п уравнений гиперболического типа имеет п различных характеристик. Характеристики системы непосредственно связаны с волновым характером движения, которое возникает в системах, описываемых этими уравнениями. Они указывают на конечную скорость распространения начальных или граничных возмущений вдоль координат. В области 1 (см. рис. 1) выбирается шаг по пространственной координате х = Их.

Для того чтобы узнать, в какой момент времени ^ = И1 граничное возмущение достигло

прямой х = Их, необходимо из (9) найти уравнение характеристик, выходящих из точки О

= 0, х = 0) в рассматриваемую область (V > 0, х > 0). Так как вектор и(^, х) уже найден,

то можно найти уравнение характеристик [5] и найти кг, соответствующий шагу Их. Для

области 2 (см. рис. 1) все делается аналогично с заменой друг на друга переменных V и х, и Их соответственно.

Рис. 2 поясняет процедуру выбора шага расчета и определения текущей точки расчета. На рис. 2 изображена одна характеристика, направленная в область (V > 0, х > 0). В случае, когда таких характеристик несколько, необходимо обращаться к системе урав-

t t

x

O

Рис. 2

нений, записанной в инвариантах Римана, и по ним определять, какие решения uk (t, x) двигаются по той или иной характеристике. Значение ht для области 1 находится как h t = min (h j ),

где hj - момент времени, соответствующий hx, на j-й характеристике, направленной в рассматриваемую область; j = 1,..., m, m < n . Для области 2 все аналогично, но с заменой ht на

hx.

Таким образом, решение uk (t, x) на прямой x = hx имеет вид

NRk,, (t) h\

uk

(t, hx) =

i=0

i!

t > ht.

Решение uk (t, x) на прямой x = ht имеет вид

uk (ht, x) = x

NRki (x)hlt

i=0

i!

при x > hx .

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

и Ык (г, hx) более простыми функциями иапп к (ь г, х), иапп к (г^ х) (например, тригонометрическими, степенными полиномами и др.).

После этого аппроксимирующие функции Ыапп к ^ г, х^ иапп к

) становятся начальными и граничными условиями на следующем шаге расчета, а начало координат 0(0,0) переносится в точку 0'(^, ^).

Пример 1. В качестве тестового примера рассмотрим уравнение

(д и/д г) + и (д и/д г) = 0.

Начальное условие u (0, x ) = f (x ) =

sin(x), 0 < x< л ;

[0 - иначе.

Это уравнение имеет известное точное решение в виде неявной функции и = / (х - ги) [5].

Нахождение аналитического решения. Имеем д и/д г = -и(д и/д х). Здесь А(Бг ,Бх ) = Бг, А(Бг ,Бх ) = 0, Н (г,х, и, Г) = -и (д и/д х).

22

Р ( ч ^ R (x)ti дu ^ R'(x)ti

Решение ищем в виде ряда u[t,x) = X—, тогда — = У w

п i! д x „ i!

i=0 i=0

, , ч ® R (x) ti ® Ri (x)ti ® T (x)ti H (t,x, u, F ) = -u ( д u¡ д x) = - X , = E"

i! i! i!

i=0 i=0 i=0

Преобразуем исходное уравнение по Лапласу с учетом выражения для н(/,х, и, Г):

да - ( )

О ^ Р,°х ^ q, ри(р, q) = и(0, q) + X.

I=0 Р 1

В ( ) ( ) и(0, q) Т 0 ^) Т 1(q) % (q)

Выражаем м(р, ^ через м(р, ^ = 4 ' + 0" + 1 + • - + ^ ' + • • ■.

Р р2 р3 pN+2

После обратного преобразования Лапласа получим искомое решение:

( ) (п ) T0(x)t T1{x)t2 TN{x)tN+1 <-)<i

<(t,x)=u (0,x)++- + и/л, + -=E

1! 2! (N +1)!

i=0

i!

где R (x) = u(0,x), Ri(x) = -0 (x) = R (x)R (x), R2 (x) = -i(x) = R (x)R (x) + Ri(x)R (x), ..., т. е.

Rq (x) = sin(x), Rl(x) =-COs(x)sin(x), R2 (x) = sin(x)sin2 (x) + cos2 (x)] + cos2 (x)sin(x) при 0 < x < n; Rq (x) = Ri (x) = R2 (x) = 0 при других значениях x.

Из аналитичности начального условия u (0, x) и вида нелинейного члена в H(t,x,u,F) следует, что все коэффициенты Ri (x), i = q, ..., да являются аналитическими функциями аргумента x.

Выбор шага расчета. Уравнение характеристик dx¡dt = u (t, x). Следовательно, характеристика, выходящая из начала координат, имеет вид x = u(0,0) t. Выберем шаг по t равным ht. Тогда шаг по x имеет вид hx = u (0,0) h t = 0 . Следовательно, начальная точка следующего шага расчета есть o(ht ,о).

/ ч NRi (х) ht

Решение u\ht, х) = > -;-аппроксимируется интерполяционными полиномами.

i=0 i!

На рис. 3 приведены графики приближенного и точного решений для разных моментов времени. Приближенное решение (кривые 1) получено аналитически-численным методом, точное (кривые 2) - с помощью системы компьютерной алгебры Maple 6. На рис. 3, а кривые 1 и 2 совпадают.

В процессе расчета шаг ht изменялся в пределах от 0.1 до 0.5 с; порядок полинома

N выбирался равным 6 или 7.

Известия вузов России. Радиоэлектроника. 2003. Вып. 1 и, о.е. _____u, о.е.

0.8 0.6 0.4 0.2

0.8 0.6 0.4 0.2

0 0.5 1.0 1.5 2.0 2.5 3.0 x, о.е. 0 0.5 1.0 1.5 2.0 2.5 3.0 x ое а б

Рис. 3

Пример 2. Исследуемое устройство описывается уравнениями

Г- (du/аx) = Сд (di/dt);!

I- (di/dx) = L 0 (du/dt)

Сд = Сд {и).

Эта система уравнений возникает при рассмотрении процессов в коронирующей длинной линии [7]. Здесь и(г, х), ¡(г, х) - напряжение и ток в линии; Сд, Ц - погонные

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

Пусть в момент времени г = 0 обесточенная линия подключается к источнику напряжения, находящемуся в точке х = 0. Тогда начальные условия и(0, х) = 0, ¡(0, х) = 0;

, Я1 I-

граничные условия и(г,0) = gl(г), ¡(г,0) = g2 (г) = (1Д/Ц0) [д/Сд (и) ди .

0

При таких начальных и граничных условиях данная система уравнений имеет известное точное решение

x) = g1[t - ^]10Сд(u)], i(t, x) = (Ул/^0 Сд(u) du .

0

Нахождение аналитического решения. Аппроксимируем нелинейность: Сд(u) =

= С0 + Qu . Запишем исходное уравнение в виде (4):

A ( Dt ,D x ) u (t,x ) = G ( Dt ,D x ) F (t,x ) + H (t,x, u, F )

H (t,x, u, F ) =

- 2J u,

^-Qu Dt U 0

где u (t,x) = [ub U2 ]T = [u, i]T, A (D¿,Dx ) =

f D x L0D Л

V QD,

0t D

, G(Dt,Dx) = 0,

x

V

У

/ \ ^ Rki (t) x

Решение ищем в виде ряда u^[t,x) = ^ —'—-

, k = 1, 2, тогда

i =0

D u = £

Ri (t) x

i=o

i!

H (t,x, u, F) = -Qu Dtu = -Q Е Е

R i (t) xi ® R i (t) xi = ® Txi (t) x

i ! - i !

i=o i=o

= Е-

i=o

i!

Преобразуем полученные уравнения по Лапласу, с учетом выражения для Н (/, х, и, Г):

Dt ^ p, Dx ^ q, A (p, q) u (p, q) = H (p, q) + Q (p, q),

где A (p, q) =

q Чр

cop q

Выразим u(p,q) и i (p,q) :

H (p, q ) =

Y TmM

^ J+1

i=o q o

, Q(p,q) = [u(p,o), i(p,o)]T.

u

(p,q)=

í (p,q)

® T ( p)

qu ( p,o ) - Lo pi ( p,o ) + Lo p X ^^

_i=o q

2 2 2 q - c p

® T ( p)

qi(p,o) -Copu(p,o)- qX ^^

_i=o q

2 2 2 , q - c p

где с = Ч0С0.

Ограничим ряд, входящий в эти выражения, степенью N. С помощью очевидных преобразований получим:

N-1

x В (р) 4

щ(р,q)=-, к = 1, 2.

x А (р )4

I=0

В результате деления имеем:

и (р 4) = и (^0) - Ь0рр (Р,0) + ь0р-0 (р) + с2р1™ (Р,{°) + ь0р— (р) + с2р2Ч0рр (Р, 0) +

q

q

з

q

4

q

i (p,q) =

i (p, o) Co pu (p, o) + To (p ) -Ti (p) + c p i (p, o)

+

q q q

+ -T2 (p) - c2p2 [Cop u (p, o) + To (p)] +

q

Известия вузов России. Радиоэлектроника. 2003. Вып. 1======================================

После обратных преобразований Лапласа получим искомые решения для параметров

u и i

u

(t,x ) = u (t,0 ) +

+

L0 Dt i (t,0) x

1!

„2^2

- +

L)DtT0 (t) + c2D2 u (t,0)

x2

2!

- +

L0DtT1 (t) + c D2 L0 Dt i (t,0 )

x

3!

-+...=Z

i =0

N Jii (t) xi i!

i (t,x ) = i (t,0 ) + J

[C0 Dt u (t,0) + T0 (t)] x +

1!

2!

- +

+

-T2 (t)-c2Dt2 [CQDtu(tfi) + T0 (t)]}x3 + = ZZ R2i (t)xi

3!

i =0

i!

Из аналитичности граничных условий и(г,0), ¡(г,0) и вида нелинейного члена в матрице Н(г,х,и,Г) следует, что все коэффициентыК^-(г), к = 1, 2; г = 0, ..., да , являются

аналитическими функциями аргумента г .

Выбор шага расчета. Для получения уравнений характеристик приведем исходное уравнение к характеристическому виду:

дг/дх + ^Ц0 Сд(и)(дг/дг) + ^Сд/Ь0 ди/дх + ^Ь0 Сд(и)(ди/дг)] = 0;

дг/дхЦ0Сд(и)(дг/дг)Сд/Ь0 ди/дх -Л\Ь0 Сд(и)(ди/дг)] = 0. Отсюда уравнение характеристики, направленной в рассматриваемую область:

{ё г/ёх) = Л1 ¿0 Сд{и) = л!¿0Со + С\и).

Следовательно, характеристика, выходящая из начала координат имеет вид г - л!¿0[С0 + С1и(0,0)] х. Шаг по х выбирается равным hx. Тогда шаг по г имеет вид

ht = ^ ¿0[С0 + С1и(0,0)] hx, а начальная точка следующего шага расчета есть 0(кг,hx). / ч (г) ^х

Решения ик\г,hx]- ^—'—-, к = 1, 2, при г > ht аппроксимируются тригономет-

i!

рическим полиномом с использованием для этого быстрого преобразования Фурье. Задачи решали при следующих численных значениях параметров системы:

C0 = 8.46998-10-12 [Ф/м] ; С1

0.104-10-15 [ф/(В • м)], u > 0; -0.104•Ю-15 [ф/(В• м)], u < 0

L0 = 0.131182 •Ю-5 [Гн/м];

u

(t,0) = g1(t) = 240 000 [1 + cos(314t)] [b] .

2

x

Рис. 4 Рис. 5

В процессе расчета шаг Их уменьшался от 280 км в начале расчета до 28 км в его

конце; порядок полинома N варьировался от 10 до 8.

Графики точных (кривая 1) и приближенных (кривая 2) решений и((,х) и х) для различных х приведены на рис. 4 и 5 соответственно. Здесь ^' - момент времени, в который сигнал достигает точки х. Так как решения являются периодическими функциями аргумента I, то графики их представлены только для первого периода. На рис. 4, а и 5, а графики точных и приближенных решений совпадают.

Таким образом, рассмотренный метод получения приблизительных аналитических решений систем квазилинейных уравнений в частных производных гиперболического типа расширяет возможности аналитически-численного метода расчета динамики нелинейных систем [1], [2], снимая наложенные на него ограничения в виде требования сосредоточенности параметров системы.

Библиографический список

1. Бычков Ю. А. Аналитически-численный расчет динамики нелинейных систем / ИПЦ СПбГЭТУ "ЛЭТИ". СПб., 1997. 368 с.

2. Бычков Ю. А., Щербаков С. В. Аналитически-численный метод расчета динамических систем. СПб.: Энергоатомиздат, 2001. 344 с.

3. Уизем Дж. Линейные и нелинейные волны. М.: Мир, 1977. 622 с.

Известия вузов России. Радиоэлектроника. 2003. Вып. 1======================================

4. Рождественский Б. Л., Яненко Н. Н. Системы квазилинейных уравнений и их приложения к газовой динамике. М.: Наука, 1978. 687 с.

5. Рождественский Б. Л., Яненко Н. Н. Системы квазилинейных уравнений и их приложения к газовой динамике. М.: Наука, 1968. 592 с.

6. Масленникова В. Н. Дифференциальные уравнения в частных производных. М.: Изд-во РУДН, 1997. 447 с.

7. Караев Р. И. Переходные процессы в линиях большой протяженности. М.: Энергия, 1978. 191 с.

A. V. Prikota

The Saint-Petersburg state electrotechnical university "LETI"

Analytical-Numerical Account of Dynamics of the Chosen Class of Nonlinear Distributed Parameter Systems Models

The account method of dynamics of nonlinear distributed parameter systems models circumscribed by systems of quasilinear differential partial equations of a hyperbolic type is offered. The method is based on an analytical-numerical account method of dynamic systems, and is its addition.

Systems of the quasilinear equations, analytical-numerical account, distributed parameters

Статья поступила в редакцию 21 октября 2002 г.