CC BY
15Математические заметки СВФУ Июль—сентябрь, 2017. Том 24, № 3
УДК 519.633
ЧИСЛЕННОЕ ВОССТАНОВЛЕНИЕ СТАРШЕГО КОЭФФИЦИЕНТА НЕЛИНЕЙНОГО ПАРАБОЛИЧЕСКОГО УРАВНЕНИЯ Д. Х. Иванов
Аннотация. Численно восстанавливается старший коэффициент параболического уравнения в многомерной области. Рассматривается случай, когда старший коэффициент зависит только от решения, и в качестве дополнительного условия берутся данные в некоторых внутренних точках области. Для численного решения задачи используется метод конечных элементов и библиотека РЕшС8. Приведены примеры идентификации старшего коэффициента двумерного параболического уравнения.
Б01 10.25587/8УРИ.2018.3Л0892
Ключевые слова: задача восстановления коэффициента, параболическое уравнение, метод конечных элементов, библиотека РЕшОЭ.
1. Введение
Работа посвящена численному моделированию коэффициентных обратных задач для нелинейных параболических уравнений. Обратные задачи часто относятся к некорректно поставленным задачам по Адамару, когда нарушено какое-либо из условий: 1) существование решения; 2) единственность решения; 3) непрерывная зависимость решения от входных данных. В математической физике среди обратных задач для уравнений с частными производными выделяют геометрические (неизвестны геометрические размеры объекта), ретроспективные (неизвестно начальное условие) и коэффициентные (неизвестны отдельные части, входящие в определяющее уравнение) [1, 2]. В данной статье рассматривается коэффициентная обратная задача — восстановление старшего коэффициента параболического уравнения в предположении, что этот коэффициент зависит только от самого решения.
В обратных задачах требуются дополнительные условия для точного определения неизвестных величин. Тип условия переопределения может быть разным, так, например, в работе [3] это значения потока в некоторой точке границы, в [4] известны данные на одном из концов отрезка, а в [5] размещался внутри области эталонный образец с известными теплофизическими свойствами и считались известными данные на его границе.
Работа выполнена при финансовой поддержке Российского фонда фундаментальных исследований (код проекта 17-01—00689).
© 2017 Иванов Д. Х.
В настоящей работе в качестве условия переопределения взяты наблюдения за решением во внутренних точках области на всем промежутке времени. Такая задача в одномерном случае была исследована в [2], в [6] рассматривался итерационный численный метод решения коэффициентной обратной задачи теплопроводности, где искомый коэффициент аппроксимировался кубическими сплайн-функциями.
Вычислительный алгоритм решения задачи реализуется с помощью метода конечных элементов на основе вычислительной библиотеки FEniCS [7] и модуля dolfin-adjoint [8, 9], входящего в ее состав.
2. Постановка задачи
Пусть дан прямоугольник О:
б?т=Пх[0,Т], £1 = (0, h) х (0,12), 0 <t<T, T = dn = TDUTN,
где Гд — часть границы с условиями Дирихле, IN — с условиями Неймана. Рассмотрим начально-краевую задачу по определению решения u(x, t) из уравнения
ut — V(k(u)Vu) = f (x,t), x = (xi,x2) € О, 0 <t < T. (1)
Задаются граничные условия Дирихле и Неймана:
u(x, t) = g(x, t), x € ГD, 0 < t < T, (2)
du
—(x,t) = о, жеГдг, о <t<T, (3)
on
а также однородные начальные условия:
u(x, 0) = 0, x € О, (4)
где f (x,t),g(x,t) — заданные функции. Будем считать, что
д(х,0) = 0, ^(x,t)> 0, 0 <t<T,
т. е. g(x, t) — монотонно возрастающая функция по второму аргументу.
В (1) зависимость k(u) неизвестна и ее необходимо найти, считая известными данные во внутренних точках области О:
u(xm,t) = dm(t), 0 < t < T, m = 1, 2,...,M, (5)
dm(t) — точные значения решения при заданном коэффициенте. Принимая во внимание экспериментальную неточность и добавляя искусственный шум, дополнительное условие (5) примем в следующем виде:
u(xm,t) « dm(t)= dm(t) + S(t), 0 < t < T, m = 1, 2,..., M, (6)
где 5(t) € t/f (0, ст) — ошибка, описываемая нормальным распределением с сред-неквадратическим отклонением <г = const.
Задача (1)—(4) нелинейна относительно нахождения старшего коэффициента и, как в большинстве обратных задач, основная трудность ее приближенного решения заключается в возможной неединственности решения и неустойчивости решения относительно малых возмущений данных наблюдения (5). Мы не будем приводить теоретические соображения корректности поставленной задачи, а лишь потребуем от старшего коэффициента k(u) ограниченность и достаточную гладкость. Предположим, что существуют такие числа а и b, что 0 < а < k(u) < b для любого u £ [umin, umax], где
umin = 0, umax = max g(x,t) = max g(x,T).
xerD xerD
te[ü,T ]
Последнее следует из принципа максимума и монотонности функции g(x,t) по времени [2].
Для нахождения коэффициента воспользуемся методом минимизации функционала невязки [1,10]
м t
J(u; k)=Y, I (u(xm, t; k) - dm(t))2 dt (7)
m=10
с учетом погрешности измерений при его минимизации.
3. Алгоритм
Зависимость k(u) будем искать на интервале [umin,umax]. Введем однородную сетку
Wu = {un, n =0,...,N | un = umin + nh, h = (umax - umin) / N}.
Искомый коэффициент будем приближать двумя способом: кусочно постоянными либо кусочно линейными функциями на сетке wu. В первом случае используются кусочно постоянные базисные функции
( 1, если u £ (u,_i,u,), <^(u)= ' lh I = 1, 2,..., N,
0 иначе,
и коэффициент k(u) представляется в виде
N
k(u) = ^2 Cn^n(u). (8)
n=1
Во втором случае используем кусочно линейные базисные функции
f , если и £ (un, ui),
0 иначе,
u — Ui — 2 ( \
--если и £ (щ-2,Щ-1),
Ui — 1 — Ui — 2 '
Uj —u Ui — Ui—1 '
= \ если ъщ), г = 2,3,..., Ж,
иначе,
. " UN 1 . если и £ (идг_1,идг), <pN+1(u) = { 1 W
0 иначе,
так что
N+1
к(и) = ^ с„^„(и). (9)
п=1
Обозначим через N количество базисных функций: N = N в (8), N = N + 1 в (9), и введем вектор управления с, состоящий из весовых множителей базисных функций,
с = (С1,С2,...,С^ )Т.
Таким образом, решение системы (1)—(4) в предположении, что к(и) имеет форму (8) или (9), неявно зависит от вектора с, и функционал (7) можно записать в виде
м т
-/(с) = -(и(х, с); с) = ^ / (и(®т, с) - о!т(4))2 <И. (10)
т=1 о
Следовательно, аргумент минимума J по всем с даст искомый коэффициент к в виде (8) или (9).
4. Вычислительная реализация
Задачу (1)-(4), (6) с теми или иными базисными функциями для к(и) будем численно решать методом конечных элементов на основе программного обеспечения ЕЕшС8 [7]. Минимизацию функционала невязки проведем с помощью дополнительного пакета (АоШп-афоШ [8, 9] с использованием метода минимизации Ь-БЕСЭ-Б [11] из библиотеки 8шРу.
Принимая во внимание конечно-элементную аппроксимацию приближенного решения, от точечных данных перейдем к интегральным. Функционал невязки записывается с учетом (10) в виде
м т г _ _
-"(с) = ^^ / / <5р(х — хт)(и(х,£; с) — вт(£))2 вх <И,
где 5p(x — xm) — гауссово приближение дельта-функции Дирака в окрестности m-й точки наблюдения:
~ Р2
Sp(x - xm) = — exp (-p2 ((xi - жтД)2 + (ж2 - xm^)2)), (11)
так что Sp(x — xm) ^ 5(x — xm) при p ^ ж.
Обратная задача по восстановлению коэффициента диффузии в (1) преобразуется к задаче
min J(c)
c
на решениях
F(u, c) = 0.
m=10 о
Здесь Т соответствует задаче нахождения приближенного решения, которая записана в вариационной формулировке. При использовании неявной схемы для уравнения (10) с шагом т имеем
Г ик _ ик -1 г
/ -ус1х+ / 22сп<Рп(ик~1)\7ик -\7vdx
n=1
^У fvdx Vv G V, k = 1, 2, .. ., T/т (12)
о
где У = {г> (Е Нх(0) : г> = 0, г> (Е Го}, с учетом граничных условий Дирихле и Неймана
ий(ж) = д(ж,кт), х е Гд, к = 1, 2,..., Т/т, (13)
дий
— (ж) = 0, жеГлг, к = 1,2,..., Т/т, (14)
дп
и начального условия
и0(х) = 0, х е О. (15)
5. Примеры
Рассматривается модельная задача, в которой область О является единичным квадратом (11 = 12 = 1) и части границ с условиями Дирихле и Неймана определены следующим образом:
Гд = {(Ж1,Ж2) : Х1 е {0, 1}}, Г№ = {(Ж1,Ж2) : е {0, 1}}.
Функция источника выбирается в виде /(ж,£) = 0, Т = 1, а также
если х1 = 1,
g(x,t) I п п
[ U, если x1 = U.
В расчетах в качестве точного решения берется следующая функция:
k(u) = 0.5 + u. (16)
Будем варьировать число N разбиения интервала [umin,umax], который в этом случае представляет собой единичный интервал, так как
Umin = 0, Umax =max g(x,T) = 1.
xerD
Наряду с этим будем менять число точек наблюдения M.
На рис. 1 показаны данные наблюдения с шумом <г = 0.01 для двух внутренних точек x1 = (0.5, 0.5) и x2 = (0.25, 0.25). Эти данные берутся из решения (12)—(15), но только с известным старшим коэффициентом (16), на регулярной расчетной сетке 33 х 33, введенной на области О, и с шагом т = 2.5 • 10 -3 по времени. Приближенное представление дельта-функции Дирака (11) берется с параметром p = 50. На рис. 2, 3 представлены результаты численных расчетов для различного числа разбиения N. В табл. 1 приведены относительные ошибки в Т2-норме.
2.00 1.75 1.50 1.25 3 1.00 0.75 0.50 0.25
Рис. 1. Данные наблюдения во внутренних точках: 1 — при Х1 = (0.5, 0.5); 2 — при Х2 = (0.25, 0.25)
0.00
2.00 1.75 1.50 1.25 £ 1.00 0.75 0.500.25
0.0
0.2
0.4
0.6
0.8
1.0
0.00
0.0
0.2
0.4
0.6
0.8
1.0
(а)
(б)
Рис. 2. Расчеты при N = 4: (а) — кусочно постоянные базисные функцие; (б) — кусочно линейные базисные функции
(а)
(б)
Рис. 3. Расчеты при N = 10: (а) — кусочно постоянные базисные функции; (б) — кусочно линейные базисные функции
и
и
Таблица 1. Относительная ошибка численного решения: слева кусочно постоянные базисные функции, справа кусочно линейные базисные функции
N M error, %
N M error, %
2 1 15,0
2 2 16,1
4 1 17,1
4 2 16,7
10 1 15, 0
10 2 14, 2
2 1 9, 0
2 2 10,7
4 1 10, 7
4 2 11, 2
10 1 10,8
10 2 15,8
Вблизи правого края имеется падение точности численного приближения искомого коэффициента, что связано с инертностью рассматриваемой эволюционной системы. При увеличении числа разбиения N относительная ошибка приближения кусочно постоянными базисами уменьшается, более того, при увеличении числа точек наблюдения, кроме случая N = 2, точнее определяется к(и). Иначе обстоит дело при кусочно линейной аппроксимации. Так, при увеличении числа N и добавлении еще одной точки наблюдения точность ухудшается. Но относительные ошибки с такими базисами почти всегда меньше, чем с соответствующими кусочно постоянными базисами.
Автор выражает благодарность научному руководителю Вабищевичу Петру Николаевичу за ценные советы при планировании исследования и рекомендации по оформлению статьи.
1. Алифанов О. М., Артюхин Е. А., Румянцев С. В. Экстремальные методы решения некорректных задач и их приложения к обратным задачам теплообмена. М.: Наука, 1988.
2. Самарский А. А., Вабищевич П. Н. Численные методы решения обратных задач математической физики. М.: Едиториал УРСС, 2004.
3. Искендеров А. Д. Об одной обратной задаче для квазилинейных параболических уравнений // Дифференц. уравнения. 1974. Т. 10, № 5. С. 890-898.
4. Музылев Н. В. Теоремы единственности для некоторых обратных задач теплопроводности // Журн. вычисл. математики и мат. физики. 1980. Т. 20, № 2. С. 388-400.
5. Музылев Н. В. О единственности одновременного определения коэффициентов теплопроводности и объемной теплоемкости // Журн. вычисл. математики и мат. физики. 1983. Т. 23, № 1. С. 102-108.
6. Артюхин Е. А. Восстановление температурной зависимости коэффициента теплопроводности из решения обратной задачи // Теплофизика высоких температур. 1981. Т. 19, № 5.
7. Logg A., Mardal K. A., Wells G. (ed.). Automated solution of differential equations by the finite element method: The FEniCS book. New York: Springer Sci. & Business Media, 2012.
8. Farrell P. E., Ham D. A., Funke S. F., Rognes M. E. Automated derivation of the adjoint of high-level transient finite element programs // SIAM J. Sci. Comput. 2013. V. 35, N 4. P. C369-C393.
9. Funke S. W., Farrell P. E. A framework for automated PDE-constrained optimisation // arXiv preprint arXiv:1302.3894. 2013.
10. Васильев Ф. П. Методы решения экстремальных задач. М.: Наука, 1981.
ЛИТЕРАТУРА
С. 963-967.
11. Byrd R. H., Lu P., Nocedal J., Zhu C. A limited memory algorithm for bound constrained optimization // SIAM J. Sci. Comput. 1995. V. 16, N 5. P. 1190-1208.
Статья поступила 26 июня 2017 г. Иванов Дьулус Харлампьевич
Северо-Восточный федеральный университет им. М. К. Аммосова, Институт математики и информатики, ул. Кулаковского, 48, Якутск 677000 i.am.djoos@gmail•com
Математические заметки СВФУ Июль—сентябрь, 2017. Том 24, № 3
UDC 519.633
NUMERICAL RECOVERING OF LEADING COEFFICIENT OF NONLINEAR PARABOLIC EQUATION D. Kh. Ivanov
Abstract: We numerically recover the leading coefficient of one parabolic equation in a multidimensional region. We consider the case when the leading coefficient depends only on solution itself and observations are taken in some interior points of the region as an additional condition. Finite element method implemented by FEniCS library is used for numerical solution of the problem. Several examples of identification of the leading coefficient of a two-dimensional parabolic equation are given.
DOI 10.25587/SVFU.2018.3.10892
Keywords: inverse coefficient problem, parabolic equation, finite element method, FEniCS.
REFERENCES
1. Alifanov O. M., Artyukhin E. A., and Rumyantsev S. V., Extreme Methods for Solving Ill-Posed Problems and Their Applications to Inverse Problems of Heat Transfer [in Russian], Nauka, Moscow (1988).
2. Samarskii A. A. and Vabishchevich P. N., Numerical Methods for Solving Inverse Problems of Mathematical Physics, Walter de Gruyter (2007).
3. Iskenderov A. D., "On an inverse problem for quasilinear parabolic equations [in Russian]", Differ. Equ., 10, No. 5, 890-898 (1974).
4. Muzylev N. V., "Uniqueness theorems for some inverse heat conduction problems [in Russian]", U.S.S.R. Comput. Math. Math. Phys., 20, No. 2, 388-400 (1980).
5. Muzylev N. V., "On the uniqueness of the simultaneous determination of the coefficients of thermal conductivity and the volume heat capacity [in Russian]", U.S.S.R. Comput. Math. Math. Phys., 23, No 1, 102-108 (1983).
6. Artyukhin E. A., "The recovering of the temperature dependence of the thermal conductivity coefficient from the solution of the inverse problem [in Russian]", High Temperature, 19, No 5, 963-967 (1981).
7. Logg A., Mardal K. A., and Wells G. (ed.), Automated Solution of Differential Equations by the Finite Element Method: The FEniCS book, 84, Springer (2012).
8. Farrell P. E., Ham D. A., Funke S. F., and Rognes M. E., "Automated derivation of the adjoint of high-level transient finite element programs", SIAM J. Sci. Comput., 35, No. 4, 369-393 (2013).
9. Funke S. W. and Farrell P. E., "A framework for automated PDE-constrained optimization", arXiv preprint arXiv:1302.3894. (2013).
10. Vasiliev F. P., Methods for Solving Extreme Problems [in Russian], Nauka, Moscow (1981).
The work was supported by RFBR (grant 17-01-00689).
© 2017 D. Kh. Ivanov
11. Byrd R. H., Lu P., Nocedal J., and Zhu C., "A limited memory algorithm for bound constrained optimization", SIAM J. Sci. Comput., 16, No. 5, 1190-1208 (1995).
Submitted June 26, 2017 Dulus Kh. Ivanov
M. K. Ammosov North-Eastern Federal University Institute of mathematics and Informatics 48 Kulakovsky Street, Yakutsk 677891, Russia i.am.djoos@gmail•com
