1. На главную
  2. Электронные версии
  3. Журнал вычислительной математики и математической физики
  4. T. 62, № 3, 2022

Журнал вычислительной математики и математической физики, 2022, T. 62, № 3, стр. 451-461

Особенности численного восстановления граничного условия в обратной задаче для уравнения типа реакция-диффузия-адвекция с данными о положении фронта реакции

Р. Л. Аргун 1, А. В. Горбачев 1, Д. В. Лукьяненко 1, 2, *, М. А. Шишленин 3, 4

1 МГУ им. М.В. Ломоносова, физический факультет, кафедра математики
119992 Москва, Ленинские горы, 1, Россия

2 Московский центр фундаментальной и прикладной математики
119234 Москва, Ленинские горы, 1, Россия

3 Институт вычислительной математики и математической геофизики СО РАН
630090 Новосибирск, пр-т Акад. Лаврентьева, 6, Россия

4 Новосибирский государственный университет
630090 Новосибирск, ул. Пирогова, 1, Россия

* E-mail: lukyanenko@physics.msu.ru

Поступила в редакцию 31.03.2021
После доработки 08.04.2021
Принята к публикации 20.05.2021

Полный текст (PDF)

Аннотация

Предлагается новый подход к восстановлению граничного условия в обратной задаче для нелинейного сингулярно возмущенного уравнения типа реакция-диффузия-адвекция с данными о положении фронта реакции. Для решения задачи используется градиентный метод минимизации целевого функционала с выбором начального приближения, поиск которого основан на применении методов асимптотического анализа. Численные эксперименты демонстрируют эффективность предложенного подхода. Библ. 55. Фиг. 5.

Ключевые слова: обратная задача с данными о положении фронта реакции, обратная краевая задача, уравнение типа реакция-диффузия-адвекция.

1. ВВЕДЕНИЕ

В статье рассматриваются особенности численного восстановления граничного условия в обратной задаче для нелинейного сингулярно возмущенного уравнения типа реакция-диффузия-адвекция. В качестве данных обратной задачи используется информация о динамике движения фронта реакции. Задачи для уравнений такого типа возникают в газовой динамике [1], химической кинетике [2]–[8], нелинейной теории волн [9], биофизике [10]–[14], медицине [15]–[18], экологии [19]–[22], финансовой математике [23] и других областях науки [24]. Если в задачах присутствуют разномасштабные процессы, то математические модели таких задач содержат нелинейные параболические уравнения с малыми параметрами при некоторых производных. Как следствие, решения этих задач могут содержать узкие пограничные и/или внутренние слои (стационарные и/или движущиеся фронты).

Часто в постановках обратных задач для уравнений в частных производных используют дополнительную информацию о решении на части границы области (см., например, [25]–[34]). В частности, в работах [35], [36] были рассмотрены особенности численного решения обратных задач для уравнений типа реакция-диффузия-адвекция с данными в финальный момент времени. Однако одной из возможных постановок обратных задач для уравнений такого типа является постановка с дополнительной информацией о динамике движения фронта реакции (см., например, [37]–[39]), который наиболее естественно наблюдать в эксперименте – положение фронта ударной волны, фронта реакции или горения являются легкоразличимыми контрастными структурами.

Данная работа является продолжением работы [40]. В работе [40] был рассмотрен вопрос о возможности применения методов асимптотического анализа для восстановления граничного условия в обратной задаче для нелинейного сингулярно возмущенного уравнения типа реакция-диффузия-адвекция с данными о положении фронта реакции. В указанной работе была использована важная особенность применения методов асимптотического анализа к исследованию нелинейных дифференциальных уравнений с малыми параметрами при старших производных, которая заключается в том, что асимптотический анализ при определенных условиях позволяет свести исходную задачу для нелинейного сингулярно возмущенного уравнения в частных производных к гораздо более простой задаче. Такая задача не содержит малых параметров и имеет меньшую размерность (а иногда и вовсе содержит не дифференциальные, а алгебраические уравнения). В результате применения такого подхода в [40] была получена редуцированная постановка обратной задачи, которая связала явным образом граничное условие, которое необходимо восстановить при решении обратной задачи, с положением движущегося фронта реакции, информация о движении которого наблюдается экспериментально. Было показано, что в случае достаточно малых значений “малого параметра”, входящего в исходное уравнение в частных производных, редуцированная постановка дает решение, близкое к решению обратной задачи в полной постановке. Однако по итогам указанной работы остался открытым вопрос о том, как быть в случае, когда значение “малого параметра” не является достаточно малым. При решении реальных прикладных задач зачастую заранее не известно, можно ли считать “малый параметр”, используемый в постановке задачи, достаточно малым для применения подходов из [40]. Может ли решение, полученное методами из [40], быть уточнено? Другими словами, может ли полученное методами асимптотического анализа приближение решения быть использовано в качестве хорошего начального приближения решения с последующим его уточнением градиентным методом минимизации целевого функционала при решении обратной задачи в полной постановке? Ответу на этот вопрос и посвящена данная работа.

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

2. ПОСТАНОВКА ОБРАТНОЙ ЗАДАЧИ И МЕТОД ЕЕ РЕШЕНИЯ, ОСНОВАННЫЙ НА ПРИМЕНЕНИИ МЕТОДОВ АСИМПТОТИЧЕСКОГО АНАЛИЗА

Рассмотрим прямую задачу для сингулярно возмущенного уравнения типа реакция-диффузия-адвекция

(1)
$\begin{gathered} \varepsilon \left( {\frac{{{{\partial }^{2}}u}}{{\partial {{x}^{2}}}} - \frac{{\partial u}}{{\partial t}}} \right) = A(u,x,t)\frac{{\partial u}}{{\partial x}} + B(u,x,t),\quad x \in (0,1),\quad t \in {{\mathbb{R}}^{1}}, \\ u(0,t) = {{u}_{{{\text{left}}}}}(t),\quad u(1,t) = {{u}_{{{\text{right}}}}}(t) \equiv q(t),\quad t \in {{\mathbb{R}}^{1}}, \\ u(x,t) = u(x,t + T),\quad x \in (0,1),\quad t \in {{\mathbb{R}}^{1}}, \\ \end{gathered} $
где $0 < \varepsilon < 1$ – “малый параметр”, а функции $A(u,x,t)$, $B(u,x,t)$, ${{u}_{{{\text{left}}}}}(t)$ и $q(t)$ являются достаточно гладкими.

Известно [40], что при определенных условиях задача (1) имеет решение в виде движущегося фронта, т.е. решение задачи близко к двум различным функциям ${{\varphi }^{l}}(x,t)$ и ${{\varphi }^{r}}(x,t)$ слева и справа от некоторой точки ${{x}_{{{\text{t}}{\text{.p}}}}}(t)$, а в малой (порядка $\varepsilon \left| {\ln \varepsilon } \right|$) окрестности этой точки наблюдается узкий внутренний переходный слой (см. фиг. 1). Функция $x = {{x}_{{{\text{t}}{\text{.p}}}}}(t)$, описывающая положение фронта реакции, может быть найдена как решение, например, следующего функционального уравнения [40]:

$u(x,t) = \varphi (x,t) \equiv \frac{1}{2}({{\varphi }^{l}}(x,t) + {{\varphi }^{r}}(x,t)),\quad t \in {{\mathbb{R}}^{1}}.$
Фиг. 1.

Типичный вид решения типа движущегося фронта в задаче (1) в фиксированный момент времени $t$.

Таким образом, мы можем сопоставить граничному условию $q(t)$ задачи (1), имеющей решение типа движущегося фронта, положение этого фронта реакции ${{x}_{{{\text{t}}{\text{.p}}}}}(t) \equiv {{f}_{1}}(t)$ в любой момент времени $t \in {{\mathbb{R}}^{1}}$ (см. фиг. 1).

Обратная задача, сформулированная в [40], состояла в определении граничного условия $q(t)$, $t \in [0,T]$, в (1) по известной дополнительной информации о положении фронта реакции

(2)
${{x}_{{{\text{t}}{\text{.p}}}}}(t) = {{f}_{1}}(t),\quad t \in [0,T].$

Отметим, что на практике вместо точных данных ${{f}_{1}}(t)$ известны их приближенные значения ${{f}_{1}}_{{{{\delta }_{1}}}}(t)$, измеренные экспериментально, такие что

${{\left\| {{{f}_{1}} - {{f}_{1}}_{{{{\delta }_{1}}}}} \right\|}_{{{{L}_{2}}}}} \leqslant {{\delta }_{1}}.$

Результатом работы [40] явилась редуцированная постановка обратной задачи (1), (2) в предположении достаточной малости “малого параметра” $\varepsilon $:

(3)
$\int\limits_{{{\varphi }^{l}}({{f}_{1}}(t),t)}^{{{\varphi }^{r}}({{f}_{1}}(t),t)} A(u,{{f}_{1}}(t),t)du = 0,\quad t \in [0,T],$
где функции ${{\varphi }^{l}}(x,t)$ и ${{\varphi }^{r}}(x,t)$ определяются как решения задач
(4)
$\begin{gathered} A({{\varphi }^{l}},x,t)\frac{{d{{\varphi }^{l}}}}{{dx}} + B({{\varphi }^{l}},x,t) = 0,\quad x \in (0,1], \\ {{\varphi }^{l}}(0,t) = {{u}_{{{\text{left}}}}}(t), \\ \end{gathered} $
и

(5)
$\begin{gathered} A({{\varphi }^{r}},x,t)\frac{{d{{\varphi }^{r}}}}{{dx}} + B({{\varphi }^{r}},x,t) = 0,\quad x \in [0,1), \\ {{\varphi }^{r}}(1,t) = q(t). \\ \end{gathered} $

Таким образом, в [40] была получена более простая связь между искомой функцией $q(t)$ и данными обратной задачи, определяющими положение фронта реакции ${{f}_{1}}(t)$. Обратная задача (1), (2) для уравнения в частных производных была заменена обратной задачей для алгебраического относительно восстанавливаемого коэффициента $q(t)$ уравнения (3) и двух обыкновенных дифференциальных уравнений (4) и (5) относительно функций ${{\varphi }^{l}}(x,t)$ и ${{\varphi }^{r}}(x,t)$, которые входят в уравнение (3). Отметим, что редуцированная постановка задачи не содержит “малого параметра”. При этом в [40] было показано, что решение задачи в такой редуцированной постановке дает достаточно точное качественное и количественное описание истинного решения обратной задачи при достаточно малых значениях параметра $\varepsilon $.

2.1. Пример численных расчетов в случае “неудачных” параметров задачи

Рассмотрим эффективность работы предложенного в [40] алгоритма, основанного на численном поиске решения редуцированной постановки задачи (3)–(5), на примере решения обратной задачи (1), (2) для следующего набора модельных параметров:

(6)
$\begin{gathered} \varepsilon {{ = 10}^{{ - 1.5}}},\quad A(u,x,t) = - u,\quad B(u,x,t) = u,\quad T = 2\pi , \\ {{u}_{{{\text{left}}}}}(t) = - 3 + 0.5{{e}^{{ - 4{{{\left( {t - \frac{{2\pi }}{3}} \right)}}^{2}}}}} + 0.5{{e}^{{ - {{{(t - 1.1\pi )}}^{2}}}}}, \\ q(t) \equiv {{q}^{{{\text{model}}}}}(t) = 3 + 0.5{{e}^{{ - 4{{{\left( {t - \frac{{2\pi }}{3}} \right)}}^{2}}}}} + 0.5{{e}^{{ - {{{(t - 1.1\pi )}}^{2}}}}}. \\ \end{gathered} $
Для указанного набора модельных параметров с помощью решения задачи в полной постановке была симулирована функция ${{f}_{1}}(t)$, которая затем была зашумлена с некоторой ошибкой ${{\delta }_{1}}$ (см. фиг. 2). Алгоритм симуляции входных данных ${{f}_{1}}_{{{{\delta }_{1}}}}(t)$ подробно описан в работе [40]. Полученная функция ${{f}_{1}}_{{{{\delta }_{1}}}}(t)$ была использована в алгоритме (описан в [40]), реализующем поиск функции ${{q}^{{{\text{inv}}}}}(t) \equiv q(t)$, $t \in [0,T]$, удовлетворяющей системе (3)–(5) при ${{f}_{1}}: = {{f}_{1}}_{{{{\delta }_{1}}}}$.

Фиг. 2.

Точная модельная функция ${{f}_{1}}(t)$ и зашумленная функция ${{f}_{1}}_{{{{\delta }_{1}}}}(t)$, симулированные для набора модельных параметров (6) с ошибкой ${{\delta }_{1}} = 1.2 \times {{10}^{{ - 3}}}$ (а) и ${{\delta }_{1}} = 3.2 \times {{10}^{{ - 2}}}$ (б), что соответствует уровню относительных ошибок $ \sim 5$ и $ \sim 20\% $ соответственно.

На фиг. 3 представлен результат восстановления функции $q(t)$ по симулированным данным ${{f}_{1}}_{{{{\delta }_{1}}}}(t)$, заданными с относительными ошибками на уровнe 5 и $20\% $. Из результатов расчетов видно, что восстановленная функция ${{q}^{{{\text{inv}}}}}(t)$ достаточно сильно отличается от модельного решения ${{q}^{{{\text{model}}}}}(t)$ как количественно, так и качественно (например, на левом рисунке количество явно различимых локальных экстремумов функций ${{q}^{{{\text{model}}}}}(t)$ и ${{q}^{{{\text{inv}}}}}(t)$ различно). Это связано с тем, что “малый параметр” $\varepsilon $ в случае рассматриваемой задачи для набора параметров (6) не удовлетворяет условиям применимости алгоритма решения, предложенного в [40]. Однако полученное с помощью этого алгоритма решение возможно использовать в качестве хорошего начального приближения при реализации градиентного метода минимизации целевого функционала при решении задачи в полной постановке.

Фиг. 3.

Точное модельное решение (функция ${{q}^{{{\text{model}}}}}(t)$) и результат его восстановления (функция ${{q}^{{{\text{inv}}}}}(t)$) по входным данным обратной задачи, определенных функцией ${{f}_{1}}_{{{{\delta }_{1}}}}(t)$, симулированной для набора модельных параметров (6) с ошибкой ${{\delta }_{1}} = 1.2 \times {{10}^{{ - 3}}}$ (а) и ${{\delta }_{1}} = 3.2 \times {{10}^{{ - 2}}}$ (б).

3. ВИДОИЗМЕНЕННАЯ ПОСТАНОВКА ОБРАТНОЙ ЗАДАЧИ И МЕТОД ЕЕ РЕШЕНИЯ, ОСНОВАННЫЙ НА МИНИМИЗАЦИИ ЦЕЛЕВОГО ФУНКЦИОНАЛА

Рассмотрим ситуацию, в которой значение “малого параметра” $\varepsilon $ в задаче (1) не удовлетворяет условиям применимости методов асимптотического анализа [41] и не позволяет искать решение обратной задачи (1), (2) в виде решения редуцированной постановки задачи (3)–(5). В этом случае видоизменим рассмотренную в [40] постановку обратной задачи следующим образом. Кроме входной информации о положении фронта реакции потребуем наличие и информации о значении функции $u$ на кривой, определяемой траекторией движения фронта. Эта дополнительная информация дает возможность свести решение обратной задачи в полной постановке к поиску решения градиентным методом минимизации целевого функционала. Таким образом, обратная задача с уточненными входными данными может быть сформулирована следующим образом.

Обратная задача состоит в определении граничного условия $q(t)$, $t \in [0,T]$, в (1) по известной дополнительной информации о положении фронта реакции и значении функции $u$ на этом положении фронта (см. фиг. 1)

(7)
${{x}_{{{\text{t}}{\text{.p}}}}}(t) = {{f}_{1}}(t),\quad u({{x}_{{{\text{t}}{\text{.p}}}}}(t),t) = {{f}_{2}}(t),\quad t \in [0,T].$

Отметим, что на практике вместо точных данных ${{f}_{1}}(t)$ и ${{f}_{2}}(t)$ известны их приближенные значения ${{f}_{1}}_{{{{\delta }_{1}}}}$ и ${{f}_{2}}_{{{{\delta }_{2}}}}$, такие что

${{\left\| {{{f}_{1}} - {{f}_{1}}_{{{{\delta }_{1}}}}} \right\|}_{{{{L}_{2}}}}} \leqslant {{\delta }_{1}},\quad {{\left\| {{{f}_{2}} - {{f}_{2}}_{{{{\delta }_{2}}}}} \right\|}_{{{{L}_{2}}}}} \leqslant {{\delta }_{2}}.$

Замечание. Отметим, что априорно известно, что ${{\delta }_{2}} \leqslant {{\left\| {{{\varphi }^{r}}(x,t) - {{\varphi }^{l}}(x,t)} \right\|}_{{{{L}_{2}}}}}$ (см. фиг. 1).

Решение обратной задачи (1), (7) может быть найдено как элемент ${{q}^{{{\text{inv}}}}}(t)$, реализующий минимум функционала А.Н. Тихонова [42]

(8)
$J[q] = \int\limits_0^T {{(u({{f}_{1}}(t),t;q) - {{f}_{2}}(t))}^{2}}dt + \alpha \Omega [q].$
Здесь $\Omega [q]$ – сглаживающий функционал. В данной работе мы используем $\Omega [q] \equiv \left\| q \right\|_{{W_{2}^{2}}}^{2}$:
$\Omega [q] = \int\limits_0^T ({{q}^{2}}(t) + q{\kern 1pt} {{'}^{2}}(t))dt.$
Также в (8): $u(x,t;q)$ – решение прямой задачи (1) для заданной функции $q(t)$, $\alpha $ – параметр регуляризации, который может быть выбран, например, по обобщенному принципу невязки [42].

Замечание. В численных экспериментах можно положить $\alpha = 0$ и использовать номер итерации $s$ как параметр регуляризации [43], [29]. Такой подход будет эквивалентен использованию сглаживающего функционала вида $\Omega [q] \equiv \left\| q \right\|_{{{{L}_{2}}}}^{2}$.

3.1. Алгоритм численного решения обратной задачи (1), (7)

Шаг 1. Зададим $s: = 0$ и ${{q}^{{(0)}}}(t)$, $t \in [0,T]$, в качестве начального приближения.

Шаг 2. Найдем решение ${{g}^{{(s)}}}(x)$ вспомогательной стационарной задачи:

(9)
$\begin{gathered} \varepsilon \frac{{{{d}^{2}}{{g}^{{(s)}}}(x)}}{{d{{x}^{2}}}} = A({{g}^{{(s)}}}(x),x,0)\frac{{d{{g}^{{(s)}}}(x)}}{{dx}} + B({{g}^{{(s)}}}(x),x,0),\quad x \in (0,1), \\ {{g}^{{(s)}}}(0) = {{\left. {{{u}_{{{\text{left}}}}}(t)} \right|}_{{t = 0}}},\quad {{g}^{{(s)}}}(1) = {{q}^{{(s)}}}(0). \\ \end{gathered} $

Шаг 3. Найдем решение ${{u}^{{(s)}}}(x,t)$ прямой задачи:

$\begin{gathered} \varepsilon \left( {\frac{{{{\partial }^{2}}{{u}^{{(s)}}}}}{{\partial {{x}^{2}}}} - \frac{{\partial {{u}^{{(s)}}}}}{{\partial t}}} \right) = A({{u}^{{(s)}}},x,t)\frac{{\partial {{u}^{{(s)}}}}}{{\partial x}} + B({{u}^{{(s)}}},x,t),\quad x \in (0,1),\quad t \in (0,T], \\ {{u}^{{(s)}}}(0,t) = {{u}_{{{\text{left}}}}}(t),\quad {{u}^{{(s)}}}(1,t) = {{q}^{{(s)}}}(t),\quad t \in (0,T], \\ {{u}^{{(s)}}}(x,0) = {{g}^{{(s)}}}(x),\quad x \in [0,1]. \\ \end{gathered} $

Шаг 4. Найдем решение ${{\psi }^{{(s)}}}(x,t)$ сопряженной задачи:

(10)
$\begin{gathered} \varepsilon \frac{{{{\partial }^{2}}{{\psi }^{{(s)}}}}}{{\partial {{x}^{2}}}} + \varepsilon \frac{{\partial {{\psi }^{{(s)}}}}}{{\partial t}} + {{(A({{u}^{{(s)}}},x,t){{\psi }^{{(s)}}})}_{x}} - ({{A}_{u}}({{u}^{{(s)}}},x,t)\frac{{\partial {{u}^{{(s)}}}}}{{\partial x}} + {{B}_{u}}({{u}^{{(s)}}},x,t)){{\psi }^{{(s)}}} + \\ \, + 2\delta (x - {{f}_{1}}(t))({{u}^{{(s)}}}(x,t) - {{f}_{2}}(t)) = 0,\quad x \in (0,1),\quad t \in [0,T), \\ {{\psi }^{{(s)}}}(0,t) = 0,\quad {{\psi }^{{(s)}}}(1,t) = 0,\quad t \in [0,T),\quad {{\psi }^{{(s)}}}(x,T) = 0,\quad x \in [0,1]. \\ \end{gathered} $
Здесь $\delta (x)$ – дельта-функция Дирака.

Шаг 5. Найдем решение ${{w}^{{(s)}}}(x)$ вспомогательной стационарной задачи

(11)
$\begin{gathered} \varepsilon {{\psi }^{{(s)}}}(x,0) + \varepsilon \frac{{{{d}^{2}}{{w}^{{(s)}}}}}{{d{{x}^{2}}}} + {{(A({{g}^{{(s)}}}(x),x,0){{w}^{{(s)}}})}_{x}} - \\ - \;({{A}_{u}}({{g}^{{(s)}}}(x),x,0)\frac{{d{{g}^{{(s)}}}(x)}}{{dx}} + {{B}_{u}}({{g}^{{(s)}}}(x),x,0)){{w}^{{(s)}}} = 0,\quad x \in (0,1), \\ {{w}^{{(s)}}}(0) = 0,\quad {{w}^{{(s)}}}(1) = 0. \\ \end{gathered} $

Шаг 6. Найдем градиент функционала (8):

$J{\kern 1pt} '{\kern 1pt} [{{q}^{{(s)}}}](t) = \left\{ \begin{gathered} - \varepsilon \frac{{d{{w}^{{(s)}}}}}{{dx}}(1) + \alpha \Omega {\kern 1pt} '{\kern 1pt} [{{q}^{{(s)}}}](t)\quad {\text{для}}\quad t = 0, \hfill \\ - \varepsilon \frac{{\partial {{\psi }^{{(s)}}}}}{{\partial x}}(1,t) + \alpha \Omega {\kern 1pt} '{\kern 1pt} [{{q}^{{(s)}}}](t)\quad {\text{для}}\quad t \in (0,T]. \hfill \\ \end{gathered} \right.$
Производная $\Omega {\kern 1pt} '{\kern 1pt} [{{q}^{{(s)}}}](t)$ сглаживающего функционала может быть вычислена либо аналитически, либо численно.

Шаг 7. Найдем приближенное решение на следующем шаге итерации:

${{q}^{{(s + 1)}}}(t) = {{q}^{{(s)}}}(t) - {{\beta }_{s}}J{\kern 1pt} '{\kern 1pt} [{{q}^{{(s)}}}](t),$
где ${{\beta }_{s}}$ – параметр спуска.

Шаг 8. Проверим условие остановки итерационного процесса. Если оно выполняется, положим ${{q}^{{{\text{inv}}}}}(t): = {{q}^{{(s + 1)}}}(t)$ решением обратной задачи. В противном случае, зададим $s: = s + 1$ и перейдем к ш. 2.

(a) В случае данных, измеренных в эксперименте с ошибками ${{\delta }_{1}}$ и ${{\delta }_{2}}$, критерий остановки:

$\int\limits_0^T [({{x}_{{{\text{t}}{\text{.p}}}}}(t;{{q}^{{(s + 1)}}}(t)) - {{f}_{1}}_{{{{\delta }_{1}}}}(t{{))}^{2}}dt + {{(u({{f}_{1}}_{{{{\delta }_{1}}}}(t),t;{{q}^{{(s + 1)}}}(t)) - {{f}_{2}}_{{{{\delta }_{2}}}}(t))}^{2}}]dt \leqslant \delta _{1}^{2} + \delta _{2}^{2}.$
Здесь ${{x}_{{{\text{t}}{\text{.p}}}}}(t;q)$ – положение фронта реакции, определяемого при решении прямой задачи (1) для заданной функции $q(t)$.

(б) В случае точных входных данных итерационный процесс останавливается, когда $J[{{q}^{{(s)}}}]$ меньше оценки ошибки конечно-разностной аппроксимации.

Наличие априорной информации об искомом решении позволяет существенно снизить число итераций при реализации алгоритма [44]. В качестве такой априорной информации будем использовать функцию ${{q}^{{(0)}}}(t)$, полученную как решение обратной задачи в редуцированной постановке (3)–(5).

3.2. Примеры численных расчетов

Рассмотрим эффективность работы предложенного алгоритма, основанного на минимизации целевого функционала (8), на примере решения обратной задачи (1), (7) для модельного набора параметров (6) и функции ${{f}_{1}}_{{{{\delta }_{1}}}}(t)$, заданной с таким же уровнем шума, как и в примере разд. 2 (см. фиг. 2). В качестве начального приближения ${{q}^{{(0)}}}(t)$ будем использовать решение задачи в редуцированной постановке (3)–(5).

Для решения прямой задачи (1) (ш. 3 алгоритма) и сопряженной задачи (10) (ш. 4 алгоритма) использовался жесткий метод прямых SMOL [45], позволяющий свести систему для уравнения в частных производных к системе обыкновенных дифференциальных уравнений. Полученные системы дифференциальных уравнений решались с помощью одностадийной схемы Розенброка с комплексным коэффициентом CROS1 [46], [47]. Подробное описание деталей численного решения подобных задач может быть найдено, например, в [35], [40]. В частности, в [35] описывается построение соответствующих численных схем с использованием динамически адаптированных сеток. Вспомогательные стационарные задачи (9) и (11) (ш. 2 и ш. 5 алгоритма) решались методом счета на установление. При реализации градиентного метода использовались сетки с $N = 300$ и $M = 300$ интервалами по пространственной и временной переменным соответственно, параметр спуска ${{\beta }_{s}} = 5 \times {{10}^{{ - 6}}}$, параметр регуляризации $\alpha = 20$. При решении сопряженной задачи (10) (ш. 3 алгоритма) использовалась следующая формула для аппроксимации $\delta $-функции [48]:

$\delta (x) = \left\{ \begin{gathered} \frac{1}{{2\omega }}\left( {1 + \cos \left( {\frac{{\pi x}}{\omega }} \right)} \right),\quad \left| {\frac{x}{\omega }} \right| \leqslant 1, \hfill \\ 0,\quad \left| {\frac{x}{\omega }} \right| > 1, \hfill \\ \end{gathered} \right.$
с размером носителя дискретной $\delta $-функции: $\omega {{ = 10}^{{ - 2}}}$. Для вычисления интеграла (ш. 8 алгоритма) использовалась формула трапеций.

На фиг. 4 представлен результат восстановления функции ${{q}^{{{\text{inv}}}}}(x)$ в случае старта с “хорошего” начального приближения ${{q}^{{(0)}}}(t)$. Полученное решение качественно и количественно уже не сильно отличается от истинного, особенно с учетом того, что оно было восстановлено по данным ${{f}_{1}}_{{{{\delta }_{1}}}}(t)$, заданным с ошибками.

Фиг. 4.

Результат восстановления функции ${{q}^{{{\text{inv}}}}}(x)$ с помощью градиентного метода минимизации целевого функционала (8). В качестве входных данных обратной задачи использовались функции ${{f}_{1}}_{{{{\delta }_{1}}}}(t)$ и ${{f}_{2}}_{{{{\delta }_{2}}}}(t)$, симулированные для набора модельных параметров (6) с ошибками $\{ {{\delta }_{1}},{{\delta }_{2}}\} = \{ 1.2 \times {{10}^{{ - 3}}},0\} $ (а) и $\{ {{\delta }_{1}},{{\delta }_{2}}\} = \{ 3.2 \times {{10}^{{ - 2}}},0\} $ (б). Старт градиентного метода осуществлялся с “хорошего” начального приближения ${{q}^{{(0)}}}(t)$ (совпадающего в случае выбранных ошибок с изображенными на фиг. 2).

Далее проведем набор расчетов для модельного набора параметров (6), но для различных значений $\varepsilon $ и разного уровня погрешности задания входных данных ${{\delta }_{1}}$. В связи с тем, что численные эксперименты продемонстрировали малую чувствительность алгоритма к ошибкам ${{\delta }_{2}}$ функции ${{f}_{2}}_{{{{\delta }_{2}}}}(t)$, во всех экспериментах мы используем ${{\delta }_{2}} = 0$. На фиг. 4 представлена зависимость ${{\left\| {{{q}^{{{\text{inv}}}}}(t) - {{q}^{{(0)}}}(t)} \right\|}_{{{{L}_{2}}}}}$ от относительной погрешности задания входных данных, определяемой по формуле ${{\delta }_{1}}{\text{/}}\left\| {{{f}_{1}}} \right\| \times 100\% $, для различных значений $\varepsilon $. Здесь функция ${{q}^{{(0)}}}(t)$ является решением задачи в редуцированной постановке (3)–(5), а функция ${{q}^{{{\text{inv}}}}}(t)$ – решением обратной задачи в полной постановке (1), (7) градиентным методом минимизации целевого функционала (8). Смысл каждого из графиков – величина расхождения между начальным приближением, выбранным исходя из методов асимптотического анализа, и решением обратной задачи в полной постановке, а именно, чем меньше величина соответствующей нормы разности, тем ближе начальное приближение лежит к истинному решению. Из полученных результатов, представленных на фигуре, можно сделать вывод, что при достаточно малых значениях $\varepsilon $ (например, для $\varepsilon {{ = 10}^{{ - 2.5}}}$) решение, найденное из редуцированной постановки задачи, достаточно близко к решению рассматриваемой обратной задачи в полной постановке (что, в частности, подтверждает основной результат работы [40]). Однако с увеличением $\varepsilon $ метод, основанный на минимизации целевого функционала (8), позволяет существенным образом уточнить начальное приближение, которое является лишь грубым приближением истинного решения. Более того, с увеличением $\varepsilon $ уменьшается зависимость этой разницы от погрешности ${{\delta }_{1}}$ задания входных данных. В случае достаточно малого значения $\varepsilon $, когда в случае данных без ошибки алгоритм из [40] дает решение, мало отличающееся от точного, при наличии ошибок в экспериментальных данных градиентный метод реализует процедуру сглаживания.

Фиг. 5.

Зависимость ${{\left\| {{{q}^{{{\text{inv}}}}} - {{q}^{{(0)}}}} \right\|}_{{{{L}_{2}}}}}$ от погрешности ${{\delta }_{1}}$ задания входных данных ${{f}_{1}}_{{{{\delta }_{1}}}}$ для различных значений малого параметра $\varepsilon $.

4. ЗАКЛЮЧЕНИЕ

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

1. При построении “хорошего” начального приближения ${{q}^{{(0)}}}(t)$ с помощью алгоритма, предложенного в работе [40], ошибка, определяющая разность между этим приближением и истинным решением, оценивалась лишь численно. Вопрос о возможности выполнения строгих аналитических оценок остается открытым и представляет существенный интерес как тема отдельной работы.

2. Вопрос о строгих условиях, при которых “хорошее” начальное приближение гарантированно лежит в окрестности глобального минимума целевого функционала, является открытым. Соответствующий вопрос также представляет существенный интерес и может являться темой отдельной работы, посвященной свойствам глобальной сходимости [49]–[51]. Отметим, что под глобально сходящимся алгоритмом подразумевается алгоритм, позволяющий найти глобальный минимум вне зависимости от выбора начального приближения. В том случае, если “хорошее” начальное приближение, выбираемое автоматически с помощью алгоритма из [40], лежит в окрестности глобального минимума, алгоритм будет обладать свойствами глобальной сходимости.

3. Для повышения эффективности численного счета при реализации градиентного метода минимизации целевого функционала, при решении прямой, сопряженной и вспомогательных стационарных задач можно использовать так называемые динамически адаптированные сетки. Техника построения таких сеток подробно описана в работе [35]. В связи с тем, что использование таких сеток не влияет на результаты этой работы, мы опустили описание их возможного применения.

4. В качестве перспектив развития предложенного метода следует отметить реализацию методов выполнения апостериорной оценки точности полученного решения [52]–[55].

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

  1. Danilov V.G., Maslov V.P., Volosov K.A. Mathematical modelling of heat and mass transfer processes. Kluwer, Dordrecht, 1995.

  2. Butuzov V.F., Vasil’eva A.B. Singularly perturbed problems with boundary and interior layers: theory and applications // Advances in Chemical Physics. 1997. V. 97. P. 47–179.

  3. Liu Z., Liu Q., Lin H.-C., Schwartz C. S., Lee Y.-H., Wang T. Three-dimensional variational assimilation of MODIS aerosol optical depth: implementation and application to a dust storm over East Asia // Journal of Geophysical Research: Atmospheres. 2010. V. 116. № D23.

  4. Egger H., Fellner K., Pietschmann J.-F., Tang B.Q. Analysis and numerical solution of coupled volume-surface reaction-diffusion systems with application to cell biology // Applied Mathematics and Computation. 2018. V. 336. P. 351–367.

  5. Yaparova N.M. Method for determining particle growth dynamics in a two-component alloy // Steel in Translation. 2020. V. 50. № 2. P. 95–99.

  6. Wu X., Ni M. Existence and stability of periodic contrast structure in reaction-advection-diffusion equation with discontinuous reactive and convective terms // Communications in Nonlinear Science and Numerical Simulation. 2020. V. 91. P. 105457.

  7. Lin G., Zhang Y., Cheng X., Gulliksson M., Forssen P., Fornstedt T. A regularizing Kohn–Vogelius formulation for the model-free adsorption isotherm estimation problem in chromatography // Applicable Analysis. 2018. V. 97. № 1. P. 13–40.

  8. Zhang Y., Lin G., Gulliksson M., Forssen P., Fornstedt T., Cheng X. An adjoint method in inverse problems of chromatography // Inverse Problems in Science and Engineering. 2017. V. 25. № 8. P. 1112–1137.

  9. Volpert A.I., Volpert V.A., Volpert Vl.A. Traveling wave solutions of parabolic systems. American Mathematical Society, 2000.

  10. Meinhardt H. Models of biological pattern formation. London: Academic Press, 1982.

  11. FitzHugh R. Impulses and physiological states in theoretical model of nerve membrane // Biophysical Journal. 1961. V. 1. № 1. P. 445–466.

  12. Murray J.D. Mathematical biology. I. An introduction. Springer, New York, 2002.

  13. Egger H., Pietschmann J.-F., Schlottbom M. Identification of nonlinear heat conduction laws // Journal of Inverse and Ill-Posed Problems. 2015. V. 23. № 5. P. 429–437.

  14. Gholami A., Mang A., Biros G. An inverse problem formulation for parameter estimation of a reaction-diffusion model of low grade gliomas // Journal of Mathematical Biology. 2016. V. 72. № 1–2. P. 409–433.

  15. Aliev R.R., Panfilov A.V. A simple two-variable model of cardiac excitation // Chaos, Solitons and Fractals. 1996. V. 7. № 3. P. 293–301.

  16. Generalov E.A., Levashova N.T., Sidorova A.E., Chumankov P.M., Yakovenko L.V. An autowave model of the bifurcation behavior of transformed cells in response to polysaccharide // Biophysics. 2017. V. 62. № 5. P. 876–881.

  17. Mang A., Gholami A., Davatzikos C., Biros G. PDE-constrained optimization in medical image analysis // Optimization and Engineering. 2018. V. 19. № 3. P. 765–812.

  18. Kabanikhin S.I., Shishlenin M.A. Recovering a Time-Dependent Diffusion Coefficient from Nonlocal Data // Numerical Analysis and Applications. 2018. V. 11. P. 38–44.

  19. Mamkin V., Kurbatova J., Avilov V., Mukhartova Yu., Krupenko A., Ivanov D., Levashova N., Olchev A. Changes in net ecosystem exchange of CO2, latent and sensible heat fluxes in a recently clear-cut spruce forest in western Russia: results from an experimental and modeling analysis // Environmental Research Letters. 2016. V. 11. № 12. P. 125012.

  20. Levashova N.T., Muhartova J.V., Olchev A.V. Two approaches to describe the turbulent exchange within the atmospheric surface layer // Mathematical Models and Computer Simulations. 2017. V. 9. № 6. P. 697–707.

  21. Levashova N., Sidorova A., Semina A., Ni M. A spatio-temporal autowave model of shanghai territory development // Sustainability. 2019. V. 11. № 13. P. 3658.

  22. Zakharova S.A., Davydova M.A., Lukyanenko D.V. A spatio-temporal autowave model of shanghai territory development // Inverse Problems in Science and Engineering. 2020. V. 29. № 3. P. 365–377.

  23. Isakov V.M., Kabanikhin S.I., Shananin A.A., Shishlenin M.A., Zhang S. Algorithm for determining the volatility function in the Black-Scholes model // Computational Mathematics and Mathematical Physics. 2019. V. 59. № 10. P. 1753–1758.

  24. Kadalbajoo M.K., Gupta V. A brief survey on numerical methods for solving singularly perturbed problems // Applied Mathematics and Computation. 2010. V. 217. № 18. P. 3641–3716.

  25. Cannon J.R., DuChateau P. An Inverse problem for a nonlinear diffusion equation // SIAM Journal on Applied Mathematics. 1980. V. 39. № 2. P. 272–289.

  26. DuChateau P., Rundel W. Unicity in an inverse problem for an unknown reaction term in a reaction-diffusion equation // Journal of differential equations. 1985. V. 59. P. 155–165.

  27. Pilant M.S., Rundell W. An inverse problem for a nonlinear parabolic equation // Communications in Partial Differential Equations. 1986. V. 11. № 4. P. 445–457.

  28. Kabanikhin S.I. Definitions and examples of inverse and ill-posed problems // Journal of Inverse and Ill-Posed Problems. 2008. V. 16. № 4. P. 317–357.

  29. Kabanikhin S.I. Inverse and Ill-posed Problems Theory and Applications. De Gruyter, 2011.

  30. Jin B., Rundell W. A tutorial on inverse problems for anomalous diffusion processes // Inverse Problems. 2015. V. 31. P. 035003.

  31. Belonosov A., Shishlenin M. Regularization methods of the continuation problem for the parabolic equation // Lecture Notes in Computer Science. 2017. V. 10187. P. 220–226.

  32. Kaltenbacher B., Rundell W. On the identification of a nonlinear term in a reaction-diffusion equation // Inverse Problems. 2019. V. 35. P. 115007.

  33. Belonosov A., Shishlenin M., Klyuchinskiy D. A comparative analysis of numerical methods of solving the continuation problem for 1D parabolic equation with the data given on the part of the boundary // Advances in Computational Mathematics. 2019. V. 45. № 2. P. 735–755.

  34. Kaltenbacher B., Rundell W. The inverse problem of reconstructing reaction-diffusion systems // Inverse Problems. 2020. V. 36. № 065011.

  35. Lukyanenko D.V., Shishlenin M.A., Volkov V.T. Solving of the coefficient inverse problems for a nonlinear singularly perturbed reaction-diffusion-advection equation with the final time data // Communications in Nonlinear Science and Numerical Simulation. 2018. V. 54. P. 233–247.

  36. Lukyanenko D.V., Prigorniy I.V., Shishlenin M.A. Some features of solving an inverse backward problem for a generalized Burgers’ equation // Journal of Inverse and Ill-posed Problems. 2020. V. 28. № 5. P. 641–649.

  37. Lukyanenko D.V., Borzunov A.A., Shishlenin M.A. Solving coefficient inverse problems for nonlinear singularly perturbed equations of the reaction-diffusion-advection type with data on the position of a reaction front // Communications in Nonlinear Science and Numerical Simulation. 2021. V. 99. P. 105824.

  38. Lukyanenko D., Yeleskina T., Prigorniy I., Isaev T., Borzunov A., Shishlenin M. Inverse problem of recovering the initial condition for a nonlinear equation of the reaction-diffusion-advection type by data given on the position of a reaction front with a time delay // Mathematics. 2021. V. 9. № 4. P. 342.

  39. Lukyanenko D.V., Grigorev V.B., Volkov V.T., Shishlenin M.A. Solving of the coefficient inverse problem for a nonlinear singularly perturbed two-dimensional reaction-diffusion equation with the location of moving front data // Computers and Mathematics with Applications. 2019. V. 77. № 5. P. 1245–1254.

  40. Lukyanenko D.V., Shishlenin M.A., Volkov V.T. Asymptotic analysis of solving an inverse boundary value problem for a nonlinear singularly perturbed time-periodic reaction-diffusion-advection equation // Journal of Inverse and Ill-Posed Problems. 2019. V. 27. № 5. P. 745–758.

  41. Vasil’eva A.B., Butuzov V.F., Nefedov N.N. Singularly perturbed problems with boundary and internal layers // Proceedings of the Steklov Institute of Mathematics. 2010. V. 268. P. 258–273.

  42. Tikhonov A.N., Goncharsky A.V., Stepanov V.V., Yagola A.G. Numerical methods for the solution of ill-posed problems. Dordrecht, Kluwer Academic Publishers, 1995.

  43. Alifanov O.M., Artuhin E.A., Rumyantsev S.V. Extreme methods for the solution of ill-posed problems. M.: Nauka, 1988.

  44. Kabanikhin S.I., Shishlenin M.A. Quasi-solution in inverse coefficient problems // Journal of Inverse and Ill-Posed Problems. 2008. V. 16. № 7. P. 705–713.

  45. Hairer E., Wanner G. Solving ordinary differential equations II. stiff and differential-algebraic problems. Springer-Verlag Berlin Heidelberg, 1996.

  46. Rosenbrock H.H. Some general implicit processes for the numerical solution of differential equations // The Computer Journal. 1963. V. 5. № 4. P. 329–330.

  47. Alshin A., Alshina E., Kalitkin N., Koryagina A. Rosenbrock schemes with complex coefficients for stiff and differential algebraic systems // Computational Mathematics and Mathematical Physics. 2006. V. 46. P. 1320–1340.

  48. Wen X. High order numerical methods to a type of delta function integrals // Journal of Computational Physics. 2007. V. 226. № 2. P. 1952–1967.

  49. Egger H., Engl H.W., Klibanov M.V. Global uniqueness and Holder stability for recovering a nonlinear source term in a parabolic equation // Inverse Problems. 2005. V. 21. № 1. P. 271–290.

  50. Belina L., Klibanov M.V. A globally convergent numerical method for a coefficient inverse problem // SIAM Journal on Scientific Computing. 2008. V. 31. № 1. P. 478–509.

  51. Klibanov M.V., Fiddy M.A., Beilina L., Pantong N., Schenk J. Picosecond scale experimental verification of a globally convergent algorithm for a coefficient inverse problem // Inverse Problems. 2010. V. 26. № 4. P. 045003.

  52. Chaikovskii D., Zhang Y. Convergence analysis for forward and inverse problems in singularly perturbed time-dependent reaction-advection-diffusion equations // arXiv. 2021. arXiv:2106.15249 [math.NA]

  53. Yagola A.G., Leonov A.S., Titatenko V.N. Data errors and an error estimation for ill-posed problems // Inverse Problems in Engineering. 2002. V. 10. № 2. P. 117–129.

  54. Dorofeev K.Y., Titatenko V.N., Yagola A.G. Algorithms for constructing a posteriori errors of solutions to ill-posed problems // Computational Mathematics and Mathematical Physics. 2003. V. 43. № 1. P. 10–23.

  55. Leonov A.S. A posteriori accuracy estimations of solutions to ill-posed inverse problems and extra-optimal regularizing algorithms for their solution // Numerical Analysis and Applications. 2012. V. 5. № 1. P. 68–83.

Дополнительные материалы отсутствуют.

Инструменты

следующая статья выпуска предыдущая статья выпуска содержание выпуска

Журнал вычислительной математики и математической физики

Архивы выпусков Информация о журнале Отправить рукопись в журнал