CC BY
8УДК 541.12.0122: 539.213.3
В. А. Сиренек1
Введение
Объектом исследования является класс квазилинейных (линейных относительно вторых производных c"t и Схх) гиперболических уравнений вида:
ci't = W2 ■ с%х + /О, t, с, c't, 4), (1)
где x — пространственная координата; t — время; с (x, t) — потенциал переносимого физического поля; w — скорость распространения возмущений поля. Такие уравнения в качестве релаксационных моделей массо- и теп-лопереноса эффективны в тех случаях, когда характерное время химической реакции (теплового источника) меньше характерного времени переноса массы (тепла) [1-3]. Например, при описании турбулентного конвективного мас-со- и теплопереноса в трубчатом реакторе, осложнённого химической реакцией (быстрой полимеризацией) с выделением (поглощением) тепла [4]. Показательным случаем применения (1) является использование его в качестве волновой модели продольного перемешивания в трубчатом реакторе с нелинейной химической кинетикой [5].
При численном решении гиперболических уравнений универсальным методом конечных разностей сталкиваются с определенными трудностями, связанными с существенными искажениями крутых линий раздела фаз — эффекты диссипации (изменение наклона фронта профиля) и дисперсии («выбросы» концентрации на фронте), возникающими при переходе от дифференциальных уравнений к разностным [6]. Эффективным средством избежать подобных проблем служат основанные на методе Монте-Карло численно-вероятностные способы решения дифференциальных уравнений, не обладающие погрешностью, обусловленной дискретизацией задачи. Решение в этом случае сводится к осреднению некоторого функционала от реализаций моделируемого случайного процес-
ЧИСЛЕННО-ВЕРОЯТНОСТНЫЙ МЕТОД РЕШЕНИЯ РЕЛАКСАЦИОННЫХ УРАВНЕНИЙ МАССОПЕРЕНОСА НА ОСНОВЕ ВЫДЕЛЕНИЯ ИХ ВОЛНОВЫХ СВОЙСТВ
Санкт- Петербургский государственный технологический институт (технический университет), Санкт-Петербург, Московский пр., 26
Метод основан на обращении дифференциального оператора. Квазилинейные гиперболические уравнения массопе-реноса (релаксационные модели) сводятся к интегральным уравнениям. Ядра интегральных уравнений выражены через эволюционный оператор волнового уравнения. Получены выражения эволюционного оператора волнового уравнения для различных краевых задач. Приведен алгоритм решения интегральных уравнений методом Монте-Карло.
Ключевые слова: гиперболические уравнения массоперено-са, интегральные уравнения, случайный процесс, точные решения волнового уравнения, метод Монте-Карло.
са — имитационного или фиктивного (специально построенного) [7]. Представляющие собой модели «случайного блуждания» и лежащие в основе вывода дифференциальных уравнений массопереноса, имитационные случайные процессы позволяют не только решать эти уравнения методом Монте-Карло, но и изучать вероятностную природу характеристик массопереноса [3].
Выигрышным в вычислительном отношении при решении какой-либо задачи является использование решений более простых задач. В [8] проведено обобщение этого приема для ряда преобразований одной задачи в другую, известных в различных областях математической физики. Применительно к гиперболическим уравнениям это — использование точных решений волнового уравнения. Для линейных гиперболических уравнений этим приемом решались как задача Коши [9-14], так и краевые задачи [15-17], в [10-12, 15-17] — одномерные, а в [11, 13, 14] — многомерные. При этом использовались имитационные случайные процессы. Решения задачи в этом случае представлены в виде математического ожидания от «рандомизированных» (согласно имитационному случайному процессу) точных решений волнового уравнения. Численная реализация таких представлений осуществлена методом Монте-Карло.
Настоящая работа посвящена численно-вероятностному методу решения уравнений (1), в том числе нелинейных, основанному на обращении дифференциального оператора. С позиций обобщённого решения абстрактного эволюционного уравнения краевая задача для (1) сведена к задаче Коши для системы интегральных уравнений, ядра в которых определяются эволюционным оператором для волнового уравнения. Введен случайный (фиктивный) процесс и функционал от него, с помощью которого полученная система интегральных уравнений решена методом Монте-Карло. Теоретические
1 Сиренек Валерий Анатольевич, канд. техн. наук, доцент, каф. системного анализа, email: wasirenek@gmail.com Дата поступления - 04 апреля 2014 года
основы предлагаемого метода изложены в [18]. Из этой работы приведены сведения, необходимые для понимания и практического использования метода, в том числе алгоритм моделирования. Кроме того, в настоящей работе дается иллюстрации используемых математических конструкций и выражения эволюционного оператора волнового уравнения («рабочие» формулы) для различных краевых задач.
Обращение дифференциального оператора
Начально-краевая задача для уравнения (1) в общем виде рассмотрена как эволюционное уравнение относительно неизвестной функции V (t) в некотором функциональном (банаховом) пространстве E [19, гл.5]:
^ = A(t)V(t) + B(t)V(ty, ПО) = V0] V{t) G Zt, Vt, (2)
оператор A (t) — линейный, B (t) — непрерывный, вообще говоря, нелинейный. Ограничимся случаем линейных неоднородных (аффинных) граничных условий, т. е. подпространство Zt с E — аффинное и является результатом параллельного переноса линейного (направляющего) подпространства функций Zt а Е, удовлетворяющих линейным однородным граничным условиям [20]. Назовем обобщенным решением задачи (2) решение уравнения [19]:
НО = t/(t, 0)V0 + /о Ü (t, s)B(s)V(s)ds, (3)
где U (t, s): Zs Zt — эволюционный оператор для «невозмущенной» задачи:
™ = A(t)V(t); ПО) = ПО G Zt, Vt (4)
Оператор U (t, s) —_аффинный, непрерывный, т. е. для любых VeZs, WeZ5 имеет^ место соотношение U(V+W) = U(V) + U (W), где U (t, s)\ Zs^Zt — линейный оператор, называемый ассоциированным с U (t, s). (Иллюстративный пример аффинных подпространств и действия аффинного оператора приведен на рисунке 1). Предполагается, что для любой непрерывной V (t) функция B (t) V (t) интегрируема (в смысле Бох-нера [20, гл.2]) и для неё выполняется локальное условие Липшица. При этих условиях существование и единственность решения уравнения (3) устанавливается методом последовательных приближений.
Рисунок 1. Аффинные подпространства (и 12) и аффинный оператор (и) в евклидовом пространстве Я3.
Относительно неизвестной W (t) = U (0, t) V (t) уравнение (3) принимает вид:
w{t) = W0 + ff U(0,s)B(s)U(s, О)W(s)ds, (5)
где Wo = № (0) = V (0) = Уо. Уравнение (5) эквивалентно системе уравнений:
1^(0 = УУ0+ I 0(0,5)14^2
■'о
Ж2(0 = = и (г, 0Ж(О (6)
Через решение системы (6) определяется решение задачи (2): V (г) = и (г,0) №1(1).
Для уравнения (1) рассмотрим три типа началь-
но-краевых задач:
c\t=0 = Со О), c't |t=0 = ct0 О), X G (о,*); (7)
а) с|*=0 = 7iW), с'х\х=х = у2(t), t > 0; (8)
б) с|*=0 = 7i(О, с|м = 7з(t), t > 0; (9)
в) = УМ, с'х\х=х = у2(t), t > 0 (10)
Пусть coeC^O, X], сюеС[0, X], yiEC1[0, T], y2eC [0, T], узеС^О, T], Y4EC [0, T], и выполнены условия согласования начальных и граничных условий. Предполагается, что функция f из уравнения (1) непрерывна по х, интегрируема по t, по остальным переменным выполняется локальное условие Липшица. Введено пространство Е функций ф(х. Г). хе[0, X], /=1,2,3, непрерывных по х и таких, что ср(х,3) = (р'х (х, 1). В полученном функциональном (банаховом) пространстве Е определены операторы А и B (^применительно к (2):
Аф, 1) = ф,2), А(р(х,2) = м?ц>'х(х,3) = w2Ц)Хх(х,1), Аф,3) = ц>'х(х,2),
B(t)ф(хД) = 0, B(t)y(x,2) =J[x,t, ф(х,1), ср(х,2), ф(х,3)), 5(ОфМ = 0
При этом уравнение (1) эквивалентно абстрактному (2) для вектор-функции V (t)(x, i)eE, которая связана с решением с (x, t) уравнения (1) соотношениями:
V(t)(x,l) = c(x,t), V{t){x,2) = c't{x,t), V(t)(x,3)=c'x(x,t)
F(0)(x,l) =Vo(x, 1) = c0(x), V(0)(x,2) = Vo(x,2) = сю(х),
V(0)(x,3) = Vo(x,3) = dco(x)/dx,
Решение уравнения (1) понимается как решение уравнения (3). U (t, s) — эволюционный оператор задачи (4), который можно выписать явно, пользуясь точными решениями волнового уравнения и производными от них.
Эволюционный оператор для волнового уравнения
На основе интегрального уравнения колебаний [21]: fR (ftdx + w*d/xdt) + Д. foc,Odxdt (11)
где область G ограничена кусочно-гладкой кривой R, получены аналитические решения различных краевых задач для неоднородного волнового уравнения с начальным и граничными условиями, совпадающими с условиями (7)-(10):
cl't = w2c"+f(x,t), (12)
c|t=o = cQ(x), c't|t=0 = ctQ(x), x £ (0,xy, (13)
а) c\x=0 = 7i(t), c'x\x=x = 72(t)/ t > 0; (14)
б) c\x=0 = 7i(t), c\x=x = 73(t), t > 0; (15)
в) c'\x=0 = vAt), c'\x=x = y2{t), t > 0 (16)
При z = t - s, где s фиксировано, t > s, T =X/w точки прямоугольной области □ = {0 < z < T; 0 < x < X} имеют не более одного отражения характеристик уравнения (12) от границ отрезка [0, X] (см. рисунок в работе [17]). Диагоналями x — wz = 0 и x + wz = X область □ разбивается на 4 зоны:
I. {0 < z < X/2w, wz < x < X — wz}
II. {0 < z < X/2w, 0 < x < wz}u{X/2w < z < X/w,
0 < x < X — wz}
III. {0 < z < X/2w, X — wz < x < X}u{X/2w < z < X/w, wz < x < X}
IV. {X/2 w < z < X/w, X — wz < x < wz}
Зона I характеризуется отсутствием отражений характеристик волнового уравнения (12) от границ области, зона II — отражением от границы x = 0, зона III — отражением от границы x = X, зона IV — отражением от обеих границ.
Обозначим точки зон □: A1 (x, z)e I, A11 (x, z)e II, A111 (x, z)e III, AIV (x, z)e IV
B (xi, s), B'(— xi, s), C (x2, s), C'(x3, s), D (0, ti), D'(X, t2), E (0, s), E'(X, s),
где xi = x - wz, x2 = x + wz, x3 = 2X - x - wz, ti = z - x/w. t2 = z - (X - x)/w (i7)
Выражения c(x,t) для всех трех типов граничных условий (i4)-(i6) в зависимости от зон области □ приведены в [i7] при решении линейных гиперболических уравнений с использованием имитационных случайных процессов, при этом принято f(x,t) = 0. В нашем случае для построения эволюционного оператора U(t,s) используются как выражения d(x,t), так и полученные из них дифференцированием выражения ci(x,t), и Cx(pc,t),. Приведем результат действия оператора U (t, s) на функцию ф^, i), (I = i,2,3) с учётом обозначений (i7). При решении уравнения (i) с граничными условиями (i4) при обозначениях y(i,l) = yi(0, y(t,2) = у2(/), y(i,3) = (b(\{t)ldt, для точек зоны I имеем: u(t, s)q>0,1) = (l/2)fo>(*1( 1) + cp(x2,1) + (1/2w) Q ф(у, 2)dy;
U(t,5)ф(х, 2) = (1/2)2) + <p(x2,2)] + (w/2)[V(x2i 3) - <p(Xl, 3)] U(t,s)<p(x, 3) = (1/2) [фС^, 3) + ф(*2,3)] + (1/21у)[ф(*2,2)-ф(*1,2)];
(i8)
для точек зоны II:
U(t, s)<f,(x, 1) = (1/2) [tpfe, 1) - <p(-Xl, 1)] + (1/2w) [ 2 <p(y, 2)dy + y(tv 1)
U(t,s)<p(x,2) = (w/2)[—cp(—3) + <p(x2,3)]+ +(1/2)[ф(*2,2) - <p(-x1,2)] + Y(«i, 3); U(t,s)q>(x,3) = (1/2) [cpi-*!, 3) + ф(х2,3)] +(l/2w)M~Xi,2) + <Pfe.2)] - y(t1;3)/w;
(19)
для точек зоны III:
г* 3
U(t, 5)ф(ж, 1) = ф(ж3, l)/2 + ф(ж1; l)/2 + l/2w ф (y, 2)dy +
+l/w J* <p(y,2)dy + w fst2 у(y, 2)dy; U(t, s)<p(x, 2) = -(w/2) Mxt, 3) + 3)] + (1/w) [<p(x3,2) + ф(%, 2)] +
WY(t2,2),
U(t,s)V(*,3) = —(1/2)[ф(ж3,3) — ф(ж!, 3)] + +(l/2w)fo>(*3,2) - ф(*1,2)] + y(t2,2);
для точек зоны IV:
i/(t,sM*, 1) = Ф(х3,1)/2 - Ф(-х1( 1)/2 +
гХ з rX rt2
+ l/2w I ф(у, 2)dy + 1/w ip(y, 2)dy + y(a, 1) + w I у(у, 2)dy J-Xi -Ix-i Js
U(t, х)ф(х, 2) = - w[(p{-x1,3) + <p(x з, 3)]/2 +
ф(*3,2)/2 - 2)/2 + y(ti, 3) + wy(t2,2)
i/(t, я)ф(х, 3) = [ф(-х1( 3) - <p(x3,3)]/2 +
[фС-*!, 2) + ф(*3,2)]/2w - 1/w Y(i!, 3) + y(t2,2);
(2i)
При решении уравнения (i) с граничными условиями (i5) при
у(i,l) = yi(0, y(i,5) = уз(0, y(i,6) = yi(t)/dt, y(t,7) = J3(t)/dt, для точек зоны III:
U(t,s) = -(1/2)[V(-Xll) - Ф(х31)] + (l/2w) ф(у, 2)dy + Y(t25);
(22)
U (t, s) ф(^2) = (w/2) [ф (-xi,3) - ф (x3,3)] + (i/2) [ф(xз,2) -ф (-xi,2)]+Y (t2,7);
U (t, s) ф(^3) = (i/2) №(-xi,3) - ф(xз,3)] - (i/2w) [ф^3,2) + ф (-xi,2)]+ (i/w) Y(t2,7);
для точек зоны IV:
t/(i,i'M*,l) = y(ii,l) + yfe5) - (1/2)[ф(дя,1) + Ф(ЛГ2,1)]+ (V2-w)Q ф(у, 2)dy; Щ^УФЛ) = у('1,6) + y(fc,7) - (w/2)[9(«,3) - <p(»,3)] - (1/2)[р(*1,2) + ф(«,2)]; М[иЖ*.3) = (l/w)[ -y(fi,6) + rtft,7)] + (1/2)[р(л,3) + ф(»,3)] -
-(1/2№)[(р(Х1,2)-ф(й,2)];
(23)
При решении уравнения (i) с граничными условиями (i6) при обозначениях Y(t,8) = Y4(t), Y(t,2) = Y2(t) для точек зоны II имеем:
Е/Мф(х,1) = (1/2)[ф(х1,1) + ф(х2,1)]+ (1/2W)/*2 ф (у, 2)dy + +(1/2) /0Ж1 ф(у, 2)dy - wf^yiy, 8)dy
U (t, sMx,2) = (w/2)[ф(Xl,3) + ф(X2,3)] + (1/2)[ф(Xl,2) + ф(x2,2)] - wY(ti,8);
U (t, s)ф(x,3) = (1/2)[ф^2,3) - ф(Xl,3)] + (1/2w)[ф(X2,2) -ф(Xl,2)] + Y(ti,8);
(24)
для точек зоны IV:
= (1/2)[ф(х,,1) + Ф(Х2,1)]+ (l/w)/f2 ф(у, 2 )dy + +(l/w) j;1 ф(у, 2)dy + (l/2w) Q ф(у, 2)dy-
-^C y(y> 8)dy + w /st2 y(y, 2)dy
(/(^)ф(х,2) = (н,/2)[ф(дс1,3) - Ф(х2,3)] + (1/2)[ф(хь2)
+ ф(х:2,2)] - wy(ii,8) + wy(i2,2);
(/(^)ф(х,3) = - (1/2)[ф(хьЗ) + Ф(х2,3)] + (1/2)[ф(хь2)
+ ф(х2,2)] + y(ii,8) + y(t2,2);
Функция U (t, s) ф (x, i) в любой точке (x, i) выражается в виде линейной комбинации значений функций Ф и y в определённых точках и интегралов от этих функций по определенным промежуткам. При х = (x,i), t={t,i) получаем х Е -ДНО,*] х {1,2,3}, t е Г=[0,Г] х {1,2,3}. Формулы типа (18)-(25) принимают вид:
i/(t,s)<p(x) = 4 QtiS(x-, dyipiy) + iT С (*; dt)y(i),
где ядро Qt,s не учитывает границы, а ядро Q^ps - учитывает. Линейный оператор Ü(t,s), ассоциированный с U (t, s), выражается интегралом с ядром Qt, s.
Применение метода Монте-Карло
Решение уравнения (3) реализовано методом Монте-Карло через решение системы (6), где интегралы сведены к интегралам по вероятностным мерам:
w1(3c,t)) = w„(;E)+ [ PBp{t -» ds}aBp(t,s) f P0,s{x -» dy}a0,s(x;y)w2(y,s)
J 0 Jx
w2(x, t) = f [x, t, w3(x, 1, t),w3(x, 2, t), w3(x, 3, t)] при x = (x, 2); 0 при x =
(x, 1), x = (x, 3),
w3(x,t) =
(26)
P{-^} — вероятностные меры; a — «веса» («вр» и «гр»
— время и граница). Такое представление не единственно — следует учитывать его влияние на точность и трудоемкость расчёта методом Монте-Карло. Для построения мер Р(х -> ■■■] достаточно заменить коэффициенты линейных комбинаций (18)-(25) на положительные числа, сумма которых равна 1. Предполагается разложение f в ряд: f(_x,t,a1,a2,a3) = ,
где A,p,ve {0,1,2,...}, сходящийся при любых а\,аг,аъ; Рл,м,у а 0 — вероятности; ¿Pa^.v = 1; ß*,M,v — «веса». Задаются вероятности pi (р,£) и ргр ^t) выхода моделируемого в методе Монте-Карло случайного процесса на начальные и граничные условия, а также величины qH = 1 - pH(*.t). qrp = 1 - рФ<*-£). Теперь все подготовлено для введения марковского процесса с ветвлением [22]. Состояние случайного процесса — это конечный набор G точек («частиц»). Реализация процесса
— это последовательность состояний {G„}, n = 0,1,2.....
где n играет роль времени. Моделирование процесса начинается с 3-го уравнения системы (26) и осуществляется согласно выбранным вероятностям. При нелинейной (степенной) зависимости функции f от своих аргументов c,c't,c'x моделируется ветвящийся случайный процесс. Использована обработка дерева ветвлений не «по поколениям», а «по ветвям» [23]. «По ходу» моделирования вычисляется значение мультипликативного функционала R — произведения всех «весов», соответствующих переходам «частиц» в уравнениях (26) и ветвлениям, а также значений функций w0 и Y, соответствующих исчезновениям «частиц» при выходе случайного процесса на начальные условия или на границы отрезка [0,X] [24] (подробнее см. приведенный алгоритм, а также рисунок 2).
im '
т .(О
t д) ^
/ т
3 Yit)
0 X
/
Граничные условия
Начальные условия
Рисунок 2. Фазовое пространство моделируемого случайного процесса
Алгоритм моделирования одной реализации случайного процесса и вычисление мультипликативного функционала R
1. R: = 1.
2. Изъять из набора G одну точку; ее координаты присвоить переменным х, i, г.
3. Моделировать случайное событие вероятности ргр (х, i, г); если это событие наступило, то R: = R/pгр(х, i, О и перейти к п.9, иначе Я: = Я/9Ф(х, /, О-
4. Моделировать у с распределением - ^3'); Д := Я -щ:=у7;.
5. Моделировать случайное событие вероятности рн (х, i, г); если это событие наступило, то R: = R/pн (х, i, г) и перейти к п.10, иначе R: = R/qн (х, i, г).
6. Моделировать se(0, г) с распределением Р^Г^ = Р-авр(Г, э).
7. Моделировать у с распределением
РоЛ* = Лчи),1 (х; у).; X ■= у; V. = 5.
8. Моделировать (Л,р,у) с распределением р\м»(х, г); к набору G добавить Л точек (х,1, г), р точек (х,2, г) и V точек (х,3, г); R: = R• рЛ,м,„; перейти к п.11.
9. Моделировать t с
• dt}-, R~R- aTtp0(x; t) ■ y(t); перейти к п.11.
10. Я == Я -»у0(х)
11. Если набор в не пуст, то перейти к п.2, иначе — закончить моделирование.
Перед началом работы алгоритма задается набор G из одной точки (х, 0 = (*,£,£) Результатом конечного числа шагов алгоритма является случайная величина
R, математическим ожиданием которого является решение исходной задачи, т. е. с (х, 0 = Е^}.
В [18] в целях состоятельности оценок решения доказана конечность Е{|Я[] и ст2{|Я|}|} при I < +со. Предложены приемы для уменьшения ст2{Я}}, растущей с ростом I. Проведено исследование влияния интенсивности ветвления моделируемого процесса (при решении нелинейных уравнений) на среднее время расчета одной реализации процесса йр. Приведены результаты апробации метода при решении различных краевых задачи для уравнений вида (1), в том числе нелинейных. В качестве контрольных использованы известные решения, в соответствие с которыми задавались начальные и граничные условия. Проведено М = 104 реализаций процесса. Относительная погрешность — менее 5 %.
Заключение
Рассмотренный численно-вероятностный метод решения гиперболических уравнений массоперено-са вида (1) с коэффициентами (кроме коэффициентов при старших производных, имеющих постоянные значения), зависящими от пространства и времени (в том числе быстро меняющимися), применим для решения задач диффузионно-химической кинетики с учётом релаксационного эффекта и нелинейных источников (стоков) потенциала переноса. Численная реализация осуществляется методом Монте-Карло. Приводится алгоритм метода и «рабочие» формулы для решения для различных краевых задач.
Литература
1. Лыков А. В. Тепломассообмен. Справочник. М.: Энергия, 1972. 560 с.
2. Таганов И. Н. Моделирование процессов мас-со -и энергопереноса. Нелинейные системы. Л. Химия, 1979. 204 с.
3. Таганов И. Н., Сиренек В. А. Волновая диффузия. НИИХ СПбГУ. 2000. 209 с.
4. ХолпановЛ. П., Поляков Ю. С. Численный анализ гиперболической модели турбулентного теплообмена при химических превращениях // Труды XX Междунар научно-техн. конф. «Математические методы в технике и технологиях".Ярославль. 28-31 мая 2007 г. Ярославль: ЯГТУ, 2007. Т. 1. 202-204.
5. Вестертерп К. Р., Дильман В. В., Кронберг А. Е., БеннекерА. Волновая модель продольного перемешивания //Теор. основы хим. технологии. 1995. Т. 29. №6. С. 510-587.
6. Абиев Р. Ш. Вычислительная гидродинамика и тепломассообмен. СПб.: НИИХимии СПбГУ, 2002. 576 с.
7. Ермаков С. М., Некруткин В. В., Сипин А. С. Случайные процессы для решения классических уравнений математической физики. М.: Наука, 1984. 205 с.
8. Hersh R. The method of transmutations // Lect. Notes Math. 1975. V. 446. P. 264-282.
9. Кац М. Несколько вероятностных задач физики и математики.М.: Наука,1967. 176 с.
10. Hersh R. Stochastic solutions of hyperbolic equations // Lect. Notes in Math. 1975. № 446. Р. 283-300.
11. Kisynski J. On M.Kac's probabilistic formula for the solution of the telegraphist's equations // Annales polonici mathematici. 1974. V. 29. P. 259-271.
12. Сиренек В. А., Кузичкин Н. В. Вероятностное решение гиперболических моделей гетерогенного катализа // Каталитические процессы и катализаторы. Л.: ЛТИ им. Ленсовета, 1987. С. 122-128.
13. Кабанихина И. И. Численная реализация вероятностного представления решения задачи Коши для телеграфного уравнения. Сб. науч. тр. Методы статистического моделирования. Новосибирск: ВЦ СО АН СССР. 1986. С. 86-90.
14. Крючков А. Ф., Сиренек В. А. О численно-вероятностных методах решения трехмерного гиперболического уравнения диффузии // Журн. вычисл. матем. и матем. физики. 1998. Т. 38. № 1. С. 107-114.
15. Сиренек В. А., Крючков А. Ф. Вероятностное решение начально-краевой задачи для гиперболического уравнения массопереноса // Матем. моделирование. 1998. Т. 10. № 6. С. 107-117.
16. Сиренек В. А., Сидоров В. А. [и др]. Вероятностный подход к исследованию волновой модели продольного перемешивания // Теор. основы хим. технологии. 1999. Т. 33. N 5. С. 539-546.
17. Сиренек В. А. Связь релаксационных уравнений массопереноса с моделями «случайного блуждания» и способы решения диффузионных задач на их основе // Известия СПбГТИ (ТУ). 2014. N 23(49). С. 93-96.
18. Сиренек В. А. Вероятностное решение квазилинейных гиперболических уравнений на основе обращения дифференциального оператора // Журн. вычисл. матем. и матем. физики. 2000. Т. 40. № 3. С. 417-427.
19. Функциональный анализ (серия «Справ.-ма-тем. библ.») / под ред. С. Г. Крейна). М.: Наука, 1972. 544 с.
20. Бурбаки Н. Алгебра. Алгебраические структуры. Линейная и полилинейная алгебра. М.: Физмат-гиз, 1962. 516 с.
21. Тихонов А. Н., Самарский А. А. Уравнения математической физики. М.: Наука, 1977. 735 с.
22. Гихман И. И., Скороход А. В. Теория случайных процессов. М.: Наука, 1973.Т. 2. 639 с.
23. Соболь И. М. Численные методы Монте-Карло. М.: Наука, 1978. 311 с.
24. Ермаков С. М., Михайлов Г. А. Курс статистического моделирования. М.: Наука, 1976. 319 с.
