Метод решения неоднородного дифференциального уравнения n-го порядка для задач схемотехнического проектирования электронных устройств Текст научной статьи по специальности «Физика»

А.Б. ШЕИН

МЕТОД РЕШЕНИЯ НЕОДНОРОДНОГО ДИФФЕРЕНЦИАЛЬНОГО УРАВНЕНИЯ N-ГО ПОРЯДКА ДЛЯ ЗАДАЧ СХЕМОТЕХНИЧЕСКОГО ПРОЕКТИРОВАНИЯ ЭЛЕКТРОННЫХ УСТРОЙСТВ

Ключевые слова: метод, проектирование, неоднородное дифференциальное уравнение, начальные условия, временной шаг, выбор, решение.

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

A.B. SHEYIN

THE METHOD OF THE DECISION OF THE N ORDER NON-UNIFORM DIFFERENTIAL EQUATION FOR PROBLEMS SCHEMATIC DESIGNING OF ELECTRONIC DEVICES

Key words: method, designing, non-uniform differential equation, entry conditions, time step, choice, decision.

It is offered a method of the decision of the n order non-uniform differential equations, based on replacement of the derivatives containing in the equation, them it is finite-difference relations.

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

( dn dn-1 dn-2 d2 d Л

an---------+ an ,--------------------------------------- + an 2- +... + a2—- + a,-+ a0

K dtn dtn-1 dtn-2 2 dt2 dt \

x(t ) =

^ dm L dm-1 L dm-2 ,d2 ,d L ^ /ч

bm----------+ bm 1--------^ + bm 2 ---------T" + ... + 2 ----7* + b,------+ bo v(t ),

v m dtm dtm-1 2 dtm-2 2 dt2 1 dt 0 J W

(1)

где х() - переменная состояния устройства во времени, v(t) - входное воздействие на объект, заданное в виде какой-либо функции времени; а0, а\, а2,...,ап_1, ап и Ь0, Ь1, Ъ2, ..., Ът-1, Ът - коэффициенты, определяемые структурой и параметрами компонентов устройства (обычно это комбинация постоянных времени или коэффициентов демпфирования отдельных цепей в схемах замещения реального устройства) и видом входного воздействия на объект.

Решение неоднородного дифференциального уравнения (1) наиболее просто находится из его линейной формы представления, которую можно получить заменой производных, содержащихся в уравнении, их конечноразностными отношениями [2]:

dx(t) _ х((к + \)к) - х((к + 0)к) dt к ......................... (2)

^(- 1)Г (—Т7кпх(к + п - г)к)

dt г=о (п - г).г\к

Для входного воздействия v(^) имеем:

) _ v((k + 1)к)- v((k + 0)к)

& к ............................ (3)

Ми- _ ^(- 1)г -—7“кгхкк+п - г)к)•

& г=о (п - г)!г!к где к - шаг дискретности переходного процесса, равный интервалу между моментами времени ґі+1 = (к + 1)к и ґі = кк : ґі+1 - ґі = (к + 1)к - кк = к, і и к = 0, 1, 2, ....

Так как все производные от х(ґ) и у(ґ) выражены через конечноразностные отношения (2) и (3), то подстановка их в неоднородное дифференциальное уравнение (1) позволяет представить его в следующем виде:

ап £ (- 1)г ---+ п - г)к) +

г=0 -п - г)г!к

+ап-1 £:(- 1)г-—,(п -1Х—+п-1 -г)к)+

г=0 (п -1 - г).г!кп 1

+ ап - 2 £ (- 1)г ----2(П )кп-2 Х(к + П - 2 - г )к) +

г=0 -п - 2 - г)!г!кп 2

+ а1 £ — 1)г к-г1)!г!к х((к +1 -г)к) + а00Х— + 0)к) =

= Ьт £ (- 1)г ----”) ^{к + т - г)к)+

г=0 -т - г}.г!кт

+ ьт-1 £ (- 1)г -------т -1 V-- + т -1 - г)к) +

г=0 (т -1 - г)!г!кт 1

+ Ьт-2 I? (- 1)г -----2т ) -2 + т - 2 - г)к) +

(т - 2 - г)г!кт 2

т-

г=0

+ Ъ2 ^ (- 1)Г ( 2) ,,2 ^ + 2 - Г )к) + г=0 (2 - г)г.к

+Ъ\ е (- 1)г кТ^Тк^+1 - г ^к ^+Ъ0 ^(к+0)). (4)

Раскрывая уравнение (4) и группируя получившиеся слагаемые, содержащие х(к + п - г)) и v((k + т - г)к) при г = 0,1, 2, 3, ..., п -1, п , получаем равенство, из которого можно выразить искомую переменную состояния х(к + п)) посредством ее значений и значений функции входного воздействия на всех шагах к сетки решения. Опуская достаточно громоздкие матема-

тические выкладки, сразу запишем формулу для нахождения значений переменной состояния х((к + п )к ):

х(к + п )к ) = -£ X (- 1)г —(п, /У)'—ч—п-^ ку х(к + п - г )к)+

г=^=0 (п - г)(г -V) ап (5)

+ ^к(п - т)]Г £ (- 1)г-V ( ^т.-(г')! ) Ът-^к"у(к + т - г)к).

ап г=0v=0 (т - г )(г -^! Ът

Как правило, для решения неоднородного дифференциального уравнения (1) задаются начальные значения как самой переменной состояния, так и ее

(Л (0) &х(0) (1)

производных до (п -1) -го порядка включительно: х(0) = х0 , —^- = х0 ,

&

С 2 х(0) ( 2) &п-2 х(0) ( п - 2) &п-1х (0) (п - 1) „

—= х0 % •••, —= х0 ^ -У- = х0 . При решении урав-

С2 Лп 2 С1п-1 0

нения (1) по формуле (5) необходимо выполнить пересчет этих начальных условий. Например, в соответствии с первым из конечно-разностных отношений (2), записанным для к = 0, имеем равенство

Сх(0) (1) х(+1))-х(0) // \ \ (0) (1),

—— = х0 =--------------, из которого находим х(0 + 1)) = х0 + х0 к .

& к

Продолжив процесс определения начальных условий, необходимых для реализации формулы (5), получим:

х((0 + 2)к) = х00^ + 2х0к + х02^к2,

х(0 + 3)к ) = х00 ^ + 3х0к + 3х02 ^ к2 + х03) к3,

х(0 + п)к)= £ ( п! V х0г}кг, п = 1, 2,..., п -1. (6)

г=0( - г)г!

Формула (6) используется также для нахождения значений функции входного воздействия у((0 + ц)к) в формуле (5):

4(0 + ц)к)= £( ^V ,^0г}кг , ц = 1, 2, ..., т . (7)

г=0(1- г )г!

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

Решение уравнения (1) по формуле (5) дает жесткое объединение частей траектории движения переменной состояния х(^) = х ((к + п )к ) в дискретные моменты времени ^ =(к + п - г)к на сетке решения этого уравнения и позволяет получать устойчивые и точные схемы расчета переходных процессов в электронных устройствах при правильно выбранном временном шаге к.

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

££ (- •Гг-ТТ' ^ к = 0. (8)

г=1 v=о (п - г)(г -V)! ап

Если раскрыть двойной ряд приведенного уравнения, то можно получить формулу для определения критического значения шага ккр, превышение которого приводит к численной неустойчивости вычислительного процесса по формуле (5):

1

. (9)

ккр =

ап

V а0 У

Формула (5) пригодна и для решения уравнения (1) в случае устройств с переменными по времени параметрами, так как каждый параметр того или иного компонента схемы принимает присущее ему конкретное значение для рассматриваемого момента времени развития процесса и для этого момента времени исходное уравнение (1) решается по формуле (5) с постоянными коэффициентами. Иными словами, функция изменения параметров компонентов схемы устройства определяется для рассматриваемого момента времени предварительно, т.е. решение уравнения (1) по формуле (5) выполняется с фиксированными для этого момента времени значениями коэффициентов.

Из формулы (5) видно, что она справедлива для любых значений п и т , что в корне меняет сложившееся представление о всегда присутствующем ограничении для выполнения равенства (1) в виде условия п > т, т.е. можно рассчитывать переходные процессы даже в гипотетических объектах, когда т > п .

Пример 1. Требуется решить неоднородное дифференциальное уравне-

С3/, С2 Л С/1 Я0 Я

ние а3—— + а2—— + а1-----------+ а0 = Ъ0 V, где а3 = 1, а 2 =--------1—,

С С & Ь0 Ь

Я0 Я 1 1 Я0 1 Я 1 1

а1 ----------------------------------------1-1-, а0 =-1-, Ъ0 =-, описывающее пе-

1 Ь0 Ь Ь0С ЬС 0 Ь0 ЬС Ь Ь0С 0 Ь0ЬС

реходный процесс в схеме параллельного инвертора [3] для переменной состояния (х), при следующих начальных условиях: (0 ) = /1(00 ^= 0,

(0)= .(1) = Е С2(0) = л(2)= Я0 Е

" - Л Л - , ~ 1

Ю ~ Т ’ 10

00

3 .

С 10 Ь0 ’ сх2 Ь0 Ь0

Исходные данные для расчета: Я0 -1-10 3 Ом, Я -1-10 1 Ом,

Ь0 -1-10-4Гн, Ь -1-10-2Гн, С -1-10-4 Ф, Е - 1-102 В.

Имеем: а3 -1, а2 -20, а1 -1,01-108, а0 -1,01-109, Ъ0 - 1-1010.

Так как п - 3 и т - 0, то согласно формуле (5) можно записать выражение:

(

/1 (к + 3)к) =

Г

1 -

а3 а

\

3 —^ И

аз у

(к + 2)И) -

3 - 2^ И + -1 И

1 7„2

'1 (к + 1)И)

3 У

+

1 - 02 И + 02 И 2 - 00 И3

Л

\

3 У

^ (к + 0)И )+—0- И3 у(к + 0)И ).

1 а3

Максимальный шаг сетки решения И определяем по формуле (9):

И=

1

^ а3 Л3 V ао У

1

Л 3

V 1,01-10 у

= (9,90099 -10-10) = 9,96688 -10-

с.

Принимаем шаг И = 1-10- с .

Полагая к = 0 , по формуле (6) находим начальные условия для реализации формулы (10):

*1(0 + 1)И ) =

.(0) .(1) Е 1-102

*10 + 71П И = — И =----------------

10 10 Т0 1-10- 4

-1-10- 4 = 1-102Л,

/1 (0 + 2)И) = /40} + 2/1(01)И + }И2 = 2ЕИ - ^ЕИ2 =

Т0 Т0 Т0

Е_

ЬП

(

Я0

\

2- 0 И

Т0 У

(

= 1-102

2 - 1-10 4 -1-10- 4 1-10 - 4

Л

= 199,9Л.

Подставляя эти начальные условия в формулу (10), находим:

( „ Л ( „ „ Л

,(0 + 3)И ) =

\

3 - 02 и

а3 у

,(0 + 2) )-

3 - 2-^И + -2И2

а3 а3 у

(0 + 1)И)-

+

^ - о!и + 01. И’- - а!и3\

V а3

а

а

1 (0 + 0)И) + -0-ИЗу(0 + 0)И) = 199,7001Л .

3 У

Если для решения задачи задаются все начальные условия, включая п-ю производную от переменной состояния /1 ():

с 3/1(0) = * (3) =

с3 10

*0 *0

1

\

т т тс

V ^0 ^0 ^0^ У

тогда этап определения 11 (0 + 3)И ) , приведенный выше, можно опустить, а вместо него для нахождения значения /1 ((0 + 3)И) нужно использовать формулу (6):

/1 (0 + 3) ) = 4°} + 3/1(01) И + 3/1(02} И2 + /1(03) И3 = 3 ЕИ - 3^ ЕИ2 +

Т0

Т0 Т0

+

1

*0 *0

1-Еи > Е , - 3 И +

= — И 3

У Т0 Т0 ч

К0 К0

1

= 199,7001Л.

В результате, применяя формулу (10), можно продолжить решение исходного неоднородного дифференциального уравнения, полагая к = 1,2,3,..., вплоть до установления переходного процесса в инверторе. При к = 1 имеем:

а

3

1

И

3

i1 (1 + 3)h) = 2,998 • i1 ((1 + 2))- 4,006 • i1 ((l + l)h)+2,00699 • i1 ((l + 0)h)+

+110 -2 v((l + 0 )h) = 2,998 • 199,7001 - 4,006 • 199,9 + 2,00699 • 100 +

+110-2 • 100 = 598,7009 - 800,7994 + 200,699 +1 = -1,3995 A.

Для k = 2 получим:

i1 (2 + 3)) = 2,998 • i1 ((2 + 2)) - 4,006 • i1 ((2 + 1)h) + 2,00699 • i1 ((2 + 0)h)+

+110-2 v((2 + 0)h) = 2,998 • (-1,3995)- 4,006 • 199,7001 + 2,00699 • 199,9 +1 =

= -4,195701 - 799,9986 + 401,1973 +1 = -401,997A.

Из результатов расчета видно, что при принятом для расчета шаге h получен расходящийся по времени процесс решения неоднородного дифференциального уравнения и необходимо уменьшить шаг сетки решения h.

Принимаем шаг h = 1-10-5 с .

Полагая k = 0 , по формуле (6) вновь находим начальные условия:

i1 (0 + 0)) = 4°} = 0A ,

i1 (0 +1)) = 4° ^ + 4^ h = 10 A ,

i1 (0 + 2)) = 40) + 2i1(01)h + i1(02} h2 = 19,999 A ,

i1 (0 + 3)h ) = 4°} + 3/\(01} h + 342} h2 + 43) h3 = 29,897 A.

Подставляя начальные условия в формулу (10) и делая первый шаг, т.е. присваивая k значение, равное 1 (k = 1), находим:

i1 (1 + 3)h) = (3 - a2h)/\ ((1 + 2)h)- (- 3 - 2a2 h + a1h2 ((1 + 1)h)+

+ (1 - a2h + a1h2 - a0h3 )i1 ((1 + 0)h)+b0h3v((1 + 0)h) = 2,9998 • 29,897 -

- 3,0097 • 19,999 +1,0098989 • 10 +1 • 10-5 • 1 • 102 =

= 89,68502 - 60,19099 +10,098989 + 0,001 = 39,594019 .

Для k = 2 получим: i1 (2 + 3)) = 2,9998^ i1 ((2 + 2)h)- 3,0097^ i1 ((2 + 1)h)+1,0098989- i1 ((2 + 0)h)+0,001= = 2,9998 • 39,594019 - 3,0097 • 29,897 +1,0098989 • 19,999 + 0,001 =

= 118,77413 - 89,981 + 20,196968 + 0,001 = 48,991098 .

Значения i1 ((k + 3)h) для k = 3, 4, 5, ... представлены в табл. 1.

Для проверки правильности выбора шага сетки решения h = 1 • 10 5 с, уменьшим значение этого шага еще на порядок. Если результат расчета на первом шаге при h = 1 • 10-5 с совпадет с результатом расчета для первых десяти шагов с шагом h = 1 • 10-6, то расчет с шагом h = 1 • 10-5 выполнен верно.

Результаты расчета i1 ((k + n )h) для n = 1,3 и k = 1,7 с шагом h = 1 • 10-6 с приведены в табл. 2.

к 11 (к + 3)) к 11 ((к + 3)И)

3 *1 ((3 + 3)И) 57,991324 26 ((26 + 3)И) 29,543458

4 *1 ((4 + 3)И) 66,500822 27 ((27 + 3)И) 18,520542

5 *1 (5 + 3)) 74,429733 28 \ ((28 + 3)И) 7,229091

6 *1(6 + 3)) 81,693163 29 ((29 + 3)И) -4,218542

7 *1 ((7 + 3)И) 88,212090 30 ((30 + 3)И) -15,707302

8 *1 ((8 + 3)И) 93,914220 31 *1 ((31 + 3)И) -27,120567

9 *1 ((9 + 3)И) 98,734785 32 \ ((32 + 3)И) -38,342310

10 (10 + 3)И) 102,617273 33 ((33 + 3)И) -49,256277

11 *1 (11 + 3)И) 105,514080 34 *1 ((34 + 3)И) -59,74816

12 (12 + 3)И) 107,387081 35 *1 ((35 + 3)И) -69,706776

13 *1 (13 + 3)И) 108,208110 36 ((36 + 3)и) -79,02521

14 *1 (14 + 3)И) 107,959344 37 *1 ((37 + 3)И) -87,601941

15 *1 ((15 + 3)И) 106,633587 38 ((38 + 3)И) -95,341926

16 ((16 + 3)И) 104,234448 39 *1 ((39 + 3)И) -102,15762

17 *1 ((17 + 3)И) 100,776413 40 \ ((40 + 3)И) -107,96994

18 ((18 + 3)И) 96,284807 41 *1 ((41 + 3)И) -112,70915

19 *1 (19 + 3)И) 90,795648 42 ((42 + 3)И) -116,31566

20 \ (20 + 3)И) 84,355389 43 *1 ((43 + 3)И) -118,74074

21 (21 + 3)И) 77,020554 44 *1 ((44 + 3)И) -119,94708

22 ((22 + 3)И) 68,857268 45 *1 ((45 + 3)И) -119,90932

23 (23 + 3)) 59,940685 46 ((46 + 3)И) -118,61439

24 \ (24 + 3)И) 50,354319 47 *1 ((47 + 3)И) -116,06178

25 1, ((25 + 3)И) 40,189286 48

Таблица 2

к 11 (к + п )) к 11 ((к + п )и)

0 *1 (0 +1)) 1 3 *1 ((3 + 3)И) 5,9978499

0 *1 (0 + 2 )) 1,99999 4 *1 ((4 + 3)И) 6,9962889

0 *1 (0 + 3)И) 2,99987 5 *1 ((5 + 3)И) 7,994118

1 *1 (1 + 3)) 3,99954 6 + 8,9912377

2 *1 (2 + 3))) 4,9989001 7 *1 ((7 + 3)И) 9,9875486

Видно, что значение контурного тока ((0 + 1)И), полученное с шагом И = 1 -10-5с, совпадает с его значением г1 ((7 + 3)И), рассчитанным с шагом И = 1 -10-6 с. Следовательно, для дальнейших расчетов, шаг сетки решения И = 1 -10-5 с можно принять за окончательный.

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

В целом предложенный метод прост, нагляден и удобен для реализации на ЭВМ.

Литература

1. Иващенко Н.Н. Автоматическое регулирование. Теория и элементы систем / Н.Н. Иващенко. М.: Машиностроение, 1978. 736 с.

2. Шеин А.Б. Решение уравнений состояний электрических и электронных цепей по ряду Тейлора /А.Б. Шеин // Изв. вузов. Электромеханика. 1998. № 2-3. С. 79-80.

3. Чиженко И.М. Основы преобразовательной техники / И.М. Чиженко, В. С. Руденко, В.И. Сенько. М.: Высш. школа, 1974. 430 с.

ШЕИН АЛЕКСАНДР БОРИСОВИЧ - кандидат технических наук, доцент кафедры промышленной электроники, Чувашский государственный университет, Россия, Чебоксары ([email protected]).

SHEYIN ALEKSANDR BORISOVICH - candidate of tehnical sciences, assistant professor of Industrial Electronics Department, Chuvash State University, Russia, Cheboksary.