Компьютерное моделирование динамики загрязнения бассейна реки с учетом запаздывания и случайных факторов Текст научной статьи по специальности «Математика»

Вычислительные технологии

Том 10, № 1, 2005

КОМПЬЮТЕРНОЕ МОДЕЛИРОВАНИЕ ДИНАМИКИ ЗАГРЯЗНЕНИЯ БАССЕЙНА РЕКИ С УЧЕТОМ ЗАПАЗДЫВАНИЯ И СЛУЧАЙНЫХ ФАКТОРОВ*

И. Е. Полосков Пермский государственный университет, Россия e-mail: Igor.Poloskov@psu.ru

A problem of estimation of the characteristics of river water pollution on the basis of stochastic differential equations with multiple delays is considered. Equations for the first moments of random pollution characteristics are derived and solved. Analytical and numerical computations as well as the presentation of results are implementedusing the Mathematica package.

Введение

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

* Работа выполнена при финансовой поддержке Российского фонда фундаментальных исследований (грант № 04-01-96031) и Министерства образования (грант № E02-1.0-151).

© Институт вычислительных технологий Сибирского отделения Российской академии наук, 2005.

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

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

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

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

1. Стохастические дифференциальные уравнения с кратными запаздываниями

Рассматривается система стохастических дифференциально-разностных уравнений вида

Х(£) = /V(ж(£),ж(£ — т),х(£ — 2т), ...,х(£ — vт),£) +

+Си(х(£), х(£ — т),х(£ — 2т),...,х(£ — ут),Ь) £(£), £ > = £0 + ^т. (1.1)

Здесь т — постоянное запаздывание; V > 0 — целое; х € — фазовый вектор; £ €

— вектор независимых гауссовых шумов с единичными интенсивностями; /V(•, •) = {/иг(•, 0}Т : х [£о, го) ^ — детерминированная вектор-функция; GV(•, •) = {д^(•, •)} :

х [£0, го) ^ х Мт — детерминированная матрица-функция; Т — символ транспонирования.

Предположим, что на полуинтервалах (£0,£1], (£1,£2], ..., (£и-1,] фазовый вектор х(£), характеризуемый одноточечной плотностью вероятности р(х, £), удовлетворяет системам

х(£) = /0(х(£),£) + Со (х(Ь),Ь)£(Ь); (1.2о)

х"(£) = ^(х^),^^ — т), £) + С1(х(£), х(£ — т), £) £(£); (1.21)

х(£) = 1и-1(х:(£),х(1 — т), х(£ — 2т),...,х(£ — (V — 1) т), £) +

+Си-1(х(£), х(£ — т), х(£ — 2т),..., х(Ь — (V — 1), т), £) £(£), (1.2и-1)

Л = {/дг}Т, Сд = {дду}, 9 = 1 2,...,V — 1

стохастических дифференциальных уравнений (СДУ). Предположим, что при £ = £0 плотность вероятности х равна р(х, £0) = р0(:г).

Для систем уравнений (1.1) и (1.2д) мы определяем их специальные случаи как линейные системы:

£(£) = Д)(£)ж(£) + с0(£) + Дю(£) £(£),

:г(£) = Р1(£)х(£) + ^п(£)ж(£ — т) + с1(£) + Д1(£) £(£),

••• (1.3)

£(£) = Ри (£)х(£) + ^1(£)х(£ — т) + фи2(£)х(г — 2т) + ...+

(£)х(£ — vт) + Си (£) + Д (£) £(£),

где Рд (£), ^дг (£), Дд (£), сд (£) (д = 0,1, 2,...^, г = 1, 2, ...,д) — любые заданные матрицы-функции и вектор-функции.

Если посмотреть на уравнения (1.1)—(1.3) с точки зрения общей теории случайных процессов, можно сделать вывод о том, что вследствие запаздывания случайные векторы х — решения этих уравнений — не являются марковскими векторными случайными процессами, а следовательно, для получения статистических характеристик векторов х (математических ожиданий, ковариаций и др.) не может быть применен хорошо известный аналитический аппарат теории марковских процессов [9], основанный на использовании уравнения Фоккера — Планка — Колмогорова (ФПК-уравнения) для (переходной) плотности вероятности распределения фазового вектора нелинейной динамической системы.

С другой стороны, естественными выглядят попытки построения на основе уравнений (1.1)—(1.3) иных математических моделей явлений, описываемых данными уравнениями, моделей, более удобных для дальнейших исследований.

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

5 € [0,т], £д = £0 + д • т, д = 0,1, 2,..., хд(в) = ж(з + £д),

Л(в) = £(в + £д), р(ж,,в)= р(жд,5 + £д), Р0(^0, 0) = р0(х0),

Ад = [£д-1,£д ], ^0 = ^0 , ^ = 001(^0, Я^, Л2 = 001(^0 , х, х),..., (1.4)

Уд = хд (0) = хд-1(т), Лд (0) = Лд-1(т), Рд (хд , 0)= Рд-1(хд ,т),

0о1(х0, х1, х2, ...) = {х01, х02, ..., х0п, хц , х12, ..., х1п, х21, х22, ..., х2п, — }Т

и рассматриваем последовательность отрезков Ад.

1. Начнем с А0. Определенный на отрезке А0 случайный вектор х0(в) удовлетворяет системе СДУ (точкой здесь и далее обозначена производная по переменной в)

х0(5) = /0(х0(з), 5 + £0) + С(хэ(5), 5 + ^^(в)

или

х 0(5) = Р)(в + £0)х0(в) + С0(в + £0) + Д (в + ^^(в).

2. Проанализируем поведение системы на отрезках До и Ді. Стохастические дифференциальные уравнения для вычисления вектора со1(Х0,Х1) можно представить в следующем виде:

Хо(5) = /оЫ^), 5 + £о) + Со(хо(^), 5 + ^0)Со(5))

X 1(5) = /і(Хі(з),Хо(з), 5 + іі) + С(Хі(з),Хо(з), 5 + іі)^)

или

X о(5) = Ро(5 + іо)хо(з) + Со(5 + іо) + До (5 + іо)/о(5),

X і(з) = Рі(5 + іі)жі(з) + ^іі(в + іі )жо(з) + Сі(5 + іі) + Лі(з + іі)/і (^).

3. Рассмотрим временные отрезки До, Ді,... , Д^. Построим систему СДУ для вектора /, в виде

Хо(в) = /о(хо(5), 5 + іо) + Со(хо(5), 5 + іо)Со(5),

X і (5) = ,/і(Хі(5),Хо(5),5 + іі ) + Сі(Хі(з),Жо(5),5 + Іі )/і(в),

;^^-і(«) = Л^-і^-і(5),XV-2(5), ...,Хо(«), 5 + і^-і) +

+ -і (XV-1 (в) , XV-2 (5) , ...,Хо(5),5 + ^-і)/^-^),

X* (5) = /V (XV (5) , XV-і (5), ...,Жо(5), 5 + ^ ) + Gv (XV (з),^-і(5), ...,Жо (^), 5 + ^ )/v («)

или

X о(5) = Ро(5 + іо^о^) + Со(5 + іо) + До (5 + іо)/о (5), а; і(й) = Рі(5 + іі)^^) + ^511 (5 + іі ^(й) + Сі(5 + іі) + Ді(5 + іі)/і (5),

X V-1(5) = Р/-1 (5 + Іv-l)Xv-1(5) + ^-іД^ + Іv -1 )Xv—2 (5) + ... + Qv-l)v-l(s + ^-і)^^^

+^-1(5 + Іv—і) + ^-1(5 + ^-і)^-^^,

X V (5) = Р (5 + ^ )/v (5) + ^1(5 + Іv )/v-l(s) + ... + Qvv (5 + ^ )/о (5) +

+Cv (5 + ^ ) + Д (5 + ^ )/v (5).

4. Рассмотрим временные отрезки До, Ді,... , Д^. Построим систему СДУ для вектора /V в виде

;^о(«) = /о^)^ + іо) + Со^о^), 5 + іо)Со(5),

X 1 (5) = /(/l(s),/о(s), 5 + іі) + G(/l(s),/о(s), 5 + іі)/і(5),

Xv-l(s) = Д-і^-і (5),XV-2(5), ..., Xо(s), 5 + Іv-l) +

+ Gv -і (XV-1 (5) , XV-2 (5) , ...^о^), 5 + ^-і)/^-^),

XV(5) = /V(XV(5),XV-1(5), ...^(й), 5 + Іv) + Gv(XV(s),/v-l(s), ...^(й), 5 + Іv^(«), Xv+l(s) = ^^+1 (5), XV(5), ...^(й), 5 + Іv+l) +

или

+ С>^і^)^(5), ...,/1(5), 5 + ^і)/^^),

/V-1(5) = /V ^ -1 ^—2 (5) , ..., ^—V— 1, 5 + іМ -1) +

+ ^(/V-1 (з), /V—2 (5) , ..., /V—V— 1(5) , 5 + ІV-О^-і^^

XV (5) = /V (XV (й),/N-1(5),...,XV-V (5), 5 + ^ ) +

+ С^(5), XV-1(5)..., XV-V(5), 5 + ім)См(5)

Xо(5) = Ро(5 + іо)/о (5) + Со(5 + іо) + До (5 + іо)/о(5),

X і(5) = Рі(5 + іі)/і(5) + 5іі(5 + іі )/о(5) + Сі(5 + іі) + Ді(5 + іі)/і(5),

/v-1(s) = Р/-1(5 + Іv-1)Xv-1(5) + ^-іД^ + і V-і^^^ + ... + Qv-1,v-1(s + і V-1)/о(5) +

+Cv-l(s + іv—і) + ^-1(5 + ^— і)^—1(5),

;с V (5) Р; (5 + іv К (5) + ^і^ + ^ )xv—1(5) + ... + Qvv (5 + іv )/о (5) +

+Cv (5 + Іv ) + Д (5 + Іv )/v (5), v+1(s) = Р/(5 + і V+1)xv+1(s) + ^і^ + іv+1)xv (5) + ... + Qvv (5 + ^+1)/1(3) +

+/v (5 + ^+0 + Д (5 + іv+1)/v+1(s),

х N (в) = Ру (в + £у )хМ (в) + ^и1(в + £у )хМ-1(5) + ... + (в + £у )ху—V (в) +

+Си (5 + £у ) + Ди (в + £У )ЛМ (в).

На основе этой схемы можно построить цепочку ФПК-уравнений для плотностей вероятности векторов Л0, Л1, Л2, ..., Лу, ..., которые принадлежат семейству вложенных фазовых пространств Мп С К2п С К3п С ... С Мп(у+1) С ..., как это было сделано в работе [10], где рассматривались СДУ с одним запаздыванием.

Но в случае линейных СДУ нет необходимости использования данного пути, так как имеется возможность построить уравнения для первых моментов фазовых векторов Лд прямо из исходных СДУ. Так, используя введенные обозначения, определим следующий алгоритм, являющийся стохастическим аналогом классического метода шагов [11].

Рассмотрим линейную систему (1.3) с гауссовой начальной плотностью

р (X) = N(X; т , V ) =

1

л/(2п)5 VI

ехр

—1 (X — то)т (Vе1) 1(/ — то)

Известно, что при таких условиях все плотности рк распределения векторов /о, /1, ..., /V, ... и г+ = со1(/о,Уо), = со1(/1,Уо), ..., /+ = co1(/V,/о), ... при любом 5 будут нормальными

гауссовыми и будут полностью определяться векторами математических ожиданий

гпк (5) = М [/*.]:

/к Рк (Xо,Xl, ...,/к, 5) с/ = Со1(?Пхо ,ГПХ1 ,...,ГПХк )

К" К"

т + (5) = М [/+] = со1(-т,к (5),-то)

и матрицами ковариаций

Vk (5) = М

(/к — тк) (/к — тк

т

V,,

V,,

V V

■^ХхХо -^ХхХх

V V

^Хк Хо ^Хк Х1

V

V.

ХоХк

Х1Хк

V

Хк Хк

где к = 0,1, 2,

Г V ХоХо V Хо Х1 V . ^ХоХк V 1 ^Хоуо

М 3 1 1 к+ т = V Х1 Хо V Х1 Х1 V . ^Х1Хк V ^Х1уо

1 J V ^ХкХо V ^ХкХ1 . V . ^Хк Хк V ^ХкУо

V ■^уохо V ^уоХ1 . V . ^Уо Хк V ^уоуо _

В связи с тем, что вектор гпк (5) и матрица Vk(5) являются блоками вектора т+ (5) и

+к

к

матрицы (в) соответственно, достаточно вычислить только последние, а затем выбрать

их необходимые элементы.

Шаг 0. Числовые характеристики вектора /+ удовлетворяют следующим системам обыкновенных дифференциальных уравнений:

т+ 00 = Ро*+ т+ («) + со+,

Х>+(5) = Р*+ »„+(*) + [Р„*+ ®0+(8)

т

+Д*+ Д*+Т,

т +(0)

т/

т/0

V+ (0)

V0 V0 V0 V0

где

р*+ ____

Ро =

Р* 0 0 0

с*+

со

с*

со

0

Д*+

Д* 0 00

Р*(5) = Ро(5 + іо), с*(5) = Со(5 + іо), Д*(5) = До (5 + іо).

Шаг 1. Вектор средних и матрица дисперсий вектора І+ являются решениями систем обыкновенных дифференциальных уравнений (ОДУ)

т/

+ (5) = Р*+ т +(5) + с1+,

V +(5) = Р*+ V+(s) + [Р*+ V+(s)]т +Д*+ Д*+т

т/ о Vе ^0Х0 (т) Vе

т+(0) = т хо(т ) , V+(0) = ^оуо (т ) Vxоxо (т) VХ0У0 (т)

т/ о Vе VУ0Х0 (т) Vе

где

Р*+

Р* 0 00

—* +

Д*

*+

Д* 0 00

Р*(5) =

Р0(5 + іо) 0

ф 11(5 + і 1) Р1(5 + і 1)

/* (5)

/0(5 + іо) Сі(3 + іі)

1

Я * (в)

Яо(в + £о) 0

0 Я і (в + і і)

Шаг V. Числовые характеристики вектора 5+ найдем из следующих систем ОДУ:

т +(з) = К + т +(5) +К+ >

,т

г>+м

Р„* +

о

; +

+Я* + Я+т,

т +(0) = со1(гп0,?тхо(т),гпЛ1 (т),. . . ,тл„—і (т),то) У+(0) =

1 У о УУожо (т) УУожі (т) . . ^—і (т) У О 1

Ужоуо (т) ^зджо (т) ^жохі (т) . . ^жох^-і (т) Ужоуо (т)

Дхіуо (т) Ужіжо (т) Дхіхі (т) . . £>*і*„-і (т) Ужіуо (т)

^х^-іуо (т) ^ж^-іхо (т) £>**—і*і (т) . . ^ —і— і (т) ®х„—іуо (т)

Уо УУожо (т) УУожі (т) . . ^—і (т) Уо

где

р *+

р р 0 ; с* + = Г с* 1 ; Я* + = о* Я 0

0 0 0 0 0

ро(в + іо) 0 0. .. 0 0

Фіі(в + іі) Рі (в + іі) 0. .. 0 0

Р*(в) = ф22(в + і2) ф2 і (в + і2) Р2(в + і2) . .. 0 0

+ і V ) ф V, V— і (в + і V ) ф V, V—2(в + ^ ) . . Фг/ і(в + і V ) (в + і V )

К (в) = со1 (со (в + іо), Сі (в + і 1), с2(в + І2),. . 4 (в + і V ));

Я* (в) = diag (Яо (в + іо), Я і (в + і і),. . Я^ (в + іг,)).

Шаг V+1. Числовые характеристики вектора і++ і найдем из следующих систем ОДУ:

т

+

(в)

£>++ і(в)

+

Р*+і у++ і (в

р*+іт ++ і(в

с*+ сг/+і,

;і ®++ і («

т

+ я * + я *+т

+Я+і Ягу+ і ,

т ++ і(0) = со1(гто,?тхо(т),гпЛ1 (т), ...,т^(т),то),

У++ і (0)

Уо Ууожо (т) УУожі (т) . . ^Уо^ (т) Уо

Ужоуо (т) Ужохо (т) Ужохі (т) . . Ужох^ (т) Ужоуо (т)

Ужіуо (т) Ужіжо (т) Ужіжі (т) . . Ужіж^ (т) Ужіуо (т)

(т) .о жо (т) жі (т) . . XV (т) (т) .о £

Уо Ууожо (т) УУожі (т) . . ^Уо^ (т) Уо

где

р

+

+

Р* р +

0

0

с*+

сг^+і

С*

+

0

Я + Я +

Р*

Я +

0

0

0

0

P0(s + tQ) О. .О О О

Q l l(s + t і ) Pl(s + tl) . .О О О

Q22 (s + t2) Q2l (s + t2) . .О О О

Qv v (s + t v ) Qv,V— l (s + tv) . . Qv v(s + tv ) Pv (s + i,v) О

О Qv v(s + t v+ l) . . Q v2(s + t v+ l) Q v l(s + ^+l) Pv (s + t v+l)

p;+ і (s)

^+l (s) = co1 (c0 (s + tQ)) cl(s + t l)) c2(s + t2), ...) O,(s + t v ), Cv (s + t l ^;

R*+і(s) = diag(Ro(s + to), Rl(s + t і),..., R(s + t^), R(s + tv+і)).

Шаг N. Числовые характеристики вектора найдем из следующих систем ОДУ:

mN(s) = PN+m N(s) +

где

V+ (s) = PN+ v+ (s) + [PN+ D+ (s)] T+RN+ RN+T,

m

(0) = col(rm0,?mXQ(т),-m,xl (т), ...,mxn_l(т),mQ),

V+ (0)

VQ DyQXQ (т) VyoXl (т) . V . Vyo XN — l (т) VQ

DXQyQ (т) VXQXQ (т ) Vxqxi (т) . V . VXQXN—l (т) Vxoyo (т)

VXlyo (т) VXlXQ (т) DXlXl (т) . V . VXlXN—l (т) Vxlyo (т)

DXN — l yo (т) DXN — lXQ (т) DXN — lXl (т) . V . VXN — lXN — -l (т) VXN—lyo (т)

VQ DyQXQ (т ) DyoXl (т) . V . Vyo XN — l (т) VQ

PN+

PN о

ОО

c*+

CN

C*

N

0

RN+

RN 0

ОО

P0(s + tQ) О О . 0 О

Q l l(s + t l) Pl (s + t l) О . 0 О

PN (s) = Q22 (s + t2) Q 2 l(s + t2) P2 (s + t2) . . 0 О

О О О . Qv l(s + tN ) Pv (s + tN)

CN(s) = col(cQ(s + to), Cl(s + t l),..., cv(s + tv),..., Cv(s + tN));

RN (s) = diag( Rq (s + to), R і (s + t і),..., R v (s + tN),..., Rv (s + tN )).

Далее рассмотренная методология демонстрируется на примере задачи математической экологии.

2. Динамика водного загрязнения

В последние годы значительный интерес проявляется к проблеме оценки и регулирования уровня загрязнения рек. При этом среди многих показателей качества водного потока можно выделить два тесно связанных:

— меру скорости поглощения кислорода водой при разложении органических материалов (BOD);

— количество растворенного в единице объема воды кислорода (DO) [12].

Правомерность выбора последнего показателя заключается в том, что если величина

DO падает ниже определенного уровня, то рыба в реке гибнет.

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

Рассмотрим конечную по длине часть водного пространства реки, в которую из установок для очистки сточных вод осуществляется управляемый сброс очищенных вод. Представим это пространство каскадом проточных водоемов, расположенных вдоль течения реки. Это приводит к тому, что значения BOD и DO водоема, лежащего выше по течению реки, влияют (с некоторым запаздыванием) на величины BOD и DO водоемов, расположенных ниже по течению. Качество воды в i-м водоеме будем определять уровнями BOD (fj, мг/л) и DO (qi, мг/л) в некоторой средней точке этого водоема, который рассматривается как реактор с идеально перемешиваемым содержимым, в результате чего параметры потока и текущие переменные имеют одни и те же значения по всему водоему, а значения Г и qi на выходе водоема эквивалентны аналогичным показателям внутри водоема.

Тогда из рассмотрения баланса масс следуют уравнения [12, 13], определяющие соотношение между BOD и DO в выбранной точке i-го водоема:

где V — объем воды в водоеме, мл; — скорость потока, поступающего в водоем, мл/сут.; К! г — скорость уменьшения BOD за сутки в водоеме; К2г — скорость изменения DO за сутки в водоеме путем реаэрации; — скорость прохождения воды через водоем, мл/сут.; 9/ — уровень насыщения DO для г-го водоема, мг/л; Пг/^г — изменение DO в связи с образованием осадков, мг/(л-сут.); тг — концентрация BOD в поступающем в водоем потоке, мг/л; г и 9 — концентрации BOD и DO, мг/л, определяемые характеристиками воды в лежащих выше по течению водоемах.

М.Б. Бек [13] определил, что для одного из участков р. Кем около г. Кембридж в Англии допустимы следующие значения коэффициентов в уравнениях (2.1), (2.2):

(2.1)

(2.2)

Ки = 0.32, K2i = 0.29, Vi = 1.0, = 4.19,

Vi Vi

t

^ = 10, -тт- = 0.1, — = 0.9.

Таким образом, для г-го водоема уравнения (2.1), (2.2) принимают вид

„.!• |

d ' ri(t) '

dt , qi(t) ,

-1.32 0

-0.32 -1.29

ri(t)

5i(t)

+

0.1

0

mi (t) +

0.9f (t) + 4.19 0.9q(t) + 1.9

(2.3)

Данная модель описывает лишь один водоем р. Кем и предусматривает только одну точку слива очищенных вод. Если таких мест на реке много, а необходимых дополнительных данных нет, то можно считать, что в каскаде многих водоемов каждый из них обладает свойствами модели М.Б. Бека. Именно такую схему моделирования динамики изменения BOD и DO предложил Х. Тамура [14]. Он учел влияние распределения источников загрязнений по реке введением в выражения для fi и f членов с задержками. Согласно модели Х. Тамуры, части a (j = 1, 2,..., s) уровней BOD и DO, характеризующих в момент t — Tj времени (i — 1)-й водоем, появляются в i-м водоеме в момент времени t, т. е. транспортные задержки распределены во времени между ti и Ts. Таким образом, для r(t), q(t) можно записать

fi(t) = a ri- 1 (t - Tj);

j=i

причем

£■

j=i

qi(t) = a qi- i(t - Tj), j=i

1, Tj > 0, Ti < T2 < ... < Ts.

(2.4)

(2.5)

Анализируя эмпирические данные, Х. Тамура для каждой из распределенных задержек нашел следующие значения переменных: в = 3, Т1 = 0, т2 = 0.5 сут., т3 = 1 сут., Го = 0, 5о = 10, а1 = 0.15, а2 = 0.70, а3 = 0.15, которые можно использовать в рассматриваемой модели.

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

Г^1(£) = — 1 .32г1 (£) + 0.1т 1^) + 0.9го(£) + 4.19 + д^^),

(?1(^) = —0.32г1(^) — 1.29^1 (£) + 0.9д0(£) + 1.9,

Г2(^) = 0.9 £ а,п(* — г,-) — 1.32,г:2(£) + 0.1ш2(£) + 4.19 + #26^),

,=1

<&(£) = 0.9 £ а, ^1 (^ — т,) — 0.32г2(£) — 1.29^2(£) + 1.9,

,=1

г^з(£) = 0.9 ^ а, Г2(г — т,) — 1.32гз(^) + 0.Шз(г) + 4.19 + дз£з(£),

,=1

s

s

j

<?з(t) = 0-9^ Оqa(t - Tj) - 0.32гэ(^) - 1.29^э(£) + 1.9 j=i

(g^ = const), причем фазовый вектор системы имеет вид

X = col(ri(t),qi(t),r2(t),q2(t),r3(t),q3(t)) = [xj(t)]T, i = 1, 6.

Для вычисления первых моментов случайных функций Xi(t) был разработан многошаговый численно-аналитический алгоритм, реализованный с помощью пакета Mathematica [15] и состоящий из следующих основных блоков:

— расширение фазового вектора;

— формирование списка новых первых моментов;

— символьное построение системы ОДУ для вектора математических ожиданий и матрицы ковариаций;

— подготовка начальных условий;

— численное интегрирование построенной системы ОДУ;

— сохранение результатов вычислений в конце каждого шага длиной т для дальнейших вычислений и для визуализации результатов.

Расчеты производились со следующими значениями параметров:

m0 = col(l0.0, 7.00, 5.94, 6.0, 5.237, 4.69), D0 = diag(0), gi = g2 = дз = 0.5,

m(t) = col(28.9,12.9,10.4), т = 0.5, v = 2, N =14.

Результаты этих расчетов приведены на рис. 1 (математические ожидания mj(t) функций Xj(t), i = 1, 2, 3,4, 5, 6), на рис. 2 (дисперсии D11(t), D33(t) и D55(t)) и на рис. 3 (D22(t), D44(t) и D66(t)).

Заметим, что на последнем шаге численно интегрировалась система 4656 линейных ОДУ с разреженной матрицей. Общее время вычислений составило около 2 ч на персональном компьютере, имеющем центральный процессор с частотой 1400 МГц. После этих расчетов была предпринята попытка сокращения времени вычислений с помощью генерирования соответствующих Fortran-программ и дальнейшего численного интегрирования ОДУ в среде Compaq Fortran (шаг 5 алгоритма). Время расчетов при этом уменьшилось на 40 %, что важно при многократном использовании созданных программ. Наряду с этим

D

0.12

0.1

0.08

0.06

0.04

0.02

5 ' '

3

1

6 t

6 t

0

1

2

3

4

5

0

1

2

3

4

5

Рис. 1.

Рис. 2.

D

Рис. З.

отметим и тот факт, что все действия по подготовке и проведению численных расчетов (генерирование РоЛгап-программ, компилирование, редактирование, запуск ехе-модуля, считывание из текстового файла результатов расчетов) проводились в автоматическом режиме без вмешательства человека.

Список литературы

[1] Марчук Г.И. Математическое моделирование в проблеме окружающей среды. М.: Наука, 1982.

[2] Знаменский В.А. Гидрологические процессы и их роль в формировании качества воды. Л.: Гидрометеоиздат, 1981.

[3] ГринвАльд Д.И. Турбулентность русловых потоков. Л.: Гидрометеоиздат, 1974.

[4] Картвешвили Н.А. Стохастическая гидрология. Л.: Гидрометеоиздат, 1975.

[5] Гиляров Н.П. Моделирование речных потоков. Л.: Гидрометеоиздат, 1973.

[6] Найденов В.И., Швейкина В.И. Нелинейные модели колебаний речного стока // Водные ресурсы. 2002. Т. 29, № 1. С. 62-67.

[7] Найденов В.И., Крутова Н.М. Исследование многолетних колебаний уровня Каспийского моря на основе теории стохастических дифференциальных уравнений // Там же. 2002. Т. 29, № 3. С. 299-310.

[8] Диментверг М.Ф. Нелинейные стохастические задачи механических колебаний. М.: Наука, 1980.

[9] МАЛАнин В.В., Полосков И.Е. Случайные процессы в нелинейных динамических системах. Аналитические и численные методы исследования. Ижевск: НИЦ “Регулярная и хаотическая динамика”, 2001.

[10] Полосков И.Е. Расширение фазового пространства в задачах анализа дифференциальноразностных систем со случайным входом // Автоматика и телемеханика. 2002. № 9. С. 58-73.

[11] Эльсгольц Л.Э., Норкин С.Б. Введение в теорию дифференциальных уравнений с отклоняющим аргументом. М.: Наука, 1971.

[12] Системы: декомпозиция, оптимизация и управление / Сост. М. Сингх, А. Титли. М.: Машиностроение, 1986.

[13] Beck M.B. The application of control and systems theory to problems of river pollution control / Ph.D. Thesis. Cambridge University, 1974.

[14] TAMURA H. A discrete dynamical model with distributed transport delays and its hierarhical optimization to preserve stream quality // IEEE Trans. SMC. 1974. N 4. P. 424-429.

[15] Wolfram S. The Mathematica Book. 4th ed. Cambridge: University Press, 1999.

Поступила в редакцию 25 декабря 2003 г., в переработанном виде —16 июня 2004 г.