Multigrid Method for Solving Inverse Problems for Heat Equation

Abstract

In this paper, the inverse problems for the boundary value and initial value in a heat equation are posed and solved. It is well known that those problems are ill posed. The problems are reformulated as integral equations of the first kind by using the separation-of-variables method. The discretization of the integral equation allowed us to reduce the integral equation to a system of linear algebraic equations or a linear operator equation of the first kind on Hilbert spaces. The Landweber-type iterative method was used in order to find an approximation solution. The V-cycle multigrid method is used to obtain more frequent and fast convergence for iteration. The numerical computation examples are presented to verify the accuracy and fast computing of the approximation solution.

1. Introduction

In the early twentieth century, Hadamard [1] labeled some well-posed problems, stating that a problem is well posed when it fulfills the following points:

• Solution exists.
• Uniqueness—this solution is unique.
• Stability (the given data are continuously dependent on the solution).

If at least one of the above points or conditions is not fulfilled, the problem is considered an ill-posed problem. The violations of 1 and 2 can often be improved with a small re-formulation of the problem. Violations of stability are much harder to remedy because they imply that a small disturbance in the data leads to a large disturbance in the estimated solution [2,3,4,5].

The inverse problem under the study of the heat equation can be solved by many methods. For example, the method of regularization Tikhonov A.N. [6], by the method of Lavrentiev M.M. [7], by the method of quasi-solutions Ivanov V.K. [8] and many others. Various methods and algorithms for solving IP “inverse problems” were explained and used in [9,10,11,12,13,14,15,16,17,18,19]. The success of these methods and algorithms is largely based on understanding and analyzing the mathematical problems related to the declarations of the properties in these IP “inverse problems” and identifying specific difficulties in solving them [20,21,22,23,24,25].

One of the many known methods of solving a linear operator equation is the iteration method, such us the Landweber iteration method. This method has several issues, such as initial guess, complexity, and convergence. The discretization process for the integral equation primes to a huge, sparse, and highly structured system of linear equations or linear operator equations. Inappropriately, theoretic convergence studies and mathematical experiments show that the convergence rate is slow for these iterative algorithms, leading to an amplified cost of iteration [26]. The direct iteration method suffers from some restricting boundaries. Multigrid methods advanced from efforts to overcome these limitations. Multigrid settings are largely successful when used with relaxation or iteration methods and they lead to fast and direct point-to-estimated solutions to solve the IP “inverse problems”.

Multigrid methods [27,28,29,30] are often used to accelerate the convergence rate of iterative or relaxing methods. This is an effective method for solving large systems of linear algebra equations, which results from the discretization of integral equations or PDE. Multigrid methods have been applied to least-squares wave front reconstruction [31]. This method has also been used to solve inverse problems with linear parabolic PDE constraints [32] and ODT problems [29].

In this work, we provide a new algorithm for numerical solutions to the integral equations of the first kind by using the V-cycle multigrid method and comparing it with the Landweber iteration method. In this article, we mixed between the Landweber iteration method and the multigrid method to obtain accelerated solutions

We consider the following two forms of the integral equation of the first kind that come from two inverse ill-posed heat equation problems.

$$Au(t)={\displaystyle \underset{\mathrm{a}}{\overset{b}{\int }}K(t,x)u\left(t\right)dx=f(t)}$$

(1)

where $t\in [c,d]$ and the kernel $K(x,t),{K^{\prime }}_t(x,t)\in C([a,b]\times [c,d]),$ and $f(t)\in L^2[c,d],$ the (1) can represent the initial value problem for heat equation.

$$Au(t)={\displaystyle \underset{0}{\overset{t}{\int }}K(t,\tau )u\left(\tau \right)d\tau =f(t)}$$

(2)

where $t\in \left[0,T\right]$ and $K(t,\tau )$ corresponding the kernel function for the integral Equation (2), $K(t,\tau ),{K^{\prime }}_t(t,\tau ),{K^{″}}_t(t,\tau )\in C([a,b]\times [0,T]),$ $u(t)\in L^2[0,T]$. This type of integral equation can represent the boundary value problem for the heat equation.

Problems (1) and (2) can be reduced to the form of the leaner operator equation after applying the discretization algorithm.

$$Au=f$$

(3)

The solution vector $u\in R^n$ should reduce the following minimizing problem.

$$\underset{u\in R^n}{\min }||Au-f|{|}_2$$

(4)

where $A\in R^{n\times n}$ is a large matrix and the vector $f\in R^n$ represents known error-contaminated data and can be expressed as the following.

$$f=f_{true}+\eta$$

(5)

where $f_{true}$ is an unknown error-free vector associated with $f$ and $\eta$ represents the error in $f$. The noise may stem from measurement and/or discretization errors.

$$||\eta ||\le \delta$$

(6)

where $\delta$ is given and known as error level.

The approaches proposed in this work are the Landweber-type iterative method to obtain a good approximation solution for (3). Then, the multigrid method was used to obtain a more complete and fast solution by using the Landweber type as relaxation iterative in post-smoothing and pre-smoothing parts.

All these steps are organized through the sections in this paper. Section 2 defines the linear partial differential equations and describes the solution as an integral equation of the first kind for the boundary value problem and the initial value problem of the heat equation. Section 3 describes the general notations and the iterative Landweber-type method. Section 4 describes and implements the V-cycle multigrid method. After that, in Section 5, the numerical examples are presented to explain our analysis. Finally, the explanation of the suggested method is summarized in Section 6.

2. Problem Statement

In this section, we briefly present two-class problems: the boundary value and initial value problems. Both problems can be reduced to an integral equation by using the separation-of-variables method. Then, by implementing the discretization process, we obtain the linear operator equation. The inverse ill-posed problems have been studied and solved in many works with a different method [11,12,13,14].

2.1. Inverse Boundary Value Problem

We consider the inverse boundary value problem given in [33,34] with time interval $\left[0,T\right]$.

$$\frac{\partial u\left(x,t\right)}{\partial t}=\frac{{\partial }^{2}u\left(x,t\right)}{\partial x^2};\hspace{1em}x\in (0,x),t\in (0,T]$$

(7)

$$u\left(x,0\right)=0; x\in [0,1]$$

(8)

$$\frac{\partial u\left(0,t\right)}{\partial x}=0; t\in (0,T]$$

(9)

$$u\left(1,t\right)=u\left(t\right); t\in (0,T]$$

(10)

suppose that the $u\left(t\right)\in H^4\left[0,T\right],$ is a function such that

$$u\left(0\right)=u^{\prime }\left(0\right)=u^{″}\left(0\right)=u^{″}\left(T\right)=u^{″}\left(T\right)=0$$

(11)

where $\underset{0}{\overset{T}{{\displaystyle \int }}}{|u^{‴}\left(t\right)|}^{2}dt\le r^2,$ $r-$known number. By using the separation-of-variables method, we obtained the following solution.

$$u\left(x_0,t\right)={\displaystyle \sum }_{n=0}^{\infty }{C}_n\left(t\right)\cos ((n+0.5)\pi x_0)+u\left(t\right)$$

(12)

where the $x_0\in \left(0,1\right),t\in \left[0,T\right]$

$$C_n\left(t\right)=\frac{2e^{-{\left(n+0.5\right)}^2{\pi }^{2}t}}{\left(n+0.5\right)\pi }\underset{0}{\overset{t}{{\displaystyle \int }}}{u}^{\prime }\left(\tau \right)e^{{\left(n+0.5\right)}^2{\pi }^2\tau }d\tau$$

(13)

integrating by parts the right-hand side of Formula (13) twice, we obtain

$$\begin{array}{c}{C}_n\left(t\right)=\frac{2}{\left(n+0.5\right)\pi }\underset{0}{\overset{t}{{\displaystyle \int }}}{u}^{‴}\left(\tau \right)e^{-{\left(n+0.5\right)}^2{\pi }^2\left(t-\tau \right)}d\tau +\frac{2}{{\left(n+0.5\right)}^3{\pi }^3}{u}^{\prime }\left(t\right)\\ -\left[\frac{2}{{\left(n+0.5\right)}^3{\pi }^3}{u}^{\prime }\left(0\right)+\frac{2}{{\left(n+0.5\right)}^5{\pi }^5}{u}^{″}\left(t\right)+\frac{2}{{\left(n+0.5\right)}^5{\pi }^5}{u}^{″}\left(0\right)\right]\end{array}$$

(14)

$$\begin{array}{c}{C}_n\left(t\right)=\frac{2}{{\left(n+0.5\right)}^5{\pi }^5}\underset{0}{\overset{t}{{\displaystyle \int }}}{u}^{‴}\left(\tau \right)e^{-{\left(n+0.5\right)}^2{\pi }^2\left(t-\tau \right)}d\tau \\ -\left[\frac{2}{{\left(n+0.5\right)}^3{\pi }^3}{u}^{\prime }\left(0\right)+\frac{2}{{\left(n+0.5\right)}^5{\pi }^5}{u}^{″}\left(0\right)\right]\end{array}$$

(15)

where $u^{\prime }\left(0\right)=u^{″}\left(0\right)=0$ by (11)

$$C_n\left(t\right)=\frac{2}{{\left(n+0.5\right)}^5{\pi }^5}\underset{0}{\overset{t}{{\displaystyle \int }}}{u}^{‴}\left(\tau \right)e^{-{\left(n+0.5\right)}^2{\pi }^2\left(t-\tau \right)}d\tau$$

(16)

From (12) and (16) we have the following integral equation of the first kind

$$f\left(t\right)=u\left(x_0,t\right)=\underset{0}{\overset{t}{{\displaystyle \int }}}\left[u\left(t\right)+{\displaystyle \sum }_{n=0}^{\infty }\frac{2\cos \left(n+0.5\right)\pi x_0}{{\left(n+0.5\right)}^5{\pi }^5}{e}^{-{\left(n+0.5\right)}^2{\pi }^2\left(t-\tau \right)}\right]g(\tau )d\tau$$

(17)

where $u^{‴}\left(\tau \right)=g(\tau )$. We need to define an operator $B:L^2[a,b]\to L^2[a,b]$ by the following formula.

$$Bg\left(\tau \right)=u\left(t\right)=\underset{0}{\overset{t}{{\displaystyle \int }}}\frac{{\left(t-\tau \right)}^2}{2}g\left(\tau \right)d\tau$$

(18)

where $Bg\left(\tau \right)\in L^2\left[0,T\right].$

We used the algorithm in [6] with a specified derivative term for the kernel and obtained the following.

$$Au(t)=\underset{0}{\overset{t}{{\displaystyle \int }}}K\left(t,\tau \right)g\left(\tau \right)d\tau =f\left(t\right)$$

(19)

where $t\in \left[0,T\right]$.

Where $A$ is finite-dimensional operator. For the next step, we need to apply the discretization algorithm [5,6] on the integral Equation (19)

$${\overline{K}}_i\left(t\right)=K\left(t,{\tau }_i\right)$$

(20)

$${\overline{K}}_i\left(t\right);{\tau }_i\le \tau \le {\tau }_{i+1}, t\in \left[0,T\right],i=0,1,\dots ,n-1$$

(21)

$${\overline{K}}_i\left(t_j\right);{\tau }_i\le \tau \le {\tau }_{i+1},t_j\le t\le t_{j+1},i=0,1,\dots ,n-1,j=0,1,\dots ,n-1$$

(22)

$$K_n\left(t,\tau \right)={\overline{K}}_i\left(t_j\right)$$

(23)

$$A\left[u\left(t\right)\right]=\underset{0}{\overset{t}{{\displaystyle \int }}}{K}_n\left(t,\tau \right)g\left(\tau \right)d\tau =f\left(t\right),t\in \left[0,T\right]$$

(24)

the value of kernel depended on the following condition.

$${\overline{K}}_i\left(t_j\right)=\{\begin{array}{c}K\left(t_j,{\tau }_i\right) i\le j\\ 0? i>j\end{array}$$

(25)

From above, we convert the boundary value problem to linear operator equation $Au=f,$ where $u=Bg,$ by reducing the integral equation to a system of linear algebra equations.

$$AB\left[\begin{array}{c}g\left({\tau }_0\right)\\ g\left({\tau }_1\right)\\ ⋮\\ g\left({\tau }_{n-1}\right)\end{array}\right]=\left[\begin{array}{c}f\left(t_0\right)\\ f\left(t_2\right)\\ ⋮\\ f\left(t_{n-1}\right)\end{array}\right]$$

(26)

where

$$A=\frac{1}{n}\left[\begin{array}{cc}\begin{array}{cc}K\left({\tau }_0,t_0\right)& 0\\ K\left({\tau }_0,t_1\right)& K\left({\tau }_1,t_1\right)\end{array}& \begin{array}{cc}\dots & 0\\ \dots & 0\end{array}\\ \begin{array}{cc}⋮& ⋮\\ K\left({\tau }_0,t_{n-1}\right)& K\left({\tau }_1,t_{n-1}\right)\end{array}& \begin{array}{cc}\dots & ⋮\\ \dots & K\left({\tau }_{n-1},t_{n-1}\right)\end{array}\end{array}\right]$$

(27)

The operator $A$ in (27) is a triangular operator.

2.2. Initial Value Problem

The second inverse problem is the initial value problem for the heat equation has been described by the following linear partial differential equation.

$$\frac{\partial u\left(x,t\right)}{\partial t}=\frac{{\partial }^{2}u\left(x,t\right)}{\partial x^2} ,\hspace{1em}x\in [0,l],t\in \left(0,T\right]$$

(28)

$$u\left(0,t\right)=0,t\in \left(0,T\right]$$

(29)

$$u\left(l,t\right)=0,t\in \left(0,T\right]$$

(30)

$$u\left(x,0\right)=u\left(x\right),x\in [0,l]$$

(31)

where $u\left(0,t\right)$ and $u\left(l,t\right)$ are boundary conditions, $u_0\left(x\right)$ is an initial condition, which needs to be found. This problem was solved in [35] by using Tikhonov’s regularization inversion method and it was solved by Picard’s method in [2]. In this work, we apply another numerical algorithm to obtain a more accurate solution and fast implementation with a high-scale problem. The integral form for this equation will be

$$Au\left(x\right)=\underset{0}{\overset{l}{{\displaystyle \int }}}K\left(x,y\right)u\left(y\right)dy=f\left(x\right)$$

(32)

where the kernel $K\left(x,y\right)$ is an infinite series and we cannot handle the infinite sum, so we need to make the sum of series finite to 10-times

$$K\left(x,y\right)=\frac{2}{l}{\displaystyle \sum }_{n=1}^{10}{e}^{\frac{-{\left(n\pi \right)}^{2}T}{l^2}}\sin (\frac{n\pi x}{l})\sin (\frac{n\pi y}{l}), T>0$$

(33)

For giving the approximate solution to $u\left(x\right)$, we can rewrite the problem as linear algebraic equations where the form will be $Au=f$.

$$A\left[\begin{array}{c}u\left(y_0\right)\\ u\left(y_2\right)\\ ⋮\\ u\left(y_{n-1}\right)\end{array}\right]=\left[\begin{array}{c}f\left(x_0\right)\\ f\left(x_2\right)\\ ⋮\\ f\left(x_{n-1}\right)\end{array}\right]$$

(34)

where

$$A=\frac{l-0}{n}[\begin{array}{cccc}K\left(x_0,y_0\right)& K\left(x_0,y_1\right)& \dots & K\left(x_0,y_{n-1}\right)\\ K\left(x_1,y_0\right)& K\left(x_1,y_1\right)& \dots & K\left(x_1,y_{n-1}\right)\\ ⋮& ⋮& \dots & ⋮\\ K\left(x_{n-1},y_0\right)& K\left(x_{n-1},y_1\right)& \dots & K\left(x_{n-1},y_{n-1}\right)\end{array}]$$

(35)

The bounded operator $A$ is ill conditioned and any numerical attempt to directly solve (35) will fail.

3. Iterative Method

We first establish the notations for the linear operator Equation (3). We used the vector $u$ to denote the true solution and used the vector $v$ to denote the estimate solution; perhaps $v$ is generated by the iterative method. There are two significant measures; the first one is the error vector $e,$

$$e=u-v$$

(36)

Unfortunately, the error is unknown because the vector $u$ is not given. The second measure is residual $r$, defining how well $v$ approximates $u$; it is given by

$$r=f-Av$$

(37)

The residual is simply the amount by which the approximation $v$ fails to satisfy the original problem (3). There is an equation defining the relationship between the error $e$ and residual $r$

$$e=Ar$$

(38)

In order to improve the estimation of the vector $v$, we solve Equation (38). After that, we compute the new estimation by using the definition of the error equation

$$u=e+v$$

(39)

Landweber-Type Method

The Landweber-type iteration method will be used in this work.

$$v^{k+1}=v^k+{\lambda }_k{A}^{T}M(f-A{v}^k),k=0,1,\dots ,$$

(40)

where ${\lambda }_k$ is the relaxation parameter and $M$ is the symmetric positive definite operator. We used the classical Landweber method, where the operator $M$ is defined as identity matrix [36]. The Landweber-type method has a self-regularization property. This property is very important to reduce the errors in the regularized solution by controlling the mechanism of convergence. In order to enhance convergence, we need to select ${\lambda }_k$ relaxation parameter. In [37], there is an important theorem, which proves how the value of the relaxation parameter make iterations of (40) converge to a solution of $\underset{u\in R^n}{\min }||Au-f|{|}_2$.

Theorem 1.

Let ${\lambda }_k=\lambda$ for $k\ge 0$. Then, the iterates of (40) converge to a solution ($v$) of (3) if and only if $0<\lambda <2/{\sigma }_1^2$ with ${\sigma }_1$ the largest singular value of $\sqrt{M}A$. If, in addition, $v^0\in R(A^T)$ then $v$ is the unique solution with minimal norm.

This theorem was proved in [37] with a constant relaxation parameter ${\lambda }_k=\lambda$. To reduce error in the iterative method (40), the following relaxation parameters are studied by [38]:

$${\lambda }_k=\frac{||w^k|{|}^2}{||\sqrt{M}A{r}^k|{|}^2}$$

(41)

where $r^k=f-A{v}^k$ and $w^k=A^{T}M{r}^k$ for $k=0,1,\dots$.

The following Algorithm 1 explains how to write the program code for the Landweber-type iteration method with a variable relaxation parameter.

Algorithm 1 L.T.1

For loop: k: =1, 2, 3,…

$n=size(u,1);$ % determine the size of domain

$v=zeros(n-1,1);$ % initial guess (zeros vector)

$M=one(n-1,n-1);$ % create Identity Matrix

$\begin{array}{l}r=f-(A{v}^k);\\ w=A^{T}Ir;\\ \lambda =\frac{||w|{|}^2}{||IAr|{|}^2};\\ v=v+(\lambda A^{T}M)(f-Av);\end{array}$

End loop

The next Algorithm 2 explains the Landweber-type iteration method with a constant relaxation parameter.

Algorithm 2 L.T.2

For loop: k: =1, 2, 3,…

$n=size(u,1);$ % determine the size of domain

$v=zeros(n-1,1);$ % initial guess (zeros vector)

$M=one(n-1,n-1);$ % create Identity Matrix

$\lambda =number;$ % determine the relaxation parameter

$v=v+(\lambda A^{T}M)(f-Av);$

End loop

Those algorithms can be used in next section as a relax function.

4. Multigrid Method Algorithms

The goal of the multigrid method is to recover the limitations in the iteration method, at least in its primary stages, by using a good initial guess. A well-known technique for obtaining an enhanced initial guess is by performing some preliminary iterations on a coarse grid. Iterations on a coarse grid are less expensive because there are fewer unknowns to be updated. We noted the grid name as with points. Add the new notations to (3).

$$A_h{u}_h=f_h, \mathrm{on} {\Omega }^h$$

(42)

the coarse grid ${\Omega }^{2h}$ with $(\frac{n}{2}-1)$ points noted as

$$A_{2h}{u}_{2h}=f_{2h}.\mathrm{on} {\Omega }^{2h}$$

(43)

The multigrid method consists of two main steps. One is the smoothing operators and the other is the grid transfer operators. We used the Landweber method as a smoothing method in the algorithm. We named this step “relaxing”. In the grid transfer operator step, we have two types of operators: the prolongation and the interpolation operator. We used a simple method to define those operators, which is explained in [30]. The prolongation operator $P$ takes the vector in a coarse grid and produces the fine-grid vector $P_{2h}^h:{\Omega }^{2h}\to {\Omega }^h$. The interpolation operator $I$ transfer a vector from a fine grid to a coarse grid $I_h^{2h}:{\Omega }^h\to {\Omega }^{2h}$ as shown in Figure 1.

The prolongation operator $P_{2h}^h:{\Omega }^{2h}\to {\Omega }^h$ according to the rule $P_{2h}^h{v}^{2h}=v^h,$ where

$$v_{2j}^h=v_j^{2h}$$

(44)

$$v_{2j+1}^h=\frac{1}{2}(v_j^{2h}+v_{j+1}^{2h}),0\le j\le \frac{n}{2}-1$$

(45)

It is clear that the $P_{2h}^h$ operator is mapping from $R^{\frac{n}{2}-1}$ to $R^{n-1}$. It has full-rank null space $N=\left\{0\right\}$. This operator has the general form

$$P_{2h}^h{v}^{2h}=\frac{1}{2}{\left[\begin{array}{ccccc}1& & & & \\ 2& 1& & & \\ 1& 2& & & \\ & 1& \ddots & & \\ & & \ddots & 1& \\ & & & 2& 1\\ & & & 1& 2\\ & & & & 1\end{array}\right]}_{(n-1)\times (\frac{n}{2}-1)}\left[\begin{array}{c}{v}_1\\ v_2\\ ⋮\\ ⋮\\ v_{\frac{n}{2}-1}\end{array}\right]=\left[\begin{array}{c}{v}_1\\ v_2\\ v_3\\ \\ ⋮\\ ⋮\\ \\ v_{n-1}\end{array}\right]=v^h$$

(46)

The interpolation operator $I_h^{2h}$ is mapping $R^{n-1}$ to $R^{\frac{n}{2}-1}$,

$$I_h^{2h}{v}^h=\frac{1}{4}{\left[\begin{array}{cccccccc}1& 2& 1& & & & & \\ & 1& 2& 1& & & & \\ & & \ddots & \ddots & \ddots & & & \\ & & & 1& 2& 1& & \\ & & & & & 1& 2& 1\end{array}\right]}_{(\frac{n}{2}-1)\times (n-1)}\left[\begin{array}{c}{v}_1\\ v_2\\ v_3\\ \\ ⋮\\ ⋮\\ \\ v_{n-1}\end{array}\right]=\left[\begin{array}{c}{v}_1\\ v_2\\ ⋮\\ ⋮\\ v_{\frac{n}{2}-1}\end{array}\right]=v^{2h}$$

(47)

In this work, we will use the V-cycle multigrid method. There are two smoothing parts: the pre-smooth and the post-smooth. We need to apply relax (iterative computing method) to each node to obtain a solution vector. We transfer the vector from fine grid to coarse grid by using an interpolation operator $I$ in the pre-smoothing part. In the post-smoothing part, the vector transfers from coarse grid to fine grid by using the prolongation operator $P$, see Figure 2.

If we apply the V-cycle method twice, we get a W-cycle multigrid. The relax action was implemented by a Landweber-type iteration (40). The relax action reduces the high frequency of error and the result of the relax stage is noted as vector $v_{Lh},L=1,2,4,8,16,32,64,128,\dots$.

Pre-smoothing

• Relax $A_h{u}_h=f$, initial guess $v_h=0$
• Compute residual $r_{2h}=I_h^{2h}{r}_h,r_h=f_h-A_h{v}_h$
• Relax $A_{2h}{e}_{2h}=r_{2h}$, initial guess $v_{2h}=0$
• Compute residual ${\widehat{r}}_{4h}=I_{2h}^{4h}{\widehat{r}}_{2h},{\widehat{r}}_{2h}=r_{2h}-A_{2h}{v}_{2h}$
• Relax $A_{4h}{e}_{4h}={\widehat{r}}_{4h}$, initial guess $v_{4h}=0$
• Compute residual ${\stackrel{⌣}{r}}_{8h}=I_{4h}^{8h}{\stackrel{⌣}{r}}_{4h},{\stackrel{⌣}{r}}_{4h}={\widehat{r}}_{4h}-A_{4h}{v}_{4h}$

Lowest Level

• Solve $A_{Lh}{e}_{Lh}=r_{Lh}$

Post-smoothing

• correct $v_{4h}\leftarrow v_{4h}+P_{8h}^{4h}{v}_{8h}$
• Relax $A_{4h}{e}_{4h}={\widehat{r}}_{4h}$, initial guess $v_{4h}$
• correct $v_{2h}\leftarrow v_{2h}+P_{4h}^{2h}{v}_{4h}$
• Relax $A_{2h}{e}_{2h}=r_{2h}$, initial guess $v_{2h}$
• correct $v_h\leftarrow v_h+P_{2h}^h{v}_{2h}$
• Relax $A_h{u}_h=f$, initial guess $v_h$
  where $A_{2h}=I_h^{2h}{A}_h{P}_{2h}^h$ and $A_{4h}=I_{2h}^{4h}{A}_{2h}{P}_{4h}^{2h}$.

5. Numerical Results

We will compare between the Landweber-type method and V-cycle multigrid method by considering the time and error level.

5.1. Boundary Value Problem

Considering the inverse boundary value problem (7)–(11) for the heat equation, we need to find $u\left(t\right)\in H^4\left[0,T\right]$. The exact solution will be $u\left(t\right),$ and the input function for inverse problem $u\left(x_0,t\right)=f\left(t\right)$, where the $x_0\in \left(0,1\right),t\in \left[0, T\right]$, $x_0=0.5$ and $T=5,$ as shown in Figure 3.

We tested the two algorithms in different domains ${\Omega }^h$ with different noise levels. The results are shown in

Table 1. Results for boundary value problem in different ${\Omega }^h$.

${\mathbf{\Omega }}^{\mathit{h}}$	Algorithm	$\mathit{\delta }$	CPU Time Seconds	$||\mathit{A}\mathit{v}-\mathit{f}|{|}_2$	No. of Iterations
64	Algorithm 1 L.T.1	0.01	0.042	0.001	400
		0.04	0.041	0.02
	Algorithm 3 V.M	0.01	0.039	0.015
		0.04	0.035	0.024
128	Algorithm 1 L.T.1	0.02	0.132	0.016	300
		0.05	0.147	0.029
	Algorithm 3 V.M	0.02	0.122	0.019
		0.05	0.103	0.032
256	Algorithm 1 L.T.1	0.02	0.46	0.02	200
		0.08	0.481	0.04
	Algorithm 3 V.M	0.02	0.593	0.026
		0.08	0.566	0.05
512	Algorithm 1 L.T.1	0.01	10.98	0.023	500
		0.06	11.03	0.039
	Algorithm 3 V.M	0.01	4.914	0.027
		0.06	4.838	0.03
1024	Algorithm 1 L.T.1	0.02	63.636	0.034	500
		0.06	66.03	0.042
	Algorithm 3 V.M	0.02	31.38	0.033
		0.06	32.04	0.041

. When we called the Call Algorithm 1 L.T.1 or Algorithm 3 V.M algorithm, we used Call Algorithm 1 L.T.1 or Algorithm 1 L.T.1 in Relax ($Av=f$) step with 100 iterations.

Algorithm 3 V.M

$grid(A,v,f)$   % define function grid with three input

$n=size(v,1);$ % determine the size of domain

$\mathrm{Relax} (Av=f$) % Call Algorithm 1 L.T.1 or Algorithm 2 L.T.2 with low iteration times

If (n > size of lowest level)

$r=f-Av;$ % Compute Residual

for i = 1:n % Create prolongation matrix

P(2*i − 1, i) = 1;

P(2*i, i)= 2;

P(2*i + 1, i) = 1;

End

$I=P^T;$ % Create interpolation matrix

$r\_to\_c=0.25*Ir;$% from fine-grid to coarse-grid

$A\_n=0.25*I*A*0.5*P;$$\%\mathrm{Create} A_{Lh}$ interpolation matrix

$v=zeros\left(\mathrm{size}\right(r\_to\_c,1),1);$ % initial guess (zeros vector)

$e\_to\_c=grid(A\_n,v,r\_to\_c);$ % used recursion function

$e=0.5*P*e\_to\_c;$ %from coarse-grid to fine-grid

$v=v+e;$ % correct

end

$\mathrm{Relax} (Av=f$) % Call Algorithm 1 L.T.1 or Algorithm 2 L.T.2 with low iteration times with low iteration times

$output=v;$

$End\_grid;$

The above table shows the Algorithm 3 V.M has speeded the computation time compared with the L.T.1 algorithm. The solutions for two algorithms in ${\Omega }^h$ are shown in Figure 4.

5.2. Initial Value Problem

Considering the inverse initial value problem (28)–(31) for the heat equation, we need to find $u\left(x\right)\in L_2\left[0,1\right]$. The exact solution will be $u\left(x\right)=\sin \pi x,0\le x\le 1$. We created the input function for the inverse problem, as shown in Figure 5. $u\left(x,T\right)=f\left(x\right),0\le x\le 1$ and $T=0.01$.

In this numerical example, we will use the Landweber-type method with constant relaxation parameter ${\lambda }_k=\lambda$; also, we will use the same iterative method in V-cycle multigrid with constant relaxation parameter. We compared the two algorithms with a high scale of domains ${\Omega }^h$. The results are shown in

Table 2. Results for initial value problem in different ${\Omega }^h$.

${\mathbf{\Omega }}^{\mathit{h}}$	Algorithm	$\mathit{\delta }$	CPU Time Seconds	No. of Iterations	$||\mathit{A}\mathit{v}-\mathit{f}|{|}_2$
512	Algorithm 2 L.T.2	0.064	0.760	77	0.041
	Algorithm 3 V.M		0.127
1024	Algorithm 2 L.T.2	0.09	4.405	84	0.0567
	Algorithm 3 V.M		0.684
2048	Algorithm 2 L.T.2	0.13	28.718	91	0.0816
	Algorithm 3 V.M		4.026
4096	Algorithm 2 L.T.2	0.18	252.61	102	0.115
	Algorithm 3 V.M		28.381

. When we called the V.M. algorithm, we used L.T.1 in Relax ($Av=f$) step with 10 iterations.

From the above results, we see the Algorithm 3 V.M Algorithm successfully applied to obtain a good approximating solution with low time in high-scale domains ${\Omega }^h$, see Figure 6.

6. Conclusions

This work deals with algorithms for solving the boundary value problem and the initial value problem for the heat equation and some results were collected. The forward problem was solved by the separation-of-variables method and then the PED is converted to an integral equation of the first kind. By using the discretization method, we converted the integral equation to a linear operator equation for the first kind. The Landweber iteration method was successfully applied to obtain an estimated solution. The V-cycle multigrid method is applied to make the computation cost for the iterative method cheap and more accurate. It is clear that the V.M. algorithm has a recursion function, named grid. W used a low number of iterations in each relaxing step inside the grid function. The results show the CPU time for solving each equation in a specific domain size is accelerated by using a low iteration number in grid function through all the low domains in V.M. algorithm. All results in this numerical example show that the multigrid method is fast, more convergent and accurate. In future work, parallel computing will be used to make the cost of computing cheaper.