On the determination of points sources in inverse heat and mass transfer problems

Толық мәтін

ВВЕДЕНИЕ

Мы рассматриваем обратные задачи об определении правой части в уравнении

$$\begin{gathered} Mu=u_t-Lu=f_0+\sum _{i=1}^r {\alpha }_i\left(t\right)\delta (x-d_i) \\ a(t,x)=\left(a_1\right(t,x),\dots ,a_n(x,t){)}^T,\nabla u={\left(\frac{дu}{дx_1},...,\frac{дu}{дx_n}\right)}^T,n=2,3, \end{gathered}$$ (1)

где Lu=div(c(x,t)∇u)-a(x,t)∇u-a₀(x,t)u, (x,t)∈Q=(0,T)×G, G – область в ℝⁿ(n=2,3) с границей Γ,c(t,x) – диагональная матрица с элементами c_i на главной диагонали, δ(x-d_i) – дельта-функция Дирака. Уравнение (1) дополняется начальными и граничными условиями:

Bu|_S = g, u|_t=0 = u₀(x), S = (0,T)×Γ, (2)

где $Bu=\frac{\partial u}{\partial N}+\sigma u-$ – производная по конормали или Bu=u и условиями переопределения

u(t,y_j )=ψ_j(t), j=1,2,…,r. (3)

В теории тепломассопереноса функция u – концентрация переносимого вещества, а правая часть характеризует источники (стоки) [1]. В самой общей постановке задачи (1)–(3) определению подлежат как сами мощности точечных источников N_i(t), так и их местоположение x_i и их число m. Описание моделей такого сорта можно найти, например, в [1]. Обратным задачам такого вида посвящено очень большое количество работ, однако основные результаты связаны с методами численного решения подобных задач, причем многие из них далеко не всегда обоснованы. Можно строить примеры, когда постановки оказываются некорректными в том смысле, что имеет место несуществование решений или их неединственность. Как правило, численные методы основаны на сведении задачи к задаче оптимального управления и минимизации соответствующего функционала, что требует больших вычислительных возможностей и не всегда приводит к желаемому результату [2, 4, 5, 6, 7]. Некоторые теоретические результаты по исследованию задачи (1)–(3) имеются в работах [8–12]. В работе [11] рассматривается стационарный случай, а в работе [10] – и нестационарный, и граничные условия переопределения (данные Коши), что позволяет, используя наборы тестовых функций и алгоритмы типа Прони, полностью решить задачу определения числа источников, их местоположения и интенсивности. В работе [12] рассматривается модельная задача (1)–(3); с помощью явного представления решений прямой задачи и использования вспомогательной вариационной задачи авторы смогли определить величины ∑_iN_ir^l_ij (здесь N_i(t)=const для всех i и r_ij=|x_i-y_j|), что позволило решить задачу при помощи алгоритма из работы [11] (см. теорему 2). Однако, как оказалось, можно решить задачу и при помощи асимптотических представлений решений стационарных задач [14]. В одномерном случае некоторые подобные результаты на эту тему приведены в [13] (асимптотическая формула определения источника и численный алгоритм). Каких-либо общих теорем существования решений задачи (1)–(3), в отличие от распределенного источника (см. [15]), фактически в литературе не имеется (см. одномерный случай в [13]). В настоящей работе мы опишем некоторые условия, гарантирующие существование и единственность решений задачи, и численный алгоритм, в какой-то степени основанный на теоретических рассмотрениях. Мы будем рассматривать случай, когда неизвестными являются только мощности источников, а координаты предполагаемых источников считаются известными. Алгоритм основан на методе конечных элементов по пространственным переменным и методе конечных разностей по времени. В конце работы приведены результаты численных экспериментов, которые показали, что решение устойчиво при случайных возмущениях данных задачи.

Определения и вспомогательные утверждения

Пусть G – область в ℝⁿ. Через L_p(G) и W_p^s(G)(1≤p≤∞) обозначим Лебега и Соболева пространства соответственно [18]. Пусть E – банахово пространство. Через L_p (G;E) (G – область в ℝⁿ) обозначим пространство сильно измеримых функций, определенных на G, со значениями в E и конечной нормой ∥∥u(x)∥_E∥_Lp(G) [18]. Мы также используем пространства C^k(G‾;E), состоящие из функций со значениями в E, имеющих в G все производные до порядка k включительно, непрерывные в G и допускающие непрерывное продолжение на замыкание G‾. Определение пространств Соболева W_p^s(G;E) также стандартное (см. [18, 19]). Для данного интервала J=(0,T) и цилиндра Q=G×J и S=Γ×J положим W_p^s,r(Q)=W_p^s(J;L_p(G))∩L_p(J;W_p^r(G)), W_p^r,s(S)=W_p^s (J;L_p(Γ))∩L_p(J;W_p^r(Γ)). Если Γ,S – некоторые множества, то символ ρ(Γ,S) обозначает далее расстояние между этими множествами. Обозначим через D(L) область определения оператора L. Символ B_r(x₀) обозначает шар радиуса r с центром в точке x₀.

Мы рассматриваем случай: G – область в ℝⁿ с компактной границей класса C², или область G совпадает с цилиндрической областью вида Ω×(0,l), где Ω имеет границу класса C² или совпадает с прямоугольником (0,X)×(0,Y). В случае неограниченной области G решение, которое мы ищем, принадлежит некоторому пространству Лебега, т. е. фактически предполагаем убывание решения при |x|→∞. Пусть a=(a₁, a₂) для n=2 и a=(a₁, a₂, a₃) для n=3. Скобками ( ⋅,⋅ ) обозначим скалярное произведение в ℝⁿ. Мы приведем известный теоретический результат из работы [16] в случае, когда c_i=1 (i=1,2,3 или i=1,2 в двумерном случае) и коэффициенты не зависят от времени. В этом случае производная по конормали в граничном условии совпадает с производной по нормали к границе. Предположим, что

u₀(x)∈W₂¹(G), u₀(x)|_Γ=g(x,0) если Bu=u. (4)

Пусть T=∞. Зафиксируем параметр λ≥0 и предположим, что

$\begin{array}{c}\begin{array}{rr}& e^{-\lambda t}g\in W_2^{1/\mathrm{4,1}/2}\left(S\right),\text{ если }Bu=\frac{\partial u}{\partial \nu }+\sigma u \left(\sigma \in C^1\left(\text{Γ}\right)\right),\\ & \end{array} \left(5\right)\end{array}$

e^-λt g∈W₂^3/4,3/2(S), если Bu=u, f₀e^-λt∈L₂(G). (6)

Введем функции

${\phi }_j\left(x\right)=\frac{-1}{2}{\int }_0^1\left(\overrightarrow{a}\left(y_j+\tau \left(x-y_j\right)\right),\left(x-y_j\right)\right)d\tau$

и предположим, что

a_i∈W_∞²(G)(i=1,…,n), ∇φ_j, Δφ_j(j=1,…,s), a₀∈L_∞ (G), σ∈C¹(Г). (7)

Решая вспомогательную задачу

w_t – Lw = f₀(t,x), Bw|_S = g, w|_t=0 = u₀(x), (8)

можем найти ее решение, и такое, что e^-λt w₀∈W₂^1,2(Q) при достаточно большом параметре λ₀ и λ≥λ₀. Заметим, что в случае цилиндрической области или параллелепипеда могут потребоваться дополнительные интегральные условия согласования данных, их можно найти, например, в работе [16]. Рассмотрим обратную задачу (1)–(3). После замены переменных w=u-w₀ мы придем к более простой задаче

$\begin{array}{c}{w}_t+Lw=\underset{i=1}{\overset{m}{}}{N}_i\left(t\right)\delta \left(x-x_i\right)\end{array}$

$\begin{array}{c}{\left.Bw\right|}_S=\mathrm{0,}{\left.w\right|}_{t=0}=0\end{array}$

$\begin{array}{c}w\left(y_j,t\right)={\psi }_j\left(t\right)-w_0\left(t,y_j\right)={\tilde{\psi }}_j\left(t\right),j=\mathrm{1,2,}\dots ,s.\end{array}$

Мы предполагаем, что справедливо представление

$\begin{array}{c}{\tilde{\psi }}_j\left(t\right)={\int }_0^t{V}_{{\delta }_j}\left(t-\tau \right){\psi }_{0j}\left(\tau \right)d\tau ,{\psi }_{0j}{e}^{-\lambda t}\in L_2\left(\mathrm{0,}T\right),\end{array}$ (9)

где V_γ(t) определяется своим преобразованием Лапласа:

${\hat{V}}_{\gamma }\left(\lambda \right)=e^{-\sqrt{\lambda }\gamma },n=3;{\hat{V}}_{\gamma }\left(\lambda \right)=\frac{i}{4}{H}_0^{\left(1\right)}\left(i\sqrt{\lambda }\gamma \right),n=2.$

Здесь H₀⁽¹⁾ – это функция Ханкеля и √λ=|λ|^1/2 e^iargλ/2 – ветвь корня аналитическая в плоскости с разрезом argλ=π. Не так трудно установить, что $V_{\gamma }\left(t\right)=\frac{e^{-{\gamma }^2/4t}}{4\pi t}$ при n=2 и при $n=\mathrm{3,}{V}_{\gamma }=\frac{\gamma e^{-{\gamma }^2/4t}}{2\sqrt{{\pi }^3/2}}$

Введем множества K={y∈G:ρ(y,∪^m_i=1d_i)≤ρ(y,Γ)}, Q_ε= {(t,x)∈Q:|x-d_i|>ε ∀i=1,2,…,r}. Пусть δ_j=min_ir_ij, j= 1,2,…,s, где r_ij=|d_i-y_j|. Введем матрицу A₀ с элементами a_ji= eφ_j(d_i), если |d_i-y_j|=δ_j, и a_ji=0 – в противном случае. Условие корректности записывается в виде

det A₀ ≠ 0. (10)

Опишем способ построения матрицы A₀. Матрица есть нормированная главная часть асимптотики по параметру λ (см. [16]) матрицы {v_i(y_j)}, где v_i есть решение задачи

λv_i-Lv_i=δ(x-d_i), Bv|_r=0.

Фиксируем p∈(1,n/(n-1)). Тогда справедлива теорема [16].

Теорема 1. Пусть T=∞ и выполнены условия (4), (7), (10) и y_i∈К для всех i=1,2,…,r. Тогда найдется λ₀≥0 такое, что при λ≥ λ₀ и выполнении условий (5), (6), (9) существует единственное решение задачи (1)–(3) такое, что u=w₀+w,w₀ есть решение вспомогательной задачи (8), e^-λtw₀∈W₂^1,2(Q), e^-λtN⃗∈L₂(0,∞), e^-λtw∈L₂(0,∞;W_p¹(G)), e^-λtw_t∈L₂(0,∞;W^-1_p,B(G)), e^-λtw∈W₂^1,2(Q_ε) для всех ε>0.

Здесь пространство W^-1_p,B(G)) есть двойственное пространство к W₂¹(G) в случае условия Дирихле и к W₂¹(G) в случае условий Неймана или Робина.

Теорема останется справедливой и в случае конечного интервала (0,T) (см. [16]). Несмотря на то, что мы рассматриваем случай не зависящих от t коэффициентов, расчеты показывают, что на нее можно опираться при рассмотрении вопросов единственности решений. В общем случае, в зависимости от расположения точек замеров {b_j}, может иметь место неединственность решений, и примеры легко строятся.

Описание алгоритма решения задачи

Опишем численный алгоритм. Считаем, что пространственная область имеет вид G=Ω×(0,Z) в случае n=3 и G – прямоугольник в случае n=2. Положим Γ=∂G,S=(0,T)×Γ. Опишем метод и результаты численных экспериментов в случае n=2, G – круг радиуса 2 с центром в 0. Естественным образом считаем, что функция с_i строго положительна. Начально-краевые условия имеют вид:

$\begin{array}{c}{\left.u\right|}_{t=0}=u_0\left(x\right),{\left.\frac{\partial u}{\partial N}\right|}_S=\mathrm{0,}\end{array}$

где производная есть производная по нормали. Для численного решения используем метод конечных элементов. Определению подлежат решение u и неизвестные функции α_i в правой части уравнения (1). Опишем метод решения прямой задачи. Задана триангуляция области G и соответствующий базис метода конечных элементов {φ_i}^N_i=1(φ_i|_Γ=0). Узлы сетки обозначим через {b_i}. Без ограничения общности считаем, что точки b_j(j=1,2,…,r) совпадают с точками источников, а точки {y_i=(y_i¹,y_i²)} – с точками b_N-r+1⋅ b_N-r+2,…,b_N. Пусть b_i=(x₁ⁱ, x₂ⁱ).

Прежде чем приводить метод решения, приведем дополнительную замену

$\begin{array}{c}v=u-\Phi ,\Phi =\underset{i=N-r+1}{\overset{N}{}}{\psi }_i\left(t\right)\frac{{\left(x_1^2+x_1^2-4\right)}^2{{\displaystyle \prod }}_{j\ne i}{\left|x-b_j\right|}^2}{{({\left(x_1^i\right)}^2+{\left(x_2^i\right)}^2-4)}^2{\left(x_1^i-X\right)}^2{{\displaystyle \prod }}_{j\ne i}{\left|b_i-b_j\right|}^2}. \#\end{array}$

Функция v есть решение задачи

$\begin{array}{c}Mv=f-M\left(\text{Φ}\right)=f_0\end{array}$

$\begin{array}{c}v\left(y_i,t\right)=\mathrm{0,}i=\mathrm{1,2,}\dots ,r,y_i\in G\end{array}$

$\begin{array}{c}{\left.v\right|}_{t=0}=u_0\left(x\right)-\Phi \left(\mathrm{0,}x\right)=v_0\left(x\right),{\left.\frac{\partial u}{\partial N}\right|}_S=0.\end{array}$

Ищем приближенное решение в виде v=∑^N_i=1C_i(t)φ_i.

Для удобства далее считаем, что точки y_i совпадают с узлами сетки b_N-r+1,…,b_N. Функции C_i определяем из системы

R₀C⃗_t+R₁(t)C⃗=F⃗+f⃗₀, C⃗=(C₁,C₂,…,C_N )^T, (11)

где F⃗ имеет координаты F_j=α_j(t) для j≤r и F_j=0 при j>r, и координаты f⃗₀ имеют вид

$f_j=\left(f_0\left(t,x\right),{\phi }_j\right)={\int }_G{f}_0\left(t,x\right){\phi }_{j}dx,$

R₀ – матрица с элементами r_ij=(φ_i,φ_j)= ∫_Gφ_i (x)φ_j(x)dx,R₁ – матрица с элементами

$\begin{array}{c}\begin{array}{rr}{R}_{jk}=\left(c_1\left(t,x\right){\phi }_{k{x}_1},{\phi }_{j{x}_1}\right)+(& \left(c_2\left(t,x\right){\phi }_{k{x}_2},{\phi }_{j{x}_2}\right)+\\ & \left(b\left(t,x\right)\nabla {\phi }_k,{\phi }_j\right)+\left(a\left(t,x\right){\phi }_k,{\phi }_j\right)\end{array}\end{array}$ (12)

Имеем, что C⃗(0)=(v₀(b₁),v₀(b₂),…,v₀(b_N))^T. Решение системы ищем методом конечных разностей. Пусть τ= T/M – шаг по времени. Заменим уравнение (11) системой

$\begin{array}{c}{R}_0\frac{{\overrightarrow{C}}_{i+1}-{\overrightarrow{C}}_i}{\tau }+A_{i+1}{\overrightarrow{C}}_{i+1}={\overrightarrow{F}}_{i+1}+{\overrightarrow{f}}_{i+1},{\overrightarrow{C}}_i={\left(C_i^1,\dots ,C_i^N\right)}^T,i=\mathrm{0,1,2,}\dots ,M-\mathrm{1,}\#\end{array}$

где C_i^k≈C_k(τi), F⃗_i≈F⃗(τi), f⃗_i=f⃗₀(τi), A_i=R₁(τi). Систему (12) можно записать в виде

R_i+1C⃗_i+1=τF⃗_i+1+τf⃗_i+1+R₀C⃗_i, C_i^k=C_k(τi),C⃗_i=(C_i¹,…,C_i^N )^T, i=0,1,2,…,M-1 (13),

где R_i+1=R₀+τA_i+1. Опишем алгоритм решения обратной задачи.

Пусть α⃗(t)=(α₁(t)…α_r(t))^T, α⃗_i=(α_i¹……,α₎^r)≈α⃗(τi).

В силу условий переопределения должно быть выполнено, что C_k^N-r+i=0. Определим начальные данные. Имеем C₀^k=v₀(b_k) при k=1,…,N.

Предположим, что мы нашли векторы α⃗_i, C⃗_i. Ищем величину C⃗_i+1 из системы

$\begin{array}{c}{R}_{i+1}{\overrightarrow{C}}_{i+1}=\tau {\overrightarrow{F}}_{i+1}+\tau {\overrightarrow{f}}_{i+1}+R_0{\overrightarrow{C}}_i.\end{array}$ (14)

Вектор α⃗_i+1 определяем из системы

$\begin{array}{c}\tau B_{i+1}{\overrightarrow{\alpha }}_{i+1}=-\tau {\text{Φ}}_0{R}_{i+1}^{-1}{\overrightarrow{f}}_{i+1}-{\text{Φ}}_0{R}_{i+1}^{-1}{R}_0{\overrightarrow{C}}_i, \end{array}$(15)

где r×r матрица имеет вид B_i+1=Φ₀R^-1_i+1Φ₁, где Φ₀ – матрица r×N, в первых N-r столбцах которой стоят нули, а в последних r столбцах стоит единичная матрица размера r×r,Φ₁ матрица N×r, в последних N-r строках стоят нули, а в первых r строках стоит единичная матрица размера r×r. Матрица B_i может быть сингулярной (с малыми элементами). Однако для корректности задачи мы считаем, что detB_i≠0 для всех i. В ряде случаев это действительно так. Вместо системы рассмотрим ее регуляризацию. Матрица B_i+1 малой размерности. Регуляризация имеет вид

$\begin{array}{c}\tau \left(B_{i+1}^{\text{*}}{B}_{i+1}+\epsilon \right){\overrightarrow{\alpha }}_{i+1}=-\tau B_{i+1}^{\text{*}}{\text{Φ}}_0{R}_{i+1}^{-1}{\overrightarrow{f}}_{i+1}-B_{i+1}^{\text{*}}{\text{Φ}}_0{R}_{i+1}^{-1}{R}_0{\overrightarrow{C}}_i,\end{array}$ (16)

где B^*_i+1 – сопряженная матрица. Определив вектор α_i+1, найдем вектор C_i+1 из (13).

Опишем некоторые свойства алгоритма. Как видно, если не использовать регуляризацию, то задача решается единственным образом (уравнение (15) имеет единственное решение), если detB_i≠0. Предположим, что это условие выполнено. Пусть P₀ – проектор в ℝ^N, сопоставляющий вектору v⃗ вектор, у которого первые N-r нулевые, а остальные совпадают с координатами v⃗. Его можно задать в виде матрицы, в которой p⁰_ij=δ_ij при i,j=N-r+1,…,N, а остальные элементы p⁰_ij равны нулю. Пусть также P₁ – аналогичный проектор, использующий первые r координат данного вектора: p¹_ij=δ_ij при i,j=1,…,r, а остальные элементы p¹_ij равны нулю. Тогда систему (15) можно записать в виде

$\begin{array}{c}\tau P_0{R}_{i+1}^{-1}{P}_1{\overrightarrow{F}}_{i+1}=\tau P_0{R}_{i+1}^{-1}{\overrightarrow{f}}_{i+1}+P_0{R}_{i+1}^{-1}{R}_0{\overrightarrow{C}}_i,B_{i+1}=P_0{R}_{i+1}^{-1}{P}_1.\end{array}$

Чтобы не загромождать изложение, мы сохранили то же обозначение для новой матрицы B_i+1, которая здесь уже отлична от определенной в (15). Имеем, что detB_i+1=0, но в силу наших предположений kerB_i+1={C⃗:P₁ C⃗=0}. В качестве скалярного произведения в ℝ^N можно рассматривать обычное скалярное произведение (⋅,⋅), а можно и скалярное произведение вида <C⃗,V⃗>=(R₀C⃗,V⃗). В первом случае B^*_i+1=P₁ (R^*_i+1)^-1P₀, во втором случае B^*_i+1=(P₁)^*(R^*_i+1)^-1(P₀)^*. По определению, в первом случае P_i=(P_i)^* и kerB^*_i+1={C⃗:P₀C⃗=0} во втором случае и kerB^*_i+1={C⃗:(P₀)^*C⃗=0} во втором. Оператор B_i+1 отображает H₀=(kerB_i+1)^⊥ на H₁=(kerB^*_i+1)^⊥. Здесь kerB_i+1={C⃗:P₁C⃗=0}, соответственно, (kerB_i+1)^⊥= {v⃗: P₁^*v⃗=v⃗}. Действительно, <v⃗,C⃗⟩=⟨P₁^* v⃗,C⃗⟩=0. Обратно, если ⟨v⃗,C⃗⟩=0 для всех C⃗∈kerB, то C⃗=(I-P₁)g⃗ и ⟨v⃗,(I-P₁)g⃗⟩=0. Тогда (I-P₁^*)v⃗=0. Поскольку матрица B_i+1 сингулярная, берем в качестве решения системы (15) (оно находится неоднозначно) вектор

$\begin{array}{c}\tau {\overrightarrow{\alpha }}_{i+1}=-{\left(B_{i+1}^{\text{*}}{B}_{i+1}\right)}^{-1}{B}_{i+1}^{\text{*}}\left(\tau P_0{R}_{i+1}^{-1}{\overrightarrow{f}}_{i+1}+P_0{R}_{i+1}^{-1}{R}_0{\overrightarrow{C}}_i\right),\end{array}$

который выбирается исходя из включения α⃗_i+1∈H₀, т. е. P₁α⃗_i+1=α⃗_i+1. Подставляя это решение в систему (14), получим

${\overrightarrow{C}}_{i+1}=\tau (R_{i+1}^{-1} {\overrightarrow{f}}_{i+1}-R_{i+1}^{-1}{\left(B_{i+1}^{*}{B}_{i+1}^{}\right)}^{-1}{B}_{i+1}^{*}{R}_{i+1}^{-1}{\overrightarrow{f}}_{i+1})+R_{i+1}^{-1}{R}_0{\overrightarrow{C}}_i-R_{i+1}^{-1}{\left(B_{i+1}^{*}{B}_{i+1}^{}\right)}^{-1}{B}_{i+1}^{*}{P}_o{R}_{i+1}^{-1}{R}_0{\overrightarrow{C}}_i$

Пусть P=R^-1_i+1(B^*_i+1B_i+1)^-1B^*_i+1P₀. Тогда последнее равенство можно записать в виде

$\begin{array}{c}{\overrightarrow{C}}_{i+1}=\tau \left(I-P\right)R_{i+1}^{-1}{\overrightarrow{f}}_{i+1}+\left(I-P\right)R_{i+1}^{-1}{R}_0{\overrightarrow{C}}_i.\end{array}$ (17)

Отметим, что оператор I-P есть проектор, т. е. P²=P. Более того, имеет место равенство P₀ (I-P)=0, т. е. R(I-P) ⊂ kerP₀. Действительно,

$P^2=R_{i+1}^{-1}{\left(B_{i+1}^{*}{B}_{i+1}^{}\right)}^{-1} B_{i+1}^{*}{P}_o{R}_{i+1}^{-1}{\left(B_{i+1}^{*}{B}_{i+1}^{}\right)}^{-1}{B}_{i+1}^{*}{P}_0=R_{i+1}^{-1}{\left(B_{i+1}^{*}{B}_{i+1}^{}\right)}^{-1}{B_{i+1}^{*}}^{}{B}_{i+1}^{}{\left(B_{i+1}^{*}{B}_{i+1}^{}\right)}^{-1}{B}_{i+1}^{*}{P}_0=P.$

Приведенные выше рассуждения показывают (см. равенство (17) Лемма 1. Oператор перехода с одного временного слоя на другой, определенный равенством (17), отличается от стандартного для прямой задачи лишь наличием конечномерного проектора P, размерность которого не увеличивается с уменьшением параметров сетки и параметра τ.

Опишем условия невырожденности матрицы B_i. Пусть функция v_i – решение эллиптической задачи с параметром

$\begin{array}{c}M{v}_i=\frac{1}{\tau }-L\left(t,x\right)v_i=\delta \left(x-d_i\right),{\left.v\right|}_{\text{Γ}}=0. \#\end{array}$

При решении этой задачи методом конечных элементов ищем решение в виде v˜_i=∑^NC_j^j (t)φ_i, где векторы C⃗_i определяются как C⃗_i=τR^-1F_i, где F_i совпадает с базисным вектором i, координата которого совпадает с 1, а остальные координаты – нули. Матрица R зависит от t как от параметра и совпадает в точках t_k=τk с матрицей R_k, определенной выше. В случае регулярного семейства конечных элементов можно показать, что будет иметь место сходимость элементов {v˜_i(t,y_j)}= C_i^j(j=N-r+1,…,N,i=1,2,…,r) к элементам матрицы {v_i(t,y_j)} при каждом t, в том числе и в точках t=t_k (см. теоремы 3.1.5,3.1.6 в [19]). В случае L=Δ, построив матрицу A₀, которая в данном случае будет зависеть от t, а роль параметра λ будет играть роль параметр 1/τ, мы придем к условию:

при всех τ≤τ₀ с некоторым τ₀>0 и всех t∈[0,T],detA₀≠0.

В общем случае возникает условие: при каждом t∈[0,T] и всех τ≤τ₀ с некоторым τ₀>0

$\begin{array}{c}\det {\tilde{A}}_0\ne \mathrm{0,}{\tilde{A}}_0=\left\{v_i\left(y_j\right)\right\},\left(j=N-r+\mathrm{1,}\dots ,N,i=\mathrm{1,2,}\dots ,r\right). \#\end{array}$

Это условие есть достаточное условие разрешимости системы (15). Однако даже при выполнении этого условия матрица B_i все же может быть плохо обусловленной. Поэтому мы используем регуляризацию (см. равенство (16)).

Запись алгоритма в случае граничного условия Неймана или Робина имеет тот же вид.

РЕЗУЛЬТАТЫ И ОБСУЖДЕНИЕ

В этом разделе проанализируем результаты численных экспериментов, полученных в рамках решения прямой и обратной задачи, в рамках экспериментов на ЭВМ. Характеристики испытуемой ЭВМ: процессор: Intel i5 9500; ОЗУ: 16,00 ГБ; тип системы: Windows 10. Численные эксперименты были проведены в программном комплексе MATLAB 2024a.

Возьмем в качестве области G круг |x|=2, примерная сетка указана на рис. 1. Возьмем Т=5 и рассмотрим однородное условие Дирихле. В качестве коэффициентов a₀=(t+1)(x²+1), a₁=2/(t²+1), a₂=(t+1)(x+1), с₁=(t-1)² с₂=2t+1.

 

Рисунок 1. Количество узлов сетки 530

 

Первая группа вычислений

Для первой группы экспериментов берем u=(x²+y²-4)² (t²+ t), ${\alpha }_2=\frac{10+cos\left(3t\right)}{exp\left(t\right)}$. Точки d_i имеют координаты: (0;-0,4) и (-0,4;0).

$f_0=t\left(x^2+y^2\right)-1-(\left(x\right)\left(\left(t^2+2\right)\right)-\left(\left(y\right)\left({\left(t-2\right)}^3\right)\right)$.

Возьмем в качестве параметра ε_a=max_i (max_j |α^j(iτ)-α_i^j|) – ошибка вычисления функции α_i в граничном условии, где α_i^j≈α^j(τi) – результаты вычисления. Процесс вычисления прекращается после того, как ε_α становится меньше ε₀=10^-3.

Результат восстановления искомых функций представлен на рисунке 2. Зеленым обозначены найденные функции, красным – эталонные.

 

Рисунок 2. Результат вычисления первой группы

 

Далее были добавлены источники генерации шума – в возмущение значений «данных с датчиков» добавляется случайный шум noise ∈(-δ;δ), где δ – процентное выражение отклонения. Итоговая функция принимает вид ψ_i(jτ)=ψ_i(jτ)(1+ noise(jτ)), j=1,…,M. Мы проводим расчеты с различным уровнем шума (добавляемый шум строится при помощи функции random пакета «Матлаб»). Результат вычислений с шумом 3 % (δ=0.03) представлен на рисунке 3 и с шумом 10 % (δ=0.1) – на рисунке 4.

 

Рисунок 3. Результат вычисления первой группы с добавлением шума 3 %

 

Рисунок 4. Результат вычисления первой группы с добавлением шума 10 %

 

Вторая группа вычислений

Для второй группы экспериментов u=cosπ(x²+y²)t построение ведется на той же сетке (рисунок 1), возьмем те же самые начальные условия и параметры уравнения, функции α_i имеют вид: α₁=(t-1)², α₂=t²+1. Коэффициенты те же самые, что и в первом случае. f₀=(x+1)(y+1)(t+1)-((x)(t²+2)-((y)((t-2)³)).

Для повышения точности было решено увеличить точность расчетов и взять значение ε₀=10^-4. Результат построения представлен на рисунке 5, зеленым обозначены найденные функции, красным – эталонные. Как можно заметить, графики сходятся намного лучше.

 

Рисунок 5. Результат вычисления второй группы

 

Рисунок 6. Результат вычисления второй группы с добавлением шума 3 %

 

Рисунок 7. Результат вычисления второй группы с добавлением шума 10 %

 

Для данных 2 группы было проведено 15 экспериментов, по результатам которых собрана следующая статистика: средняя ошибка равна 0,0035, среднее время выполнения – 2,2 с; при добавлении шума 3 % показатели такие: средняя ошибка равна 0,0045, среднее время выполнения – 3,5 с; при добавлении шума 10 % средняя ошибка равна 0,0037, время выполнения – 5,1 с.

ЗАКЛЮЧЕНИЕ И ВЫВОДЫ

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

Авторлар туралы

Sergey Shmorgan

Yugra State University

Хат алмасуға жауапты Автор.
Email: s_shmorgan@ugrasu.ru

Postgraduate Student, Engineering School of Digital Technologies

Ресей, Khanty-Mansiysk

Lyubov Neustroeva

Yugra State University

Email: l_neustroeva@ugrasu.ru

Candidate of Physical and Mathematical Sciences, Researcher

Ресей, Khanty-Mansiysk

Әдебиет тізімі

Қосымша файлдар