Вычислительные технологии Том 18, № 2, 2013
Численное моделирование динамики теплообмена модифицированным квадратичным полиномом Вольтерры*
С. В. СолодушА
Институт систем энергетики им. Л.А. Мелентьева СО РАН, Иркутск, Россия
e-mail: [email protected]
Доказана теорема сушествования решения некоторого парного двумерного интегрального уравнения Вольтерры I рода. Рассмотрен вопрос численного решения полиномиального уравнения Вольтерры I рода второй степени. Приводятся иллюстративные расчёты на примере эталонной динамической системы, описывающей процессы теплообмена.
Ключевые слова: интегральное уравнение Вольтерры I рода, моделирование, нелинейная динамическая система, автоматическое управление, теплообмен.
Введение
В теории математического моделирования нелинейных динамических систем хорошо известен универсальный аппарат функциональных рядов Вольтерры [1]. Полином Вольтерры Ж-й степени, дающий представление отклика системы типа вход-выход на внешнее возмущение х(£), имеет вид
N
y w = £ £ ./..■/
V =1 1<ii<... <i„<p П n
K
i1,...,il
(t, Si,..., Sv) Xii(si)... Xiv(sv)dsi,..., dsv,
t е [0,Т], (1)
где х(£) = (х^), х2(^), ...,хр(¿)) — вектор-функция времени, у(£) — скалярная функция времени, причём у(0) = 0, у(£) е С^Т]. В (1) ядра Вольтерры К^,..,^ симметричны по тем переменным, которые соответствуют совпадающим индексам ¿1, ...,г^. Построить математическую модель в виде (1) — значит идентифицировать ядра Вольтерры по информации об откликах системы на те или иные тестовые входные сигналы.
Если система стационарна в том смысле, что её динамические характеристики остаются неизменными за время Т переходного процесса, то вместо (1) используется модель
N t t
(t) = £ £ /.■/
V=1 1< ii< ...<iv< pJn i
y
0 0
Kil,...,iv (S1, ..., Sv)Xii (t - S1)...Xiv (t - Sv)dsb ...,dSv,
t e [0,T].
* Работа выполнена при финансовой поддержке РФФИ (грант № 12-01-00722а).
Проблема применения аппарата рядов Вольтерры для моделирования технических (в частности, электротехнических и теплотехнических) систем рассмотрена в работах [2 - 5]. Из зарубежных источников отметим имеющие обширную библиографию монографии [6-8].
Работы Института систем энергетики им. Л.А. Мелентьева (ИСЭМ) СО РАН в области идентификации ядер Вольтерры начались в 90-х годах прошлого века. В [9, 10] была предложена оригинальная методика идентификации ядер Вольтерры, базирующаяся на задании специальных многопараметрических семейств кусочно-постоянных тестовых входных сигналов. При этом задача идентификации сводится к решению линейных многомерных уравнений Вольтерры I рода с переменными верхними и нижними пределами (что отличает их от классических уравнений типа уравнений Вольтерры). Специфика таких уравнений позволяет находить решения в явном виде, а их сеточные аналоги допускают построение высокоэффективных вычислительных процедур, обладающих саморегуляризующим свойством. Важный для приложений случай (2) при N = 2 исследовался в [11, 12]. В работах [13-15] изложено применение данного подхода для моделирования процессов теплообмена на установке ВТК (высокотемпературный контур) ИСЭМ СО РАН.
Статья продолжает исследования, начатые в [12-15]. В первом разделе рассмотрена проблема идентификации несимметричных ядер Вольтерры, во втором дано краткое описание соответствующего программного обеспечения и приводятся иллюстративные расчёты на примере моделирования процессов теплообмена, третий раздел посвящён численному решению полиномиальных уравнений Вольтерры I рода, связанных с задачей об определении входного сигнала ж(£), которому соответствует заданный выход у(£).
1. Об одном способе идентификации несимметричных ядер Вольтерры
В работе [16] предложен алгоритм повышения точности моделирования нелинейной динамической системы полиномами Вольтерры в скалярном случае. В его основе — комбинация линейной нестационарной и билинейной стационарной составляющих. В случае N = 2 комбинация (1), (2) дает уравнение
р * р * *
0 М=1 0 0
Р М-1 * *
+ / / К^Дзь - - 52)^51^52, £ € [0,Т].
^=2 ь-=1
00
Наша задача — получить необходимые условия существования К^^, 52) в классе непрерывных на П2 = |в1, 52/0 < 51, 52 < Т} несимметричных функций.
Отметим, что процедуру разделения отклика моделируемой системы на составляющие у"г, г = 1, 2, обусловленные индивидуальным влиянием ^-й компоненты вектора ж(£), реализуют тестовые сигналы
(£) = (£) - I(£ - ш)), жА(£) = 0, 0 < ш < £ < Т, (4)
где I(С) — функция Хевисайда, Л = 1,р, Л = Вещественные числа аг,* = 0, г = 1, 2, = а2,*, ^ = 1,р, характеризуют высоту возмущающих воздействий по входу х*. Подстановка (4) в (3) приводит к системе двух линейных относительно К* и К** интегральных уравнений Вольтерры I рода с переменными верхними и нижними пределами интегрирования, при этом справедливы формулы обращения [16]
V и \ д/*(С,ш)
Кп(г,ш) = ---, 0 < ш < с < Т ,
дш
с (С С ш) Vд2/**(С ,ш) + /МЛ 0< ,,< + < Т
где
ш) - <*уа2'и_ «1,*^ (с, ш) - (С, ш)
УмС^Н —-7-!-ч-, / \ ,
«1,*«2,*(а2,* - «1,*) «1,*а2,*(а2,* - «1,*)
причём из необходимого условия существования К**^,^) в классе непрерывных на П2 симметричных функций [11]
следует равенство
«1,* + «2,* = 0.
Постановка и решение экстремальной задачи выбора а* = ±аг,* в (4) для идентификации ядер К*, КС** в (3) приведены в [17].
Остановимся на проблеме восстановления на квадрате П2 ядра К^*^, з2). По аналогии с [12] рассмотрим две группы тестовых сигналов
х^ = вV(I(С) - I(с - ш)), х^ = в*1 (С - ш), хА(С) = 0, 0 < ш < С < Т, (5) хв^ = вvI(С - ш), х*и = в*(1 (С) - I(С - ш)), ха(с) = 0, 0 < ш < С < Т, (6)
где в^ = 0, в* = 0 — амплитуды сигналов по входам xV, х*, ^ = 2,р, V = 1, ^ - 1, Л = 1,р, V = ^ = Л. Подстановка (5), (6) в (3) приводит к парному двумерному интегральному уравнению Вольтерры I рода относительно несимметричного на П2 ядра К,*
4 4—ш
I У К^1, 32)^2 = (С,ш), 0 < ш < С < Т, (7)
4—ш 0
4—ш 4
J У К^(31, = (С,ш), 0 < ш < С < Т, (8)
0 4—ш
где
в^ви Ш 4
Лгви (С, ш) = ^ в;ш) - Ц К (С, - в- I к* (с,
V , оу юО I 1 * V
0
* * ш ш
- 1 1 К™ (51,52)^51^52 - ^ ! j 52)^51^52,
"*—ш*—ш О О
* ш
Р^ 2д \
•А вд
^вд (£,ш) = - К(М)^ - I К"(М)^-
0
4—ш 4—ш * *
- 1 1 К™ (51,52)^51^52 - вв^ ^ J 52)^1^2. (10)
Здесь у2^2д (£,ш), у2^2д(£,ш) — отклики системы на сигналы вида (5) и (6) соответственно.
Теорема. Условия
+ д!^) € Сд, г =1,2, Д= ш/0 < ш < £ < Т}, (11)
■ув"2"((,0) = ув"в"(М) = в-("2-- - ^> у1,(М) + &- ) у».,(и), (12)
ув^в"(£, 0) = у^2"(М) = в"(а2." - в"\ уа1,д (£,£) + в"(в" - а1."\ уа2,д (£,£), (13)
а1."(«2." - «1.") а2."(«2." - «1.")
/V (2^ - /V) ( д V1-(£,ш) + дУ1- (£,ш) N / д2ув-2д (£,ш) + д2ув-2д (£,ш)') +
(а2^ - )\ д£дш дш2 /ш=0 V д£дш дш2
ш=0 \ ^^^ ^^ / ш=0
+ /V(/V - 2а1.у) (дУ2- (£,ш) + дУ2>- (£,ш) «2^ («2^ - А д£дш дш2 , ш=0
= в"(2а2." - в") (д2у"1,Д (£,ш) + д2уа1,д (£,ш) \ - (д^вд (£,ш) + д2у2в^вд (£,ш) N + а1."(«2." - «1."Д д£дш дш2 /ш=0 \ д£дш дш2 /ш=0
+ /"(/" - 2а1.") (д2уа2,Д(£,ш) + д2уа2,Д(£,ш) \ (14)
а2."(«2." - «1."Д д£дш дш2 /ш=0
необходимы и достаточны для существования решения парного уравнения (7), (8) в классе Сл2. Пара
К ) вд (£,ш) + д вд (¿,ш)\ Т (15)
- ш) = -( д£дш +-дЩ2-Ь 0 < ш < ' < Т, (15)
i^v"(t-ш,£) = + , 0 <ш <£<Т, (16)
определяет решение во всей области П2.
Доказательство. Воспользуемся результатами [11], согласно которым формулы обращения (15), (16) справедливы тогда и только тогда, когда
(^ + € < =1,2, Д = {£,ш/0 < ш < £ < Т} (17)
gfv/3M(t, 0) = 0, gfv^(t,t) = 0, i = 1,2, t G [0,T], (18)
"(t, W) + d 2gf"(t,^) \ = / d 2g2e "(Ы + d2 g2e v (t,^ (1Q)
5ш2 У,, п \ 5ш2 ,,, п
/ ш—и V / ш—и
С учётом (9), (10) из (17)—(19) следуют (11)-(14). □
2. Программное обеспечение для моделирования нелинейной динамики
Для построения и тестирования интегральных моделей нелинейной динамики теплообмена был создан программно-вычислительный комплекс [18], использующий эталонную модель. Эталоном служило описание процесса теплообмена в элементе теплообменно-го аппарата (теплообменника) с независимым подводом тепла, представленное в [19]. Согласно [19], отклонение энтальпии на выходе Дг(£) (кДж/кг) при произвольных законах возмущений расхода вещества ДТ>(£) (кг/с) и тепловой нагрузки Д2(£) (кВт) описывается зависимостью
t t t
Ai Ao С f On \ / -Л1 Г -Л2 f D(s)ds\
Aiet(t)= лТ-Т^У (AQ(n) (e n - e n j dn, t G [0,T ]. (20) о
Здесь t — время, Ai и A2 — некоторые константы, индексами "0" обозначены параметры начального стационарного режима, D0 = 0.16 (кг/с), Q0 = 100 (кВт), A — приращение, например D(t) = D0 + AD(t). Числовые характеристики, входящие в (20), принимались соответствующими реальной установке ВТК.
Вычислительный комплекс создан в объектно-ориентированной среде программирования C+—+ Builder и основан на функционально-модульном принципе. Он включает в себя блоки настройки входных параметров, идентификации и моделирования. Соотношение (20) используется в блоке идентификации для получения необходимого набора откликов на тестовые сигналы, а в блоке моделирования — в качестве эталона для оценки точности интегральных моделей. Процедуры идентификации базируются на разностных аналогах формул обращения, что обеспечивает быстродействие в режиме "on-line". Процедуры вычисления откликов интегральных моделей основаны на использовании методов средних прямоугольников и интегрирования произведения. Расчёт выходных значений эталонной модели (20) проводится с помощью метода трапеций. Реализованы функции ввода и оцифровки входных воздействий, отображения результатов вычислительного эксперимента в графической форме. Вся выходная информация хранится в соответствующих файлах на диске и может быть использована для подробного ознакомления и анализа.
Проведём тестирование квадратичной модели (3), где x(t) = (AD(t), AQ(t)). В ходе вычислительного эксперимента будем учитывать результаты расчётов, полученные с помощью (2), где N = 2. С учётом специфики (20) вместо (2) и (3) для p = 2 используем соответственно уравнения
t t t Aiimod(t) = y KKi(s)AD(t - s)ds + J j Kii(si,S2)AD(t - si)AD(t - s2)dsids2+ 0 0 0
Ь Ь Ь
+ у КЭДДе^ - з)^ + 1 У Х^ь^Д^ - - 82)^1^2, * € [0, Т], (21)
0 0 0
и
Ь Ь Ь
Д22тоа (*) = / (М)ДР(в)^ + 1 У Хц^ь^Д^ - 51 )Д£>(* - 82)^1^2 + 0 0 0 Ь Ь Ь
+ У + 1 у Кг 12(81 ,в2)Д^(* - в1)Дд(* - ^м^, * € [0,Т]. (22)
0 0 0
Введём сетку узлов г = 1,п, пЛ = Т. Пусть Л, = 1 (с) (такой выбор шага связан с реальными данными, полученными в ходе эксперимента на ВТК ИСЭМ СО РАН). С помощью вычислительного комплекса [18] выполняем построение квадратичных моделей (21), (22). В данном программном обеспечении используются алгоритмы численного решения уравнений Вольтерры I рода, основанного на саморегуляризующем свойстве [20, 21] процедуры дискретизации. Отметим, что в работе [22] проведено подробное исследование применимости простейших кубатурных формул (правых и средних прямоугольников) и методов типа метода Рунге — Кутты для устойчивого решения соответствующих многомерных интегральных уравнений Вольтерры I рода. Следуя [22], в качестве "базовой" используем формулу средних прямоугольников. Такой выбор обусловлен относительной простотой реализации и получением приближенного решения с погрешностью порядка 0(Л2) в случае отсутствия возмущений исходных данных. Считаем далее в качестве допустимого множество X(В,Т) входных сигналов вида
Д^(*) = №(/(*) - 21 (* - и) - 1(* - Т)), Д^(*) = 0, (23)
где в € [0,В], 0 < и < * < Т, В = 35, 30, 25.
Во многих физических приложениях важно значение выходного сигнала динамического объекта в конце переходного процесса, поэтому в качестве критерия точности моделирования выбираем значения невязок между Дг^Ь(^) и Дг^ а(*г) (Дг2 й(*г)) при
и = Т то то
= В± |Дг2'(ТД- ^(Т" 100%, ; = 1,2. (24)
в 1
Для сравнения эффективности применения моделей вида (21), (22) были найдены значения е1(Т), е2(Т) по (24) с точностью 8 = 10-4 для 1 < и < 20, Т = 30. Рисунок иллюстрирует результаты расчётов. Следует заметить, что области предпочтительности использования того или иного алгоритма не являются универсальными, так как выявлены на основе анализа результатов вычислительных экспериментов, проведённых с использованием одной эталонной модели (20).
Вычисление откликов Дг2 ¿(¿г), а(*г) на заданное входное возмущение (23) при 0.0016 < в^0 < 0.0016В проведено по разностным аналогам (21), (22), ядра Вольтерры в которых настроены на тестовые сигналы с амплитудами 0.032 < |а|Т>0 < 0.0016В, от выбора которых существенно зависит точность моделирования исследуемого нелинейного процесса. Естественно, что для обеспечения е1(Т) = е2 (Т) = 0% требуется согласование амплитуд пробных сигналов с величиной действующих возмущений. Приведём
■■■■
Результаты вычислительных экспериментов для В = 35 (а), 30 (б), 25 (в)
значения невязок для наиболее "неприятных" входных сигналов из множества допустимых. Было получено, что для В = 25 £\(Т) < 5.62%, £2(Т) < 5.18%, для В = 30 £1 (Т) < 6.89%, £2(Т) < 6.38%, для В = 35 £1 (Т) < 8.41%, £2(Т) < 7.93%.
3. Полиномиальные уравнения Вольтерры в задаче автоматического управления
До сих пор центральной проблемой в работе было восстановление ядер Вольтерры в (3), (22). Допустим, что с помощью того или иного подхода эта проблема решена. В качестве следующего этапа математического моделирования можно рассмотреть типичную задачу автоматического управления динамическим объектом при отсутствии обратной связи — нахождение входного сигнала х(£), которому соответствует заданный отклик у(£). Продолжая исследование, начатое в [23], выберем в качестве управляющего воздействия сигнал ж1(^) = ДР(£). Для определённости входное возмущение х2(£) = ДQ(t), £ € [0,Т], считаем заданным и перепишем (22) в следующем виде:
г г г
! К1(£,Й)ДР(Й)^5 + У J Кп (Й1 ,з2)ДХ>(£ - 51 )Д£>(£ -
0 0 0 г г
+ // К12(«1,52)Д^(£ - 51 - 52)^1 ^2 =
00
Д^тоа(и) - К2(и,8)ДС(в)Ж, и € [0,Т]. (25)
2mod
0
При известных К1, К2, КС11, К 12 и у(и) = Дг2тс^ (¿) (25) является полиномиальным интегральным уравнением относительно ДР(и). При малом и численное решение (25) заведомо существует [24]. Найдем его кубатурным методом средних прямоугольников. Обоснование сходимости численного решения полиномиального уравнения Вольтерры I рода второй степени, основанного на методе средних прямоугольников, дано в [25].
1
Введём сетку узлов и = гЛ, 1 = (г - -1 Л, г = 1, п, пЛ = Т. Аппроксимируем интегралы в (25) суммами. Для вычисления ДР(и) в - ^ -м узле получим квадратное
уравнение относительно ДТ2 1
г 2
2 / г
Л2К111,1 (Д^г- 2) + 1 ,_1+2Л2 Е ц^-1,1 ДТ ^+
+Л2У К121 1Д22 1 ДТ2 1 = Дг2 „ - 1 1Д22 1 - (26
г=1 7 ¿=1
где
г- 1 г г
' КГ,,_ ,,_ 2 _2 2
¿=1 .7=2 1=2
V 1,1_ 2 ДТ ^1 Д0-1+1, г = 1,п.
.7=2 1=1
Выбор нужного корня уравнения (26) определяется условием
Дг'2 (0)
ДТ1 ^ дт (0) = при Л ^ 0. (27)
2 К1 (0, 0)
Интересующее нас численное решение уравнения (26) обозначим через ДТ^ 1 Рассуждая по аналогии, рассмотрим вместо (21) уравнение
Ь Ь Ь
У^^ДТСи - ^ + У У КК11(51,52)ДТ(* - 51)ДТ(* - +
0 0 0
Ь Ь
+ УУ 12(51, 52)ДТ(* - 51)Дд(* - 52)^51^52 = 00
Ь
= Дг^(и) - У ^2(5)Да(и - в)Ж, и € [0,Т], (28)
0
Ь
Таблица 1. Значения 71тах, при которых существует решение Др"^), 0 < и < Т - Н
Т, с |а|, %
40 45 50 55 60
30 54.98 55.56 56.68 58.43 61.00
35 53.86 54.55 55.82 57.77 60.63
40 53.12 53.89 55.27 57.37 60.45
45 52.59 53.43 54.90 57.10 60.32
50 52.22 53.11 54.63 56.93 60.26
Таблица 2. Значения 72тах, при которых существует решение ), 0 < < Т - Н
Т, с |а|, %
40 45 50 55 60
30 55.59 56.12 57.17 58.83 61.28
35 54.25 54.98 56.12 58.00 60.76
40 53.38 54.13 55.47 57.51 60.50
45 52.78 53.60 55.03 57.20 60.36
50 52.35 53.23 54.73 56.99 60.28
где К1, К2, К11, К12, х2(£) = ДQ(í) и у(£) = Дг1тоа(£) известны. Аппроксимируем (28) методом средних прямоугольников и, используя условие типа (27), выполним построение Др (расчётная схема приведена в [23]).
Предположим, что допустимые входные возмущения ДQ(í) = 7Q0/(£), 7 < 0.75, £ € [0,Т]. Вычисление управляющих сигналов Др 1) , Др , обеспечиваю-
щих отклики системы Дг^^ = 0, = 0 соответственно, проводилось с помощью
уравнений (28), (25), ядра Вольтерры в которых были настроены на тестовые сигналы с амплитудами |а|Т0 < 0.096 (кг/с) и |а^0 < 60 (кВт). В табл. 1, 2 приведены некоторые результаты вычислительных экспериментов для г0 = 434 (кДж/кг), к = 1 (с), Т = 30 ^ 50 (с).
Основная специфика рассматриваемых полиномиальных уравнений Вольтерры I рода (25), (28), как и в скалярном случае, связана с локальностью области существования вещественных непрерывных решений. При вычислении сеточных аппроксимаций ДТ1(£), ДТ2(£), 0 < £ < Т, получено, что специфика (28), (25) проявляется при возмущающих воздействиях ДQ(í) с 7 > 71тах и 7 > 72тах (см. табл. 1, 2), что означает возможную потерю управляемости изучаемого процесса теплообмена. Вычислительные эксперименты показали, что на 71тах, 72тах влияет выбор а, используемых для восстановления ядер Вольтерры. Кроме того, из таблиц 1, 2 видно, что при фиксированных значениях Т и а имеет место неравенство 71тах < 72тах.
Таким образом, в работе теоретически обоснованы необходимые условия существования К^*, ^ = 2,р, V = 1,^ - 1, в классе непрерывных на квадрате П2 несимметричных функций. Данный результат получен в терминах амплитуд тестовых сигналов, что в сочетании с серией специальных экстремальных постановок [17] позволит в дальнейшем снять произвол в выборе амплитуд при построении интегральных моделей динамики теплообмена. Рассмотрен алгоритм численного решения полиномиального интегрально-
го уравнения Вольтерры I рода второй степени, возникающего в задаче автоматического регулирования нелинейной динамической системы. Приведены результаты тестовых расчётов для эталонной модели теплообмена. Вычислительные эксперименты проводились с помощью авторского программного обеспечения [18, 26].
Список литературы
[1] Вольтерра В. Теория функционалов, интегральных и интегро-дифференциальных уравнений. М.: Наука, 1982. 302 с.
[2] Веников В.А., Суханов О.А. Кибернетические модели электрических систем. М.: Энергоиздат, 1982. 327 с.
[3] Данилов Л.В., Матханов Л.Н., Филиппов В.С. Теория нелинейных динамических цепей. М.: Энергоиздат, 1990. 252 с.
[4] Дейч А.М. Методы идентификации динамических систем. М.: Энергия, 1979. 240 с.
[5] Пупков К.А., Капалин В.И., Ющенко А.С. Функциональные ряды в теории нелинейных систем. М.: Наука, 1976. 448 с.
[6] Doyle III F., Pearson R., Ogunnaike B. Identification and Control Using Volterra Models. Springer-Verlag, 2002.
[7] Rugh W.J. Nonlinear System Theory: The Volterra—Wiener Approach. Johns Hopkins Univ. Press, 1981.
[8] Ogunfunmi T. Adaptive Nonlinear System Identification. The Volterra and Wiener Approaches. Santa Clara Univ., Springer, 2007.
[9] Apartsyn A.S. Mathematical modelling of the dynamic systems and objects with the help of the Volterra integral series // EPRI-SEI Joint Seminar. Beijing, China, 1991. P. 117-132.
[10] Apartsyn A.S. On some identification method for nonlinear dynamic systems// ISEMA-92. Shenzhen, China, 1992. P. 288-292.
[11] Апарцин А.С. Теоремы существования и единственности решений уравнений Вольтерра I рода, связанных с идентификацией нелинейных динамических систем (векторный случай). Иркутск: СЭИ СО РАН, 1996 (Препр. № 8).
[12] СолодушА С.В. Построение интегральных моделей нелинейных динамических систем с помощью рядов Вольтерра. Автореф. дис. ... канд. физ.-мат. наук. Иркутск: Иркутский гос. ун-т, 1996.
[13] Апарцин А.С., Таиров Э.А., СолодушА С.В., Худяков Д.В. Применение интегро-степенных рядов Вольтерра к моделированию динамики теплообменников // Изв. РАН. Энергетика. 1994. № 3. С. 138-145.
[14] Апарцин А.С., СолодушА С.В., Таиров Э.А. Математические модели нелинейной динамики на базе рядов Вольтерра и их приложения // Изв. Академии естественных наук. Математика, матем. моделирование, информатика и управление. 1997. Т. 1, № 2. С. 115-125.
[15] Apartsyn A.S., Solodusha S.V. Mathematical simulation of nonlinear dynamic systems by Volterra series // Eng. Simulation. 2000. Vol. 17, No. 2. P. 143-153.
[16] Апарцин А.С. О повышении точности моделирования нелинейных динамических систем полиномами Вольтерра // Электронное моделирование. 2001. Т. 23, № 6. С. 3-12.
[17] Аплрцин А.С., СолодушА С.В. Об оптимизации амплитуд тестовых сигналов при идентификации ядер Вольтерра // Автоматика и телемеханика. 2004. № 3. С. 116-124.
[18] СолодушА С.В. Программно-вычислительный комплекс для моделирования нелинейной динамики теплообмена на базе квадратичных полиномов Вольтерра: Свид-во о гос. регистрации программ для ЭВМ № 2012614246. 12. 05. 2012.
[19] Таиров Э.А. Нелинейное моделирование динамики теплообмена в канале с однофазным теплоносителем // Изв. АН СССР. Энергетика и транспорт. 1989. № 1. С. 150-156.
[20] Апарцин А.С., Бакушинский А.Б. Приближенное решение интегральных уравнений Вольтерра I рода методом квадратурных сумм // Дифференциальные и интегральные уравнения. Иркутск: Иркутский гос. ун-т. 1972. Вып. I. С. 248-258.
[21] Апарцин А.С. Дискретизационные методы регуляризации некоторых интегральных уравнений I рода // Методы численного анализа и оптимизации. Новосибирск: Наука, Сиб. отд-ние. 1987. С. 263-297.
[22] Апарцин А.С. Неклассические уравнения Вольтерра I рода в интегральных моделях динамических систем: Теория, численные методы, приложение. Дис. ... д-ра физ.-мат. наук. Иркутск: Иркутский гос. ун-т, 2000.
[23] СолодушА С.В. Приложение нелинейных уравнений Вольтерра I рода к задаче управления динамикой теплообмена // Автоматика и телемеханика. 2011. № 6. С. 133-140.
[24] Апарцин А.С. Полиномиальные интегральные уравнения Вольтерра I рода и функция Ламберта // Труды Института математики и механики УрО РАН. 2012. Т. 18, № 1. С. 69-81.
[25] Апарцин А.С. О сходимости численных методов решения билинейного уравнения Вольтерра I рода // Журн. вычисл. математики и матем. физики. 2007. Т. 47, № 8. С. 1378-1386.
[26] СолодушА С.В. Программное обеспечение и алгоритмы для моделирования нелинейной динамики полиномами Вольтерра // Программные продукты и системы. 2012. № 4. С. 137-141.
Поступила в 'редакцию 24 декабря 2012 г., с доработки — 6 февраля 2013 г.