Работа выполнена при поддержке АВЦП Развитие научного потенциала высшей школы (2009Ц2010 годы), проект № 4761 Разработка методов исследования немарковских систем массового обслуживания и их применение к сложным экономическим системам и компьютерным сетям связи.
РАСЧЕТ ДИСПЕРСИЙ ЧИСЛЕННОСТИ ВОЗРАСТНЫХ ГРУПП В МОДЕЛИ ДВИЖЕНИЯ НАСЕЛЕНИЯ ТЕРРИТОРИИ И. В. Новикова Филиал Кемеровского государственного университета в г. Анжеро-Судженске Рассматривается модель динамики возрастного состава населения [1], которая представлена в векторно-матричном виде следующим образом:
x(k +1) = A(k)x(k), x(0) = x0, (1) где x(k) - m -мерный вектор составов населения территории; A(k) - матрица естественного движения населения.
Обозначим M{x0} = x0; D(0) = D(x0) = D0.
Предположим, что коэффициенты рождаемости являются случайными последовательностями. Поэтому матрица A имеет вид:
~ A(k) = A + A(k), ( 2) где A - матрица естественного движения населения со средними значе~ ниями коэффициентов, а матрица A имеет вид:
0 L 0 al al+1 L al+m 0 L 0 L 0 0 0 L 0 0 L ~, A = (3) L L L L L L L L L L 0 L 0 0 0 L 0 0 L где ai (i = l,l + m) - случайные отклонения от средних значений, которые описываются стохастическими разностными уравнениями вида:
a (k +1) = g a (k) + r wj (k), (4) j j j j где a (0) ~ N(0;s2 ); w(k) - дискретный гауссовский белый шум с харакj a j теристиками wj ~ N(0;1); wj ^ wi "i j.
Введем обозначения:
1) m(k) = M{x(k)};
~(k) 2) x = x(k) - m(k);
3) D(k) = M{~(k)~T (k)};
x x ~ ~ T 4) R(k) = M{A(k)[x(k) - m(k)] } = M{A(k)x(k)}; (5) Уравнение для математического ожидания вектора общей численности населения x(k) для момента времени k + 1, учитывая введенные обозначения, имеет вид m(k + 1) = Am(k) + R(k); m(0) = M{x(0)}. (6) Уравнение для ковариационной матрицы D(k) с учетом того, что в ~ ~ ~ силу специфики матрицы A ABAT = 0 для любой квадратной матрицы B таких же размеров, что и матрица A, имеет вид ~ ~ D(k +1) = AD(k)AT + AM{~(k)mT (k)AT (k)}+ AM{~(k)~T (k)AT (k)} + x x x ~ ~ + M{A(k)m(k)~T (k)}AT + M{A(k)~(k)~T (k)}AT - R(k)RT (k);
x x x D(0) = D0. (7) ~ ~ Матрицы AM{~(k)~T (k)AT (k)}, M{A(k)~(k)~T (k)}AT можно выx x x x ~ числить в предположении нормальности вектора x.
~ Если дисперсия x много меньше m2, то слагаемыми, в которых T ~ ~ присутствуют одновременно x и x, можно пренебречь, и для D получается уравнение:
~ D(k +1) AD(k)AT + AM{~(k)mT (k)AT (k)}+ x (8) ~ x + M{A(k)m(k)~T (k)}AT - R(k)RT (k), где R - вектор-столбец с ненулевым элементом на первом месте; RRT - матрица с одним ненулевым элементом на месте 1,1.
Для более удобного представления выкладок рассмотрим отдельно элемент матрицы, входящей в состав второго слагаемого в (8):
l+m ~ l+m ~ ~ (M{~mT AT }) = M x msa = M{~ia }, (9) x x ij ~i js ms js s=l s=l ~ где a = 0 " j >1.
js Следовательно, l+m ~ ~ (M{~mT AT }) = M{~ia1s}d1 j.
x x ij ms s=l Аналогично получаем уравнение для вычисления элементов матрицы, входящей в состав третьего слагаемого в (8):
l+m l+m ~ ~ ~ (M{Am~T }) = M{~jais} = M{~ja1s}d1i, (10) x x x ij ms ms s=l s=l где dij - символ Кронекера.
~ Так как a1s = as, то выражения (9) и (10) примут вид:
l+m ~ (M{~mT AT }) = M{~ias}d1 j ; (11) x x ij ms s=l l+m ~ (M{Am~T }) = M{~jas}d1i. (12) x x ij ms s=l Введем обозначение:
gij (k) = M{~i (k)a (k)}. (13) x j Далее получим для gij разностное уравнение gij (k + 1) = M{~i (k +1)a (k +1)}. (14) x j Рассмотрим поочередно составляющие выражения (14):
l+m ~ ~ ~ ~ xi (k +1) = xs (k) + ais (k)ms (k) + ais (k)~s (k)]- Ri (k), (15) x [Ais s=l где l+m Ri (k) = (k)d1i. (16) gss s=l Второй множитель в выражения (14) имеет вид (4). Таким образом, (14) запишется следующим образом:
l+m gij (k + 1) = g gsj (k) + g ms (k)M{as (k)a (k)}d1i + [Ais j j j (17) s=l g M{~s (k)as (k)a (k)}d1i ].
x j j Введем обозначение:
qsj (k) = M{as (k)a (k)}. (18) j Предположим, что " s j M{as (k)a (k)} = 0, поэтому j qsj (k) = s2 dsj. (19) a j ~ Если дисперсия x много меньше m2, то последним слагаемым в (17) можно пренебречь, и получаем формулу l+m gij (k + 1) g gsj (k) + g m (k)s2 d1i;
jAis j j a j s=l gij (0) = 0. (20) Таким образом, система разностных уравнений для m(k), D(k), gij (k) позволяет рекуррентно вычислять m(k) = M{x(k)} и D(k) = M{[x(k) - m(k)][x(k) - m(k)]T }. Диагональные элементы матрицы D являются дисперсиями возрастных составов населения территории.
итература 1. Тихомиров Н. П. Демография. Методы анализа и прогнозирования. - М.: Экзамен, 2005. - 256 с.
ЧИСЛЕННОЕ МОДЕЛИРОВАНИЕ ЭЛЕКТРОФИЗИЧЕСКИХ ХАРАКТЕРИСТИК МАТЕРИАЛА ПРОВОДНИКА С МАЛЫМИ ДИЭЛЕКТРИЧЕСКИМИ ПОРАМИ Ю. А. Попкова, В. Д. Власенко Тихоокеанский государственный университет, г. Хабаровск Вычислительный центр ДВО РАН, г. Хабаровск Численное решение задач с неоднородностями разностными методами требует специальных сеток для разрешения особенностей. Для этого необходимо использовать сетки, сгущающиеся в окрестности особенностей или достаточно мелкие сетки с шагом h и большим количеством точек. Первый вариант требует применения специальных алгоритмов, второй - соответствующей памяти компьютеров. В то же время проявления особенностей зачастую являются локальными, сосредоточенными в мелкомасштабных подобластях. Наличие областей сосредоточения неоднородностей позволяет ввести сетку с характерным размером H >> h, узлы и ребра которой проходят по участкам относительной гладкости решения.
При этом сетка размером H не позволит разрешить особенности при использовании обычных численных методов, но зато число ее узлов достаточно мало.
Для решения таких задач на сетках размером H предложен метод конечных суперэлементов (МКСЭ) [1 - 3].
Работа посвящена применению МКСЭ для расчета распределений электрического потенциала и плотности тока в проводящих объектах. Материал таких объектов содержит малые диэлектрические поры. Область, занятая проводником, составляет значительную величину по сравнению с областью, занятой порами.
Пространственная область представляет собой куб с гранями, параллельными координатным плоскостям. Выделим, - подобласти, занятые проводящим материалом и порами в нем соответственно. Через, обозначим границы этих областей, b, i - основания куба, j - его боковую границу.
Целью решаемой задачи является определение эффективных характеристик: удельного электрического сопротивления и удельной электрической проводимости пористого материала. Пусть проводящий материал, занимающий область, обладает заданной проводимостью 0 и удельным сопротивлением 0=1/0.
Плотность электрического тока j в материале задается уравнением div j = 0, (1) выражающим закон сохранения полного заряда.
Электрическое поле потенциально и удовлетворяет уравнению rot E = 0. (2) Связь плотности тока j и напряженности поля E определяется свойствами вещества и выражается законом Ома j= 0E. (3) В однородном проводнике 0 = const, откуда divE = 0. Поэтому в нем потенциал электрического поля удовлетворяет уравнению Лапласа - 2j(x) = 0, x W, (4) где E(x)= - j(x).
(5) На границах раздела пор и проводника нормальные компоненты плотности тока и напряженности обращаются в нуль (из подобласти ):
jn=0 или En=0, x w. (6) Поставим задачу в области. Применим МКСЭ к решению уравнения для определения потенциала :
- 2j(x) = 0 W R3, в (7) с граничными условиями j Wj j = 0 = 0 W w, Wi Wb на, на = 1 на, = 0 на. (8) n n Wj При этом считаем, что ток через боковую границу куба не выWi текает. На нижней по x3 границе считаем заданным единичный поWb тенциал, на верхней по х3 границе - нулевой.
В области потенциал электрического поля также удовлетворяет уравнению Лапласа (7):
- 2j(x) = 0 xw.
, Для однозначного определения потенциала в диэлектрической подобласти необходимо добавить условие непрерывности потенциала при переходе через границу поры. Потенциал на границе проводника определяется из решения задачи (7Ц8).
В соответствии с МКСЭ расчетная область разбита на KE подобластей - суперэлементов. Подобласти пор в материале расположены строго внутри суперэлементов. Потенциал внутри суперэлемента может быть найден в виде линейной комбинации некоторого числа базисных функций Фi:
n j (х) = Fi(x), a i i=где ai - неизвестные коэффициенты разложения.
Базисные функции являются решениями рассматриваемой системы уравнений внутри суперэлемента, а на его внешней границе совпадают с ji(х):
некоторыми функциями - F (x) = 0,"x W, F ( x) = ji ( x).
i k i xW xW k k С узлов на каждую грань суперэлемента граничные базисные функции продолжены одним из следующих способов:
Х линейная конечно-элементная интерполяция на треугольниках;
Х квадратичная конечно-элементная интерполяция на треугольниках;
Х сокращенная кубическая интерполяция на треугольниках.
Соответствующее число и расположение узлов на каждой из граней показано на рис. 1 - 3. Число таких треугольников может быть различно.
Рис. 1. Линейная интерпо- Рис. 2. Квадратичная Рис. 3. Кубическая инляция на 8 треугольниках интерполяция на 8 тре- терполяция на 8 трена грани суперэлемента угольниках на грани супер- угольниках на грани суэлемента перэлемента Узловые значения ai определяются из уравнений метода БубноваГалеркина в определенном пространстве функций, удовлетворяющих заданной системе уравнений в каждом из суперэлементов в отдельности i j j j (x) F (x)dW = 0,"F, W соответствующая система линейных алгебраических уравнений имеет вид n (Fi,F ) = 0,"F ai j j.
i=Были выполнены расчеты для различных распределений и размеров пор в области. При этом были использованы различные разбиения области на суперэлементы, шаги расчетной сетки, типы граничной интерполяции МКСЭ, разные значения проводимости материала.
Результаты расчетов показывают, что при изменении способа граничной интерполяции МКСЭ, а также числа треугольников, на которых она проведена величина также практически не меняется. Это связано с достаточной гладкостью решения в окрестностях границ суперэлементного разбиения.
Различные способы граничной интерполяции МКСЭ требуют конечного числа суперэлементных узлов в области. Их число совпадает с числом искомых параметров решения.
итература 1. Страховская Л. Г., Федоренко Р. П. Об одном варианте метода конечных элементов // Журнал вычислительной математики и математической физики. - 1979. - Т. 19, № 4. - С. 950Ц960.
2. Галанин М. П., Савенков Е. Б. Метод конечных суперэлементов в задачах математической физики в неоднородных областях // Информационные технологии и вычислительные системы. - 2005. - № 3. - С. 34Ц49.
3. Галанин М. П., Савенков Е. Б. Совместное использование метода конечных элементов и метода конечных суперэлементов // Журнал вычислительной математики и математической физики. - 2006. - Т. 46, № 2. - С. 270Ц283.
ЧИСЛЕННОЕ РЕШЕНИЕ ОДНОМЕРНЫХ УРАВНЕНИЙ БАРОТРОПНЫХ ТЕЧЕНИЙ ДВУХКОМПОНЕНТНЫХ СМЕСЕЙ ВЯЗКИХ СЖИМАЕМЫХ ЖИДКОСТЕЙ Д. А. Прокудин, О. С. Трофимова Кемеровский государственный университет Рассмотрим систему уравнений одномерного движения двухкомпонентных смесей вязких сжимаемых жидкостей [1, 2]:
ri + riui = 0, (x,t) [0;T ], i =1,2, (1) ( ) t x ui ui pi 2ui i+ri + riui + = ni + (-1 a u2 - u1, (x,t) [0;T ], i =1,2, ) ( ) t x x xg где ni, a, T = const > 0, pi =ri, g = const >1, i =1,2, W = x R : 0 < x <1.
{ } Для системы (1) поставим задачу Коши: в цилиндре W 0;T ( ) ищется решение r,u, r = r1,r2, u = u1,u2, удовлетворяющее системе ( ) ( ) ( ) (1), периодическое (с единичным периодом) по переменной x и удовлетворяющее следующим начальным данным:
0 i ri |t=0= ri x, ui |t=0= u0 x, xW, i =1,2, (2) ( ) ( ) где r0 x, i =1,2 - строго положительные и ограниченные функции и ( ) i i r0 x, u0 x, i =1,2 - периодические (с единичным периодом) по про( ) ( ) i странственной переменной x функции.
Целью данной работы является численное решение сформулированной выше задачи.
Алгоритм численного решения задачи (1)Ц(2) основывается на идее расщепления (слабой аппроксимации) [3] уравнений (1) по физическим процессам, т.е. расчет каждого целого шага производится в три этапа. Сначала учитываются те члены системы уравнений (1), которые определяют вязкостные свойства среды, а затем на втором и третьем этапах происходит перерасчет полученных параметров, обусловленных влиянием давлений i -х составляющих смеси и обмена импульсом между компонентами смеси соответственно. Другими словами, на каждом целом шаге [nt,nt+t], n = 0,..., N -1, Nt= T (N - целое) системе уравнений (1) сопоставляются следующие системы:
1 Qi Qi 3 t + ui x = 0, i =1,2, 1 ui ui 2ui 3 R(Qi) t + R(Qi)ui x = niR1(Qi) x2, i =1,2, (3) t nt< t nt+, 1 Qi ui 1 Qi 3 t + = 0, i =1,2, 3 t = 0, i =1,2, x 1 ui Qi ui R(Qi) + = 0, i =1,2, R(Qi ) = R1(Qi )(-1)i+1a(u2 - u1), i =1,2, t x t 3 t 2t 2t nt+ 3 < t nt+ 3, nt + 3 < t nt + t.
Здесь pi ii Qi, ui, Qi = lnri, R(Qi) = c-2 eQ, c2(ri) =, R1(Qi) = R(Qi)e-Q, i =1,2.
( ) ri Поскольку расщепленная задача доставляет приближенное решение задачи (1), то мы линеаризуем системы (3) в пределах каждого расчетного шага [nt,nt+t], сопоставив тем самым задаче (1) следующую задачу (индекс t у функций опущен):
1 Qi Qi 3 t + uin = 0, i =1,2, x 1 ui ui 2ui 3 R(Qin) t + R(Qin)uin = niR1(Qin) x2, i =1,2, (4) x t nt< t nt+, 1 Qi 1 Qi ui 3 t = 0, i =1,2, 3 t + = 0, i =1,2, x 1 ui n+ n+ 1 ui Qi 3 R(Qi ) t = R(Qi 3 ) + = 0, i =1,2, t x 3 n+ t 2t = R1(Qi 3 )(-1)i+1a(u2 - u1 ), i =1,2, nt+ 3 < t nt+ 3, 2t nt+ < t nt+t, Qi |t=0= lnr0, ui |t=0= ui0, i =1,2.
i Здесь Qin = Qi |t=nt, uin = ui |t=nt, i =1,2.
Задаче (1) сопоставим следующую разностную задачу, ассоциированную с расщепленной задачей (4).
В области Wh Wt, Wh = xi = ih :i =1,...,M -1, hM =1, { } t W = tn = nt : n =1,..., N -1, Nt= T ищутся сеточные функции { } Qi = Qih,t, ui = uih,t, i =1,2 такие, что:
1 11 j,n+ j+1,n+ j-1,n+ j+1,n+ j,n+ j-1,n+ 1 Qi 3 - Qij,n Qi 3 - Qi 3 a Qi 3 - 2Qi 3 + Qi =, i = 1,2, 3 D + uij,n 2hh 1 j,n+ j+1,n+ j-1,n+ ui 3 - uij,n ui 3 - ui = 3 R(Qij,n) D + R(Qij,n)uij,n 2h j+1,n+ j,n+ j-1,n+ = niR1(Qij,n) ui - 2ui + ui, i =1,2, h j =1,..., M -1, 2 1 1 j,n+ j,n+ j,n+ j,n+ j+1,n+ j,n+ j-1,n+ 1 Qi 3 - Qi 3 ui 3 - ui 3 g Qi 3 - 2Qi 3 + Qi 3, i =1,2, = 3 + D h 2h 2 1 1 j,n+ uij,n+ 3 - uij,n+3 Qij,n+ 3 - Qij,n+ += 3 R(Qi ) D 2h j+1,n+ j,n+ j-1,n+ j,n+ = g R(Qi 3) ui - 2ui + ui, i =1,2, 3 h j =1,..., M -1, j,n+ 1 Qij,n+1 Qi 3 0, i =1,2, = D j,n+ j,n+ j,n+ 1 uij,n+1 - ui 3 R(Qi ) D = R1(Qi )(-1)i+1a(u2j,n+1 u1j,n+1), i =1,2, j =1,..., M -1, n = 0,..., N -1.
К этим уравнениям присоединим начальные условия:
Pages: | 1 | ... | 7 | 8 | 9 | 10 | 11 | ... | 27 | Книги по разным темам