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

Математические заметки СВФУ Апрель—июнь, 2015. Том 22, №2

УДК 519.6

ИССЛЕДОВАНИЕ И ЧИСЛЕННОЕ РЕШЕНИЕ ОДНОЙ ОБРАТНОЙ ЗАДАЧИ МОДЕЛИРОВАНИЯ ЦИРКУЛЯЦИИ В АКВАТОРИЯХ С «ЖИДКИМИ» ГРАНИЦАМИ

В. И. Агошков, Д. С. Гребенников, Т. О. Шелопут

Аннотация. В геофизической гидродинамике известна проблема моделирования физических процессов в водоемах с «жидкими» границами. Одним из подходов к ее решению является использование теории оптимального управления и ассимиляции данных наблюдений. В настоящей работе рассматривается задача о вычислении неизвестной функции в граничном условии для уравнений, относящихся к типу «уравнений мелкой воды». С использованием известных методов исследования и решения обратных задач и задач оптимального управления предложен итерационный алгоритм решения, получены условия однозначной и плотной разрешимости задачи и сходимости итерационного алгоритма. Приведены результаты численной реализации данного алгоритма применительно к акватории Балтийского моря. Ключевые слова: обратные задачи, открытые («жидкие») границы, некорректно поставленные задачи, итерационный алгоритм, уравнения мелкой воды.

V. I. Agoshkov., D. S. Grebennikov, T. O. Sheloput. Investigation and the numerical solution of the problem of modeling of the circulation in water areas with "liquid" boundaries.

Abstract: In geophysical fluid dynamics a problem of modeling physical processes in water areas with "liquid" boundaries is known. One way to solve this problem is to apply the optimal control theory and data assimilation methods. In this paper the problem of unknown function calculation in the boundary condition of the system of "shallow-water equations" is considered. By using the known methods for inverse problems investigation and optimal control theory an iteration algorithm is proposed. Conditions for the unique and dense solvability of the considered problem and conditions for the iteration algorithm convergence are obtained as well. The results of numerical implementation of this algorithm related to Baltic Sea area are also demonstrated in this work.

Keywords: inverse problem, "liquid" (open) boundary, ill-posed problem, iteration process, shallow-water equations.

1. Введение

В задачах геофизической гидродинамики часто возникает проблема моделирования физических процессов в водоемах (морях, океанах, реках и т. д.) с «жидкими» границами. Например, «жидкими» (открытыми) являются южные

Работа выполнена при частичной поддержке РНФ (грант 14—11—00609, в рамках которого был сформулирован итерационный процесс и проведены численные эксперименты) и РФФИ (грант 16—01—00548, в рамках которого задача была сформулирована и исследована).

© 2015 Агошков В. И., Гребенников Д. С., Шелопут Т. О.

границы Индийского океана, северные границы Баренцева, Карского морей, границы, проходящие по проливам, устьям рек и т. д. Данная работа посвящена задаче уточнения граничных функций на «жидких» границах.

Существуют различные приближения, которые можно применить для задания граничных условий на «жидких» границах. Иногда используется приближение «материальной» границы — жидкая граница считается подвижной, и на ней задается условие «непротекания» [1, с. 82-141]. Такое приближение удобно тогда, когда деформация моделируемой области не слишком велика. Однако в данном случае граница является дополнительным неизвестным задачи [2], и использование многих современных численных методов, алгоритмов и инструментов, а также теоретическое исследование, затруднены. Еще одним распространенным приемом является использование осредненных данных о потоках через открытую границу [3]. Иногда возможно провести предварительный расчет по всей акватории Мирового океана на грубой сетке и использовать полученные данные в качестве граничных условий на «жидкой» границе. Вероятно, наиболее перспективным является использование комбинации одного из этих методов с ассимиляцией данных наблюдений.

Идея использования теории оптимального управления и ассимиляции данных наблюдений для решения проблемы «жидких» границ изучена, например, в работах [4-6]. В частности, в [5] предложен и исследован итерационный алгоритм восстановления по данным наблюдений неизвестной граничной функции, учитывающей влияние океана на открытую границу расчетной области, где в качестве модели, описывающей физические процессы в акватории, выбрана система уравнений динамики приливов. Отметим, что итерационные алгоритмы, сформулированные в [4, 5], необходимо реализовывать на каждом шаге по времени.

В настоящей работе проведено исследование вопросов существования и единственности решения обратной задачи о вычислении неизвестной функции в граничном условии для уравнений, относящихся к типу «уравнений мелкой воды», применяющихся для моделирования некоторых типов циркуляций жидкости в акваториях. Суть подхода к исследованию этих вопросов состоит к сведении их к аналогичным вопросам о «граничной функции» для волнового уравнения, к которому (при известных ограничениях) сводится исходная система. При этом само граничное условие для волнового уравнения также строится на основе рассматриваемых уравнений мелкой воды. Также ниже предложен итерационный алгоритм, который в дальнейшем применяется к акватории Балтийского моря. Проведен тестовый эксперимент, в котором «жидкая» граница проходит в районе шведского города Треллеборг и отделяет Северное море от Балтийского. На примере данного эксперимента в работе проиллюстрировано, что предложенный алгоритм является достаточно точным.

§ 2. Постановка задачи

1. Введем следующие обозначения: (ж, у, г) — прямоугольная система координат, (ж, у) € О, где О — ограниченная область в М2, Ь — переменная времени, Ь € [0, Т], Т < то, Qт = О х (0, Т) — цилиндр над О, Г = дО — «кусочно-

гладкая» класса C2 граница области О, удовлетворяющая условию Липшица, Гт = Г х (0, T) — боковая поверхность Qt, Гст = Гс х (0, T), где Гс — «жидкая» граница. Пусть также u, v — скорости жидкости соответственно по осям Ож, Оу, —£(x,y,t) < z < H(ж, у), где z = £(ж, y,t) — уравнение свободной поверхности океана, z = H(ж, у) — уравнение рельефа дна (считаем для простоты H(ж, у) гладкой функцией), g = const — ускорение свободного падения, р — плотность жидкости, pa — атмосферное давление, ti , т2 — напряжения трения ветра, l — параметр Кориолиса.

Запишем осредненную по глубине (по координате z) систему уравнений гидродинамики [7, с. 47]:

f -,v+4x = -lPf*+m,Ti bQt• (1)

dV 11

— +lU + g^ = — Pl + jr- r2 bQt, (2)

dt dy р o y Hpo

6 + (UH)x + (VH)y = 0 в Qt, (3)

где U, V — функции средних скоростей жидкости соответственно по осям Ож, Оу (далее - «скорости»):

H н

11

U = — u dz, V = — v dz;

HH oo

В дальнейшем силой Кориолиса пренебрегаем, т. е. полагаем I = 0. Система дополняется начальными и граничными условиями:

и (ж, у, 0) = Цо(ж, у), V (ж, у, 0) = И>(ж,у), £(ж,у, 0) = £(ж,у) в О, (4)

(И, п) = тсис на Гст, (5)

где и = (и, V)т — вектор скорости, п — вектор внешней нормали, тс — характеристическая функция границы Гст. Введем следующие обозачения:

Н 1 Н 1

/1 =--Р°х + -Т1, /2 =--Рау + -Т2, ^ = ( /1 > /2 ) ^ ■

Ро Ро Ро Ро

Умножая уравнения (1), (2) на Н, дифференцируя первое уравнение по ж, второе — по у, а уравнение (3) по комбинируя их, получаем

Первые два уравнения можно переписать в векторном виде:

в Qт.

[д£/дх\ (fA_d_(UH\ 9 \дЩду) \h) dt\VH)

Рассматривая эти два уравнения на Г и умножая их на вектор внешней к О нормали п, получаем граничное условие вида

дС д

Таким образом, задачу (1)—(5) для уравнений мелкой воды можно переформулировать как следующую задачу для волнового уравнения:

д2С

- сИу (дНЧО = - сНу { в С}т,

С^=о = Со в О дЬ

¿=о дх ду

(6)

дС

дН-^ = ({-п) на (Г\ГС) х (0,Т),

дН= ({ ■ п) - = ({ ■ п) + тсис на Гс х (О, Г),

где ис = —тсИдпс/дг. Предполагается, что выполнены необходимые условия гладкости и согласованности при рассмотрении классической постановки задачи типа (6).

Пусть далее ис есть дополнительное неизвестное на Гс х (0, Т), тогда введем уравнение замыкания:

тоС = то^аЪв на Гт, где то — характеристическая функция границы Гот С Гт, Гот = Г0 х (0,Т), а РоЪя - результаты наблюдений за уровнем С на границе Гот.

2. В дальнейшем будем рассматривать только вещественные переменные, функции и функциональные пространства. Введем следующие гильбертовы пространства (см., например, [8]):

т

Ь2(Ят) : (п,у)Ь2(дт) = (п,у)2,дт = J ! пуйШг,

о о

„.^ . , . ¡' ( дп дь дпдь\

п„ V »=1 /

¿шг,

Ят

W21т = {п : п е W21(Qт), п = 0 при I = Т}. Считаем, что норма в каждом из указанных пространств порождается соответствующим скалярным произведением.

Введем пространство Ио — подпространство пространства Ь2(Гт), состоящее из элементов Ь2(Гт), равных нулю на Гт \ Гот. Также введем пространство Ис — подпространство пространства следов функций из W21(Qт) на Гт, состоящее только из тех элементов, которые равны нулю на Гт \ Гст. Пусть { е (W21(Qт))2, Со е W21(Qт), С1 е ¿2^т), 0 < V < дИ(х), фoЪs е Ио. Сформулируем следующую обратную задачу: найти С е W21(Qт) на Qт и ис е Ис такие что

д2С

- сИу(дНЧО = - сИу { п. в. в С}т, (7)

д£

£|*=о = £(0), -^т =С(1) П. в. в О, (8)

от ¿=о v '

дН-£- = ({ ■ п) п. в. на (Г\ГС) х (О, Г), (9)

дп

дН^ = ({-п)+ис п. в. на Гс х (0,Т), (10)

д п

то£ = то^еЪв п. в. на Гт. (11)

Для обобщенной постановки задачи (7)—(10) умножим (7) скалярно на £ € т(От) и выполним интегрирование по частям с учетом краевых условий. В результате получаем равенство

I) = Р(!) + Ь(ис,|) VI € ж!,т(Зт), (12)

где

i(l,f) = J (-|tlt + gHV|vf) dOdt,

Qt

F(|) = J f •Vf dOdi + J |i|(x,y, 0) dO, b(Uc,|) = J Ucl drdt.

Qt О Гст

Обобщенная постановка (7)—(10) формулируется следующим образом: найти l G W^Qt), удовлетворяющую (12) и условию ||t=o = |о п. в. в O. Теперь вводим обобщенную постановку обратной задачи так: найти | G W^Qt) и Uc G Hc такие, что выполнено соотношение (12), условие ||t=o = |о п. в. в O и дополнительное условие (11).

В последующем изложении задачи типа (7)—(10) понимаются в обобщенной постановке, хотя для наглядности будем часто приводить их классические формы записи вида (7)—(10).

§ 3. Задача оптимального управления

Перейдем к обобщенной постановке задачи (7)—(11), в которой (11) понимается уже в смысле наименьших квадратов: необходимо найти | G W^Qt), Uc G Hc, удовлетворяющие условиям (7)—(10) и решающие задачу минимизации функционала Ja:

Inf Ja(Uc,|(Uc)),

UcbHc

где a > 0 и

Ja{Uc,i{Uc)) = f // mcUc dTdt + i jI ш0(С - tpobsfdTdt. (13)

Гт Гт

Нетрудно показать, что при а > 0 функционал является строго выпуклым и задача минимизации имеет единственное решение. Условие оптимальности = 0 приводит к следующему уравнению:

а JJ тсис¿т + JJ то(£ - = 0, (14)

Гт Гт

где удовлетворяют системе

д2^

д42

- =0 в Qт,

= 0,

щ т

0 в О,

(15)

¿=о

дН—^ = тс511с на Гу. дп

Преобразуем вид равенства (14). Для этого введем сопряженную задачу вида

дЬ2

— а1у(дНУд) = 0 в Qт,

дд д1

0 в О,

<=т

д|<=т = 0,

дН^- = то0(£ - (роЪз) на Гт. д п

(16)

Тогда

0 = // " <М<7ЯУд))

ло [ 9К

СЙ1-

(1С1 + / / д [ ^^ - сИу(дЯУ50 ) ¿ПсИ

Ят

=о

+ Лч (9НШ)

Гт V.-V-' Гт *-V-'

=тс

Следовательно,

то(£ — = дтс¿ис

Гт

Гт

(17)

и (14) принимает вид

а jI тсис^ис + j j дтс^ис = 0.

Гт Гт

(18)

Поскольку ¿ис — независимая вариация, условие оптимальности может быть записано так:

атс ис + тсд = 0 на Гт.

(19)

о

о

Теперь мы можем выписать полную систему вариационных уравнений и сформулировать итерационный процесс получения приближенного решения

обобщенной задачи. Система вариационных уравнений имеет вид д 2£

~ = - сЦу { в (Зт,

= 6 в П,

t=0

д£

gff-^ = (f-n) на (Г\ГС) х (О,Г), д n

д£

ffff^ = (f-n)+t/cHarcx(0,T),

d2q

-ф - dW(gHVq) = 0 в QT,

I п ^

^ = at

= 0 в п,

t=T

дЯтг- = то0(£ - </?obs) на Гт, д n

amcUc + mcq = 0 на Гт.

Прежде чем формулировать итерационный процесс, исследуем обратную задачу на разрешимость.

§ 4. Исследование разрешимости задачи

4.1. Однозначная разрешимость. Перейдем к рассмотрению проблемы однозначной разрешимости задачи (7)—(11).

Предположим, что существуют два решения задачи: £' = £'', U = U'. Тогда для £ = £' — £'', Uc = U — справедливо

д 2£

- div(gHVO = 0 в QT,

«I—»■ I

д£

дЯ^ = mcUc на Г х (О, Т), £ = 0 на Г0 х (0,T).

В случае Гс = Г0 мы можем трактовать задачу (21) как смешанную начально-краевую задачу. По теореме (5.1) из [8] получаем, что система (21) имеет единственное решение £ = 0 в классе W^Qt), следовательно, Uc также равно нулю на Гст.

В случае Гс = Г0 имеем задачу с однородными граничными условиями типа Коши по пространству и времени:

д2£

- div(gHVO = 0 в QT,

= 0 в П, t=T (21)

«I—»■ I

= 0 в П, (22)

t=T

д£

gH= £ = 0 на Г0 х (0,Т). дп

Исследование однозначной разрешимости данной задачи проведено в работе [9] (см., например, теорему 1.2.1 и следствие 1.2.5 на с. 4-10). Не будем приводить здесь условия этих теорем, однако отметим, что предложенные в них достаточные условия единственности решения задачи (22) (назовем их условиями I) включают слишком «жесткие» требования к границе Г0, часто несовместимые с практическими задачами.

4.2. Плотная разрешимость. Перейдем к рассмотрению проблемы плотной разрешимости задачи (7)—(11) (см. [10]).

Как видно из (18), при а = 0 условие оптимальности имеет вид тсд = 0 п. в. на Гст, где д — решение системы (16). В случае Гс = Г0 условия оптимальности имеют вид системы

д2д

- (Ну (дНЧд) = 0вЯт, (23)

I п д(1

=Т = ¥

= 0 в О, (24)

<=т

дНр- = тп0(£-ф0Ьз) наГ0 х (0,Г), (25)

дп

дН^- = 0 на Г\Г0 х (0,Т), (26)

дп

д = 0 на Г0 х (0,Т). (27)

Данная система имеет единственное, причем нулевое обобщенное решение, следовательно, из условия (25) имеем т0(| — ф0ъв) = 0 и минимум функционала ^ при а = 0 тоже нуль, что и означает плотную разрешимость задачи (7)—(11) (см. [10]).

В случае Гс = Г0 имеем задачу с граничными условиями типа Коши на Гс, поэтому для плотной разрешимости необходимо потребовать выполнения дополнительных условий (условия I, см. предыдущий пункт).

Итак, на основании приведенных рассуждений можно сформулировать следующие утверждения.

1. В случае Гс = Г0 задача (7)—(11) однозначно и плотно разрешима.

2. В случае Гс = Г0 и при выполнении условий I имеем однозначную либо плотную разрешимость.

§ 5. Итерационный алгоритм

Поскольку в силу плотной разрешимости задачи (7)—(11) имеем ^ = ^ 0 при а ^ 0, при достаточно малом а > 0 можно предположить, что £ = £(а), ис = ис(а) (где £(а), ис(а) — точное решение задачи минимизации функционала поэтому достаточно построить приближения к £(а), ис(а), например, подходящими итерационными алгоритмами (см. [10]).

Сформулируем метод простой итерации для системы вариационных уравнений (20), который по виду совпадает с методом градиентного спуска для

функционала

а2{к дг2

- а1у(дяУ£к) = - {в ,

= £1 в О,

к

С к=о =

к I

и. гч = /~Г\ —

т

¿=0

д2дк

д£к

дЯ^ = (Г • п) на (Г\ГС) х (0,Т), д п

д£2

- а1у(дяу9к) = о в ,

= 0 в О,

ддк ~д¥

¿=т к

ддк

дН— = то0(£к - </?оЬг!) на Гт,

д п

ик+1 = ик - Тк(аЦк + дк), на Гс х (0, Т).

Здесь Тк — параметр итерационного процесса. От выбора Тк и параметра регуляризации а > 0 зависит сходимость приближенных решений £к(а), (а) к решениям Ц задачи (7)—(11). Так, например, из результатов работы [10] следует, что при любом а > 0 и достаточно малом т = т? итерационный алгоритм (28) сходится.

Заметим также, что согласно теории экстремальных задач параметры итерационного процесса можно выбрать следующим образом [11]:

_ ) - Л* Тк =

,Л>к)112 '

где Ш = Л*. Из плотной разрешимости следует Л* « 0, и можно принять [10]:

,7а(ук) = \\шо(е ~ ФоЬЯЩ2{Гт) П 11</>к)112 4||шс^||22(Гт) 1 ;

— оптимальный набор параметров итерационного процесса в данной задаче.

Как показано выше, на каждой итерации необходимо решать прямую и сопряженную задачи. Для численной реализации этих задач можно, к примеру, использовать проекционно-сеточные методы или методы конечных разностей.

§ 5. Численные эксперименты

Приведем результаты реализации итерационного алгоритма (28) применительно к акватории Балтийского моря. В рамках данной работы были проведены две серии экспериментов: первая — на тестовых функциях (расчет в реальной акватории), вторая — на близких к реальным данных по Балтийскому морю. Целью проведения численных экспериментов на тестовых функциях являлось тестирование и оценка эффективности разработанных программ, анализ сходимости итерационного процесса, оценка относительной ошибки полученного решения. Проведение экспериментов с близкими к реальным данными

L2-norm of the solution error

Рис. 1. Относительная ошибка решения при t = T от числа итераций.

L2-norm of the residual

Рис. 2. Относительная невязка решения при Ь = Т от числа итераций.

дало возможность оценить эффективность разработанного метода и возможность его практического применения. Во всех проведенных экспериментах в качестве «жидкой» была выбрана граница, проходящая в районе шведского города Треллеборг и отделяющая Северное море от Балтийского. Для решения прямой и сопряженной задач использовался метод конечных разностей (см. [8]). Данные о границе Балтийского моря были представлены в виде «маски» из нулей и единиц, а сама граница аппроксимировалась отрезками, параллельными осям координат.

Для проведения исследовательских испытаний программ в качестве тестовой была выбрана функция 8ш(ж/Ь) 8ш(у/Ь) 8ш(£/2Т), по ней была вычислена правая часть, начальные и граничные условия, далее эта функция и граничная функция на «жидкой» границе считались неизвестными и восстанавливались при помощи приведенного выше итерационного алгоритма. На рис. 1 приведен график зависимости относительной ошибки решения от числа итераций при разных параметрах регуляризации а, а на рис. 2 — норма невязки (иначе — квадратный корень от величины функционала Как видно из этих рисун-

ков, алгоритм сходится достаточно быстро (за 20 итераций) и монотонно при больших а, однако ошибка решения при этом велика. При меньших параметрах имеет смысл либо завершать итерационный процесс при 40—60 итерациях,

Рис. 3. Уровень (в см) при Ь = Т — после трех итераций.

1<е 12е 14е 1ве 18е зое 22е 24е 18е 28е зое

-вп -50 -до -30 -90 -1п п 1п 50 10 40

Рис. 4. Уровень (в см), полученный по результатам расчета модели [12].

либо увеличивать точность решения прямой и сопряженной задач (в частности, уменьшать шаг сетки по пространству и времени). Экспериментально показано, что относительная ошибка полученного решения распределена равномерно

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

Из полученных результатов можно сделать следующие выводы. Во-первых, выбор оптимального параметра Tk по формуле (29) обеспечивает быструю сходимость итерационного процесса (20 итераций). Во-вторых, выбор параметра регуляризации зависит от точности решения прямой и сопряженной задач (в частности, от шага по времени и пространству).

Приведем результаты моделирования гидродинамики Балтийского моря с учетом «жидких» границ при исходных данных, близких к реальным. Необходимые для экспериментов данные об атмосферных форсингах были взяты с ресурса ERA-Interim, а начальные данные и наблюдения получены из результатов расчетов трехмерной модели гидротермодинамики Балтийского моря, разработанной в ИВМ РАН [12]. В качестве начальных данных были использованы данные за 1 января 2012 года. Также были выбраны следующие параметры итерационного процесса: Tk по формуле (29), а = 10-3. При таком выборе параметров процесс сходится за 3 итерации (т. е. достаточно быстро), что согласуется с теорией, приведенной в [10]. На рис. 3 представлен уровень (в сантиметрах) в финальный момент времени, полученный на последней итерации. Как видно из рисунка, полученные результаты воспроизводят данные из модели (рис. 4) приемлемо. Колебания уровня, изображенные на рис. 3 в центральной части акватории, вызваны большими перепадами глубины Балтийского моря в данной области. Полученная норма невязки решения (норма разности между полученным и модельным уровнем на «жидкой» границе) на последней итерации составила 6 • 10-4, что позволяет судить о точности предложенных выше алгоритмов решения проблемы уточнения вида граничных условий на «жидких» границах.

Нетрудно заметить, что алгоритмы и подходы, изложенные в настоящей работе, могут быть применены при решении проблем граничных условий на «жидких» границах при рассмотрении циркуляций воды и в других акваториях.

ЛИТЕРАТУРА

1. Гилл А. Динамика атмосферы и океана. М.: Мир, 1986. Т. 1.

2. Agoshkov V. I. Investigation of a class of inverse problems on optimal boundaries // Computational Science for the 21-st Century (Ed. by M.-O. Bristeau, G. Etgen and others). Chichester; New York; Toronto: John Wiley and Sons, 1997. P. 589-598.

3. Чернов И. А., Толстиков А. В. Численное моделирование крупномасштабной динамики Белого моря // Тр. Карел. науч. центра РАН. 2014. № 4. С. 137-142.

4. Дементьева Е. В., Карепова Е. Д., Шайдуров В. В. Восстановление граничной функции по данным наблюдений для задачи распространения поверхностных волн в акватории с открытой границей // Сиб. журн. индустр. математики. 2013. Т. 16, № 1. С. 10-20.

5. Agoshkov V. I. Inverse problems of the mathematical theory of tides: boundary-function problem // Russ. J. Numer. Anal. Math. Modelling. 2005. V. 20, N 1. P. 1-18.

6. Agoshkov V. I. Application of mathematical methods for solving the problem of liquid boundary conditions in hydrodynamics // Numerical Analysis, Scientifc Computing, Computer Science, Special Volume of ZAMM (Proc. ICIAM-95). Berlin, 1996. P.337-338

7. Вольцингер Н. Е., Пясковский Р. В. Основные океанологические задачи теории мелкой воды. Л.: Гидрометеоиздат, 1968.

8. Ладыженская О. А. Краевые задачи математической физики. М.: Наука, 1973.

9. Isakov V. Inverse sourse problems. Providence, RI: Amer. Math. Soc., 1996.

1G. Агошков В. И. Методы оптимального управления и сопряженных уравнений в задачах математической физики. М.: ИВМ РАН, 2003.

11. Васильев Ф. П. Методы решения экстремальных задач. М.: Наука, 1981.

12. Zalesny V. В., Gusev A. V., Chernobay S. Yu., Aps R., Tamsalu R., Kujala P., Rytkönen J. The Baltic Sea circulation modelling and assessment of marine pollution // Russ. J. Numer. Anal. Math. Modelling. 2014. V. 29, N 2. P. 129-138.

Статья поступила 27 августа 2015 г.

Агошков Валерий Иванович Институт вычислительной математики, ул. Губкина, 8, Москва 119333 [email protected]

Гребенников Дмитрий Сергеевич, Шелопут Татьяна Олеговна Московский физико-технический институт (государственный университет),

Институтский пер., 9, г. Долгопрудный Московской обл. 141700 [email protected], [email protected]