CC BY
24Математические заметки СВФУ Апрель—июнь, 2014. Том 21, №2
УДК 517.956
ОБ ОПРЕДЕЛЕНИИ ФУНКЦИИ ИСТОЧНИКА В МАТЕМАТИЧЕСКИХ МОДЕЛЯХ КОНВЕКЦИИ-ДИФФУЗИИ С. Г. Пятков, Е. И. Сафонов
Аннотация. Рассматривается задача об определении функции источников в математических моделях массопереноса по данным дополнительных измерений в отдельных точках области. Искомая неизвестная функция является правой частью в параболическом уравнении конвекции-диффузии относительно концентрации переносимого вещества. В частности, эта функция характеризует источники загрязнения и их местоположение. В работе приведены некоторые теоретические результаты, построен численный алгоритм и описаны результаты численных экспериментов. Известные ранее алгоритмы в основном основаны на минимизации некоторого функционала, не являющегося выпуклым. Отличие построенного алгоритма от известных ранее состоит в том, что он является прямым, не требует большого количества вычислений и показывает хорошую сходимость к решению. Численная реализация основана на методе конечных элементов.
Ключевые слова: алгоритм, параболическое уравнение, обратная задача, метод конечных элементов.
Введение
Рассмотрим вопрос об определении, вместе с решением, правой части специального вида (функции источника) в параболическом уравнении. Пусть G — область в Rn с границей Г класса C2 и Q = G х (0, T). Параболическое уравнение имеет вид
ut - A(t,x,D)u = fc, (t,x) G Q, (1)
где A — эллиптический оператор второго порядка вида
n n
A(t,x, D)u ^^ aij(t,x)uXiXj + ai(t,x)uXi + ao(t,x)u.
i,j= 1 i=l
Уравнение (1) дополняется начальными и граничными условиями
u|t=o = uo, Bu|s = g(t, x), (2)
где Bu = u или Bu = f^ + a(x,t)u, или Ви = f + a(x,t)u и S = (0,T) x Г. Таким образом, рассматриваем одно из классических граничных условий Дирихле, Неймана, Робина, или условие с косой производной. Здесь n — единичный
Работа выполнена при финансовой поддержке Российского фонда фундаментальных исследований (код проекта 12-01—00260а).
(g 2014 Пятков С. Г., Сафонов Е. И.
вектор внешней нормали и I = (11(х, £),..., 1п(х, £)) — гладкое некасательное векторное поле на 5. Правая часть в (1) имеет вид
г
/с = £ + /. (3)
г=1
Неизвестными в (1), (2) являются решение и и функции qi(t) (г = 1, 2,. .., г), входящие в правую часть (1). Рассмотрим два вида условий переопределения.
В первом случае условия переопределения имеют вид
= 2 = 1, 2'...'Г' (4)
б;
где ^(х), ф^) — некоторые гладкие функции, условия на них уточним ниже, и Gi С С — некоторые области. Во втором случае рассматриваем условия вида
и(х^)= ф^), жi е С, г = 1, 2 ...,г. (5)
Задача о нахождении функций и, qi с использованием краевых условий и условий переопределения может быть сформулирована и как задача управления. Обратные задачи подобного вида возникают при описании процессов теп-ломассопереноса, диффузионных процессов, процессов фильтрации и во многих других областях [1—4]. В литературе рассматривались условия переопределения как вида (4), так и вида (5). В частности, обратные задачи об определении коэффициентов уравнения (1), зависящих от переменной t, с условием переопределения (4), где г = 1 и Gi = С, рассматривались в [5—11]. Линейные обратные задачи об определении правой части исследовались в [12,13]. Линейные и коэффициентные обратные задачи с условием переопределения (5) рассматривались в [14,15] и в [16-18] соответственно. Большое количество физических постановок и численных методов решения обратных коэффициентных задач и задач об определении функции источника с условиями переопределения вида (5) приведены в [19, 20]. Задачи определения функции источника вида (3) с условиями переопределения вида (5) рассмотрены в [21, гл. 3], где основное внимание уделено численным методам решения обратных задач. В этой монографии рассматривается задача определения произвольной функции источника С(х, исходя из произвольного количества измерений вида (5). При этом функция источника заменяется его приближением вида (3), которое и подлежит определению. Однако, как в этой монографии, так и в других работах, в основном рассматриваются модельные уравнения в случае п =1. Можно отметить работы [22, 23] одного из авторов, где рассмотрены задачи вида (1), (2), (5) в общей постановке. Сошлемся также на монографии [2,5,15,24,25], где имеется большое количество постановок обратных задач для параболических уравнений и систем. Кроме того, для численного решения задачи ранее главным образом использовались численные алгоритмы, основанные на минимизации некоторого функционала, т. е. задача сводилась к задаче оптимального управления (см. [19-21, 26-28]). Этот функционал невыпуклый, и это создавало ряд дополнительных сложностей. Случаев, когда используются другие методы, в литературе описано крайне мало. В частности, в [14] рассматривалась более простая, чем у нас, задача: определялась временная компонента q(t) функции источника q(t)ф(x,y) и задача рассматривалась при более жестких, чем у нас, условиях на правую часть, точнее, требовалось, чтобы ф|г = 0.
В данной работе мы строим численный алгоритм решения задачи (1)—(3), (5) и приводим результаты численных экспериментов. Алгоритм основан на приближении решения задачи (1)—(3), (5) решениями задачи (1)—(4) со специально выбранными функциями р^ = (¿(х, е). Обратная задача (1)—(4) сводится к некоторому интегральному уравнению типа Вольтерра 2-го рода, которое затем решается численно.
В § 1 приводятся теоретические результаты, на которых основывается построение численного алгоритма решения задачи (1)—(4). В § 2 описывается алгоритм решения задачи, в § 3 — программная реализация алгоритма, а в § 4 — результаты численных экспериментов.
§ 1. Основные предположения и вспомогательные результаты
Опишем теоретические результаты, на которых основан наш алгоритм. Обозначения функциональных пространств, используемых ниже, в частности, пространств Соболева и Лебега, стандартны (см. [29]). Фиксируем р > п + 2. Используем следующие условия на данные задачи:
ио(ж) е %2-2/р(С),
я{х,г) е w¡-1/2P'2-1/p(S) или д(х,г) е wp1/2-1/2P'1-1/p(S),
где первое из условий на д используется в случае задачи Дирихле, второе — в случае оставшихся краевых условий:
/ е ьр(<2), (7)
д(х, 0) = В(х, 0)ио(х)\дс. (8)
е Wp1(0,T), ^(0)=У ио(х)^г(х) ¿х, { = 1, 2,..., г, (9)
Gi
фг(1) е Wp1(0,T), ^(0)= ио(хг), г = 1, 2,..., г. (10)
Условия на коэффициенты операторов А и В более или менее стандартны. Для простоты выкладок будем использовать не самые точные условия на коэффициенты. Считаем, что
х) € Ьоо((3), г = 0,1,..., п, йу € С(С}), = \,... ,п,
Ц € С1^), 3 = 1,...,п, а{х,г) € С1^), ^^
Ьг(х,г) е Ьж(0,т; ЬР(О)), г = 1, 2,..., г. (12)
Пусть {О3- } — набор областей с границей класса С1, вложенных в О. Будем использовать два вида условий на весовые функции {р3- (х)}:
зиррр.сС" (1/д + 1/р=1), 3 = 1,2,...,г, (13)
вирр р3- С О3-, р3- е Ь1(О), з = 1, 2,..., г. (14)
При выполнении минимальных требований (14) нам понадобятся дополнитель-
г
ные условия на данные задачи. Пусть Оо = и О3-, <о = Оо х (0,Т). Введем
3 = 1
также вспомогательную область
<г = Ог х (0, Т), Ог = У О], О] = {х е О3 : р(х, <9О3) > 5}.
3=1
Здесь р(х, дО2) обозначает расстояние от х до дО^,
V/ е ьр(Я0), уь3 е ь^(0,т,ьр(Оо)), (15)
г,] = 1,...,п, к = 0,1,...,п.
Здесь и далее V рассматривается только по пространственным переменным. Построим матрицу В с элементами Ьц = Ь2-(£,х^) в случае задачи (1)—(3), (5) и Ьц = § Ь2(х,£)^(х) ¿х в случае задачи (1)-(4). Потребуем, чтобы существовала Gi
постоянная ¿о > 0 такая, что
| В(^)| > ¿о для п. в. I е [0, Т]. (17)
Условие эллиптичности оператора А имеет вид: существует постоянная ¿о > 0 такая, что
п
13 ^ ^ ¿о1^12> М) е ((А е Мп. (18)
Доказательство приведенных ниже теорем 1, 2, 4 может быть найдено в [30], а теоремы 3 — в [22].
Теорема 1. Пусть выполнены условия (6)-(9), (11)-(13), (17), (18). Тогда существует единственное решение (и,д1,..., дг) задачи (1)-(4) такое, что
и е "^(О, Чг(г) е ьр(о,т ), г = 1,2,..., г.
Решение удовлетворяет оценке
г
'Ни^ю) + 13 Н^Нмо.т)
г= 1
< \\/\\ьр(я) + \\я\\тч,2к0{8) + \\по\\ш2-2/Р(^ + ± ) ^ ,
где ко = 1 — 1/2р в случае условия Дирихле и ко = 1/2 — 1/2р в случае остальных краевых условий.
Теорема 2. Пусть выполнены условия (6)-(9), (11), (12), (14)-(18). Тогда существует единственное решение (и, д1,..., дг) задачи (1)-(4) такое, что
и е шр}'2((), цг(ь) е ьр(о,т), г = 1,2,...,г, Vxu е "р1'2^)
для всех ¿ > 0. При фиксированном ¿ > 0 решение удовлетворяет оценке
\\и\\^-2(д) + ) + 13 НиМЬ»(о,т) < с \\/\\ьр(д) + \^х/\\ьр(д0)
1=1 \
+ \\9\\ц?к0^0(8) + 1Ы1 2-2 2-2 + Е 11^11^(0,Т) .
УУр \УР Р(С) \УР Р(С о) р у
Теорема 3. При выполнении условий (6)-(8), (10)-(12), (15)-(18) существует единственное решение задачи (1)—(3), (5) такое, что д € Ьр(0,Т), и €
^р^), Уи € ), 5 > о.
Предположим, что выполнены условия теоремы 3. Пусть е > 0. Рассмотрим функцию р(ж) € С0ю(В1), В1 = {ж € К" : |ж| < 1}, такую, что / р(ж) ¿ж = 1
и р(ж) > 0 для всех ж. Определим
Имеем
= У 1Рзе{х) dx=-^J ж=1.
Введем
г,-® = У (ж)и0(ж) ¿ж - и0(ж3).
о
В силу теорем вложения и условий на функцию ио легко видеть, что найдется постоянная М > 0 такая, что |гje | < Ме для всех Рассмотрим задачу (1)-(4), где = , е < 51, Gj = Б$1 (жj) — шар радиуса 51 с центром в точке Жj, а в качестве фj возьмем функции фje = фj (£) + гje. По построению и в силу условий (10) имеем
Ф,-®(0) = У р3-еио(ж) ¿ж.
о
Параметр 51 выбираем настолько малым, что П G¿ = 0 при г = j и С С для всех ^ Решение этой задачи (1)-(4) (оно существует и обладает свойствами, указанными в теореме 2) обозначим через и®, д® = (д®,..., д®), а решение задачи (1)-(3), (5) — через и, д = (дь .. ., дг).
Теорема 4. Пусть выполнены условия теоремы 3. Тогда
||д® - ^мд) + 1К - и|^-2(д) ^ 0 пРи е ^
Предположим, что выполнены условия теоремы 3. Пусть и(ж, £) — решение задачи (1)-(4), где ^¿(ж) = (ж), а уравнение (1) имеет вид
и - &у(с(ж, £)Уи) + Ь(ж, £)Уи + а(ж, £)и = /, (19)
/ = + /о, Ъ(х,1) = Ых,1)ММ)т, .
¿=1 ^ ' Считаем, что в качестве граничных условий рассматриваются условия Неймана, таким образом,
ди , ,
Соответствующие изменения в рассуждениях в случае других краевых условий очевидны. Интегрируем уравнение (19) по С с весом (ж). Получаем, используя наши данные, что
фjt + У сУи • Ур ¿ж + У (ЬУи + аи)р- ¿ж — J ж)р ¿Г о о г
п
е
В силу условий | ¿е! В| > ¿о для п. в. Ь € [0, Т]. Ввиду определения матрицы В имеем
ч = в-
фи + (сУи, У<х) + (ЬУи, <х) + (аи, <х) — / д<х КГ — (/, <х)
фгъ + (сУи, У<г) + (ЬУи, <г) + (аи, <г) — / д<г КГ — (/, <г).
(20)
Здесь (и, у) обозначает скалярное произведение в пространстве Ь2(О), т. е. (и, у) = / и(х)у(х) 1х. Функция (20) в правой части — некоторый оператор а
Я(д), сопоставляющий вектор-функции д правую часть (20), где и — решение задачи (1)-(3). Тогда на уравнение (20) можно смотреть как на уравнение для нахождения функции д: д = Н(д). Можно показать, что оператор Н(д) при малых Ь сжимающий, значит, справедлива теорема о неподвижной точке. Более того, можно показать, что метод последовательных приближений сходится к решению и на произвольном промежутке времени [0,Т]. Именно это уравнение и будет использоваться при построении численного алгоритма ниже.
§ 2. Описание алгоритма
Опишем численный алгоритм в случае краевых условий Неймана, в случае остальных краевых задач изменения незначительны. Уравнение
щ — &у(с(х, Ь)Уи) + Ь(х, Ь)Уи + а(х, Ь)и = /, (21)
/ = ]Г/г(МЫ^ + /о, Ь(х,Ь) = (Ь1(х,Ь),Ь2(х,Ь))т, ,
¿=х V х 2 /
рассматривается в области Q = О х (0,Т), х = (хх,х2) € О С К2, Ь € [0,Т], с начально-краевыми условиями
ди дп
и|
t=о
: ио(х),
д.
(22)
оах(о,т)
Функции определяются на основе дополнительных измерений вида
и(хг,Ь) = ф(), I = 1, 2,...,г. (23)
В теории тепломассопереноса [1] функция и — концентрация переносимого вещества, Ь — вектор скорости течения, а — коэффициент, определяющий скорость осаждения в результате химических реакций, с — коэффициент диффузии, правая часть / — плотность функции источников.
При построении численного алгоритма обратной задачи (21)-(23) приходится находить старшие частные производные приближенных решений, что не очень удобно, поэтому заменим условие (23) (см. теорему 4) условием вида
/ <р1е (х)и(х,Ь) <1х = фг(Ь),
где функции <е обладают свойствами, приведенными перед теоремой 4, а именно, считаем, что они неотрицательны, имеют носитель, лежащий в области О и стягивающийся к соответствующей точке х1 при е ^ 0, и интеграл по О от каждой из этих функций равен 1.
х
Алгоритм итерационный и основан на методе конечных элементов (МКЭ). Пусть на к-м шаге построено приближение д к = (д^,..., д^?) вектор-функции д = (дх,..., дг), д 0 = 0. Задаем триангуляцию области С, узлы триангуляционной сетки Ж1,Ж2,..., жт и соответствующие базисные кусочно-линейные функции {^¿(ж)} (таким образом, ) = ¿¿^, где — символ Кронекера). Без ограничения общности считаем, что точки замеров жх,ж2,... ,жг являются узлами сетки. Тогда найдутся номера г = 1,2,..., г, такие, что ж^ = ж^, г = 1,2,... , г. Выберем в качестве <р¿е функции где ^ = /¿ж. Такие
1 о
функции удовлетворяют всем нашим условиям (в качестве параметра е можно брать наибольший диаметр треугольника сетки). Приближенное решение (21) ищем в виде
т
ик = Сг(£)рг(ж).
¿=1
Вектор-функция С(£) = (С1(4), С2(£),..., Ст(£)) удовлетворяет системе обыкновенных дифференциальных уравнений
MCt + KC = Fo + G + Fiqk, (24)
i=i
где M и K — матрицы с элементами
Mjj = | pipj dx, Kjj = + bVpipj + apipj) dx,
G G
Fo, G и Fj — векторы с координатами J- /pj dx, J- gpj dx и J- /¿pj dx, ¿,j =
G г G
1, 2,. .., m, соответственно. Пусть
C (t) = (Gi(i),G2(i),...,Gm (t)), C (0) = ((uo(xi),uo(x2),...,uo(xm)).
Чтобы решить систему (24), используем метод конечных разностей (МКР). Заменяем (24) конечно-разностным уравнением
MCl + KCi = Fi, Go = U\t=o, i = 1,2,...,N, At = T/N, (25)
где At — шаг по времени и Fj — значение правой части в (24) в точке ¿At. Таким образом, кусочно-постоянное приближение решения G(t) системы (24) есть кусочно-постоянная функция, равная вектору Gj на множестве ((i — 1)At, ¿At]. В соответствии с равенством (20) (подставляем в правую часть (20) функции pje и используем их определение) для нахождения следующей итерации q k+1 используем равенство
' ^itV1 + (cVufe, Vpj) + (6Vufe + auk, pj) — / gpj dr — (/, pj)'
q k+1 = B-1
где
+ (cVufe, Vpjr) + (bVufe + auk, pj) — / ffPjV dr — (/, pj)
Г (26)
(/, Pji) = (/o, Pji) + E (/j, Pji), i = 1, 2,...,
j=1
r
B
(/ъЫ ••• (/r
.(/1,<Pr) ••• (/r ,<Pr). Отметим, что поскольку x, — внутренние точки области G, а носители функций tfji стягиваются к точкам xi при измельчении сетки, без ограничения общности можно считать, что слагаемые J giPj1 dr в (26) равны нулю. При численной реа-
г
лизации (см. ниже) заменяем функции Ci(t), входящие в определение функции Uk, их кусочно-постоянными приближениями.
§ 3. Программная реализация алгоритма
Реализация алгоритма осуществлялась в среде инженерных и научных расчетов Matlab 2008, полученные функции Matlab компилировались и собирались в dll-библиотеку, которая впоследствии используется в .Net веб-приложении на движке Razor. Серверная и клиентская части программы писались на языке C# в среде разработки Visual Studio 2012 Express. Для наибольшей наглядности трехмерный график функции u и графики функций q(t) выводились при использовании JavaScript-библиотек Three.js и flot.js соответственно. Веб-приложение позволяет пользователю проще вводить данные задачи для его последующего решения.
Программная реализация алгоритма состоит из трех частей.
1. Инициализация массива, описывающего геометрию области и граничного вектора.
В этом пункте по заданным пользователем типу геометрии (эллипс или прямоугольник), параметрам геометрии (для эллипса — радиусы по горизонтали и по вертикали, для прямоугольника — стороны по горизонтали и по вертикали) и по координатам точек замера при помощи функции decsg создаем матрицу декомпозиционной геометрии, содержащую информацию о граничных участках. Используя декомпозиционную геометрию области и условия на границе, задаем граничный вектор.
2. Реализация итерационного расчета методом конечных элементов.
Для реализации МКЭ используем библиотеку Matlab, Partial Differential Equation Toolbox. На первой итерации генерируем триангуляционную сетку Делоне для ранее заданной декомпозиционной геометрии, используя функцию initmesh. Если пользователем задано количество разбиений сетки, то разбиваем сетку нужное количество раз с помощью функции refinemesh. Обозначим через X1,X2,...,Xm узлы триангуляционной сетки, которую строим таким образом, что все точки замеров {xi} являются узлами сетки. Далее определяем индексы этих точек, скажем j1,j2,... ,jr. При помощи функции assemb собираем граничный вектор правой части G, используя функцию assema и assemban, находим матрицу масс M, вектор правой части Fo и матрицу жесткости K. Далее при помощи функции assema находим векторы правой части Fj, j = 1, 2,... ,r.
3. Реализация метода конечных разностей.
Из уравнения (25) получаем
Ci = (M + At • K )-1 • (At • F, + M • Ci-1), i = 1, 2,...,N,
где вектор (Co совпадает с начальным вектором U |t=o и C, — значение кусочно-постоянного приближения вектора C(t) на ((i — 1)At,iAt].
Также на каждой итерации добавляется возмущение функции фjt случайными погрешностями в соответствии с формулой
Фj (t) = Фjt + S(2an — 1).
Здесь an — случайная функция, нормально распределенная на [0,1], а величина S задает уровень погрешности. Положим
ji(t)=(Vi$i(t),...,Vr фг (t)). Таким образом, заменяем (26) следующим равенством:
= B-1(iAt) ■ (C(iAt) + K ■ Ci — F0(iAt)), i = 1, 2,...,N,
где F 0 — вектор с координатами ((fo ,^j1),..., (fo,lPjr)) и qik+1 — значение кусочно-постоянного приближения вектора q k+1(t) на множестве ((i — 1) At, iAt].
Далее проверяем истинность неравенства
max |Cik+1 — Cik | < s,
i
где s — погрешность, задаваемая пользователем, i = 1, 2,... ,r. Если неравенство верно или число k превысит некоторую наперед заданную пользователем величину I, то произойдет выход из программы.
§ 4. Результаты численных экспериментов
В этом параграфе проанализируем результаты численных экспериментов для двух групп входных данных задачи, полученных в результате решения модельных обратных задач в рамках эксперимента на ЭВМ. Характеристики испытуемой ЭВМ следующие: процессор Intel(R) Core(TM) i7-3517U CPU @ 1.90GHz 2.40GHz, ОЗУ 4,00 Гб, 64-разрядная операционная система Windows 7 Ultimate.
В результате расчетов определялось приближенное значение решения (u, q1, q2,... ,qr) задачи (21)—(23) в узлах сетки в моменты времени ti = iAt. Чтобы не загромождать изложение, ниже приведем результаты вычислений только функций qi(t), i = 1, 2, .. . ,r.
Зададим решение и начальные условия для модельной задачи
u = (x2 + y2)(1 + t), u|t=o = x2 + y2, ut = x2 + y2.
Перейдем к обсуждению результатов решения модельных обратных задач по восстановлению зависимости плотности источника от времени в окружности единичного радиуса с центром в точке (0, 0). Будем считать, что f1 = 1, f2 = x, fs = У и дополнительная информация задается в трех точках наблюдения: (0.3, —0.3), (0.1, 0), (—0.5, 0.5). Для численного решения задачи рассматриваем две сетки для описанной области с количеством узлов N1 = 263 и N2 = 1015 (рис. 1).
Все проводимые вычислительные эксперименты разобьем на две группы в зависимости от искомых коэффициентов qi, граничных условий, коэффициентов уравнения a, b, c, d и правой части f. Входные данные для первой группы таковы:
искомые коэффициенты q1 = 1, q2 = t + 2, q3 = t2 + 4;
граничные условия Неймана ^ = 2(i + 1);
коэффициенты уравнения d = 1, c = 1, b1 = 1, b2 = 1, a = 1;
Рис. 1. Сетки: a) N = 263; b) N2 = 1015
правая часть f = (x2 + y2) ■ (t + 1) - 4t + 2x(t + 1) + 2y(t + 1) + x2 + y2 - 5 -x(t + 2) - y(t2 + 4).
Приведем результаты численных экспериментов для первой группы. Сначала сравним результаты вычислений коэффициентов qi для двух сеток без за-шумления данных, 5 = 0 и е = 10-5 (погрешность, задаваемая пользователем). Обозначим через i количество полных итераций работы алгоритма и через т — время, затраченное на вычисление, в секундах. Пусть алгоритм остановился на k-й итерации. В качестве еще одной погрешности работы алгоритма определим
величину ео = max \qik - q(iAt) I, i = 1, 2,..., N. Рассмотрим промежуток [0, T],
i
T =1, положим At = T/N — шаг по времени, N = 100.
Рис. 2. Сравнение результатов вычислений для сеток:
a) сетка N1: т = 24.18, 1 = 16, е0 = 0.109;
b) сетка М2: т = 510.01, 1 = 16, е0 = 0.027
На рис. 2(а) и 2(Ь) приведены графики функций qi (¿) (синие) и графики их приближения (красные). Как можно видеть, они практически сливаются. Далее будем изменять параметры е и 5 так, чтобы описать зависимость ошибки решения ео от уровня погрешности данных 5 и параметра е. Результаты экспериментов представлены на рис. 3.
На рис. 3 видно, что вне зависимости от уровня вводимого возмущения 5 результаты вычислений qi повторяют искомые значения (рис. 3(а,Ь)) или располагаются рядом с ними (рис. 3(е,ё)). Видим, что использование сетки N2 приводит к меньшей погрешности вычислений.
Рис. 3. Результаты вычислений для первой группы данных:
a) сетка N1, е = 10-4, 6 = 0.002, т = 23.29, 1 = 15, ео = 0.112;
b) сетка N2, е = 10-4, 6 = 0.002, т = 477.46, 1 = 14, ео = 0.05;
c) сетка N1, е = 10-5, 6 = 0.01, т = 24.59, 1 = 16, ео = 0.257; а) сетка N2, е = 10-5, 6 = 0.01, т = 530.44, 1 = 16, ео = 0.209
Рис. 4. Результаты вычислений для первой группы данных:
a) сетка N1, е = 10-4, 6 = 0.002, т = 14.22, 1 = 9, ео = 0.078;
b) сетка N, е = 10-4, 6 = 0.002, т = 292.7, 1 = 9, е0 = 0.032; с) сетка №ь е = 10-3, 6 = 0.01, т = 11.8, 1 = 7, е0 = 0.172;
а) сетка N2, е = 10-3, 6 = 0.01, т = 245.66, 1 = 7, ео = 0.177
Приведем результаты численных экспериментов для второй группы дан-
искомые коэффициенты ql = г, д2 = 2г2, дз = 3г3; граничные условия Неймана
ди
дп
= 2;
коэффициенты уравнения 1
d = 1,
правая часть
1 +t
bi =
2(1 +1)'
b2 =
2(1 +1)
1+ t
/ = 3 • (ж2 + у2) - 4 - г - 2х • г2 + 3у • г3. Так же, как и с первой группой, будем изменять параметры е и 6.
Рис. 5. Результаты вычислений для первой группы данных: а) сетка N: е = 10-3, S = 0.01, т = 118.19, i = 7, е0 = 0.27; b) сетка N2: е = 10-3, S = 0.01, т = 2546.07, i = 7, ео = 0.192
Для двух экспериментов уменьшим шаг по времени в 10 раз (At = 0.001). Как видим из результатов экспериментов, в связи с изменением данных время выполнения т для второй группы уменьшилось. Так же заметим, что при уменьшении шага по времени в 10 раз время, затраченное на просчет, увеличивается в 10 раз, а погрешность вычислений не уменьшается.
Подводя итоги, отметим, что эксперименты с использованием сетки N2 показывают лучшую точность, чем с сеткой N1, но выполнение программы с использованием сетки N2 занимает больше времени. Необходимо принять во внимание тот факт, что число узлов сетки N2 в « 3.86 больше N1, а время выполнения программы с сеткой N2 больше в среднем в « 20 раз, чем с сеткой N1. Увеличение переменной е приводит к уменьшению времени вычислений т и количества итераций i, а результаты вычислений становятся менее схожи с искомыми. Увеличение зашумленности данных 5 в 5 раз приводит к увеличению выходной погрешности между искомым решением и найденным приблизительно в 5 раз.
ЛИТЕРАТУРА
1. Алексеев Г. В. Оптимизация в стационарных задачах тепломассопереноса и магнитной гидродинамики. М.: Науч. мир, 2010.
2. Belov Yu. Ya. Inverse problems for parabolic equations. Utrecht: VSP, 2002.
1
x
У
c =
« =
3. Levandowsky M., Childress W. S., Hunter S. H., Spiegel E. A. A mathematical model of pattern formation by swimming microorganisms //J. Protozoology. 1975. V. 22. P. 296—309.
4. Capatina A., Stavre R. A control problem in biconvective flow //J. Math. Kyoto Univ. 1997. V. 37, N 4. P. 585-595.
5. Iskenderov A. D., Akhundov A. Ya. Inverse problem for a linear system of parabolic equations // Dokl. Math. 2009. V. 79, N 1. P. 73-75.
6. Ismailov M. I., Kanca F. Inverse problem of finding the time-dependent coefficient of heat equation from integral overdetermination condition data // Inverse Probl. Sci. Eng. 2012. V. 20, N 24. P. 463-476.
7. Ivanchov M. I. Inverse problem of simultaneous determination of two coefficients in a parabolic equation // Ukr. Math. J. 2000. V. 52, N 3. P. 379-387.
8. Li J., Xu Y. An inverse coefficient problem with nonlinear parabolic equation // J. Appl. Math. Comput. 2010. V. 34. P. 195-206.
9. Kamynin V. L., Franchini E. An inverse problem for a higher-order parabolic equation // Math. Notes. 1998. V. 64, N 5. P. 590-599.
10. Kerimov N. B., Ismailov M. I. An inverse coefficient problem for the heat equation in the case of nonlocal boundary conditions // J. Math. Anal. Appl. 2012. V. 396. P. 546-554.
11. Кожанов А. И. Параболические уравнения с неизвестным коэффициентом, зависящим от времени // Журн. вычисл. математики и мат. физики. 2005. Т. 45, № 12. С. 2168-2184.
12. Криксин Ю. А., Плющев С. Н., Самарская Е. А., Тишкин В. Ф. Обратная задача восстановления плотности источника для уравнения конвекции-диффузии // Мат. моделирование. 1995. Т. 7, № 11. С. 95-108.
13. Васин И. А., Камынин В. Л. Об асимптотическом поведении решений обратных задач для параболических уравнений // Сиб. мат. журн. 1997. Т. 38, № 4. С. 750-766.
14. Калинина Е. А. Численное исследование обратной задачи восстановления плотности источника двумерного нестационарного уравнения конвекции-диффузии // Дальневост. мат. журн. 2004. Т. 5, № 1. С. 89-99.
15. Prilepko A. I., Orlovsky D. G., Vasin I. A. Methods for solving inverse problems in mathematical physics. New York: Marcel Dekker, Inc., 1999.
16. Трянин А. П. Определение коэффициентов теплообмена на входе в пористое тело и внутри него из решения обратной задачи // Инженерно-физ. журн. 1987. Т. 52, № 3. С. 469-475.
17. Dehghan M., Shakeri F. Method of lines solutions of the parabolic inverse problem with an overspecification at a point // Numer. Algorithms. 2009. V. 50, N 4. P. 417-437.
18. Dehghan M. Numerical computation of a control function in a partial differential equation // Appl. Math. Comput. 2004. V. 147. P. 397-408.
19. Алифанов O. M., Артюхов E. A., Ненарокомов A. В. Обратные задачи сложного теплообмена. М.: Янус-К, 2009.
20. Alifanov O. M. Inverse heat transfer problems. Berlin; Heidelberg: Springer-Verl., 1994.
21. Ozisik M. N., Orlando H. A. B. Inverse heat transfer. New York: Taylor & Francis, 2000.
22. Pyatkov S. G., Samkov M. L. On some classes of coefficient inverse problems for parabolic systems of equations // Sib. Adv. Math. 2012. V. 22, N 4. P. 287-302.
23. Pyatkov S. G. On some classes of inverse problems for parabolic equations // J. Inverse Ill-Posed Probl. 2011. V. 18, N 8. P. 917-934.
24. Ivanchov M. Inverse problems for equations of parabolic type. Lviv: VNTL publ., 2003. (Math. Stud. Monogr. Ser.; V. 10).
25. Кабанихин С. И. Обратные и некорректные задачи. Новосибирск: Сиб. науч. изд-во, 2009.
26. Самарский А. А., Вабищевич П. Н. Вычислительная теплопередача. М.: Эдиториал УРСС, 2003.
27. El Badia A., Ha-Duong T., Hamdi A. Identification of a point source in a linear advection-dis-persion-reaction equation: application to a pollution source problem // Inverse Probl. 2005. V. 21. P. 1-17.
28. Панасенко А. Е., Старченко А. В. Численное решение некоторых обратных задач с различными типами источников атмосферного загрязнения // Вестн. ТГУ. 2008. Т. 2, № 3. С. 47-55.
29. Ладыженская О. А., Солонников В. А., Уральцева Н. Н. Линейные и квазилинейные уравнения параболического типа. М.: Наука, 1967.
30. Пятков С. Г., Сафонов Е. И. О некоторых классах линейных обратных задач для параболических систем уравнений // Науч. ведомости Белгород. гос. ун-та. Математика. Физика. 2014. Т. 35, № 7. С. 61-74.
Статья поступила 20 мая 2014 г.
Пятков Сергей Григорьевич, Сафонов Егор Иванович Югорский гос. университет, ул. Чехова 16, Ханты-Мансийск 628012 Б_ру^ко¥@т^гаБи. ги, <1с . gerz. ЬсК^таз.! . сот
