Numerical simulations for initial value inversion problem in a two-dimensional degenerate parabolic equation

1.   Introduction

In this work, we consider an inverse problem of recovering the initial value in a two-dimensional degenerate parabolic equation from terminal observation. The problem can be stated in the following form:

● Problem P: Consider the following two-dimensional degenerate parabolic equation

where $a(X)$ is a given smooth function satisfying

and $\Omega {\rm{ = }}{(0, l)^2} \subset {R^2}$ is a two-dimensional bounded plane region. $\varphi (X)$ is an unknown initial function in (1.1). Assume that an additional condition is given as:

where $T > 0$ is a fixed time, and $\psi (X)$ is the observation data. In this paper, we use $X$ to denote the two-dimensional spatial variable $(x, y).$ The inverse problem is to determine the functions $u$ and $\varphi $ satisfying (1.1) and (1.2).

In recent years, many inverse problems of degenerate parabolic equations have been proposed to meet actual needs (see ^{[2,6,16,29,33]}), which often appear in many engineering and financial fields. It is noted that the degenerate parabolic equation has attracted the attention of many scholars and has been successfully applied to porous media (see ^{[1,3,4,5,25,27]}). These results have a profound impact on the development of related fields. Eq (1.1) also belongs to degenerate parabolic equations. Compared with the traditional parabolic equation, the degenerate one can get rid of the limitation of boundary conditions in some cases. According to the famous Fichera theory, whether or not the boundary conditions of the degenerate parabolic equation should be given is determined by the corresponding symbol of the Fichera function (see ^[23]). If the sign is nonnegative, the boundary condition is not necessary; otherwise, it is the opposite. Using the Fichera's theorem, we can easily calculate that the value of Fichera function of Eq (1.1) is zero at the corresponding degenerate boundary. Therefore, Eq (1.1) does not need to give boundary conditions.

For general heat conduction problems, the initial value function usually represents the initial temperature distribution field of the object. While for the inverse heat conduction problem, the main task is to reconstruct the initial temperature with some additional observation data. This is a kind of severly ill-posed problem, that is, the small change of observation data will lead to the big change of solution. This is easy to understand. In fact, most things in nature are irreversible. If we can easily reconstruct the past temperature by using the current measurement data, then we can easily obtain the temperature of the earth 65 million years ago. Maybe we can know how dinosaurs died out, but it's impossible. (The reason why dinosaurs died out has always been a mystery. One is that 65 million years ago, the earth's climate suddenly changed, the temperature dropped sharply, dinosaurs were frozen to death. Another view is that because of the eruption of volcanoes, a large amount of carbon dioxide spewed out, resulting in the earth's rapid greenhouse effect, making food death, and eventually the extinction of dinosaurs.) However, this is a problem that must be solved in some engineering fields. For example, in steel smelting, we need to obtain the ideal temperature distribution at a certain time, so how to accurately control the initial temperature is such a problem.

For the inverse initial value problems of nondegenerate parabolic equations, there have been a lot of research work, and abundant theoretical and numerical results have been obtained (see ^{[9,10,11,12,13,14,15,17,18,19,20,31]}). In ^[17], the data assimilation problem for the soil water movement model is analysed carefully, where the numerical computation is realized through searching the minimum of the cost function iteratively by a descent gradient method.

For the case of degeneration, the literatures are quite few. These works mostly discuss the coefficient inversion of degenerate parabolic equation from the angle of theoretical analysis, such as source term coefficient or zero-order term coefficient, etc. Documents ^[7,8,22] are earlier ones involving degenerate parabolic equations. In these papers, the authors establish the Carleman estimate in the case of boundary degradation, and discuss the null controllability of related equations. Since then, these works have been widely extended and successfully applied to the inverse coefficient(s) problem of degenerate parabolic equations (see ^[2,6,29]). In ^[26], the uniqueness of the inverse source problem of one-dimensional degenerate heat equation

is proved by the extremum principle, and the numerical simulation is also given. In ^[29], the authors consider the inverse problem of determining the insolation function, where the underlying equation is the nonlinear Sellers climate model. The uniqueness and Lipschitz stability of the solution are also obtained. In ^[2], the inverse source problem for multidimensional parabolic operators of Grushin type

where ${\Omega _1}$ is a bounded open region that contains the origin and thus the principle coefficient is degenerate.

In ^[32], numerical inversion for an inverse initial value problem is studied carefully, where the mathematical model is governed by a one-dimensional degenerate parabolic equation. In ^[33], the authors consider an inverse problem of simultaneously reconstructing the initial value and source term in a degenerate parabolic equation. Uniqueness and stability of the solution are obtained by the Carleman estimate. It should be pointed out that the vast majority of literature in the previous work is about one-dimensional degenerate parabolic equation, and there are few relevant literature about higher-dimensional case. The main reason is that in the case of higher dimension, the degree and way of coefficient degeneration is much more complex than that in the case of one dimension. However, the application in high dimension is more extensive than that in one dimension. For example, in the field of computer signal processing, one-dimensional functions can only represent one-dimensional signals, but two dimensional functions can represent more complex information, such as digital images. If it is a color image or a video, the dimension of the function will be higher.

This work studies the inverse initial value problem for a two-dimensional degenerate diffusion equation. In a sense, this is an extension of the one-dimensional form, and its thermal conductivity is also extended from one-dimension to two-dimension. The change of thermal conductivity enables us to take into account both anisotropic diffusion and slow diffusion, but it also brought essential difficulty for theoretical analysis and numerical simulations. In particular, these results in the case of two-dimension can be easily applied to the field of digital image processing. The outline of the manuscript is as follows: In Section 2, based on the finite volume method, a difference scheme of the forward problem is proposed. The stability and convergence of the scheme are proved in Section 3. In Section 4, the Landweber iteration and conjugate gradient method (CGM) are used to solve the numerical solution for the inverse problem. Some typical numerical examples are presented in the last section to show the validity of the inversion method.

2.   Difference scheme of forward problem

The main purpose of this paper is to study the problem P from the perspective of numerical analysis. The conclusions on the theoretical aspects are only briefly introduced, and those proofs can be found in the references (see ^[20,23]). Let

where $\varepsilon $ is an arbitrarily small positive constant. If $a$ and $\varphi $ satisfy

then there exists a unique weak solution $u(X, t) \in {W^{k, \infty }}(\bar Q)$ to the Eq (1.1) (see ^[23]). In this paper, we always assume that the solution of the direct problem has sufficient regularity (at least belongs to ${C^{2 + \alpha, 1 + \alpha /2}}(\bar Q), {\rm{ 0 < }}\alpha < 1$). Considering the embedding property of space, it seems that $k \geqslant 5$ is enough. Moreover, the uniqueness and stability of the solution for problem P can be proved by the logarithmic convexity method (see ^[20]).

In this article, we use the finite difference method to obtain the numerical solution of Eq (1.1). For the sake of simplicity, we assume that the domain $\Omega $ is a square. For the case of rectangle, we can discuss it similarly, which is not the key of the problem.

Assumed that the domain $\Omega $ is divided into a mesh $J \times J$ with the spatial step size $h = {h_x} = {h_y} = \frac{l}{J}$ and the time step size $\tau = \frac{T}{N}, $ respectively.

Grid points $({x_i}, {y_j})$ of the layer ${t_n}({t_n} = n\tau, n = 0, 1, 2, \cdots, N)$ are defined by

The notation $u_{ij}^n$ is used for the finite difference approximation of $u(ih, jh, n\tau)$.

Remark: The region can be generalized to a general rectangular region where the spatial steps can be taken as

Using the finite volume method (see^[21,28]), one can derive the following implicit difference scheme of (1.1).

for $1 \leqslant i$, $j \leqslant J - 1$ and $0 \leqslant n \leqslant N - 1$, where${a_{ij}} = a({x_i}, {y_i}).$

In this section, we always assume that $a(X)$ at least belongs to ${C^2}(\bar \Omega)$. Note that on the boundary of $\Omega $, Eq (1.1) is degenerated into the following first order hyperbolic equations

Therefore, the boundary conditions of (1.1) can be divided into three cases according to the value of $\nabla a(X)$.

Case 1. If ${a_x}\left| {_{x = 0, l}} \right. = {a_y}\left| {_{y = 0, l}} \right. = 0$, we get the follow boundary condition from (2.2) and (1.1)

Case 2. If, the boundary conditions can be discretized as

Case 3. If at least one of the four values ${a_x}\left| {_{x = 0, l}} \right., {a_y}\left| {_{y = 0, l}} \right.$, is not zero, we can discuss it similarly.

The initial condition

can be discretized by

Without loss of generality, in the following discussions, the boundary conditions are assumed to be Case 2.

In order to write the difference equation conveniently, we introduce the following symbol

Setting

the difference scheme of (2.1) can be rewritten as

Let

Combining (2.3) and (2.5), the diference scheme can be further written as the following matrix equation:

Solving Eq (2.6), we can obtain the numerical solution of the forward problem.

3.   Stability and convergence

In order to obtain the stability and convergence of the difference equation, we introduce the following notations:

and the norms

where

Before proving the stability, we first introduce the following lemmas.

Lemma 3.1. (discrete Gronwall inequality). Let be a non-negative sequence which satisfies

where $C > 0$ and $G \geqslant 0$ are two constants. Then we have the following estimate:

The proof can be found in [28,lemma 3.8].

Lemma 3.2. Assuming that $\{ u_{ij}^n|0 \leqslant i, j \leqslant J, 0 \leqslant n \leqslant N\} $ be the solution of difference Eqs.(2.1), (2.3) and (2.4), then we have for relatively small $\tau $,

where the constant $A$ is dependent on $a(X)$.

Proof. Using the above notations the difference equations (2.1), (2.3), (2.4) can be rewritten as the following system:

Multiplying ${h^2}u_{ij}^{n + 1}$ on both sides of (3.1) and summing from $i$, $j = 1$ to $J - 1$, we have

For the first term of (3.7), we have

where we have used the following basic inequality:

From (3.7) and (3.8), we obtain

For the second term of (3.7), we have

Similarly, for the third part of (3.7), we have

where we have used the boundary informations of $a(X).$

Combining (3.8)-(3.11), we obtain

Multiplying $\frac{1}{2}{h^2}u_{0, j}^{n + 1}$ on both sides of Eq (3.3) and summing from $j = 0$ to $J$, and using the Cauchy inequality, we get

which implies

Likewise, for (3.4), (3.5) and (3.6), we have

From (3.12) to (3.16), we obtain

Noticing the smoothness of $a(X)$, we have from the Taylor's formulation (Because in the expansion, only one variable is expanded, and the other keeps unchanged. For the sake of simplicity, another variable is omitted when we take the value in the expansion.)

Similarly

From (3.17) -(3.21), we can get

By the smoothness of $a(X)$, we may assume

Using (3.22) and (3.23), we have

which implies

So, if we choose $ \mathrm{\tau } $ small enough such that

then we have

From (3.25) and (3.26), we get

By the Lemma 3.1 and (3.27), we get

This completes the proof of Lemma 3.2.

The stability of the scheme is the direct consequence of Lemma 3.2.

Remark: For a more general case:

we may have

where $C$ is dependent on $a(X)$.

The proof is similar to that of Lemma 3.2.

Theorem 3.3. Assume that the initial data $\varphi (X) \in {C^2}(\Omega)$ in (1.1), $\{ u({x_{ij}}, {y_{ij}}, {t_n})|0 \le i, j \le J, 0 \le n \le N\} $ is the solution of (1.1), and $\{ u_{ij}^n|0 \le i, j \le J, 0 \le n \le N\} $ is the solution of (3.1)-(3.6). Then, as $\tau \to 0$ and $h \to 0$, the approximate solution $u_{ij}^n$ unconditionally converges to the real solution and the convergence rate is $O\left({\tau + {h^{\frac{3}{2}}}} \right)$, that is

Proof. Let

Subtracting (1.1) from (3.1) to (3.6), we get the error equation:

Using the regularity of the exact solution $u(X, t) \in {C^{2 + \alpha, 1 + \alpha /2}}(\bar Q), {\rm{ 0 < }}\alpha < 1, $ it can be easily seen that

Combining (3.29-3.39) and using (3.28), one can easily get the conclusion.

4.   Numerical experiment

4.1.   Landweber iterative method

It can be easily seen that the Problem P is a linear problem. Let $K$ be the parameter-to-data mapping, that is

where $u$ is the weak solution of (1.1). In this article, we would like to use the Landweber iteration method to solve the operator Eq (4.1).

Lemma 4.1.1. The mapping defined in (4.1) is a linear bounded operator from to and it has the following bound

Proof. From (1.1), we have

That is

This completes the proof of Lemma 4.1.

Note that (4.1) can be rewritten as

for some $\alpha > 0$, where ${K^ * }$ is the adjoint operator of $K$. Then, we use the iteration method to solve (4.2), that is, computing

From (4.3) and the definition of the operator $K$, we have

where ${u_{m - 1}}$ is the solution of (1.1) with $\varphi = {\varphi _{m - 1}}$.

The following Lemma (4.2) illustrates the specific form of the adjoint operator${K^ * }$.

Lemma 4.1.2. For any given $\phi (X)$, let ${K^ * }\phi = v(\bullet, 0).$Then $v$ satisfies the following backward parabolic equation:

Proof. From (1.1) and (4.5) and using the Green's theorem, we have

that is

By the definition of $K$, we have

that is

This completes the proof of Lemma 4.1.2.

Based on (4.3), (4.4) and Lemma 4.1.2, the procedure of the iterative algorithm can be stated as follows:

Step 1. Choose an initial value of iteration For simplicity, we can choose

Step 2. Solve the degenerate parabolic problem (1.1) to obtain the solution, where

Step 3. Solve the adjoint problem of (1.1)

to get the solution $v(X, t)$ which is denoted by ${v_0}(X, t).$

Step 4. Let

where $\alpha > 0$, and let ${u_1}(X, t)$ be the solution of (1.1) with $\varphi = {\varphi _1}(X)$.

Step 5. Select an arbitrarily small positive constant as error bound. Compute

and compare it with $\varepsilon $.

If

then stop the iteration scheme and take $\varphi (X) = {\varphi _1}(X).$

Otherwise go to Step 3. Let ${\varphi _1}(X)$ be a new initial value of iteration and go on computing by the induction principle.

The proof of convergence for the Landweber iterative method can be found in ^[18,19].

Remark 4.1.1: The paramter is called acceleration factor, and is usually required. The iterative algorithm (including the CGM below) is also applicable to the observation data with noise, i.e., can be replaced by, where satisifes

and $\delta > 0$ is the noise level.

4.2.   Conjugate gradient method (CGM)

In this part, we introduce the CGM iteration to solve the inverse problem of (1.1). Generally speaking, the convergence rate of CGM is much faster than that of Landweber iterative method. But at the same time, the iterative procedure of conjugate gradient method is more complex than that of Landweber iterative method. Both methods have advantages and disadvantages.

The CGM iteration algorithm can be stated as follows:

Step 1. Set and choose an initial value of iteration. For simplicity, we can choose

Step2. Solve the initial value problem (1.1) to obtain the solution${u_k}(X, t)$, where$\varphi (X, t) = {\varphi _k}(X, t)$.

Step3. Solve the following conjugate equation

to obtain the solution ${V_k}(X, t)$, where ${r_k}(X): = \psi (X) - {u_k}(X, T)$ is the residual.

Step4. Calculate

where

Step5. Solve the following equation

to obtain the solution $u(\bullet, T) = {u_k}(\bullet, T): = {q_k}.$ Set

and let

Step6. Increase and go to (Step2). The stopping rule is taken as follows: Let be a fixed number. Then we stop the algorithm at the first occurrence of such that ${\left\| {{r_k}} \right\|_{{L^2}(\Omega )}} \le \lambda \delta .$

The proof of stability and convergence for the CGM can be found in ^[24].

5.   Numerical results

In this section, we try to reconstruct the initial temperature field through the terminal observation data of $t = T$, which is a serious ill-posed problem. For a given initial condition, the change of temperature will decrease exponentially with the increase of time. So the terminal moment $T$ cannot be too large. We take $T = 0.5$ in the process of numerical simulation. The other basic parameters are

and the diffusion coefficient $a(X)$ is taken as

The direct problem is numerically solved by the finite difference method. In the process of numerical computations, the spatial step size $h$ and the time step size $\tau $ are taken as:

Example 1. In the first numerical experiment, we take

the additional condition $\psi (X)$ is given by

where $u(x, y, T; \varphi (X))$ is the numerical solution of (1.1) with the given initial value $\varphi (X)$.

The noisy data is also considered in the numerical experiment and generated in the following form

5.1.   The results of Landweber iteration

As can be seen from the above Figures 1, 2, when the number of iterations reaches 1000, the reconstruction effect is remarkably good. The initial guess is taken as zero. Obviously, this initial guess is not good at all, but the convergence of the iterative algorithm is very stable and the error order is about ${10^{ - 2}}$. After adding noise, the reconstruction effect is also very good (see the Figures 3, 4).

5.2.   The results of CGM iteration

When the number of CGM iterations reaches 300, the error reaches ${10^{ - 3}}$, which shows that the reconstruction effect is satisfactory (see the Figures 5, 6). Furthermore, it shows that the convergence speed of CGM is faster than that of the Landweber iterative method. For the case of adding noise, it can be seen from the figure that the reconstruction effect is also satisfactory (see the Figures 7, 8).

Example 2. In the second numerical experiment, we consider a digital image deblurring problem. We chose two gray images as examples. For the convenience of calculation, the pixels of these two images are $ 128\times 128 $. The diffusion coefficient is taken as

and the observation time is taken as $T = 0.5$. The initial value $\varphi (X)$ is the original clear image, while the blur image $\psi (X)$ is obtained from the clear image through the above diffusion process. Our task is to restore the original clear image from the blur image. We have tested the two algorithms, and the results are similar. The convergence speed of CGM algorithm is faster than that of Landweber iterative algorithm. The Landweber iterative algorithm needs 100 iterations, while CGM only needs about 50 iterations. The reconstruction effect of the two images is shown in Figures 9-12 (since the reconstructed results are similar, we only use the results obtained by the Landweber iterative algorithm).

From these four groups of pictures, we can see that the image restoration effect is very good. Because the fuzzy part is mainly concentrated in the middle of the image (this is true of most images in reality (see ^[30])), our diffusion function can describe the evolution process quite well. After a period of iterations, the middle blurred part of the image is almost completely restored, and the edge part of the image is also well preserved.

6.   Conclusions

In this paper, the inverse initial value problem for a two-dimensional degenerate parabolic equation is discussed. The finite volume method is used to construct the difference scheme for the forward problem, and the stability and convergence of the difference scheme is also proved. Then, we use the Landweber iteration and CGM to construct the calculation process of the inverse problem and do numerical experiments. Numerical results show that our algorithm is stable and converges quickly. These methods and results can be widely applied in engineering heat conduction problems and computer graphics and image processing.

Acknowledgments

We would like to thank the anonymous referees for their valuable comments and helpful suggestions to improve the earlier version of the paper. The work is partially supported by the National Natural Science Foundation of China (Grant Nos. 61663018, 11961042), Foundation of A Hundred Youth Talents Training Program of Lanzhou Jiaotong University, and NSF of Gansu Province of China (No. 18JR3RA122).

Conflict of interest

The authors declare no conflict of interest in this paper.