Журнал вычислительной математики и математической физики, 2019, T. 59, № 9, стр. 1537-1545

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

П. Н. Вабищевич^1, 2, *

¹ ИБРАЭ РАН
115191 Москва, Б. Тульская ул., 52, Россия

² СВФУ им. М.К. Аммосова
677000 Якутск, ул. Белинского, 58, Россия

^* E-mail: vabishchevich@gmail.com

Поступила в редакцию 04.04.2019
После доработки 04.04.2019
Принята к публикации 15.05.2019

DOI: 10.1134/S0044466919090151

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

Аннотация

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

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

ВВЕДЕНИЕ

Математическое моделирование многих прикладных проблем приводит к необходимости численного решения обратных задач для уравнений с частными производными [1], [2]. При теоретическом исследовании таких задач [3], [4] рассматриваются принципиальные вопросы единственности решения, его устойчивости в тех или иных функциональных классах. При построении вычислительных алгоритмов их приближенного решения основное внимание уделяется разработке устойчивых вычислительных алгоритмов, максимальному учету особенностей обратных задач [5], [6]. В частности, многие обратные задачи относятся к классу некорректных (условно корректных) задач и для их приближенного решения необходимо использовать соответствующие регуляризирующие алгоритмы [7], [8].

Среди обратных задач для уравнений с частными производными можно выделить коэффициентные обратные задачи, которые связываются с идентификацией коэффициентов уравнения и/или правой части при задании некоторой дополнительной информации. При рассмотрении нестационарных задач обычно выделяют как самостоятельные задачи идентификации зависимости коэффициентов от времени или от пространственных переменных [3], [4]. Задачи идентификации правой части уравнения относятся к классу линейных обратных задач, что существенно облегчает их исследование.

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

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

Разработка вычислительных алгоритмов решения коэффициентных обратных задач для гиперболических уравнений может базироваться на соответствующих результатах для параболических уравнений. В частности, при рассмотрении задач нахождения зависимости правой части и коэффициентов от времени имеет смысл ориентироваться на безытерационные методы [10]–[12], которые основаны на специальной декомпозиции решения на каждом временном слое. Переход на новый слой по времени обеспечивается решением стандартных дискретных задач. Тем самым решение обратной задачи в вычислительном плане ненамного сложнее приближенного решения прямой задачи. В настоящей работе такой подход применен для задачи определения зависимости правой части многомерного гиперболического уравнения от времени по информации о решении во внутренней точке расчетной области.

1. ПОСТАНОВКА ЗАДАЧИ

Для простоты ограничимся двумерной задачей. Переход к более общим трехмерным задачам носит, во многом, редакционный характер. Пусть $x = ({{x}_{1}},{{x}_{2}})$ и $\Omega $ – ограниченный многоугольник. Прямая задача ставится следующим образом. Ищется $u(x,t)$, $0 \leqslant t \leqslant T$, $T > 0$, которая является решением гиперболического уравнения второго порядка

(1)

$\frac{{{{\partial }^{2}}u}}{{\partial {{t}^{2}}}} - \operatorname{div} (k(x)\operatorname{grad} u) = f(x,t),\quad x \in \Omega ,\quad 0 < t \leqslant T.$

Заданы также граничные и начальные условия:

(2)

$k(x)\frac{{\partial u}}{{\partial n}} = 0,\quad x \in \partial \Omega ,\quad 0 < t \leqslant T,$

(3)

$u(x,0) = {{u}_{0}}(x),\quad \frac{{\partial u}}{{\partial t}}(x,0) = {{{v}}_{0}}(x),\quad x \in \Omega ,$

где $n$ – нормаль к $\Omega $. В виде (1)–(3) формулируется прямая задача, в которой правая часть, коэффициенты уравнения и граничных условий, начальное условие заданы.

Будем рассматривать обратную задачу, когда в уравнении (1) правая часть имеет вид

(4)

$f(x,t) = {{f}_{0}}(x,t) + p(t){{f}_{1}}(x,t),\quad x \in \Omega ,\quad 0 < t \leqslant T.$

Коэффициент $p(t)$ неизвестен, а функции ${{f}_{\alpha }}(x,t)$, $\alpha = 1,2$, заданы. Дополнительное условие (переопределение) часто формулируется в виде

(5)

$\int\limits_\Omega {u(x,t)\omega (x)dx} = \varphi (t),\quad 0 < t \leqslant T,$

где $\omega (x) \geqslant 0$ – некоторая весовая функция. В частности, при выборе $\omega (x) = \delta (x - x*)$ ($x* \in \Omega $), где $\delta (x)$ есть $\delta $-функция, из (5) имеем

(6)

$u(x*,t) = \varphi (t),\quad 0 < t \leqslant T.$

Будем считать, что поставленная обратная задача по нахождению пары $u(x,t)$, $p(t)$ из условий (1)–(4) и дополнительных условий (5) или (6) является корректной. Соответствующие условия существования и однозначной разрешимости можно найти в книге [4]. В данной работе мы основное внимание уделяем проблемам численного решения поставленных обратных задач.

2. ЧИСЛЕННОЕ РЕШЕНИЕ ПРЯМОЙ ЗАДАЧИ

Для приближенного решения задачи (1)–(3) будем использовать равномерную сетку по времени: ${{t}^{n}} = n\tau $, $n = 0,1,\; \ldots ,\;N$, $\tau N = T$ и использовать обозначения ${{y}^{n}} = y({{t}^{n}})$. Применяется стандартная конечно-элементная аппроксимация по пространству [13]. В многоугольнике $\Omega $ проводится триангуляция, на этой расчетной сетке определим конечномерное пространство конечных элементов $V \subset {{H}^{1}}(\Omega )$.

Аппроксимация по времени базируется на использовании трехслойной схемы с весами [14], [15]. Приближенное решение на новом слое по времени определяется из следующей вариационно-разностной схемы:

(7)

$\begin{gathered} \int\limits_\Omega {\frac{{{{y}^{{n + 1}}} - 2{{y}^{n}} + {{y}^{{n - 1}}}}}{{{{\tau }^{2}}}}{\kern 1pt} } {v}dx + \int\limits_\Omega k \operatorname{grad} (\sigma {{y}^{{n + 1}}} + (1 - 2\sigma ){{y}^{n}} + \sigma {{y}^{{n - 1}}})\operatorname{grad} {v}dx = \\ = \;\int\limits_\Omega {{{f}^{n}}{v}dx} \quad \forall {v} \subset V,\quad n = 1,2,\; \ldots ,\;N - 1. \\ \end{gathered} $

Начальные условия (3) дают

(8)

$\int\limits_\Omega {{{y}^{0}}{v}dx} x = \int\limits_\Omega {{{u}_{0}}{v}dx} ,\quad \int\limits_\Omega {{{y}^{1}}{v}dx} = \int\limits_\Omega {{{u}_{1}}{v}dx} \quad \forall {v} \subset V.$

Принимая во внимание уравнение (1) для ${{u}_{1}} \approx u({{t}_{1}})$ положим

${{u}_{1}}(x) = {{u}_{0}}(x) + \tau {{{v}}_{0}}(x) + \frac{{{{\tau }^{2}}}}{2}(\operatorname{div} (k(x)\operatorname{grad} {{u}_{0}}) + f(x,0)),\quad x \in \Omega .$

Схема (7), (8) имеет второй порядок точности по времени на гладких решениях.

Введем скалярное произведение и норму в ${{L}_{2}}(\Omega )$:

$(u,{v}) = \int\limits_\Omega {u{v}dx} ,\quad \left\| u \right\| = {{(u,u)}^{{1/2}}}.$

В этих обозначениях схема (7), (8) записывается в виде

(9)

$\begin{gathered} \left( {\frac{{{{y}^{{n + 1}}} - 2{{y}^{n}} + {{y}^{{n - 1}}}}}{{{{\tau }^{2}}}},{v}} \right) + (k\operatorname{grad} (\sigma {{y}^{{n + 1}}} + (1 - 2\sigma ){{y}^{n}} + \sigma {{y}^{{n - 1}}}),\operatorname{grad} {v}) = \\ = \;({{f}^{n}},{v}),\quad n = 1,2,\; \ldots ,\;N - 1, \\ \end{gathered} $

(10)

$({{y}^{0}},{v}) = ({{u}_{0}},{v}),\quad ({{y}^{1}},{v}) = ({{u}_{1}},{v})\quad \forall {v} \subset V.$

Условия устойчивости трехслойных операторно-разностных схем с весами хорошо известны [15], [16]. Аналогично рассматривается и вариационно-разностная схема (7), (8) . Сформулируем типичный результат об устойчивости по начальным данным и правой части.

Теорема 1. Вариационно-разностная схема (7), (8) безусловно устойчива при $\sigma \geqslant 0.25$ и для решения при $\tau \leqslant 0.5$ верна априорная оценка

(11)

${{\mathcal{E}}^{{n + 1}}} = {{\varrho }^{n}}{{\mathcal{E}}^{1}} + \tau \sum\limits_{k = 1}^n {{{\rho }^{{n - k}}}} {{\left\| {{{f}^{k}}} \right\|}^{2}},\quad n = 1,2,\; \ldots ,\;N - 1,$

в которой $\varrho = {\text{exp}}(3\tau )$ и

${{\mathcal{E}}^{n}} = {{\left\| {\frac{{{{y}^{n}} - {{y}^{{n - 1}}}}}{\tau }} \right\|}^{2}} + \left( {\sigma - \frac{1}{4}} \right){{\tau }^{2}}{{\left\| {{{k}^{{1/2}}}\operatorname{grad} \frac{{{{y}^{n}} - {{y}^{{n - 1}}}}}{\tau }} \right\|}^{2}} + {{\left\| {{{k}^{{1/2}}}\operatorname{grad} \frac{{{{y}^{n}} + {{y}^{{n - 1}}}}}{2}} \right\|}^{2}}.$

Доказательство. Введем новые функции:

${{s}^{n}} = \frac{{{{y}^{n}} + {{y}^{{n - 1}}}}}{2},\quad {{r}^{n}} = \frac{{{{y}^{n}} - {{y}^{{n - 1}}}}}{\tau }.$

Принимая во внимание, что

${{y}^{n}} = \frac{1}{4}({{y}^{{n + 1}}} + 2{{y}^{n}} + {{y}^{{n - 1}}}) - \frac{1}{4}({{y}^{{n + 1}}} - 2{{y}^{n}} + {{y}^{{n - 1}}}),$

запишем уравнение (9) в виде

(12)

$\left( {\frac{{{{r}^{{n + 1}}} - {{r}^{n}}}}{\tau },v} \right) + \left( {\sigma - \frac{1}{4}} \right)\tau (k\operatorname{grad} ({{r}^{{n + 1}}} - {{r}^{n}}),\operatorname{grad} {v}) + \left( {k\operatorname{grad} \frac{{{{s}^{{n + 1}}} + {{s}^{n}}}}{2},\operatorname{grad} {v}} \right) = ({{f}^{n}},{v}).$

Выбирая

${v} = 2({{s}^{{n + 1}}} - {{s}^{n}}) = \tau ({{r}^{{n + 1}}} + {{r}^{n}})$

в (12), получаем

(13)

$\begin{gathered} {{\left\| {{{r}^{{n + 1}}}} \right\|}^{2}} - {{\left\| {{{r}^{n}}} \right\|}^{2}} + \left( {\sigma - \frac{1}{4}} \right){{\tau }^{2}}\left( {{{{\left\| {{{k}^{{1/2}}}\operatorname{grad} {{r}^{{n + 1}}}} \right\|}}^{2}} - {{{\left\| {{{k}^{{1/2}}}\operatorname{grad} {{r}^{n}}} \right\|}}^{2}}} \right) + \\ + \;{{\left\| {{{k}^{{1/2}}}\operatorname{grad} {{s}^{{n + 1}}}} \right\|}^{2}} - {{\left\| {{{k}^{{1/2}}}\operatorname{grad} {{s}^{n}}} \right\|}^{2}} = \tau ({{f}^{n}},{{r}^{{n + 1}}} + {{r}^{n}}). \\ \end{gathered} $

Для правой части привлекается оценка

$({{f}^{n}},{{r}^{{n + 1}}} + {{r}^{n}}) \leqslant {{\left\| {{{f}^{n}}} \right\|}^{2}} + \frac{1}{4}{{\left\| {{{r}^{{n + 1}}} + {{r}^{n}}} \right\|}^{2}} \leqslant {{\left\| {{{f}^{n}}} \right\|}^{2}} + {{\left\| {{{r}^{{n + 1}}}} \right\|}^{2}} + {{\left\| {{{r}^{n}}} \right\|}^{2}}.$

При $\tau \leqslant 0.5$ имеем

$\frac{1}{{1 - \tau }} < 1 + 2\tau < \exp (2\tau ).$

С учетом этого и введенных обозначений из (13) получим

${{\mathcal{E}}^{{n + 1}}} = \varrho {{\mathcal{E}}^{n}} + \tau {{\left\| {{{f}^{n}}} \right\|}^{2}}.$

Из такой послойной оценки следует доказываемое неравенство (11).

Отметим, что ограничения на шаг по времени не носят принципиального характера и введены для упрощения доказательства.

3. ЧИСЛЕННОЕ РЕШЕНИЕ ОБРАТНОЙ ЗАДАЧИ

При приближенном решении обратной задачи (1)–(5) аналогично (9) используется вариационно-разностная схема

(14)

$\begin{gathered} \left( {\frac{{{{y}^{{n + 1}}} - 2{{y}^{n}} + {{y}^{{n - 1}}}}}{{{{\tau }^{2}}}},{v}} \right) + (k\operatorname{grad} (\sigma {{y}^{{n + 1}}} + (1 - 2\sigma ){{y}^{n}} + \sigma {{y}^{{n - 1}}}),\operatorname{grad} {v}) = \\ = \;(f_{0}^{n},{v}) + {{g}^{n}}(f_{1}^{n},{v}),\quad n = 1,2,\; \ldots ,\;N - 1, \\ \end{gathered} $

где ${{g}^{n}} \approx p({{t}^{n}})$. Дополнительные соотношения (5) дают

(15)

$({{y}^{{n + 1}}},\omega ) = {{\varphi }^{{n + 1}}},\quad n = 1,2,\; \ldots ,\;N - 1.$

Отметим, что обратная задача определения ${{y}^{{n + 1}}}$, ${{g}^{n}}$, $n = 1,2,\; \ldots ,\;N - 1$, на каждом шаге по времени является линейной. Поэтому можно в том или ином варианте использовать линейную суперпозицию.

Для приближенного решения на новом временном слое используется следующая (см. [6], [10]) декомпозиция решения

(16)

${{y}^{{n + 1}}}(x) = {{w}^{{n + 1}}}(x) + {{g}^{n}}{{z}^{{n + 1}}}(x).$

Задачи для ${{w}^{{n + 1}}}(x)$ и ${{z}^{{n + 1}}}(x)$ формулируются после подстановки (16) в (14). Для нахождения ${{w}^{{n + 1}}}(x)$ используется уравнение

(17)

$\begin{gathered} \left( {\frac{{{{w}^{{n + 1}}} - 2{{y}^{n}} + {{y}^{{n - 1}}}}}{{{{\tau }^{4}}}},{v}} \right) + (k\operatorname{grad} (\sigma {{w}^{{n + 1}}} + (1 - 2\sigma ){{y}^{n}} + \sigma {{y}^{{n - 1}}}),\operatorname{grad} {v}) = (f_{0}^{n},{v}), \\ n = 1,2,\; \ldots ,\;N - 1. \\ \end{gathered} $

Функция ${{z}^{{n + 1}}}(x)$ определяется из

(18)

$\frac{1}{{{{\tau }^{2}}}}({{z}^{{n + 1}}},{v}) + \sigma (k\operatorname{grad} {{z}^{{n + 1}}},\operatorname{grad} {v}) = (f_{1}^{n},{v}),\quad n = 1,2,\; \ldots ,\;N - 1.$

При декомпозиции (16)–(18) уравнение (14) выполняется автоматически при любом ${{g}^{n}}$. Для определения ${{g}^{n}}$ привлекается условие (15). Подстановка (16) в (15) дает

(19)

${{g}^{n}} = \frac{1}{{({{z}^{{n + 1}}},\omega )}}({{\varphi }^{{n + 1}}} - ({{w}^{{n + 1}}},\omega )),\quad n = 1,2,\; \ldots ,\;N - 1.$

Применимость рассматриваемого вычислительного алгоритма будет обеспечена условием

(20)

$({{z}^{{n + 1}}},\omega ) \ne 0,\quad n = 1,2,\; \ldots ,\;N - 1.$

Вспомогательная функция ${{z}^{{n + 1}}}(x)$ определяется из дискретного эллиптического уравнения (18). Условие (20) выполняется, например, в случае $f_{1}^{n}(x) > 0$, $x \in \Omega $, если для дискретной задачи имеет место принцип максимума (монотонные аппроксимации), когда мы имеем ${{z}^{{n + 1}}}(x) > 0$, $x \in \Omega $.

Конечно-элементные аппроксимации краевых задач для эллиптических уравнений второго порядка являются монотонными при использовании линейных конечных элементов на сетках Делоне (см., например, [17], [18]). Дополнительные ограничения порождаются (см., например, [19], [20]) на величину коэффициента реакции. Применительно к нашей задаче (18) коэффициент реакции равен ${{\tau }^{{ - 2}}}$ и мы имеем ограничение снизу на шаг сетки по времени: $\tau \geqslant \mathcal{O}(h)$, где $h$ – шаг сетки по пространству. Ограничение на сетку, связанные с коэффициентом реакции (с шагом сетки по времени $\tau $), можно снять за счет применения lumping процедуры (см., например, [13]). Тем самым ограничение (20) для широкого класса задач не является обременительным.

Итогом нашего рассмотрения является

Теорема 2. Решение обратной задачи (10), (14), (15) представляется в виде (16), (19), где ${{w}^{{n + 1}}}(x)$ есть решение дискретного эллиптического уравнения (17), а ${{z}^{{n + 1}}}(x)$ есть решение (19) при условии, что выполнено (20).

Предложенный вычислительный алгоритм решения обратной задачи идентификации зависимости правой части от времени является эволюционным – последовательное определение решения на каждом слое по времени. Вычислительная сложность алгоритма не большая: его реализация связана с решением двух стандартных эллиптических задач на каждом слое по времени. Более того, если в (4) функция ${{f}_{1}}$ не зависит от времени, то задача для ${{w}^{n}}(x) = w(x)$ решается один раз.

4. БОЛЕЕ ОБЩИЕ ЗАДАЧИ

Отметим возможности применения описанного вычислительного алгоритма с некоторыми модификациями для решения более общих, чем (1)–(5), который базируется на декомпозиции (16).

Прежде всего отметим обратные задачи для более общих, чем (1), уравнений. Можно рассмотреть, в частности, задачи для гиперболического уравнения

$\rho (x)\frac{{{{\partial }^{2}}u}}{{\partial {{t}^{2}}}} + \gamma (x)\frac{{\partial u}}{{\partial t}} - \operatorname{div} (k(x)\operatorname{grad} u) + (v \cdot \operatorname{grad} )u + c(x)u = f(x,t).$

Вместо (2) имеет смысл сформулировать более общие краевые условия, например, смешанные краевые условия, условия III рода.

Для практических приложений важна возможность повышения точности решения обратной задачи за счет использования множественного переопределения. Пусть вместо одного наблюдения за решением (5) есть несколько:

(21)

$\int\limits_\Omega {u(x,t)} {{\omega }_{\alpha }}(x)dx = {{\varphi }_{\alpha }}(t),\quad \alpha = 1,2,\; \ldots ,\;m,\quad 0 < t \leqslant T.$

Например, проводятся измерения в m точках $x_{\alpha }^{ * }$, $\alpha = 1,2,\; \ldots ,\;m$, области $\Omega $. По отдельным наблюдениям в соответствии с (19) находятся $g_{\alpha }^{n}$. В качестве приближенного решения ${{g}^{n}}$ берется некоторое взвешенное среднее из них. Например, естественно положить

${{g}^{n}} = \frac{1}{{{{\xi }^{n}}}}\sum\limits_{\alpha = 1}^m {\xi _{\alpha }^{n}g_{\alpha }^{n}} ,\quad {{\xi }^{n}} = \sum\limits_{\alpha = 1}^m {\xi _{\alpha }^{n}} ,\quad \xi _{\alpha }^{n} = \left| {({{z}^{{n + 1}}},\omega )} \right|,\quad n = 1,2,\; \ldots ,\;N - 1.$

Множественные наблюдения позволяют определить несколько неизвестных зависимостей правой части от времени. Вместо (4) используется представление

$f({\mathbf{x}},t) = {{f}_{0}}(x,t) + \sum\limits_{\alpha = 1}^m {{{p}_{\alpha }}(t)} {{f}_{\alpha }}(x,t),\quad x \in \Omega ,\quad 0 < t \leqslant T.$

Функции ${{p}_{\alpha }}(t)$, $\alpha = 1,2,\; \ldots ,\;m$, определяются по наблюдениям (21). Приближенное решение обратной задачи ищется в несколько более общем, чем (16), виде:

${{y}^{{n + 1}}}(x) = {{w}^{{n + 1}}}(x) + \sum\limits_{\alpha = 1}^m {g_{\alpha }^{n}z_{\alpha }^{{n + 1}}} (x).$

Для $z_{\alpha }^{{n + 1}}(x)$ получим задачу

$\frac{1}{{{{\tau }^{2}}}}(z_{\alpha }^{{n + 1}},{v}) + \sigma (k\operatorname{grad} z_{\alpha }^{{n + 1}},\operatorname{grad} {v}) = (f_{\alpha }^{n},{v}),\quad \alpha = 1,2,\; \ldots ,\;m,\quad n = 1,2,\; \ldots ,\;N - 1.$

С учетом наблюдений (21) параметры $g_{\alpha }^{n}$, $\alpha = 1,2,\; \ldots ,\;m$, находятся из системы m линейных уравнений.

5. ЧИСЛЕННЫЕ ПРИМЕРЫ

Демонстрация работоспособности вычислительного алгоритма идентификации зависимости правой части гиперболического уравнения от времени выполнена на модельной двумерной задаче. Задача (1)–(5) рассматривается в случае, когда расчетная область $\Omega $ есть единичный квадрат:

$\Omega = \left\{ {x\,|\,x = ({{x}_{1}},{{x}_{2}}),\;0 < {{x}_{\beta }} < 1,\;\beta = 1,2} \right\}.$

Пусть $k(x) = {{10}^{3}}$ в (1), а начальные условия – однородны: ${{u}_{0}}(x) = {{{v}}_{0}}(x) = 0$ в (3). Правая часть задается в виде (4) при ${{f}_{0}}(x) = 0$, ${{f}_{1}}(x) = {{x}_{1}}{{x}_{2}}$.

Численные эксперименты проведены для точного решения задачи $p(t) = \sin (2\pi {{t}^{\gamma }})$ при различных значениях $\gamma $ при условии, что $T = 1$. Влияние параметра $\gamma $ на восстанавливаемую зависимость правой части от времени иллюстрируется фиг. 1. Для базового варианта выбрано $\gamma = 1$.

Фиг. 1.

Зависимость правой части от времени.

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

Прямая задача решается на равномерной треугольной сетке, когда шаг на сетке по каждому направлению есть $h = 0.02$. Используются кусочно-линейные лагранжевые конечные элементы. В рассматриваемых задачах принципиальной является аппроксимация по времени – точность приближенного решения зависит, прежде всего, от сетки по времени.

Измерения связываются с выбором весовой функции $\omega (x)$ в (5). В приведенных ниже расчетах

$\omega (x) = exp\left( { - \frac{{x - x{\text{*}}}}{{{{d}^{2}}}}} \right),$

когда $x{\text{*}}$ – есть точка наблюдения, а величину $d$ будем рассматривать как характеристику измерительной аппаратуры. В приведенных ниже результатах расчетов мы брали $d = 0.05$. По численному решению прямой задачи рассчитывались функции $\varphi (t)$ для различных точек наблюдения – см. фиг. 2. Для базового варианта положим $x{\text{*}} = (0.5,0.5)$.

Фиг. 2.

Измеряемые данные в трех точках $x{\text{*}}$.

Погрешность приближенного решения задачи идентификации будем оценивать величиной $\varepsilon ({{t}^{n}}) = p({{t}^{n}}) - {{g}^{n}}$. Результаты расчетов при использовании сгущающихся сеток по времени показаны на фиг. 3. Точность решения обратной задачи существенно зависит от точки наблюдения. Этот факт иллюстрируется фиг. 4, 5 . Хотя ошибка при использовании грубых сеток по времени большая, мы наблюдаем сходимость приближенного решения задачи идентификации.

Фиг. 3.

Погрешность решения задачи идентификации при $x{\text{*}} = (0.5,\;0.5)$.

Фиг. 4.

Погрешность решения задачи идентификации при $x{\text{*}} = (0.25,\;0.25)$.

Фиг. 5.

Погрешность решения задачи идентификации при $x{\text{*}} = (0.75,\;0.75)$.

Сходимость приближенного решения при использовании более подробных сеток иллюстрируется данными в табл. 1. Здесь показана максимальная погрешность ${{\varepsilon }_{a}}$ ($\begin{gathered} {{\varepsilon }_{a}} = max\varepsilon ({{t}^{n}}) \hfill \\ ^{n} \hfill \\ \end{gathered} $, $0 < n < N$). При использовании таких сеток по времени максимум погрешности достигается вблизи начального состояния и слабо зависит от точки наблюдения.

Таблица 1.

Погрешность ${{\varepsilon }_{a}}$

$x{\kern 1pt} *$	$N = 250$	$N = 500$	$N = 1000$
(0.25, 0.25)	0.0115248	0.0054003	0.0032672
(0.5, 0.5)	0.0096114	0.0054030	0.0032672
(0.75, 0.75)	0.0097412	0.0054042	0.0032672

Ключевой параметр вычислительного алгоритма связан со знаменателем в формуле (19). В нашем случае он не зависит от времени. В табл. 2 приведена величина $\xi = ({{z}^{{n + 1}}},w)$ на различных сетках по времени и для трех точек наблюдения.

Таблица 2.

Величина $\xi $

$x{\kern 1pt} *$	$N$	$\xi $
(0.25, 0.25)	25	0.00025392
	50	4.2621e–05
	100	6.3621e–06
(0.5, 0.5)	25	0.00031421
	50	7.8558e–05
	100	1.9641e–05
(0.75, 0.75)	25	0.00039626
	50	0.00013056
	100	4.0618e–05

Точность идентификации зависимости правой части от времени при других значениях $\gamma $ (см. фиг. 1) показана на фиг. 6, 7 . В этих расчетах точка наблюдения $x{\text{*}} = (0.5,0.5)$. Погрешность растет в области сильного изменения искомой зависимости правой части.

Фиг. 6.

Погрешность решения задачи идентификации при $\gamma = 0.5$.

Фиг. 7.

Погрешность решения задачи идентификации при $\gamma = 2$.

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

Алифанов О.М. Обратные задачи теплообмена. М.: Машиностр., 1988.
Лаврентьев М.М., Романов В.Г., Шишатский С.П. Некорректные задачи математической физики и анализа. М.: Наука, 1980.
Isakov V. Inverse Problems for Partial Differential Equations. Berlin: Springer, 2017.
Prilepko A.I., Orlovsky D.G., Vasin I.A. Methods for solving inverse problems in mathematical physics. Marcel Dekker, Inc., 2000.
Vogel C.R. Computational Methods for Inverse Problems. Society for Industrial and Applied Mathematics, 2002.
Самарский А.А., Вабищевич П.Н. Численные методы решения обратных задач математической физики. М.: УРСС, 2004.
Тихонов А.Н., Арсенин В.Я. Методы решения некорректных задач. М.: Наука, 1979.
Engl H.W., Hanke M., Neubauer A. Regularization of Inverse Problems. Kluwer Academic Publishers, 2000.
Hussein S.O., Lesnic D. Determination of forcing functions in the wave equation. Part I: the space-dependent case // J. Eng. Math. 2016. V. 96. № 1. P. 115–133.
Вабищевич П.Н., Васильев В.И., Васильева М.В. Вычислительная идентификация правой части параболического уравнения // Ж. вычисл. матем. и матем. физ. 2015. Т. 55. № 6. С. 1020–1027.
Vabishchevich P.N., Vasil’ev V.I. Computational algorithms for solving the coefficient inverse problem for parabolic equations // Inverse Problems in Science and Engineering. 2016. V. 24. № 1. P. 42–59.
Вабищевич П.Н., Клибанов М.В. Вычислительная идентификация старшего коэффициента параболического уравнения // Дифференц. ур-ния. 2016. Т. 52. № 7. С. 896–903.
Thomée V. Galerkin Finite Element Methods for Parabolic Problems. Springer Verlag, 2006.
Самарский А.А. Теория разностных схем. М.: Наука, 1989.
Самарский А.А., Гулин А.В. Устойчивость разностных схем. М.: Наука, 1973.
Samarskii A.A., Matus P.P., Vabishchevich P.N. Difference Schemes with Operator Factors. Kluwer Academic Publ., 2002.
Letniowski F.W. Three-dimensional Delaunay triangulations for finite element approximations to a second-order diffusion operator // SIAM J. Scientific and Statistical Computing. 1992. V. 13. № 3. P. 765–770.
Huang W. Discrete maximum principle and a Delaunay-type mesh condition for linear finite element approximations of two-dimensional anisotropic diffusion problems // Numerical Mathematics: Theory, Methods and Applications. 2011. V. 4. № 3. P. 319–334.
Ciarlet P.G., Raviart P.A. Maximum principle and uniform convergence for the finite element method // Computer Methods in Applied Mechanics and Engineering. 1973. V. 2. № 1. P. 17–31.
Brandts J.H., Korotov S., Křížek M. The discrete maximum principle for linear simplicial finite element approximations of a reaction–diffusion problem // Linear Algebra and Its Applications. 2008. V. 429. № 10. P. 2344–2357.

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

Инструменты

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

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

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