\documentclass[reqno]{amsart} \usepackage{hyperref} \usepackage{graphicx} \AtBeginDocument{{\noindent\small \emph{Electronic Journal of Differential Equations}, Vol. 2015 (2015), No. 267, pp. 1--14.\newline ISSN: 1072-6691. URL: http://ejde.math.txstate.edu or http://ejde.math.unt.edu \newline ftp ejde.math.txstate.edu} \thanks{\copyright 2015 Texas State University.} \vspace{9mm}} \begin{document} \title[\hfilneg EJDE-2015/267\hfil Stable algorithm] {Stable algorithm for identifying a source in the heat equation} \author[L. Chorfi, L. Alem \hfil EJDE-2015/267\hfilneg] {Lahc\`ene Chorfi, Le\"ila Alem} \address{Lahc\`ene Chorfi \newline Laboratoire de Math\'ematiques Appliqu\'ees, Universit\'e B. M. d'Annaba, B.P. 12, 23000 Annaba, Alg\'erie} \email{l\_chorfi@hotmail.com} \address{Le\"ila Alem \newline Laboratoire de Math\'ematiques Appliqu\'ees, Universit\'e B. M. d'Annaba, B.P. 12, 23000 Annaba, Alg\'erie} \email{alemleila@yahoo.fr} \thanks{Submitted June 27, 2015. Published October 16, 2015.} \subjclass[2010]{35K05, 65M32, 65T50} \keywords{Inverse problem; heat equation; fourier regularization; finite difference} \begin{abstract} We consider an inverse problem for the heat equation $u_{xx}=u_t$ in the quarter plane $\{ x>0, t>0\}$ where one wants to determine the temperature $f(t)=u(0,t)$ from the measured data $g(t)=u(1,t)$. This problem is severely ill-posed and has been studied before. It is well known that the central difference approximation in time has a regularization effect, but the backward difference scheme is not well studied in theory and in practice. In this paper, we revisit this method to provide a stable algorithm. Assuming an a priori bound on $\|f\|_{H^s}$ we derive a H\"older type stability result. We give some numerical examples to show the efficiency of the proposed method. Finally, we compare our method to one based on the central or forward differences. \end{abstract} \maketitle \numberwithin{equation}{section} \newtheorem{theorem}{Theorem}[section] \newtheorem{remark}[theorem]{Remark} \newtheorem{corollary}[theorem]{Corollary} \allowdisplaybreaks \section{Introduction}\label{sec:1} In many engineering applications, we need to determine the temperature on both sides of a thick wall, but one side is inaccessible to measurements \cite{carasso82}. This problem leads to the following parabolic equation in the quarter plane: \begin{equation}\label{P1} \begin{gathered} u_{xx}=u_{t}, \quad x>0, \; t>0,\\ u(1,t)=g(t), \quad t\geq0, \\ u(x,0)=0, \quad x\geq0. \end{gathered} \end{equation} Our purpose is to determine the boundary condition source $f(t)=u(0,t)$ from the temperature $g(t)=u(1,t)$ at the interior point $x=1$. Since the data $g$ is based on (physical) observations, we have a measured data function $g^\delta\in L^2(\mathbb{R})$ which satisfies $$ \|g^\delta-g\|\leq \delta $$ where $\|\cdot\|$ denotes the $L^2$-norm, and the constant $\delta>0$ represents the level noise. The problem of identifying the source $f$ is ill-posed in the sense that the solution (if it exists) does not depend continuously on the data $g$. This can be seen by solving \eqref{P1} in the frequency domain. Let $\hat v$ denote the Fourier transform of function $v(t)\in L^2(\mathbb{R})$ defined by $$ \hat v(\xi):=\frac{1}{\sqrt{2\pi}}\int_{-\infty}^{+\infty}v(t)e^{-i\xi t}dt,\quad \xi\in \mathbb{R}, $$ and $\|\cdot\|_s$ denote the norm in Sobolev space $H^s(\mathbb{R})$ defined by $$ \|v\|_s:=\Big(\int_{-\infty}^{+\infty}(1+\xi^2)^s|\hat v(\xi)|^2d\xi\Big)^{1/2}. $$ When $s=0$, $\|\cdot\|_0:=\|\cdot\|$ denotes the $L^2(\mathbb{R})$ norm. To use Fourier techniques, we extend the functions $u(x,t)$ and $g(t)$ to the whole real t-axis by defining them to be zero for $t < 0$. Problem \eqref{P1} can now be formulated, in the frequency space, as follows: \begin{equation}\label{FP} \begin{gathered} \hat u_{xx}(x,\xi)=i\xi{\hat u}(x,\xi), \quad x>0, \; \xi\in \mathbb{R},\\ \hat u(1,\xi)=\hat g(\xi), \quad \xi\in \mathbb{R}, \\ \hat u(.,\xi), \quad \text{is bounded when } x\to +\infty. \end{gathered} \end{equation} The (formal) solution of problem \eqref{FP} is \begin{equation}\label{Fu} \hat u(x,\xi)=e^{\sqrt{i\xi}(1-x)}\hat g(\xi) \end{equation} where \[ \sqrt{i\xi}=\begin{cases} (1+i)\sqrt{|\xi|/2}, &\xi\geq 0 , \\ (1-i)\sqrt{|\xi|/2}, &\xi< 0. \end{cases} \] and, in particular, for $x=0$, \begin{equation}\label{Ff} \hat f(\xi)=e^{\sqrt{i\xi}}\hat g(\xi). \end{equation} It is clear, from equation \eqref{Ff}, that the transform $\hat g$ must decay faster than the factor $\exp{(-\sqrt{|\xi|/2})}$. This implies that $g$ must belong to Sobolev space $H^s(\mathrm{R})$ for all $s\geq0$. However, in general, the measured data $g^\delta$ does not possess such a decay property. Thus, the numerical simulation is very difficult and a special regularization is required. Some papers have presented mathematical and effective algorithms of these problems \cite{Bern99, Elden95, Chu04, Chu05, Dinho97, Dinho01, Murio90}. In the papers \cite{Elden95, Dinho01, Murio90} the authors investigated finite difference methods but they have only considered the central (or forward ) difference scheme in time. To our knowledge the backward scheme has not been well studied numerically until now. In this paper, we examine this question with the objective of providing a stability result. Our algorithm consists of two steps. In the first one, we solve a well-posed problem in the interval $x\geq 1$ with perturbed data $g^\delta$ such that $\|u(1,.)-g^\delta\|_{L^2}\leq \delta$. In the second step, we solve a Cauchy problem for $x\in [0,1]$ with perturbed Cauchy data $(g^\delta, H_mg^\delta)$, where $H_m$ is a regularized Fourier integral operator \cite{Chu05,Chu10}. We approximate the previous problem by backward finite differences in time. Then, we solve, at each step, an initial value problem for a second order equation in the space variable. In section 2, we give some information on the forward problem. In section 3, the inverse problem is reduced to the Cauchy problem in the interval $[0,1]$. In section 4, we propose a numerical procedure and derive a stability estimation under an a priori bound on $\|f\|_s$ and we show that the regularization parameter $\tau$ (step length in time) can be chosen by a discrepancy principle. Finally, numerical results are given, in section 5, to show the efficiency of the method and compared with other schemes. \section{Forward problem} \label{sec:2} To give some numerical examples, we need the solution of the forward problem in an explicit form. The direct problem is set as follows: given the source $f$, determine $u$ which satisfies the system \begin{equation}\label{PD} \begin{gathered} u_{xx}(x,t) =u_{t}(x,t),\quad x>0, \; t>0, \\ u(x,0)=0,\quad \forall x\geq 0,\\ u(0,t)=f(t), \quad t>0, \quad \lim_{x\to+\infty}u(x,t)=0. \end{gathered} \end{equation} Applying the Fourier-sine transform with respect to the variable $x$, $$ \widehat{u}(\xi,t)=\int_0^\infty u(x,t)\sin(\xi x)dx, \quad \xi \geq 0, $$ and its inverse $$ u(x,t)=\frac{2}{\pi}\int_0^\infty \widehat{u}(\xi,t)\sin(\xi x)d\xi, \quad x \geq 0, $$ yields the following problem in frequency and time, \begin{equation} \label{FSPD} \begin{gathered} \frac{\partial\widehat{u}}{\partial t}(\xi,t)+\xi^{2}\widehat{u}(\xi,t) =\xi f(t), \quad t>0, \\ \widehat{u}(\xi,0)=0. \end{gathered} \end{equation} It is easy to see that the solution of problem \eqref{FSPD} is $$ \hat u(\xi,t)=\int_0^t \xi f(s)\exp [\xi^{2}(s-t)]ds. $$ Then the solution of \eqref{PD} is given by the integral \begin{equation}\label{u_ex} u(x,t)=\int_0^t\frac{x}{t-s} k(x,t-s)f(s)ds \end{equation} where $k(x,t)$ is the heat kernel $$ k(x,t)=\frac{1}{2\sqrt{\pi t}}\exp\big(\frac{-x^2}{4t}\big). $$ For more details see the book \cite{Cop}. The heat flux at $x=1$ is \begin{equation} \label{Flux} h(t):=u_x(1,t)=\int_0^t (1-\frac{1}{2(t-s)})\frac{1}{t-s} k(1,t-s)f(s)\,ds. \end{equation} Problem \eqref{PD} is well-posed in the following sense. \begin{theorem} \label{thm2.1} \begin{enumerate} \item The solution $t\to u(.,t)$ is unique in the space $$ \mathcal{H}=C^0([0,+\infty[,H^2(\mathbb{R}^+))\cap C^1(]0,+\infty[,L^2(\mathbb{R}^+)). $$ \item Assume that $f\in L^2(\mathbb{R})$. Then, for all $s\geq0$, $$ u(x,.)\in C^0([0,+\infty[, L^2(\mathbb{R}))\cap C^\infty(]0,+\infty[, H^s(\mathbb{R})), $$ and we have the stability estimate \begin{equation}\label{stab} \|u(x,.)\|_s\leq C(s,x_0)\|f\|\quad \text{for all } x\geq x_0>0. \end{equation} \end{enumerate} \end{theorem} \begin{proof} (1) Assume that $t\to u(.,t)$ is in the space $\mathcal H$. After multiplying in the PDE with $u$, integrating with respect to $x$, the following identity results $$ \frac{1}{2}\frac{d}{dt}\int_0^\infty|u(x,t)|^2dx +\int_0^\infty|u_x(x,t)|^2dx=-f(t)u_x(0,t). $$ If $f=0$, it follows that the energy $E(t)=\frac{1}{2}\|u(.,t)\|^2$ is a decreasing function. Since $E(0)=0$, $u$ must vanish identically. (2) Using \eqref{Fu} and \eqref{Ff} we get $\hat u(x,\xi)=e^{-x\sqrt{i\xi}}\hat f(\xi)$. Since $|e^{-x\sqrt{i\xi}}|=e^{-x\sqrt{|\xi|/2}} $, we see that $u(x,.)\in H^s(\mathbb{R})$ for all $s\geq0$. If $x\geq x_0>0$, we have the uniform bound $$ \| u(x,.)\|_s\leq \big[\sup_{\xi\geq0}(1+\xi^2)^se^{-x\sqrt{2\xi}}\big]^{1/2} \|\hat f\|\leq C(s,x_0)\|f\| $$ with $C(s,x_0,)=(5s/x_0)^{s}$. From the representation~\eqref{u_ex}, we see that $u(\cdot,t)$ is $C^\infty$ for $x>0$ and is rapidly decreasing as $t\to\infty$. On the other hand $\hat u(\cdot,\xi)$ is $C^\infty$ for $x>0$ and we have, for all $n\in \mathbf{N}$, $$ \widehat{\frac{\partial^n u}{\partial x^n}}(x,\xi) =\frac{\partial^n}{\partial x^n}\hat u(x,\xi) =(-\sqrt{i\xi})^ne^{-x\sqrt{i\xi}}\hat f(\xi). $$ Using the rapid decay of the factor $e^{-x\sqrt{|\xi|/2}}$, this proves that $\frac{\partial^nu}{\partial x^n}(x,\cdot)\in H^s(\mathbb{R})$.\\ It remains to show that the mapping $x\to u(x,\cdot)$ is continuous at $x=0$, i.e., $\lim_{x\to0} \|(1-e^{-x\sqrt{i\xi}})\hat f\|=0$. Using the inequalities \begin{gather*} | 1- e^{-x\sqrt{i\xi}}|\leq 2 \quad \text{for $x\geq 0$ and $\xi\in\mathbb{R}$},\\ |1- e^{-x\sqrt{i\xi}}|\leq x\sqrt A e^{x\sqrt A}\quad \text{for $x\geq 0$ and $|\xi|\leq A$}, \end{gather*} it follows that for all $x\in[0,1]$ and all $A\geq 1$, \[ \int_\mathbb{R}(1- e^{-x\sqrt{i\xi}})^2 \|\hat f\|^2d\xi\leq Ax^2e^{2x\sqrt A} \|\hat f\|^2+4\int_{|\xi|\geq A}|\hat f|^2d\xi. \] For any $\epsilon>0$, if we choose $Ax\leq 1$ and $A$ large enough, we can make the right hand side less than $\epsilon$. Which ends the proof. \end{proof} We remark that our inverse problem is equivalent to the following integral equation of Volterra type, \begin{equation}\label{IE} g(t)=\int_0^t\frac{k(1,t-s)}{t-s}f(s)ds, \end{equation} with a $C^\infty$-kernel, which proves again that problem \eqref{P1} is severely ill-posed. This equation is regularized by Tikhonov method in the paper \cite{Bern01}. \begin{remark}\label{causa} \rm From the integral \eqref{u_ex} we see that the direct problem $f\mapsto g=u(1,.)$ satisfies the causality principle. That is if we change $f(t)$ for $t\geq t_*$ then the solution $g(t)$ can only change for $t\geq t_*$ as well. On the other hand the integral equation \eqref{IE} shows that for inverse problem the numerical solution $f(t_*)$ depends on $g(t)$ for $t\in[0,t_*]$. Numerically this may creates more noise amplifications as has been noted by Carasso in \cite{carasso92}. To reduce this phenomenon, we must observe $g(t)=u(1,t)$ for enough long period. \end{remark} \section{Inverse problem} \label{sec:3} \subsection*{Cauchy problem} Consider the well-posed problem \begin{equation}\label{P2} \begin{gathered} u_{xx}=u_{t}, \quad x>1, \; t>0,\\ u(1,t)=g(t),\quad t>0, \\ u(x,0)=0, \quad x\geq1. \end{gathered} \end{equation} If $g\in H^s(\mathbb{R})$, $s\geq 1/2$, the actual solution is \begin{equation}\label{u_ext} u(x,t)=\int_{-\infty}^\infty \widehat{g}(\xi)\exp[\sqrt{i\xi}(1-x)]\exp(i\xi t)d\xi\quad\text{for } x\geq 1. \end{equation} Since $|\exp\left[\sqrt{i\xi}(1-x)\right]|=e^{-(x-1)\sqrt{|\xi|/2}} $ , the integral \eqref{u_ext} is convergent for $x>1$ and $u(x,.)\in H^\sigma(\mathbb{R})$ for all $\sigma>0$. The heat flux $h(t)=u_x(1,t)$ is in $H^{s-\frac{1}{2}}(\mathbb{R})$ and is given by: \begin{equation}\label{flux} h(t):=Hg(t)=-\int_{-\infty}^{+\infty}\sqrt{i\xi}\widehat{g}(\xi)\exp(i\xi t)d\xi. \end{equation} We see that \eqref{P1} is equivalent to the Cauchy problem \begin{equation}\label{cauchy} \begin{gathered} u_{xx}(x,t) =u_{t}(x,t),\quad 00, \\ u(1,t)=g(t),\quad u_{x}(1,t)=h(t)\; t>0,\\ u(x,0)=0 ,\quad 0\xi_m}|\xi||\hat g(\xi)|^2d\xi\hfill\\ &\leq \max_{\xi>\xi_m}\frac{\xi}{(1+\xi^2)^s}\| g\|_s^2 \leq M^2(\xi_m)^{1-2s}. \end{aligned} \end{equation} On the other hand \begin{equation}\label{eq3} \|h_m-h_{m,\delta}\|\leq \sqrt\xi_m \|g-g_{\delta}\|\leq \delta\sqrt \xi_m . \end{equation} If we choose $\xi_{m}= (M/\delta)^{1/s}$, then \begin{equation}\label{eq4} \|h-h_{m,\delta}\|\leq \|h-h_{m}\|+\|h_m-h_{m,\delta}\| \leq 2 M^{1/(2s)} {\delta}^{1-\frac{1}{2s}}. \end{equation} \end{proof} \begin{corollary}\label{cor1} Assume that $\hat g(\xi)=e^{-x\sqrt{i\xi}}\hat f(\xi)$ with $\|f\|\leq E$, and $g^\delta\in L^2(\mathbb{R})$ satisfying $\|g-g^\delta\|\leq \delta$. If we select $\xi_{m}=5s (E/\delta)^{1/s}$ then we get the error bound \begin{equation}\label{estim2} \|h-h_{m,\delta}\|\leq 2\sqrt {5s}{E}^{\frac{1}{2s}}{\delta}^{1-\frac{1}{2s}}, \quad \forall s\geq\frac{1}{2}. \end{equation} \end{corollary} \section{Discretization of the Cauchy problem and stability}\label{sec:4} The solution of the Cauchy problem \eqref{cauchy} is unique but is not stable with respect of the data $(g,h)$. To stabilize the problem, we propose a scheme in two steps. \subsection*{Time discretization} The problem for $0\leq t\leq T$ can be discretized by replacing the time derivative $u_{t}$ by the backward difference with step length $\tau$. Indeed, let $t_{n}=n\tau$, $n=0$ to $N$, with $\tau=\frac{T}{N}$ is the time step, then we have the approximation, for $n\geq1$: $$ u_{t}(x,t_{n})\approx\frac{u(x,t_{n})-u(x,t_{n-1})}{\tau}. $$ Furthermore if we assume that $|u(x,t)|,|u_{t}(x,t)|,|u_{tt}(x,t)| \leq M $, for all $(x,t) \in [0,1]\times[0,T]$, then $$ u_{t}(x,t_{n})=\frac{u(x,t_n)-u(x,t_{n-1})}{\tau}+\psi(x,t_{n}) \quad \text{with } \psi(x,t)=O(\tau). $$ Noticing $w_{n}(x)=u(x,t_{n})$ and $\psi_n(x)=\psi(x,t_n)$, the equation $u_{xx}=u_{t}$ becomes an ordinary differential equation in the space variable $x$: $$ w_{n}''-\theta^{2}w_{n}=-\theta^{2}w_{n-1}+\psi_{n}(x), $$ with $\theta^2=1/\tau$. Thus we consider the semi-discrete problem \begin{equation}\label{Cn} \begin{gathered} v_{n}''(x)-\theta^{2}v_{n}(x)=-\theta^{2}v_{n-1}(x),\quad \text{for } 0\leq x\leq 1 ,\; (v_{0}=0) \\ v_{n}(1)=g_{n},\quad v_{n}'(1)=h_{n}\quad (g_{n}=g(t_{n}), h_{n}=h(t_{n})) . \end{gathered} \end{equation} The solution $v_n$ has the representation \begin{equation} v_{n}(x)= g_{n}\cosh (\theta (1-x))-\frac{h_{n}}{\theta }% \sinh (\theta (1-x)) + \theta \int_{x}^{1}\sinh \theta(x-y)v_{n-1}(y)dy. \end{equation} Starting with $v_0=0$, we obtain recursively the expression \begin{equation}\label{semi-disc} v_{n}(x)=\varphi_n+S\varphi_{n-1}+S^2\varphi_{n-2}+\dots +S^{n-1}\varphi_{1}, \end{equation} with \begin{gather*} \varphi_n(x)=g_{n}\cosh (\theta (1-x))-\frac{h_{n}}{\theta } \sinh (\theta (1-x)), \\ S\varphi(x)=\theta\int_{x}^{1}\sinh \theta(x-y)\varphi(y)dy. \end{gather*} \subsection*{Space discretization} Now we discrete the interval $[0,1]$ by the sequence $x_j=jk$, $j=1,L$ ($k=1/L$) and we approximate the integral operator $S$ by quadrature by considering the matrix $S_k$: $$ S_k(i,j)=\begin{cases} k\theta\sinh(\theta(j-i)k)& \text{if } j> i,\\ 0& \text{if } j\leq i. \end{cases} $$ For a function $\Psi(x,t)$ defined on the grid ${\mathcal G}=\{x_m = mk; t_n =n\tau : m = 0,L; \, n=0,N\}$, we introduce the vector $$ \Psi_n=(\Psi(x_1,t_n), \Psi(x_2,t_n), \dots , \Psi(1,t_n)). $$ Now we define the discrete solution $\bar u(x,t)$ on the grid $\mathcal G$ by \begin{equation}\label{ubar} \bar u_n:=\varphi_n+S_k\varphi_{n-1}+S_k^2\varphi_{n-2}+\dots +S_k^{n-1}\varphi_{1} \end{equation} with $\varphi_n=(\varphi_n(x_1),\varphi_n(x_2),\dots, \varphi_n(1))$. For simplicity, we assume in the following that $T=1$. \begin{theorem}\label{thm1} Let $w_n(x)=u(x,t_n)$ be the exact solution of the Cauchy problem \eqref{cauchy} at $t=t_n$ and $v_n(x)$ be the semi-discrete solution given by \eqref{semi-disc} . Then for all $n\leq N$ and all $x\in]0,1]$, \begin{equation}\label{estim3} |w_n(x)-v_n(x)|\leq M \tau^2 2^{-n}\exp\big(\frac{n(1-x)}{\tau^{\frac{1}{2}}}\big). \end{equation} \end{theorem} \begin{proof} Let $z_n=w_n-v_n$. Then \begin{gather*} z_{n}''-\theta^{2}z_{n}=-\theta^{2}z_{n-1}+\psi_n(x),\quad (z_{0}=0) \\ z_{n}(1)=0,\quad z_{n}'(1)=0 . \end{gather*} Hence \begin{equation} z_{n}(x)={\theta } \int_x^1\sinh \theta (x-y)z_{n-1}(y)\,dy -\frac{1}{\theta } \int_{x}^{1}\sinh \theta(x-y)\psi_{n}(y)\,dy. \end{equation} It follows that \begin{equation} z_{n}(x)=-\tau[S\psi_n+S^2\psi_{n-1}+\dots +S^n\psi_{1}]. \end{equation} Since $|\psi_n(x)|\leq M\tau/2$, we have \begin{equation} |z_{n}(x)|\leq \frac{M}{2}\tau^2\sum_{i=1}^n\|S\|^i. \end{equation} But $\|S\|\leq \frac{1}{2}\exp\theta(1-x)$, which leads to \begin{equation} \begin{aligned} |z_{n}(x)| &\leq \frac{M}{2}\tau^2 [2^{-1}\exp\theta(1-x)+2^{-2}\exp 2\theta(1-x)+\dots +2^{-n}\exp n\theta(1-x)]\\ &\leq M\tau^2 2^{-n} \exp\big(\frac{n(1-x)}{\tau^{1/2}}\big). \end{aligned} \end{equation} \end{proof} \begin{remark} \label{rmk4.2} \rm Since $n\leq N=1/\tau$, the right hand of \eqref{estim3} behaves like $r(x,\tau)=\tau^2\exp\big(-\frac{\log 2}{\tau}+\frac{1-x}{\tau\sqrt\tau}\big)$. In particular, if $\exp (\frac{1}{\tau\sqrt\tau})=\frac{1}{\epsilon} \Leftrightarrow \tau=1/(\log\frac{1}{\epsilon})^{2/3}$, then $r(x,\tau)\leq \frac{1}{\epsilon^{1-x}}(\log\frac{1}{\epsilon} )^{-4/3}$ does not approach zero as $\epsilon\to0$. This means that the difference scheme is (perhaps) inconsistent. But, in the numerical test (see Figure \ref{fig1}) the algorithm converges at least when the Cauchy data are exact. It seems that the error bound \eqref{estim3} is sharp. \end{remark} \begin{remark} \rm (1) Elden \cite[Corollary 3.2]{Elden95} investigated central difference discretization in time. He approximates the heat equation by the central difference equation: $$ v_{xx}(x,t)=\frac{1}{2\tau}(v(x,t+\tau)-v(x,t-\tau)), $$ and he proved an error estimate between $u(x,t)$ and $v(x,t)$ subjected to the conditions $u(1,t)=g(t)$ and $v(1,t)=g^\delta$. More precisely, he gets asymptotically, as $\delta\to0$, $$ \|u(x,.)-v(x,.)\|\sim \frac{M}{(\log(\frac{M}{\delta}))^2} $$ with $\tau=1/2(\log(\frac{M}{\delta})^2$. (2) To prove a similar estimate for the backward difference equation $$ v_{xx}(x,t)=\frac{1}{\tau}(v(x,t)-v(x,t-\tau)), $$ is an open problem. Indeed, one of the conclusions of Elden in his paper is \it{It is important that the time difference has a substantial forward component}. \end{remark} The stability of the discrete scheme \eqref{ubar} is proved in the following theorem. \begin{theorem}\label{thm3} Let $\bar u_{n}$ be the solution defined by \eqref{ubar} associated with $(g_n,h_n)$. Then we have \begin{equation} \label{estim4} \|\bar u_{n}\|\leq 2{\tau}^{1/4}(\sqrt k\theta)^n\exp(n(1-k)\theta)\left(|g_{n}|+\sqrt\tau|h_{n}|\right),\quad \forall n\leq N, \end{equation} where $\|\cdot\|$ is the discrete $L^2$-norm, i.e. $\|\bar u_{n}\|^2=k\sum_{m=1}^L (\bar u_{n}^m)^2$. \end{theorem} \begin{proof} From \eqref{ubar} we have \begin{equation} \|\bar u_{n}\|\leq \sum_{i=0}^{n-1}\|S_k\|^i\|\varphi_{n-i}\|. \end{equation} Moreover, \begin{align*} |\varphi_{n}^m| &\leq |g_n|\cosh(\theta(1-mk))+\frac{|h_n|}{\theta}\sinh(\theta(1-mk))\\ &\leq \sqrt {2}\exp(\theta(1-mk))( |g_n|+\sqrt \tau|h_n| ), \end{align*} then $$ \|\varphi_{n}\|\leq \frac{1}{\sqrt\theta}\exp(\theta(1-k)) ( |g_n|+\sqrt \tau|h_n|). $$ On the other hand \begin{align*} \|S_{k}\|^2 &\leq k\theta^2\sum_{i=1}^L\sum_{j=i}^L\sinh^2[\theta k(j-i)] \leq\frac{1}{4}k\theta^2\sum_{i=1}^L\sum_{j=i}^L\exp[2\theta k(j-i)]\\ &\leq \frac{1}{4} k\theta^2\sum_{i=1}^L\sum_{j=0}^{L-i}\exp[2j\theta k] \\ &\leq \frac{1}{2}k\theta^2\sum_{i=1}^L\exp(2k\theta(L-i))\leq k\theta^2\exp(2\theta(1-k)); \end{align*} that is, $$ \|S_{k}\|\leq \sqrt k\theta\exp(\theta(1-k)), $$ which leads to estimate \eqref{estim4}. \end{proof} As a consequence we prove the following the main result. \begin{theorem}\label{thm4} Assume that $\|f\|\leq E$. Let $\bar u$ the discrete solution associated with $(g,h=Hg)$ and $\bar u^\delta$ the discrete solution associated with the perturbed data $(g^\delta,h_{m,\delta}=H_mg^\delta)$, where $H$ (resp. $H_m$) is the operator given by \eqref{flux} (resp. \eqref{hm}). If we select $ \xi_{m}={5s}(\frac{E}{\delta})^{1/s}$, $\tau=(2s/\log \frac{1}{\delta})^{2/3}$ and $k\leq \tau$, then we have \begin{equation} \label{estim5} \|\bar u-\bar u^\delta\|\leq 8\sqrt{10}s E^\frac{1}{2s}\big(\log \frac{1}{\delta}\big)^{-1/3}{\delta}^{1-\frac{1}{s}}, \quad \forall s\geq 1, \end{equation} where $\|\cdot\|$ is the discrete $L^2$-norm. \end{theorem} \begin{proof} Since $\sqrt k\theta\leq 1$ and $n\leq N=1/\tau$, it follows from \eqref{estim4} that \begin{equation} \label{estim6} \|\bar u-\bar u^\delta\|\leq 2\sqrt2 \exp\big(\frac{1}{\tau\sqrt\tau}\big)[\|g-g^\delta\|+\sqrt\tau\|h-h_{m,\delta} \|]. \end{equation} With \eqref{estim2}, it follows that \begin{equation} \label{estim7} \|\bar u-\bar u^\delta\|\leq 2\sqrt2 \exp\big(\frac{1}{\tau\sqrt\tau}\big)[\delta+ 2\sqrt{5s}{E}^\frac{1}{2s}{\delta}^{1-\frac{1}{2s}}\sqrt\tau]. \end{equation} If we choose $\tau=(2s/\log \frac{1}{\delta})^{2/3}$, then $\exp\big(\frac{1}{\tau\sqrt\tau}\big)=(\frac{1}{\delta})^{\frac{1}{2s}}$ and we obtain \begin{equation} \|\bar u-\bar u^\delta\|\leq 8\sqrt{10}s E^\frac{1}{2s}\delta^{1-\frac{1}{s}}\big(\log \frac{1}{\delta}\big)^{-1/3}, \end{equation} which establishes the estimate \eqref{estim5}. \end{proof} \begin{remark} \rm Assume the a priori bound $\|f\|_s\leq E$ with $s\geq 1$. Since $ \hat g(\xi)=e^{-\sqrt{i\xi}}\hat f(\xi)$, then $\|g\|_s\leq E$. As a consequence, we can establish, as in theorem \ref{thm4}, the estimate \begin{equation} \label{estim8} \|\bar u-\bar u^\delta\|\leq 8\sqrt{2}s^\frac{1}{3} E^\frac{1}{2s}\delta^{1-\frac{1}{s}}\big(\log \frac{1}{\delta}\big)^{-1/3}. \end{equation} \end{remark} \section{Algorithms and numerical examples} \label{sec:5} \subsection*{Numerical implementation of the Fourier integral} To use the Fast Fourier Transform (FFT) it is necessary to assume periodicity. Therefore we extend $g$ to the interval $[T,2T]$ (as in \cite{EBR2000}). The Fourier integral $$ h_{N}(t)=-\int_{-\xi_N}^{+\xi_N}\sqrt{i\xi}\widehat{g}(\xi)\exp(i\xi t) d\xi,\quad \xi_N=\frac{\pi N}{T}, $$ is approximated by Discrete Fourier Transform (DFT). Let $ g=\{ g_k\}_{k=1}^{2N}$ be a discrete vector and $\hat g=\{ \hat g_k\}_{k=1}^{2N}$ its DFT, then the vector $h_N(t)$ is given by the trigonometric interpolation polynomial of the form $$ h_N(t)=-\frac{1}{2T}\sum_{k=-N+1}^{N}\sqrt {-i\xi_k} {\hat g}_{N+k}e^{i\xi_kt}; \quad \xi_k=\frac{k\pi}{T}. $$ \subsection*{Algorithms} The proposed algorithm is described as follows: \smallskip \noindent\textbf{Method 1} (Backward) \begin{enumerate} \item Given an exact source $f$, we provide a data $(g=Af, h=Bf)$ by solving the forward problem (see section \ref{sec:2} where $A$ and $B$ are the integral operators \eqref{u_ex} and \eqref{Flux} respectively). We extend $g$ to $[0,2T]$ such that the extension $\tilde g$ vanishes at the end point $t=2T$ and we set $g_k=\tilde g(t_k)$, $k=1,\dots, 2N$. \item We perturb the data $g^\delta=g+\delta\sigma$, where $\sigma(t)$ is the Gaussian random function and $\delta$ is the level noise. \item Calculate $\hat g=F(g^\delta)$ and $h_m=F^*(\sqrt \Lambda\hat g)$, with $\Lambda=(i\xi_k)$,where $F$ is the Fast Fourier Transform (FFT) and $F^*$ its inverse. \item Calculate $\bar u$ by the procedure:\\ $u(1,:)=zeros(1,N)$;\\ for $i=2:N+1$;\\ $w=D*u(i-1,:)'$; \% $D(i,j)=\sinh(r*(i-j)*k)$\\ \ \ \ for $j=1:M+1$;\\ \ \ \ $u(i,j)=g(i)*\cosh(r*(1-(j-1)*k))-$\\ \ \ \ \ \ \ $h(i)*\sinh(r*(1-(j-1)*k))/r+ r*k*w(j)$;\ \ \ \ \ \ \% with ($r=\sqrt\tau$)\\ \ \ \ end;\\ end; \item The approximate solution is $f_{ap}=(u(:,1))$. \end{enumerate} To compare our algorithm (method 1), we recall the forward and central time difference schemes. Letting $w:=u_x$, the Cauchy problem \eqref{cauchy} can be rewritten as \begin{equation}\label{cauchy_2} \begin{gathered} u_x(x,t)=w(x,t),\quad w_{x}(x,t) =u_{t}(x,t),\quad 00, \\ u(1,t)=g(t),\quad w(1,t)=h(t),\quad t>0,\\ u(x,0)=0 ,\quad 0