\documentclass[10pt]{amsart}
\usepackage[cp1251]{inputenc}
\usepackage[english,russian]{babel}
\usepackage{amsmath}
\usepackage{amssymb}
\usepackage{amsfonts}
\usepackage{graphicx}
\usepackage{subfig}
%\usepackage[left=45mm, top=35mm, right=45mm, bottom=35mm]{geometry}
\graphicspath{{graf/}}
\usepackage{capt-of}

\usepackage[usenames]{color}
\usepackage{colortbl}

\def\udcs{519.6} %Здесь автор определяет УДК своей работы
\def\mscs{65M32:65R32} %Здесь автор определяет классификаторы AMS своей работы
\setcounter{page}{144}

\newtheorem{lemma}{Лемма}
\newtheorem{theorem}{Теорема}
\newtheorem{definition}{Определение}
\newtheorem{corollary}{Следствие}
\newtheorem{proposition}{Предложение}

\def\logo{{\bf\huge S\raisebox{0.2ex}{\hspace{0.55ex}\raisebox{0.05ex}e\hspace{-1.65ex}$\bigcirc$}MR}}

\def\semrtop
     {
  \vbox{
     \noindent\logo
     \hspace{80mm}\raisebox{1ex}{ISSN 1813-3304 }

     \vspace{5mm}

     \begin{center}
     {\huge СИБИРСКИЕ \ ЭЛЕКТРОННЫЕ} \\[2mm]
     {\huge МАТЕМАТИЧЕСКИЕ ИЗВЕСТИЯ} \\[2mm]
     {\large Siberian Electronic Mathematical Reports} \\[1mm]
     {\LARGE\tt{http://semr.math.nsc.ru}}\\[0.5mm]
%     {\small Институт математики им. С.Л.Соболева СО РАН}\\[-1mm]
%     {\small Sobolev Institute of Mathematics SB RAS}
     \end{center}
     \vspace{-3mm}
     \noindent
     \begin{tabular}{c}
     \hphantom{aaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaa} \\
     \hline\hline
     \end{tabular}

     \vspace{1mm}
     {\flushleft\it Том 19, стр. 144--144 (2022) \hspace{65mm}{\rm\small УДК \udcs}}
     \newline
     	{\rm\small DOI~10.33048/semi.2022.19.xxx}\hphantom{aaaaaaaaaaaaaaaaaaaaaaaaaaaaaaa}{\rm\small MSC\ \ \mscs }
  }%\vbox
}

\begin{document}
\renewcommand{\refname}{References}

\thispagestyle{empty}

\title[Итерационное решение ретроспективной обратной задачи]{Итерационное решение ретроспективной обратной задачи теплопроводности с помощью интеграла Пуассона}
\author{{В.И. Васильев, А.М. Кардашевский }}%
\address{Vasil'ev V.I., Kardashevsky A.M.  % Контактные данные всех авторов и место работы указываются только на английском языке
\newline\hphantom{iii} North-Eastern Federal University,
\newline\hphantom{iii} str. Belinskogo, 58,
\newline\hphantom{iii} 677000, Yakutsk, Russia}%
\email{vasvasil@mail.ru, kardam123@gmail.com}%

\thanks{\sc Vasil'ev V.I., Kardashevsky A.M.
Iterative solution of the retrospective inverse heat conduction problem using the Poisson integral}
\thanks{\copyright \ 2021 Васильев В.И., Кардашевский А.М}
\thanks{\rm Работа В.И.Васильева выполнена при поддержке мегагранта Правительства РФ, соглашение № 14.Y26.31.0013. Работа А.М.Кардашевского выполнена при поддержке Минобрнауки РФ, соглашение от 02.02.2022 № 075-02-2022-881.}
\thanks{\it Поступила  2021 г., опубликована  2021 г.}%


\semrtop \vspace{1cm} \maketitle {\small
\begin{quote}
\noindent{\bf Abstract.}  This paper considers the inverse problem of identification of the finite initial condition of the Cauchy problem for the homogeneous heat conduction equation using the first kind linear Fredholm integral equation. 
Its discretization is carried out with the help of the quadrature rectangular formula. For the numerical realization of the obtained system of linear algebraic equations with almost complete, symmetric, positively determined, ill-conditioned matrix it is proposed to use the method of conjugate gradients. 
Examples of reconstruction of smooth, nonsmooth and discontinuous initial conditions in one- and two-dimensional cases, including the introduction of <<noise>>, characteristic of redefinition condi- tions of inverse problems, are given.
\medskip

\noindent{\bf Keywords:} retrospective inverse heat conduction problem, Poisson inte- gral, first kind Fredholm integral equation, system of linear equations with ill-conditioned matrix, method of conjugate gradients.
\end{quote}
}



\begin{center}
{\bf 1. Введение}
\end{center}

Обратные задачи для дифференциальных уравнений с частными производными имеют многочисленные приложения в науке и технике. Они, как правило, не относятся к классу корректно поставленных задач. Типичными для них являются неединственность решения и нарушение требования непрерывной зависимости решения от входных данных. С точки зрения приложений актуальной является разработка эффективных вычислительных методов решения обратных задач, они востребованы  физиками, геологами, биологами и медиками. Неудивительно, что при таком широком наборе приложений теория обратных и некорректных задач с момента своего появления стала одной из наиболее стремительно развивающихся разделов современной вычислительной и прикладной математики. Многочисленные публикации по вопросам теоретического исследования корректности, условной корректности, разработки эффективных численных методов решения и актуальные приложения обратных задач в науке и технике обобщены в монографиях \cite{LL} -- \cite{I}. Отметим также статью С.И.Кабанихина \cite{K1}, в которой приведен обзор современных численных методов решения обратных задач, в том числе градиентных, стохастических алгоритмов, применение генетических алгоритмов и нейронных сетей.
 
Приведем краткий обзор итерационных методов решения дискретных аналогов обратных задач для уравнений математической физики. В статьях проф. Кабанихина С.И., Баканова Г.Б. \cite{Kab1} -- \cite{Kab2} предложен оптимизационный итерационный метод решения дискретной обратной задачи для гиперболического уравнения и доказано, что функционал невязки дискретной обратной задачи имеет единственную стационарную точку, совпадающую с решением обратной задачи. Проф. А.Л.Карчевским в работе \cite{Karch}  доказано, что нелинейный функционал невязки конечно-разностного аналога коэффициентной обратной задачи имеет единственную точку глобального минимума, получена оценка скорости сходимости к решению исходной обратной задачи.  
  
Особое место в обратных задачах теплопередачи занимает ретроспективная обратная задача определения начального условия $u(x,0)$ по заданному распределению температуры тела в финальный момент времени $u(x,T)$. 
В работе \cite{SVV} предложен итерационный метод численной реализации конечно-разностного аналога ретроспективной обратной задачи теплопроводности. В ней на каждой итерации решается дискретный аналог корректной прямой задачи теплопроводности.
В работе \cite{LE} предложен и теоретически обоснован итерационный метод граничных элементов.
Для численного решения ретроспективной обратной задачи для  уравнения Бюргерса в работе \cite{KT} предложен вариационный  метод, позволяющий  свести исходную задачу к серии прямых, следовательно, корректных задач.  
Бессеточная численная схема решения обратной задачи теплопроводности предложена в статье \cite{M}. В ней метод численного решения базируется на использовании в качестве базисной функции фундаментального решения уравнения теплопроводности. Численная реализация полученной системы линейных алгебраических уравнений с плохо обусловленной матрицей проведена с помощью стандартного метода регуляризации А.Н.Тихонова.
В статье \cite{JD} предложен итерационный алгоритм устойчивой идентификации как начального условия, так и стационарного источника, основанный на последовательном решении прямых задач для уравнения  теплопроводности, которые решаются на каждой итерации с использованием метода граничных элементов.  
В статье \cite{IK} представлено сравнение численных методов решения ретроспективной обратной задачи для нелинейной эволюционной среды. 

В работе \cite{ZZ} предложен модифицированный метод регуляризации А.Н.~Тихонова, в которой представлены правила апостериорного выбора параметра регуляризации. 
В работе \cite{VK} дискретизация обратной задачи для параболического уравнения выполняется по разностной схеме с неположительным весовым множителем, а в \cite{VK1}, предложен численный метод решения ретроспективной обратной задачи для многомерного параболического уравнения с помощью метода сопряженных градиентов. 

В \cite{TS} изучается алгоритм регуляризации приближенного решения интегральных уравнений Фредгольма первого рода. 
Решение прямой задачи Коши для однородного уравнения теплопроводности задается в виде интегральной формулы Пуассона (интеграла Пуассона). Ретроспективную обратную задачу теплопроводности в случае задания финального условия в виде финитной функции с помощью интеграла Пуассона можно свести к интегральному уравнению Фредгольма 1-го рода \cite{BT}, \cite{D}. 
Дискретный аналог интегрального уравнения, полученный с помощью той или иной кубатурной формулы, приводит к системе линейных алгебраических уравнений с полной матрицей. При измельчении сетки матрица быстро становится плохо обусловленной (сингулярной), и поэтому прямые методы работают только на достаточно крупной сетке. 
Авторы в работе \cite{VK2} для обхода этой ситуации, воспользовавшись высокой точностью квадратурной формулы для приближенного вычисления интеграла даже на достаточно грубой сетке, предложили метод многократного применения квадратурной формулы на множестве крупных сеток, каждая из которых решается прямым методом с последующим обьединением результатов. 

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



\medskip
\begin{center}
{\bf 2. Постановка задачи}
\end{center}

Рассмотрим ретроспективную обратную задачу Коши для уравнения теплопроводности, записанную в безразмерных переменных
\begin{equation}\label{eq:1}
 \frac{\partial u}{\partial t}= \sum_{\alpha=1}^p \frac{\partial^2 u}{\partial x_\alpha^2},  \quad {\mathbf x} \in \mathcal{R}^p ,\quad 0 \leqslant t < T.  
\end{equation} 

Для определения единственного решения искомой функции задано значение искомой функции в финальный момент времени    
\begin{equation}\label{eq:2}
  u({\bf x}, T) = \varphi({\bf x}), \quad {\mathbf x} {\in} \overline{\Omega}.  
\end{equation}

В задаче Коши (\ref{eq:1}) -- (\ref{eq:2}) подлежит определению начальное условие $u({\bf x}, 0)=f(x), \ {\mathbf x} {\in} \overline{\Omega}$.
Отметим, что уравнение (\ref{eq:1}) имеет аналитическое решение, которое в дальнейшем понадобится для оценки точности полученного численного решения
\begin{equation}\label{eq:3}
  u({\mathbf x},t)=\frac{1}{\sqrt{(\beta+4t)^p}}e^{-\frac{\|{\mathbf x}\|^2}{\beta+4t}}, \quad \beta>0, \quad {\mathbf x} \in \mathcal{R}^p ,\quad 0 \leqslant t \leqslant T. 
\end{equation}

Пусть $\Omega \in \mathcal{R}^p$ представляет собой $p$-мерный куб,
\[
\Omega = \{{\mathbf x} \ | {\mathbf x} =\left(x_1,\dots,x_p\right), \quad  0 < x_\alpha < L, \quad \alpha = 1,\dots,p\}.
\]

Предположим, что финальное условие $\varphi({\mathbf x})$ является финитной, отличной от нуля в  некоторой подобласти $\Omega \subset \mathcal{R}^p$. 
Решение прямой задачи для однородного параболического уравнения (\ref{eq:1}) с начальным условием $u({\mathbf x},0)=f({\mathbf x})$ задается в виде  интеграла Пуассона
\begin{equation}\label{eq:4}
  \varphi({\bf x})=\frac{1}{(2\sqrt{\pi T})^p}\int\limits_{\overline{\Omega}} e^{-\frac{\|{\mathbf y}-{\mathbf x}\|^2}{4T}}f({\mathbf y})d{{\mathbf y}}, \qquad {\mathbf x} {\in} \overline{\Omega}.
\end{equation}
Здесь $T>0$ -- финальное время.

В $p$ -- мерном кубе $\overline\Omega$ введем кубическую сетку $\overline\omega$ с шагом $h$
$$  
  \overline\omega = \prod_{\alpha=1}^p \overline\omega_\alpha, \qquad
  \overline\omega_\alpha=\{x_{\alpha i}=-L+i_\alpha h, \ i_\alpha=0, 1, \ldots, n; \ h=2L/n\}.
$$

Если рассматриваемая задача не имеет аналитического решения, то для приближенного определения финального решения следует использовать дискретный аналог интеграла Пуассона (\ref{eq:4}) в узлах достаточно подробной сетки
\begin{equation}\label{eq:5}
\varphi({\bf x}_j)=\frac{h^p}{(2\sqrt{\pi T})^p}\sum_{{\bf y}_i\in\overline{\omega}} e^{-\frac{\|{\bf x}_i-{\bf y}_j\|^2}{4T}}f({\bf y}_i), \quad {\bf x}_j \in\overline{\omega}.
\end{equation}

По кубатурной формуле (\ref{eq:5}) можно с необходимой точностью, достаточно быстро, находить приближенное решение прямой задачи Коши с заданным начальным условием $f({\mathbf x}), \ {\mathbf x}\in \overline\omega$ в момент времени $T>0$ на пространственной сетке $\overline\omega$. Вычислительный процесс является устойчивым, а полученное численное решение имеет достаточно высокую точность даже на грубых сетках.    

В случае ретроспективной обратной задачи (\ref{eq:1}) -- (\ref{eq:2}), когда по заданному финальному условию $\varphi ({\bf x}), \ {\bf x}\in \overline\omega$ подлежит определению начальное условие $f({\bf x}), \ {\bf x}\in \overline\omega$, для определения начального условия $f({\bf x}), \ {\bf x}\in \overline\omega$ получаем интегральное уравнение Фредгольма первого рода :
\begin{equation}\label{eq:6}
  \frac{1}{(2\sqrt{\pi T})^p}\int\limits_{\Omega} e^{-\frac{\|{\mathbf y}-{\mathbf x}\|^2}{4T}}f({\mathbf y})d{{\mathbf y}}=\varphi({\bf x}), \qquad {\mathbf x} {\in} \overline{\Omega}.
\end{equation}

Следовательно, дискретный аналог рассматриваемой обратной  задачи для определения сеточного аналога начального условия в узлах пространственной сетки приводит к необходимости решения следующей системы линейных алгебраических уравнений: 
\begin{equation}\label{eq:7}
\sum_{{\bf y}_i\in\overline{\omega}} e^{-\frac{\|{\bf x}_i-{\bf y}_j\|^2}{4T}}f({\bf y}_i)=g\varphi({\bf x}_j), \quad {\bf x}_j \in\overline{\omega}.
\end{equation}
где $g=\Big(2\sqrt{\pi T}/h\Big)^p$.

Систему линейных алгебраических уравнений (\ref{eq:7}) перепишем в векторно-матричной форме
\begin{equation}\label{eq:8}
   a_{ij}\varphi_i=gf_j, \quad j=0,1,\ldots,N, \quad N=(n+1)^p.
\end{equation}
Здесь использованы обозначения  
$$
   a_{ij}=e^{-s_{ij}},  \quad i=0,...,N, \quad f_j=f({\bf y}_j), \quad  j=0,...,N,
$$
где 
$$
  s_{ij}=\frac{1}{4T}\displaystyle{\sum_{\alpha=1}^p {(y_{\alpha j}-x_{\alpha i})^2}},  \ 
  j=\sum_{{\alpha=1}}^p (n+1)^{\alpha-1} j_\alpha, \ 
  i=\sum_{{\alpha=1}}^p (n+1)^{\alpha-1} i_\alpha.
$$

Построенная таким образом матрица $\mathcal{A}$ c элементами $a_{ij}, \ i=0,...,N; \ j=0,...,N$  является квадратной, почти плотной, симметричной, все ее элементы неотрицательные, элементы стоящие на главной диагонали равны единице, остальные элементы по мере удаления от главной диагонали убывают по экспоненциальному закону. Более того, матрица $\mathcal{A}$ на решении системы уравнений (\ref{eq:8}) является положительно определенной. Следует отметить, что она не обладает свойством диагонального преобладания, вследствие чего определитель матрицы даже на достаточно грубых сетках может обратиться в нуль.
\begin{table}[h!]
\center
\begin{tabular}{ | c | c | c | c | c | c | } \hline
N  & 20 & 40 & 100 & 200 & 500 \\ \hline
$Det(\mathcal{A})$ & 3.950e$-$52 & $-$3.123e$-$317  & --- & --- & --- \\  \hline
$Cond(\mathcal{A})$ & 2.899e+9 & 3.332e+18  & 3.126e+18 & 1.990e+19 & 6.895e+19 \\ \hline 
\end{tabular}
\bigskip
\caption{Определитель и число обусловленности матрицы 
$\mathcal{A}$ для разных значений числа неизвестных $N$.}
\label{table1}
\end{table}

В таблице 1 представлены определитель $Det(\mathcal{A})$ и число обусловленности $Cond(\mathcal{A})$ матрицы $\mathcal{A}$, полученные при дискретизации интегрального уравнения Фредгольма первого рода (\ref{eq:6}), эквивалентного ретроспективной обратной задаче для одномерного параболического уравнения (\ref{eq:1}) -- (\ref{eq:2}) на разных пространственных сетках, т.е. в зависимости от количества разбиений области определения задачи {$N$}. Из таблицы видно, что при $N>40$, прямые методы решения системы линейных алгебраических уравнений не работают, поскольку определитель матрицы $\mathcal{A}$ обращается в машинный нуль. Для решения системы уравнений с матрицей, число обусловленности которой чрезвычайно велико, а определитель равен машинному нулю, целесообразно использовать метод сопряженных градиентов. Приведенные ниже результаты вычислительного эксперимента подтвердили справедливость данного утверждения.

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

Искусственное зашумление финального условия осуществим с помощью следующей функции 
\begin{equation}\label{eq:9}
   \tilde{\varphi}({\bf x})=\varphi({\bf x})+\sigma \cdot rand(-0.5, 0.5), \quad {\bf x} \in  \omega,   
\end{equation} 
где $\sigma=0,01$ -- амплитуда задаваемой ошибки генератора случайных чисел $rand(-0.5, 0.5)$, равномерно распределенных на интервале $(-0.5,0.5)$. 

В одномерном случае для сглаживания зашумленного финального условия в конечный момент времени $\varphi(x)$ целесообразно пользоваться фильтром Савицкого-Голея \cite{SG}, осуществляющим полиномиальную аппроксимацию входного сигнала построенным с помощью метода наименьших квадратов
\begin{equation}\label{eq:10}
   \overline{\varphi}(x)=savgolfilter(\tilde{\varphi}(x), n, m), \quad x \in \omega.     
\end{equation}

В двумерном случае для сглаживания зашумленного финального условия,  используем двумерный аналог фильтра Савицкого-Голея, предложенный в работе B.Gorbi с соавторами \cite{GHS} и реализованный в свободной библиотеке программ \cite{GH} 
\begin{equation}\label{eq:11}
   \overline{\varphi}({\bf x})=sgolay2.SGolayFilter2(window{\textunderscore}size=n,poly{\textunderscore}order=m)(\tilde{\varphi}({\bf x})), \ {\bf x} \in \omega.     
\end{equation}

В формулах (\ref{eq:10}) -- (\ref{eq:11}) параметр $n$ -- задает количество точек в окрестности соответствующего узла сетки, используемых при сглаживании полиномом порядка $m$. 



\medskip
\begin{center}
{\bf 3. Результаты численного решения модельных \\ одномерных  ретроспективных обратных задач}
\end{center}

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

\medskip
{\bf Модельная задача 1:} Будем определять численное решение ретроспективной обратной задачи (\ref{eq:1}) -- (\ref{eq:2}) посредством решения системы линейных уравнений (\ref{eq:8}) -- дискретного аналога интегрального уравнения (\ref{eq:7}), являющегося эквивалентной постановкой задачи (\ref{eq:1}) -- (\ref{eq:2}). Для численной реализации системы линейных алгебраических уравнений (\ref{eq:7}) с симметричной и положительно определенной матрицей $\mathcal{A}$ используем итерационный метод сопряженных градиентов. 

Итак, ищем решение следующей ретроспективной обратной задачи 
\begin{equation}\label{eq:12}
 \frac{\partial u}{\partial t}= \frac{\partial^2 u}{\partial x^2},  \quad -L\leqslant x \leqslant L, \quad 0 \leqslant t < T,  
\end{equation}
\begin{equation}\label{eq:13}
   u(x,T)=\varphi(x),  \quad -L\leqslant x\leqslant L,
\end{equation}
где финальное условие задано в виде финитной функции 
\begin{equation}\label{eq:14}
  \varphi(x)=\frac{1}{\sqrt{1+4T}}e^{-\frac{x^2}{1+4T}}\geqslant 0 ,  \quad -L\leqslant x\leqslant L. 
\end{equation}
Ретроспективная обратная задача теплопроводности (\ref{eq:12}) -- (\ref{eq:14}) имеет точное аналитическое решение, задаваемое формулой (\ref{eq:3}) при $\beta=1$. 
Следовательно, искомое начальное условие задается формулой 
$$
  u(x,0)=f(x)=e^{-x^2}\geqslant 0 ,  \quad -L\leqslant x\leqslant L.
$$
\begin{figure}[htp]
{\includegraphics[width=1.\columnwidth]{F1.eps}}
\caption{Восстановление гладкого начального условия}
\end{figure}

Численные расчеты проводились при $L=6, \ T=1, \ N=120.$ На рис.~1 представлены результаты вычислительного эксперимента по восстановлению гладкого начального условия $f(x)$ с помощью итерационного метода сопряженных градиентов. На рис.1A показано отклонение численного решения, вычисленное по явной формуле (\ref{eq:5}), от точного аналитического решения в финальный момент времени $t=T$:
$$
  R=\frac{\|y^J-\varphi\|}{\|\varphi\|}.
$$ 

На рис.1В приведены: 1 -- график точного начального условия $f(x)$, $2$ -- график точного финального условия $\varphi(x)$ в момент времени $T$, $3$ -- график восстановленного начального условия, $4$ -- график погрешности идентификации начального условия $z=y-u0$. 
На рис.1С представлены те же графики с той лишь разницей, что на графике $3$ -- приведено восстановленное начальное условие в случае, когда зашумленное финальное условие $\varphi(x)$ получено с помощью формулы (\ref{eq:9}) с среднеквадратическим отклонением $\sigma=0.01$. 
Аналогично на рис.1D представлены те же графики, где $3$ -- график восстановленного начального условия, полученное при зашумленном финальном условии и сглаженное с помощью фильтра Савицкого-Голея (\ref{eq:10}) c параметрами $n=25, \ m=4$.

Дискретный аналог финального условия определен с чрезвычайно высокой точностью, хотя сетка была недостаточно подробной $N=120$. Вследствие этого точность восстановления начального условия также достаточно высокая, чем больше итераций, тем выше точность. На этом примере мы ограничились 70 итерациями. При вводе возмущения финального условия с помощью формулы (\ref{eq:10}) с достаточно большим $\sigma=0.01$ восстановленная функция ожидаемо сильно искажается. Применение фильтра Савицкого-Голея на возмущенную финальную функцию существенно повысило точность восстановления начального условия. В численных расчетах с зашумленными финальными условиями итерации продолжались до тех пор, пока невязка итерационного процесса убывает монотонно. Получилась так, что в обеих случаях проведено по 5 итераций. 

\medskip
{\bf Модельная задача 2}. Будем идентифицировать начальное условие $f(x)$, представляющее собой негладкую функцию с тремя точками разрыва первой производной $x=-L/3$,  $x=0$ и $x=L/3$: 
\begin{equation}\label{eq:15}
f(x)= \begin{cases}
0, & \text{при} \ \ x\in [-L,-L/3)\cup(L/3,L], \\
1+3x/L, & \text{при} \ \ -L/3 < x< 0 \\ 
1-3x/L, & \text{при} \ \ 0 < x<L/3  \\
\end{cases} 
\end{equation}


\begin{figure}[htp]
{\includegraphics[width=1.\columnwidth]{F2.eps}}
\caption{Восстановление негладкого начального условия}
\end{figure}

Расчеты проводились при $L=6, \ T=1, \ N=120.$ На рис.~2 представлены результаты вычислительного эксперимента по восстановлению негладкого начального условия $f(x)$, заданного по формуле (\ref{eq:15}). Данная задача не имеет точного аналитического решения. Численная реализация интеграла Пуассона по квадратурной формуле (\ref{eq:5}) является корректной задачей, поэтому точность определения финального условия растет по мере измельчения сетки. За эталонное финальное условие возьмем решение, полученное на сетке с минимальным шагом. На рис.~2A показана разница вычисленных по формуле (\ref{eq:5}) финальных условий на крупной и на более мелкой сетке $N=120$ и $N=240$. Обратную задачу будем решать на сетке с $N=120$. 
На рис.2В приведены: 1 -- график заданного негладкого начального условия $f(x)$, $2$ -- график финального условия $\varphi(x)$ в момент времени $T$ на сетке с $N=240$, $3$ -- график восстановленного начального условия, $4$ --  погрешность идентификации начального условия $z=y-u0.$
На рис.2С представлены графики восстановленного негладкого начального условия и погрешность идентификации, вычисленных при задании зашумленного финального условия $\tilde\varphi(x)$ с амплитудой задаваемого шума $\sigma=0.01$. 
Графики $3$ и $4$ на рис.2D, соответственно, восстановленного начального условия и погрешности идентификации в случае когда зашумленное финальное условие $\overline\varphi(x)$ с возмущением сглажено с помощью фильтра Савицкого-Голея с $n=25, m=4$.

Из представленных результатов счета видно, что разница решения в финальный момент времени на выбранных сетках имеет порядок $R=10^{-4}$. Наибольшая ошибка определения финального условия находится в центре области определения задачи. Выход из итерационного процесса восстановления начального условия осуществлено на 70--ой итерации. Из рис.2D видно, что точность восстановления приемлемая. Введенное достаточно ощутимое случайное возмущение финального условия ожидаемо сильно влияет на идентификацию начального условия, которое имеет сильно осциллирующий характер, дает ошибку порядка $\pm0.3$. Фильтр Савицкого-Голея улучшает результат кардинально, получается сильно сглаженное начальное условие и максимальное отклонение от искомого начального условия имеет величину не превышающую $\pm0.1$.


\medskip
{\bf Модельная задача 3}. Рассмотрим в качестве искомого начального условия разрывную функцию с двумя точками разрыва при $x=-L/3$ и $L/3$: 
\begin{equation}\label{eq:16}
f(x)=
\begin{cases}
0, & \text{при} \ \ x\in [-L,-L/3]\cup[L/3,L], \\
1, & \text{при} \ \ -L/3 < x< L/3. \\ 
\end{cases} 
\end{equation}

\begin{figure}[htp]
{\includegraphics[width=1.\columnwidth]{F3.eps}}
\caption{Восстановление разрывного начального условия}
\end{figure}

Расчеты проводились при $L=6, \ T=1, \ N=120.$ На рис.3 представлены результаты вычислительного эксперимента по восстановлению разрывного начального условия обратной задачи для параболического уравнения (\ref{eq:12}). Она не имеет аналитического решения. Поэтому для задания эталонного финального условия $\varphi(x)$, соответствующего начальному условию (\ref{eq:16}) решаем прямую задачу для уравнения (\ref{eq:12}) с начальным условием (\ref{eq:16}) на более мелкой сетке $N=240$. 
На рис.3A показана разница решений прямой задачи, полученных с помощью дискретного аналога интеграла Пуассона (\ref{eq:5}) на крупной сетке с $N=120$ и на более мелкой сетке с $N=240$. 
На рис.3В $1$ -- приведены графики восстановленного разрывного начального условия $f(x)$, $2$ -- график финального условия $\varphi(x)$ в момент времени $T$, $3$ -- восстановленное начальное условие. Линия $4$ -- погрешность идентификации начального условия $z=y^0-f$. 
На рис.3С представлены те же графики, здесь $3$ -- восстановленное начальное условие в случае, когда условие переопределения $\varphi(x)$ задано с возмущением $\sigma=0.01$. На рис.3D представлены аналогичные графики, при этом график $3$ -- восстановленное начальное условие когда финальное условие с <<шумом>> пропущено через фильтр Савицкого-Голея $n=25, m=4$.

Для оценки точности восстановления разрывного начального условия провели вычисление на разных сетках с $N=120$ и $N=240$, она оказалась величиной порядка $10^{-3}$, поэтому за эталонное финальное условие берем решение, полученное на сетке с $N=240$. За счет введения возмущения получается сильно осциллирующее начальное условие, а введение фильтра Савицкого-Голея дает сглаженное и не сильно осциллирующее начальное условие.
Отметим, что использование подробных сеток не приводит к существенному улучшению точности восстановления разрывного начального условия, поскольку в точках разрыва точность восстановления разрывного начального условия невозможно улучшить. 

\medskip
\begin{center}
{\bf 4. Численная апробация итерационного метода \\ на двух двумерных модельных задачах}
\end{center}

Проведем апробацию предложенного итерационного метода численного восстановления начального условия с помощью метода сопряженных градиентов (\ref{eq:8}) на двух двумерных модельных задачах (\ref{eq:17}) -- (\ref{eq:20}). 

\medskip
{\bf  Модельная задача 4.} Приведем результаты численной реализации предложенного итерационного метода решения ретроспективной обратной задачи на примере двумерной задачи с гладким начальным условием
\begin{equation}\label{eq:17}
 \frac{\partial u}{\partial t}= \frac{\partial^2 u}{\partial x^2}+\frac{\partial^2 u}{\partial y^2},  \quad -l\leqslant x \leqslant l,\quad -l\leqslant y \leqslant l, \quad 0 \leqslant t < T,  
\end{equation}

\begin{equation}\label{eq:18}
   \varphi(x,y)=\frac{1}{\sqrt{1+4T}}e^{-\frac{x^2+y^2}{1+4T}},  \quad -l\leqslant x\leqslant l, \quad 
   -l\leqslant y\leqslant l.
\end{equation}

\begin{figure}[htp]
{\includegraphics[width=0.8\columnwidth]{F4.eps}}
\caption{Восстановление гладкого начального условия}
\end{figure}

Ретроспективная обратная задача (\ref{eq:17}) -- (\ref{eq:18}) имеет точное аналитическое решение
\begin{equation}\label{eq:19}
  u(x,y,t)=\frac{1}{\sqrt{1+4t}}e^{-\frac{x^2+y^2}{1+4T}}. 
\end{equation}

На рис.4 представлены результаты численной идентификации гладкого начального условия на примере ретроспективной обратной задачи (\ref{eq:17}) -- (\ref{eq:18}). Расчеты проводились при $N=40;\ L=6; \ T=1.$ Выбрали квадратную пространственную сетку $x_i=-L+ih, \ y_i=ih, \ i=0,1,\ldots,N, h=2L/N$.

На рис.4A приведен график искомого точного начального  $f(x,y)$. 
На рис.4В представлен график восстановленного начального условия. 
На рис.4С приведен график восстановленного начального условия при введении искусственного возмущения (\ref{eq:11}) с $\sigma=0.01$ в финальное условие $\varphi(x,y)$. 
На рис.~4D представлено восстановленное начальное условие в случае, когда условие возмущенного финального условия сглажено с помощью двумерного фильтра Савицкого-Голея \cite{GHS} c параметрами $n=9, \ m=4$.
 Так как задача с таким начальным условием имеет аналитическое решение имеем возможность использовать точное финальное условие для восстановления начального условия и сравнения результатов.
 
\medskip
{\bf Модельная задача 5.} Проведем вычислительный эксперимент по итерационному восстановлению разрывного начального условия для дифференциального уравнения (\ref{eq:17}), имеющего форму толстостенного цилиндра единичной высоты с внутренним и внешними радиусами $2$ и $4$, соответственно
\begin{equation}\label{eq:20}
u(x,0) = f(x)=
\begin{cases}
1, & \text{при $4<x^2+y^2<16$,} \\
0, & \text{в остальных случаях}. \\ 
\end{cases}
\end{equation}

\begin{figure}[htp]
{\includegraphics[width=0.8\columnwidth]{F5.eps}}
\caption{Восстановление  разрывного начального условия}
\end{figure}

На рис.5 представлены результаты численного восстановления разрывного начального условия, заданного по формуле (\ref{eq:20}). Задача не имеет точного аналитического решения. Поэтому для оценки точности финального условия расчеты прямой задачи (\ref{eq:17}), (\ref{eq:20}) проводились на мелкой и крупной сетках с $N=32, 64; \ L=8, \ T=1.$ Пространственная сетка квадратная $x_i=-L+ih, \ y_i=ih, \ i=0,1,\ldots,N, h=2L/N$. Максимум разности результатов на мелкой и крупной сетках получилось порядка 0.001 вследствие этого в качестве финального условия выбрано решение при  $N=64$, как более точное. На рис.5A показана искомое начальное условие. На рис.5В приведен график восстановленного начального условия. На рис.5С представлен график восстановления когда финальное условие $\varphi(x,y)$ зашумлено с помощью формулы (\ref{eq:9}) с амплитудой задаваемого возмущения $\sigma=0.01$. На рис.5D представлен график восстановления начального условия в случае, когда финальное условие $\tilde\varphi(x,y)$ заданное с искусственным <<шумом>>  $\sigma=0.01$ сглажено с помощью двумерного фильтра Савицкого-Голея (\ref{eq:11}) с параметрами $n=9, \ m=4$.

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



\medskip
\begin{center}
{\bf 6. Заключение}
\end{center}

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

\bigskip
\begin{thebibliography}{1}

%\begin{thebibliography}{99}
\bibitem{LL} R.~Lattes and J.L.~Lions, \textit{The Method of Quasi-Reversibility}, Amer. Elsevier Publ., New York, 1969.

\bibitem{IVT} В.К.~Иванов, В.В.~Васин, В.П.~Танана,  \textit{Теория линейных некорректных задач и ее приложения}, Наука, Москва, 1978.

\bibitem{TA} А.Н.~Тихонов, В.Я.~Арсенин, \textit{Методы решения некорректных задач}, Наука, Москва, 1979.

\bibitem{BT} А.Б.~Васильева, Н.А.~Тихонов, \textit{Интегральные уравнения}, Изд--во МГУ, Москва, 1989. 

\bibitem{TG} А.Н.~Тихонов, А.В.~Гончарский, В.В.~Степанов, А.Г.~Ягола, \textit{Численные методы решения некорректных задач},  Наука, Москва, 1990.

\bibitem{D} А.М.~Денисов, \textit{Введение в теорию обратных задач}, Изд--во МГУ, Москва, 1994.

\bibitem{LR} М.М.~Лаврентьев, В.Г.~Романов, С.Т.~Шишатский, \textit{Некорректные задачи математической физики и анализа}, Наука, Москва, 1980.

\bibitem{R} В.Г.~Романов, \textit{Обратные задачи математической физики}, АН СССР, Москва, 1984.

\bibitem{SB} А.А.~Самарский, П.Н.~Вабищевич, \textit{Численные методы решения обратных задач математической физики}, Издательство ЛКИ, Москва, 2009.

\bibitem{K} С.И.~Кабанихин, \textit{Обратные и некорректные задачи}, Изд--во СО РАН, Новосибирск, 2018.

\bibitem{PO} A.I.~Prilepko, D.G.~Orlovsky, I.A.~Vasin, \textit{Methods for solving inverse problems in mathematical physics}, Dekker,  Marcel, 2000.

\bibitem{I} V.~Isakov, \textit{Inverse Problems for Partial Differential Equations}, Springer, New York, 2006.

\bibitem{K1} С. И. Кабанихин, \textit{Обратные задачи естествознания}, Ж. вычисл. матем. и матем. физ., \textbf{60}:6 (2020), 935--938. DOI:10.31857/S0044466920060046.

\bibitem{Kab1} Kabanikhin S.I., Bakanov G.B. The optimizational method for solving the discrete inverse problem for hyperbolic equation. \textit{Journal of Inverse and Ill-Posed Problems}, 1996, Vol. 4, No. 6, pp. 513--530.

\bibitem{Kab2} Kabanikhin S.I., Bakanov G.B. The optimization method for solving a discrete inverse problem for a hyperbolic equation. \textit{Doklady Mathematics,} 1998, Т. 58, No 1, pp. 32--35.

\bibitem{Karch} A.L.~Karchevsky. Finite-Difference Coefficient Inverse Problem and Properties of the Misfit Functional. \textit{ Journal of Inverse and Ill-Posed Problems}, 1998, Vol. 6, No 5,  pp. 431--452.

\bibitem{SVV} А.А.~Самарский, П.Н.~Вабищевич, В.И.~Васильев, \textit {Итерационное решение ретроспективной обратной задачи теплопроводности}, Матем. моделирование, \textbf{9}:5 (1997), 119--127.

\bibitem{LE} D.~Lesnic, L.~Elliott, D.~Ingham, \textit {An iterative boundary element method for solving the backward heat conduction problem using an elliptic approximation}, Inverse Problems in Engineering,   \textbf{6}:4 (1998), 255--279. DOI:10.1080/174159798088027679.

\bibitem{KT} А. И. Короткий, И. А. Цепелев, \textit {Решение ретроспективной обратной задачи для одной нелинейной эволюционной модели}, Тр. ИММ УрО РАН, \textbf{9}:2 (2003), 73--86.

\bibitem{M} N.S.~Mera, \textit {The method of fundamental solutions for the backward heat conduction problem}, Inverse Problems in Science and Engineering, \textbf{13}:1 (2005), 65--78. DOI:10.1080/10682760410001710141.

\bibitem{JD} B.T.~Johansson, D.~Lesnic, \textit {A procedure for determining a spacewise dependent heat source and the initial temperature}, Applicable Analysis, \textbf{87}:3 (2008), 265--276. DOI:10.1080/00036810701858193.

\bibitem{IK}  A.~Ismail-Zadeh, A.~Korotkii, G.~Schubert, I.~Tsepelev, \textit{Numerical techniques for solving the inverse retrospective problem of thermal evolution of the Earth interior}, Computers \& Structures, \textbf{87}:11 (2009), 802--811. DOI:10.1016/j.compstruc.2009.01.005.

\bibitem{ZZ} Zhenyu Zhao, Zehong Meng, \textit{A modified Tikhonov regularization method for a backward heat equation}, Inverse Problems in Science and Engineering. \textbf{19}:8 (2011), 1175--1182. DOI:10.1080/17415977.2011.605885.

\bibitem{VK} V.I.~Vasil'ev, A.M.~Kardashevsky, P.V.~Sivtsev, \textit{Computational experiment on the numerical solution of some inverse problems of mathematical physics}, IOP Conf. Series: Mater. Sci. Eng., \textbf{158} (2016), 012093. DOI:10.1088/1757-899X/158/1/012093.

\bibitem{VK1} V.I.~Vasil'ev, A.M.~Kardashevsky, \textit{Iterative Solution of the Retrospective Inverse problem for a Parabolic Equation Using the Conjugate Gradient Method}, LNCS, \textbf{10187} (2017),  698--705. https://doi.org/10.1007/978-3-319-57099-0-80.

\bibitem{TS} В.П.~Танана, А.И.~Сидикова, \textit{Об оценке погрешности приближенного решения, вызванной дискретизацией интегрального уравнения первого рода}, Труды ИММ  УрО РАН, \textbf{22}:1 (2016), 263--270, DOI:10.1134/S0081543817090231.

\bibitem{VK2} В.И.~Васильев, А.М.~Кардашевский, \textit{Численное решение ретроспективной обратной задачи теплопроводности с помощью интеграла Пуассона}, Сиб. журн. индустр. матем., \textbf{21}:3 (2018), 26--36. DOI: https://doi.org/10.17377/sibjim.2018.21.303.

\bibitem{SG} A.~Savitzky and Marcel J.E.~Golay, \textit{Smoothing and Differentiation of Data by Simplified Least Squares Procedures}, Anal. Chem., \textbf{36}:8 (1964), 1627--1639. DOI:10.1021/ac60214a047. 

\bibitem{GHS} Ganga Gowri B., Hariharan V., Thara S. etc., \textit{2D Image data approximation using Savitzky Golay filter — Smoothing and differencing}, (2013), International Mutli-Conference on Automation, Computing, Communication, Control and Compressed Sensing, 365--371. DOI:10.1109/iMac4s.2013.6526438.

\bibitem{GH} https://github.com/espdev/sgolay2/blob/master/sgolay2.py

\end{thebibliography}



%\medskip

\vspace{2mm}{\it Васильев Василий Иванович

Россия, Якутск, Северо-Восточный федеральный университет}

{\tt vasvasil@mail.ru}

\vspace{2mm}
{\it Кардашевский Анатолий Михайлович

Россия, Якутск, Северо-Восточный федеральный университет}

{\tt kardam123@gmail.com}



\end{document}