- Email: [email protected]
A multigrid–homotopy method for nonlinear inverse problems
Computers and Mathematics with Applications xxx (xxxx) xxx
Contents lists available at ScienceDirect
Computers and Mathematics with Applications journal homepage: www.elsevier.com/locate/camwa
A multigrid–homotopy method for nonlinear inverse problems Tao Liu School of Mathematics and Statistics, Northeastern University at Qinhuangdao, Qinhuangdao 066004, PR China
article
info
a b s t r a c t In the present contribution, we develop a novel method combining the multigrid idea and the homotopy technique for nonlinear inverse problems, in which the forward problems are modeled by some forms of partial differential equations. The method first attempts to use the multigrid method to decompose the original inverse problem into a sequence of sub-inverse problems which depend on the grid variables and are solved in proper order according to the grid size from the coarsest to the finest, and then carries out the inversion on the coarsest grid by the homotopy method. The strategy may give a rapidly and globally convergent method. As a practical application, this method is used to solve the nonlinear inverse problem of a nonlinear convection–diffusion equation, which is the saturation equation within the two-phase porous media flow. We demonstrate the effectiveness and merits of the multigrid–homotopy method on two actual model problems. © 2019 Elsevier Ltd. All rights reserved.
Article history: Received 2 April 2019 Received in revised form 25 September 2019 Accepted 30 September 2019 Available online xxxx Keywords: Multigrid Homotopy Regularization Nonlinear inverse problems Convection–diffusion equation
1. Introduction Generally speaking, we can formulate the nonlinear inverse problem as the nonlinear operator equation Φ(q) = f,
(1.1)
or the nonlinear least-square problem min ∥Φ(q) − f∥2 ,
(1.2)
q
where q is the parameter to be inverted, f is the measurement data, and Φ is the forward model, which can be computed from the appropriate governing partial differential equation. By using a standard procedure of Newton-type iterative method to linearize Eq. (1.1), we can get Φ′ (q)(q − qi ) = − [Φ(qi ) − f] ,
(1.3) ′
which can be used to solve Eq. (1.1). But it is general that Φ is not invertible, so we must regularize it [1]. The most famous regularization approach is the Tikhonov’s regularization approach [2]. Many researchers [3,4] have studied its analysis for the nonlinear ill-posed problems, and Bakushinskii [5] developed an iteratively regularized Gauss–Newton method qi+1 = qi − Φ′ (qi )∗ Φ′ (qi ) + βi I
[
]−1 [
i = 0, 1, 2, . . .
Φ′ (qi )∗ (Φ(qi ) − f) + βi (qi − q0 ) ,
]
(1.4)
E-mail address: [email protected]. https://doi.org/10.1016/j.camwa.2019.09.023 0898-1221/© 2019 Elsevier Ltd. All rights reserved.
Please cite this article as: T. Liu, A multigrid–homotopy method for nonlinear inverse problems, Computers and Mathematics with Applications (2019), https://doi.org/10.1016/j.camwa.2019.09.023.
2
T. Liu / Computers and Mathematics with Applications xxx (xxxx) xxx
Here, superscript ∗ refers to the adjoint of an operator. Introducing Tikhonov’s regularization to Eq. (1.3), we can obtain the Levenberg Marquardt method [6] qi+1 = qi − Φ′ (qi )∗ Φ′ (qi ) + βi I
[
]−1
Φ′ (qi )∗ [Φ(qi ) − f] ,
i = 0, 1, 2, . . .
(1.5)
Furthermore, the Landweber method is in essence another interesting regularization method. For example, Hanke et al. [7] proposed a Landweber iterative scheme to solve Eq. (1.1) qi+1 = qi − Φ′ (qi )∗ [Φ(qi ) − f] ,
i = 0, 1, 2, . . .
(1.6)
To solve inverse problems, such Newton-type methods work by performing the computations using a fixed discrete grid. While tremendous advance has been made in aspect of reducing the computational complexity of the fixed grid methods, calculation cost still deserves great attention. Perhaps more importantly, all these methods are locally convergent. In another word, good initial estimates have to be provided when these methods are used to solve Eq. (1.1). In order to accelerate the convergence rate and widen the convergence region, this paper combines the multigrid idea with the homotopy technique, and thereby proposes a multigrid–homotopy method for nonlinear inverse problems, in which the forward problems are modeled by some forms of partial differential equations. The multigrid method is known to be one of the best methods to solve the forward problems of partial differential equations [8–11]. By computing corrections on a coarser grid level and then prolongating them back to the original problem, it can reduce the computation complexities and remove the smooth error components. In recent years, this method has been extended to the study of inverse problems. Successful studies of the multigrid method include the inverse problems of the acoustic wave equation [12], the nonlinear diffusion equation [13], the nonlinear convection– diffusion equation [14], the three-dimension elliptic equation [15], the multiphase porous media flow [16], the electrical impedance tomography [17,18] and the Bayesian optical diffusion tomography [19,20]. It is shown in these papers that the multigrid method is effective on inverse problems. The homotopy method, due to its global convergence property, has been extensively used in the field of nonlinear problems. In 1998, Liao [21] employed the basic ideas of the homotopy in topology to propose a homotopy analysis method for nonlinear problems. Different from the homotopy analysis method, the homotopy perturbation method proposed by He [22,23] is actually a new perturbation technique coupled with the homotopy technique, and has been successfully used to solve nonlinear differential equations [24,25]. Based on the He’s homotopy perturbation method, Marinca and Herişanu [26] proposed an analytical technique, called the optimal homotopy perturbation method, and employed it to study a nonlinear dynamic analysis of an electrical machine rotor-bearing nonlinear system. The optimal homotopy asymptotic method established by Marinca and Herişanu [27] combines the features of the homotopy concept with an efficient computational algorithm which provides a simple and rigorous procedure to control the convergence of the solution. A series of papers by Marinca and Herişanu [28–30], Hashmi et al. [31] have proved the effectiveness, generalization and reliability of this method and obtained solutions of currently important application in science and engineering. Recently, the homotopy method has been also used to solve a wide range of inverse problems. Keller and Perozzi [32] may be the first authors who applied homotopy to the inverse problems. Then, Vasco [33] presented the homotopy method for the inverse problem and demonstrated its use for the traveltime tomography. Everett [34] transformed a discrete inverse problem of the two-dimension Helmholtz equation to a polynomial system, and used the homotopy method to solve it. Jegen et al. [35] used the homotopy method to solve the geophysical inversion of electromagnetic data, which depended on the normal equations. Liu et al. [36] used the homotopy method to identify the diffusion parameters in a nonlinear convection–diffusion equation. Han and Fu developed the homotopy method for the inverse problem of the acoustic wave equation [37] and waveform inversion constrained by well-log data [38]. The rest of the paper is organized as follows. In Sections 2 and 3, the homotopy method and the multigrid method are described respectively. Next, in Section 4 the proposed multigrid–homotopy method is presented. In Section 5, we demonstrate the merits and effectiveness of the proposed method by the inversion of a nonlinear convection–diffusion equation, which is the saturation equation within the two-phase porous media flow. Finally, Section 6 summarizes the conclusions. 2. Homotopy method If we have a good prior estimate q˜ for the solution q† of Eq. (1.1), we can introduce the Tikhonov regularization to approximate Eqs. (1.1) or (1.2), that is min ∥Φ(q) − f∥2 + ϵ∥q − q˜ ∥2 ,
[
]
q
(2.1)
where ϵ is the regularization parameter. Then, the regularized Gauss–Newton method (1.4) can be used to solve Eq. (2.1). However, in practical situation, this assumption is not realistic. In general, a good prior estimate q˜ cannot be provided for Eq. (2.1), so min ∥Φ(q) − f∥2 + ϵ∥q∥2 ,
[
q
]
(2.2)
Please cite this article as: T. Liu, A multigrid–homotopy method for nonlinear inverse problems, Computers and Mathematics with Applications (2019), https://doi.org/10.1016/j.camwa.2019.09.023.
T. Liu / Computers and Mathematics with Applications xxx (xxxx) xxx
3
is used to replace Eq. (2.1). Then, the Levenberg Marquardt method (1.5) can be used, and many trials of initial estimates q0 are made to obtain a satisfactory approximation to q† . This is a serious problem especially for the reservoir history matching. In the field of reservoir history matching, one is interested in inferring as many parameters as there are grid cells in the simulator, and this number usually exceeds 105 for a field simulation. A single field-model reservoir simulation may take several hours to complete and many such simulator runs are required to obtain an inversion result with inversion methods such as Levenberg Marquardt method, so many trials of initial estimates will result in the very high computational cost [39]. To overcome this difficulty, we use the homotopy method to solve Eq. (2.2). First, consider the normal equation of the objective functional in Eq. (2.2) Φ′ (q)∗ (Φ(q) − f) + ϵ q = 0,
(2.3)
and then construct its fixed-point homotopy equation H(q, η) = η Φ′ (q)∗ (Φ(q) − f) + ϵ q + (1 − η)(q − q0 ) = 0,
[
]
(2.4)
where η ∈ [0, 1] is a homotopy parameter, q is an optional initial estimate. When η = 0, the solution of H(q, 0) = 0 is the known initial estimate q0 . When η = 1, H(q, 1) = 0 is just Eq. (2.3). In order to solve Eq. (2.3), we divide the interval [0, 1] into 0 = η0 < η1 < · · · < ηK = 1. For η = ηk (k = 1, 2, . . . , K ), use the iterative method to solve Eq. (2.4) in sequence from k = 1 to k = K . As the solution q0 of H(q, η0 ) = 0 is known, it can be used as the initial estimate to the next equation H(q, η1 ) = 0. Assuming that the approximation solution qk−1 of H(q, ηk−1 ) = 0 has already been obtained, a successive linearization method could be used to compute the solution qk of H(q, ηk ) = 0. Suppose that the ith approximation qki to qk is known. To compute the next approximation qki+1 , in equation H(q, ηk ) = 0, we use the linear function Φ′ (qki )(q − qki ) + Φ(qki ) to replace Φ(q), then have 0
[ ] ηk Φ′ (qki )∗ (Φ′ (qki )(q − qki ) + Φ(qki ) − f) + ϵ q + (1 − ηk )(q − q0 ) = 0. The solution of Eq. (2.5) is just the (i + 1)th approximation qki+1
to q
(2.5)
k
[ ]−1 = qki − ηk Φ′ (qki )∗ Φ′ (qki ) + (ϵ + 1 − ηk )I
[ ] · ηk Φ′ (qki )∗ (Φ(qki ) − f) + ϵ qki + (1 − ηk )(qki − q0 ) , qk0
qki+1
k−1
=q
,
k
q =
qkNk +1
,
i = 0, 1, . . . , Nk ,
(2.6)
k = 1, 2, . . . , K .
Eq. (2.6) has a fast convergence rate similar to the regularized Gauss–Newton method, therefore a good approximation to qk can be obtained by only one iteration when ηk −ηk−1 is small enough. In order to save unnecessary computational cost, we can let ηk = Kk (k = 0, 1, . . . , K ) by an isometric division, and Nk ≡ 0. The solution of Eq. (2.4) will not be affected if K is appropriately chosen. Then, Eq. (2.6) can be transformed into a simplified version
) ]−1 k ϵ+1− I K K [ ( ) ] k ′ k ∗ k · Φ (q ) (Φ(qk ) − f) + ϵ qk + 1 − (qk − q0 ) ,
qk+1 = qk −
[
k
K
Φ′ (qk )∗ Φ′ (qk ) +
(
K
(2.7) k = 0, 1, . . . , K − 1.
In fact, the final iteration result qK of Eq. (2.7) is not the solution to H(q, ηK ) = 0 (also Eq. (2.3)), but it is the one iteration result of Eq. (2.6) with the initial estimate qK −1 and ηK −1 = K K−1 . qK is such a good approximation to the solution of the normal Eq. (2.3) that it can be used as the initial estimate q˜ to the regularized Gauss–Newton method (1.4). Therefore, in order to find the solution q† , we need to perform a regularized Gauss–Newton method. Then, Eqs. (2.7) and (1.4) can form a globally convergent homotopy method for the nonlinear inverse problem (1.1). For the convergence analysis of the homotopy method for nonlinear inverse problems see [40]. 3. Multigrid method In most inverse problems, it is important to choose a discrete grid spacing to describe both the forward problem and its subsequent inverse operation. Usually, a fine grid is attractive because it can reduce modeling error and increase the resolution of the solution, but these are improved at the cost of a tremendous increase in calculation amount. To overcome the problem of calculation complexity, we can use the multigrid method to solve the nonlinear inverse problems. From the beginning we introduce not only one grid but a sequence of grids {Ωl } (l = 0, 1, . . . , L). Let Ω0 be the finest discretization grid, and the grid spacing of Ωl (l = 1, . . . , L) can be determined by multiplying the grid spacing of Ω0 by 2l . The nonlinear operator equation obtained by discretizing the inverse problem on grid Ωl is ΦΩl (qΩl ) = fΩl ,
(3.1)
and the corresponding least-square problem with regularization is min JΩl (qΩl ) = ∥ΦΩl (qΩl ) − fΩl ∥2 + ϵΩl ∥qΩl ∥2 . qΩ l
(3.2)
Please cite this article as: T. Liu, A multigrid–homotopy method for nonlinear inverse problems, Computers and Mathematics with Applications (2019), https://doi.org/10.1016/j.camwa.2019.09.023.
4
T. Liu / Computers and Mathematics with Applications xxx (xxxx) xxx
old Let qold Ωl be some given approximation on grid Ωl , and a coarser resolution representation qΩl+1 can be obtained from qold Ωl by the relation l+1 old qold Ωl+1 = Il qΩl ,
(3.3)
where Ill+1 is the restriction operator. We would like to improve the finer grid approximation qold Ωl by first performing the fixed grid iteration method on the coarser grid Ωl+1 , starting with the initial estimate qold Ωl+1 w old qne Ωl+1 ←− Fixed_Grid_Iteration(qΩl+1 , JΩl (·)),
(3.4)
w where qne Ωl+1 is the iteration result, and then using this result to correct the finer grid approximation old new l old w qne Ωl ←− qΩl + Il+1 (qΩl+1 − qΩl+1 ),
(3.5)
where Ill+1 is the prolongation operator. The fixed grid iteration method can be chosen as any one of the above three common iterative methods: regularized Gauss–Newton method (1.4), Levenberg Marquardt method (1.5), and Landweber method (1.6). old w It is ideal that the new approximation qne Ωl should be better than the old approximation qΩl , specifically, w old JΩl (qne Ωl ) ≤ JΩl (qΩl ).
However, this would not be the case unless the objective functionals are consistent. In reality, if we naively choose a group of objective functionals such as JΩl (·) in Eq. (3.2), the approximation after correction will easily be far away from the optimal solution. To eliminate the problem of inconsistent objective functionals, we use the conditions introduced in [20], which guarantee that the multigrid inversion method is monotonely convergent. Firstly, the adjusted objective functionals are defined as follows adj
JΩl (qΩl ) = JΩl (qΩl ) − rΩl qΩl
(3.6)
= ∥ΦΩl (qΩl ) − fΩl ∥2 + ϵΩl ∥qΩl ∥2 − rΩl qΩl ,
where rΩl is a row vector, whose function is to adjust the functional’s gradient. On the finest discretization grid Ω0 , all adj variables choose their fine grid values and rΩ0 = 0, such that JΩ0 (qΩ0 ) = JΩ0 (qΩ0 ). Secondly, the initial error between the forward model and measurement data on the coarse grid is equal to the one on the fine grid l+1
ΦΩl+1 (Il
qΩl ) − fΩl+1 = ΦΩl (qΩl ) − fΩl .
(3.7)
This leads to the recursive expression for fΩl+1 fΩl+1 ←− fΩl − ΦΩl (qΩl ) − ΦΩl+1 (Ill+1 qΩl ) .
]
[
(3.8)
Thirdly, the initial regularization term on the coarse grid is equal to the one on the fine grid
ϵΩl+1 ∥Ill+1 qΩl ∥2 = ϵΩl ∥qΩl ∥2 .
(3.9)
This leads to the recursive expression for ϵΩl+1
ϵΩl+1 ←−
∥qΩl ∥2 ∥Ill+1 qΩl ∥2
ϵΩl .
(3.10)
Finally, the gradients of the coarse and fine adjusted objective functionals are the same at the current values of Ill+1 qΩl and qΩl adj l+1 l ∇ Jadj Ωl+1 (Il qΩl ) = ∇ JΩl (qΩl )Il+1 .
(3.11)
This leads to the recursive expression for rΩl+1 rΩl+1 ←− ∇ JΩl+1 (Ill+1 qΩl ) − ∇ JΩl (qΩl ) − rΩl Ill+1 .
(
)
(3.12)
The multigrid-V method is obtained by simply replacing the fixed grid iteration method on Ωl+1 of this two-grid method with another two-grid method recursively. After initialization of fΩ0 ← f, ϵΩ0 ← ϵ , and rΩ0 ← 0, we can solve the nonlinear inverse problem (1.1) by the multigrid-V method, which moves from the finest grid to the coarsest grid, and then moves back to the finest grid, as described in Algorithm 3.1. Please cite this article as: T. Liu, A multigrid–homotopy method for nonlinear inverse problems, Computers and Mathematics with Applications (2019), https://doi.org/10.1016/j.camwa.2019.09.023.
T. Liu / Computers and Mathematics with Applications xxx (xxxx) xxx
5
Algorithm 3.1. The multigrid algorithm. (1) Initialization. Given q0Ω0 ← q0 , fΩ0 ← f, ϵΩ0 ← ϵ , rΩ0 ← 0. (2) For l = 1, 2, . . . , L, compute q0Ωl ←− Ill−1 q0Ωl−1 ,
]
[
fΩl ←− fΩl−1 − ΦΩl−1 (q0Ωl−1 ) − ΦΩl (q0Ωl ) ,
ϵΩl ←−
∥q0Ωl−1 ∥2 ∥q0Ωl ∥2
ϵΩl−1 ,
)
(
rΩl ←− ∇ JΩl (q0Ωl ) − ∇ JΩl−1 (q0Ωl−1 ) − rΩl−1 Ill−1 , and obtain the adjusted objective functional adj
JΩl (qΩl ) = ∥ΦΩl (qΩl ) − fΩl ∥2 + ϵΩl ∥qΩl ∥2 − rΩl qΩl . †
adj
(3) For l = L, L − 1, . . . , 1, use the common iterative method to compute the solution qΩl of min JΩl (qΩl ) starting from the initial estimate q0Ωl , and carry out the correction
(
†
)
q0Ωl−1 ←− q0Ωl−1 + Ill−1 qΩl − q0Ωl . (4)
Use the common iterative method to compute the solution q† of min JΩ0 (qΩ0 ) starting from the initial estimate q0Ω0 .
From Algorithm 3.1, it can be seen that the gradients of JΩl (l = 0, 1, . . . , L) are only calculated once at the beginning of the multigrid algorithm. The gradients on the coarser grids are less costly, and the most computationally costly part is the gradient on the finest grid. The common iterative method on the finest grid Ω0 also requires the calculation of the gradient of JΩ0 at every iteration step. In fact, the computational cost of this part of the multigrid algorithm, which is part (2) in Algorithm 3.1, is less than that of one iteration of the common iterative method on the finest grid, and it does not influence how many grids are chosen. 4. Multigrid–homotopy method The multigrid inversion method is one special kind of multiresolution techniques, which decomposes the original inverse problem into a sequence of sub-inverse problems which depend on the different grid variables, from the coarsest to the finest. Since the coarse grid problems involve the reduced number of variable dimension, the multigrid inversion method can efficiently compute the solution for a desired fine grid problem. Importantly, the multigrid inversion method can greatly accelerate convergence and reduce computational cost because both the forward and inverse problems are more coarsely discretized on the coarser grids. In [20], Oh et al. gave the following monotone convergence theorem of the multigrid inversion method, where condition (1) states that the monotone convergence of the multigrid inversion method depends on the choice of the fixed grid iteration method. When the fixed grid iteration method is chosen to be the locally convergent common iterative method, the multigrid inversion method will not achieve convergence if the initial estimate q0 is far enough away from the solution q† . To enlarge the convergence region, we make the combination of the multigrid method and the homotopy method, and thereby propose a novel inversion method, named as the multigrid–homotopy method. Theorem 4.1 (Multigrid Monotone Convergence). For 0 ≤ l ≤ L − 1, define the functional ξΩl+1 as
ξΩl+1 (qΩl+1 ) = JΩl+1 (qΩl+1 ) − JΩl (qΩl + Ill+1 (qΩl+1 − Ill+1 qΩl )), where JΩl (·) and JΩl+1 (·) are continuously differentiable. Assume that the following conditions are satisfied. (1) (2) (3) (4)
The fixed grid iteration method is monotone for 0 ≤ l ≤ L.
ξΩl+1 (·) is convex for 0 ≤ l < L.
The adjustment vector rΩl+1 is given by Eq. (3.12) for 0 ≤ l < L. The iteration number on grid Ωl is greater than or equal to 1 for 0 ≤ l ≤ L.
Then, the multigrid inversion method of Algorithm 3.1 is monotone for JΩ0 (·). In the inversion process, we first discretize the inverse problem on the different discretization grids Ωl (l = 0, 1, . . . , L) to obtain the objective functionals JΩl (qΩl ) (l = 0, 1, . . . , L), and then adjust these objective functionals at the values of q0Ω0 = q0 , q0Ωl = Ill−1 q0Ωl−1 (l = 1, . . . , L) in the manner the same as that described in Section 3, thereby leading to the adj
adjusted objective functionals JΩl (qΩl ) (l = 0, 1, . . . , L). Next, we use the homotopy method described in Section 2 to adj
†
minimize the adjusted objective functional JΩL (qΩL ) starting from the initial estimate q0ΩL , and use the obtained result qΩL Please cite this article as: T. Liu, A multigrid–homotopy method for nonlinear inverse problems, Computers and Mathematics with Applications (2019), https://doi.org/10.1016/j.camwa.2019.09.023.
6
T. Liu / Computers and Mathematics with Applications xxx (xxxx) xxx
to correct the finer grid value q0ΩL−1 . On the second coarsest grid ΩL−1 , the corrected result serves as the initial estimate, adj
and the common iterative method can be used to minimize JΩL−1 (qΩL−1 ). Successively fining upward in this way succeeds in finding the solution q† of the original fine grid inverse problem. The corresponding algorithm is summarized as follows: Algorithm 4.1. The multigrid–homotopy algorithm. (1) Initialization. Given q0Ω0 ← q0 , fΩ0 ← f, ϵΩ0 ← ϵ , rΩ0 ← 0. (2) For l = 1, 2, . . . , L, compute q0Ωl ←− Ill−1 q0Ωl−1 ,
]
[
fΩl ←− fΩl−1 − ΦΩl−1 (q0Ωl−1 ) − ΦΩl (q0Ωl ) ,
ϵΩl ←−
∥q0Ωl−1 ∥2 ∥q0Ωl ∥2
ϵΩl−1 ,
)
(
rΩl ←− ∇ JΩl (q0Ωl ) − ∇ JΩl−1 (q0Ωl−1 ) − rΩl−1 Ill−1 , and obtain the adjusted objective functional adj
JΩl (qΩl ) = ∥ΦΩl (qΩl ) − fΩl ∥2 + ϵΩl ∥qΩl ∥2 − rΩl qΩl . †
adj
(3) Use the homotopy method to compute the solution qΩL of min JΩL (qΩL ) starting from the initial estimate q0ΩL , and carry out the correction
(
†
)
q0ΩL−1 ←− q0ΩL−1 + ILL−1 qΩL − q0ΩL . †
adj
(4) For l = L − 1, L − 2, . . . , 1, use the common iterative method to compute the solution qΩl of min JΩl (qΩl ) starting from the initial estimate q0Ωl , and carry out the correction
(
†
)
q0Ωl−1 ←− q0Ωl−1 + Ill−1 qΩl − q0Ωl . (5)
Use the common iterative method to compute the solution q† of min JΩ0 (qΩ0 ) starting from the initial estimate q0Ω0 .
On the coarsest grid, there is no prior information about the initial estimate q0ΩL , so the convergence may not be achieved if directly applying the locally convergent common iterative method. To overcome this difficulty, the globally convergent homotopy method is used on the coarsest grid. On the other grids, the initial estimates are all corrected, so it is better to use the common iterative method. If the homotopy method is applied on all grids, unnecessary computational cost will be introduced. 5. Inversion of a nonlinear convection–diffusion equation The saturation equation within the two-phase porous media flow can be described by the following nonlinear convection–diffusion equation ut +
∂ ∂ ϕ (u) + ψ (u) − ∇ · (q(x, y)N(u)∇ u) = s(x, y, t), ∂x ∂y
t > 0,
(5.1)
where q(x, y) is the permeability at (x, y) in the medium, s(x, y, t) is the piecewise smooth source function, N(u) is the positive nonlinear diffusion function, and ϕ (u), ψ (u) are the nonlinear S-shaped flux functions of Buckley–Leverett type in the x and y directions respectively. For simplicity, we will study the problem in the domain Ω = [0, 1] × [0, 1], with the boundary condition u(x, y, t)|(x,y)∈∂ Ω = 0,
(5.2)
and the initial condition u(x, y, 0) = φ (x, y).
(5.3)
Eqs. (5.1)–(5.3) constitute the forward problem of the nonlinear convection–diffusion equation. But in practice, the permeability q(x, y) is not known. What we know is just some measurement data of the concentration field, e.g., u(xp , yp , t) = f (xp , yp , t),
p = 1, 2, . . . , P ,
0 ≤ t ≤ T.
(5.4)
We should reconstruct the permeability from Eqs. (5.1)–(5.3) and the measurement data f , then Eqs. (5.1)–(5.3) and (5.4) constitute the inverse problem of the nonlinear convection–diffusion equation. Please cite this article as: T. Liu, A multigrid–homotopy method for nonlinear inverse problems, Computers and Mathematics with Applications (2019), https://doi.org/10.1016/j.camwa.2019.09.023.
T. Liu / Computers and Mathematics with Applications xxx (xxxx) xxx
7
Fig. 1. The exact model and inversion results on different grids by MGHT in the first example.
After the finite-difference discretization, Eqs. (5.1)–(5.4) can be written as
⎧ m m−1 ui,j −ui,j −1 m−1 m m m ⎪ ⎪ + ∇ · (ϕ (um ⎪ ∆t i,j ), ψ (ui,j )) − ∇ · (qi,j Ni,j ∇ ui,j ) = si,j , ⎪ ⎪ ⎪ i = 1, 2, . . . , I − 1; j = 1, 2, . . . , J − 1, ⎪ ⎪ ⎨ m i = 0, 1, . . . , I , ui,0 = um i,1 = 0, m m ⎪ u = u = 0 , j = 0, 1, . . . , J , ⎪ 0 , j 1 , j ⎪ ⎪ ⎪ 0 ⎪ u = φ , i = 0 , 1, . . . , I ; j = 0, 1, . . . , J , i , j ⎪ i,j ⎪ ⎩ m uxp ,yp = fxmp ,yp , p = 1, 2, . . . , P ,
(5.5)
m m where um i,j = u(i∆x, j∆y, m∆t), qi,j = q(i∆x, j∆y), si,j = s(i∆x, j∆y, m∆t), φi,j = φ (i∆x, j∆y), fxp ,yp = f (xp , yp , m∆t), ∆t,
∆x, ∆y are respectively the time and spatial step lengths, and I =
1
∆x
,J =
1
∆y
−1 m−1 m m . ∇ · (ϕ (um i,j ), ψ (ui,j )) and ∇ · (qi,j Ni,j ∇ ui,j )
are the discretization expressions of the convection term and the diffusion term, respectively, which have been described in detail in our previous paper [41]. The difference Eq. (5.5) defines a vector-valued function Φ : q → f, where q and f are the vectors composed of qi,j and fxmp ,yp q⊤ = (q1,1 , q1,2 , . . . , q1,J , q2,1 , q2,2 , . . . , q2,J , . . . , qI ,1 , qI ,2 , . . . , qI ,J ), f⊤ = (fx11 ,y1 , fx12 ,y2 , . . . , fx1P ,yP , . . . , fxM1 ,y1 , fxM2 ,y2 , . . . , fxMP ,yP ), where M =
T
∆t
. Let f¯xmp ,yp denote the measurement data and form the vector f¯ in the same sequence as f
f¯⊤ = (f¯x11 ,y1 , f¯x12 ,y2 , . . . , f¯x1P ,yP , . . . , f¯xM1 ,y1 , f¯xM2 ,y2 , . . . , f¯xMP ,yP ), then the inversion of q is reduced to the nonlinear operator equation Φ(q) = f¯.
(5.6)
Please cite this article as: T. Liu, A multigrid–homotopy method for nonlinear inverse problems, Computers and Mathematics with Applications (2019), https://doi.org/10.1016/j.camwa.2019.09.023.
8
T. Liu / Computers and Mathematics with Applications xxx (xxxx) xxx
Fig. 2. The inversion results by HT and MG in the first example.
Fig. 3. ε with respect to the CPU times with different initial estimate values in the first example.
To examine the effectiveness of the proposed multigrid–homotopy method, we perform the numerical simulations for u2 (1−5(1−u2 )) two models on a 2.6 GHz PC with 4 GB RAM in MATLAB 2012b environment under Windows 7. Take ϕ (u) = u2 +(1−u)2 ,
ψ (u) =
u2 , u2 +(1−u)2
N(u) = u2 − u + 1, s(x, y, t) = 0, φ (x, y) = sin(π x) sin(π y), T = 0.1, ∆t = 0.005, and ϵ = 10−8 . By
trial-and-error, we set K = 10 in the homotopy method. The initial estimate q0 is set as four different constant values 1 1.51, 3.89, 6.59, and 9.05. The inverse problem is discretized on three grids Ω0 , Ωl , Ω2 , that is L = 2, and ∆x = ∆y = 24 on the finest grid Ω0 . The measurements points are located at each grid points of Ω2 , and the measurement data are taken from t = 0 to t = 0.1 with time span of ∆t = 0.005. To examine the ability of noise suppression of the methods, Please cite this article as: T. Liu, A multigrid–homotopy method for nonlinear inverse problems, Computers and Mathematics with Applications (2019), https://doi.org/10.1016/j.camwa.2019.09.023.
T. Liu / Computers and Mathematics with Applications xxx (xxxx) xxx
9
Table 1 The CPU run times (seconds) by four different methods. Example number
Initial estimate
MGHT
HT
MG
GN
1
1.51 3.89 6.59 9.05
783.347 874.550 783.936 867.263
1265.565 1409.426 1297.970 1336.795
852.307 No convergence No convergence No convergence
1112.221 No convergence No convergence No convergence
2
1.51 3.89 6.59 9.05
483.956 464.332 452.011 446.015
1144.559 1037.194 984.816 933.361
479.844 471.176 463.266 No convergence
915.021 837.895 757.788 No convergence
Table 2 σ by four different methods. Example number
Initial estimate
MGHT
HT
MG
GN
1
1.51 3.89 6.59 9.05
0.0835 0.0851 0.0863 0.0863
0.0836 0.0851 0.0863 0.0864
0.0851 No convergence No convergence No convergence
0.0852 No convergence No convergence No convergence
2
1.51 3.89 6.59 9.05
0.0338 0.0333 0.0334 0.0339
0.0340 0.0338 0.0341 0.0342
0.0346 0.0346 0.0341 No convergence
0.0342 0.0339 0.0341 No convergence
we add 1% Gaussian noise to the measurement data. The common iterative method in the multigrid–homotopy method and the multigrid method is chosen to be the regularized Gauss–Newton method. The error ε of the iterated inversion result and the exact permeability model is computed by
ε = ∥qiter − qexact ∥, where qiter and qexact are respectively the iterated inversion result and exact permeability model. The relative error σ of the inversion result is computed by
σ =
∥qresu − qexact ∥ , ∥qexact ∥
where qresu is the inversion result. In the first example, the exact permeability model is a stratified medium consisting of three vertical layers, as shown in Fig. 1(a), and we respectively use the multigrid–homotopy method (MGHT), the homotopy method (HT), the multigrid method (MG), and the regularized Gauss–Newton method (GN). The typical inversion results on grids Ω0 , Ω1 , and Ω2 by MGHT are respectively shown in Figs. 1(b), 1(c), and 1(d), and the typical inversion results by HT and MG when the convergence can be achieved are shown in Fig. 2. To compare the performance of MGHT, HT, MG, and GN, the errors ε of the iterated inversion results and the exact permeability model vs. the CPU times with different initial estimate values 1.51, 3.89, 6.59, and 9.05 are plotted in Fig. 3. If the initial estimate is poor (i.e. the cases of 3.89, 6.59, and 9.05), GN and MG fail to converge, while MGHT and HT can succeed to converge. Therefore, there are no black lines of GN and green lines of MG in Figs. 3(b), 3(c), and 3(d). The reason that MG achieves no convergence is that the condition (1) in Theorem 4.1 is not satisfied. It is shown that MGHT has a larger convergence region than MG and GN. It is also shown from Fig. 3 that MGHT requires less computational cost than HT and GN. In the second example, the exact permeability model has two anomalous bodies in a homogeneous medium, as shown in Fig. 4(a), and we respectively use MGHT, HT, MG, and GN. The typical inversion results on grids Ω0 , Ω1 , and Ω2 by MGHT are respectively shown in Figs. 4(b), 4(c), and 4(d), and the typical inversion results by HT and MG when the convergence can be achieved are shown in Fig. 5. In the second example, the errors ε of the iterated inversion results and the exact permeability model vs. the CPU times with different initial estimate values 1.51, 3.89, 6.59, and 9.05 are plotted in Fig. 6. In addition, the CPU run times and the relative errors σ of the inversion results by these four methods in above two examples are respectively listed in Tables 1 and 2. These results further demonstrate that (1) the multigrid–homotopy method is stable, fast, and globally convergent; (2) the multigrid–homotopy method has a larger convergence region than the multigrid method and regularized Gauss–Newton method; (3) the multigrid–homotopy method requires less computational cost than the homotopy method and regularized Gauss–Newton method. Please cite this article as: T. Liu, A multigrid–homotopy method for nonlinear inverse problems, Computers and Mathematics with Applications (2019), https://doi.org/10.1016/j.camwa.2019.09.023.
10
T. Liu / Computers and Mathematics with Applications xxx (xxxx) xxx Table 3 σ by MGHT for several values of ϵ with the initial estimate value 1.51. Example number
ϵ = 10−5
ϵ = 10−6
ϵ = 10−7
ϵ = 10−8
1 2
0.1584 0.0408
0.1367 0.0370
0.1111 0.0355
0.0835 0.0338
Fig. 4. The exact model and inversion results on different grids by MGHT in the second example.
Remark 5.1. From Eq. (3.8), it can be noted that fΩl+1 and fΩl are of the same length. To do this, the measurements points are located at the grid points of the coarsest grid. Therefore, the choice of L is related to the locations of the measurements points. Remark 5.2. The regularization parameter ϵ is determined by trial-and-error, i.e., we test several ϵ values and check which performs best. The relative errors σ of the inversion results by MGHT for several values of ϵ with the initial estimate value 1.51 are listed in Table 3. It can be seen from Table 3 that the regularization parameter ϵ = 10−8 performs best. 6. Conclusions For the general nonlinear inverse problems, this paper proposed a multigrid–homotopy method which combines the advantages of the multigrid method and homotopy method. This method first uses the multigrid idea to decompose the original inverse problem into a sequence of sub-inverse problems which depend on the grid variables and are solved in proper order according to the grid size from the coarsest to the finest, and then carries out the inversion on the coarsest grid by the homotopy technique. As a practical application, this new method has been used successfully to solve the inverse problem of the nonlinear convection–diffusion equation, which can be seen as the saturation equation within the two-phase porous media flow. The numerical results show the effectiveness of the multigrid–homotopy method in the aspects of stability, global convergence, and fast computation. Most important of all, the comparative numerical examples show that the multigrid–homotopy method not only has a larger convergence region than the multigrid Please cite this article as: T. Liu, A multigrid–homotopy method for nonlinear inverse problems, Computers and Mathematics with Applications (2019), https://doi.org/10.1016/j.camwa.2019.09.023.
T. Liu / Computers and Mathematics with Applications xxx (xxxx) xxx
11
Fig. 5. The inversion results by HT and MG in the second example.
Fig. 6. ε with respect to the CPU times with different initial estimate values in the second example.
method and regularized Gauss–Newton method, but also requires less computational cost than the homotopy method and regularized Gauss–Newton method. The proposed multigrid–homotopy method can be well suited for other nonlinear inverse problems. Acknowledgments The author would like to thank Editor-in-Chief Daniele Boffi and anonymous referees for their valuable help. This work was supported by the National Natural Science Foundation of China (11271102), the Natural Science Foundation of He-bei Province of China (A2017501001), the Fundamental Research Funds for the Central Universities, PR China (N182304023). Please cite this article as: T. Liu, A multigrid–homotopy method for nonlinear inverse problems, Computers and Mathematics with Applications (2019), https://doi.org/10.1016/j.camwa.2019.09.023.
12
T. Liu / Computers and Mathematics with Applications xxx (xxxx) xxx
References [1] H.W. Engl, M. Hanke, A. Neubauer, Regularization of Inverse Problems, Kluwer, Dordrecht, 1996. [2] A.N. Tikhonov, V.Y. Arsenin, Solution of Ill-posed Problems, John Wiley and Sons, New York, 1977. [3] H.W. Engl, M. Hanke, A. Neubauer, Convergence rates for tikhonov regularization of nonlinear ill-posed problem, Inverse Problems 5 (1989) 523–540. [4] B. Kaltenbacher, Some Newton-type methods for the regularization of nonlinear ill-posed problems, Inverse Problems 13 (1997) 729–753. [5] A.B. Bakushinskii, The problem of the convergence of the iteratively regularized Gauss-Newton method, Comput. Math. Phys. 32 (1992) 1353–1359. [6] M. Hanke, A regularization Levenberg–Marquardt scheme, with application to inverse groundwater filtration problems, Inverse Problems 13 (1997) 79–95. [7] M. Hanke, A. Neubaner, O. Scherzer, A convergence analysis of the Landweber iteration for nonlinear ill-posed problems, Numer. Math. 72 (1995) 21–27. [8] J. Sogn, S. Takacs, Robust multigrid solvers for the biharmonic problem in isogeometric analysis, Comput. Math. Appl. 77 (2019) 105–124. [9] R. Dai, P. Lin, J. Zhang, An EXCMG accelerated multiscale multigrid computation for 3D Poisson equation, Comput. Math. Appl. 77 (2019) 2051–2060. [10] K.J. Brabazon, M.E. Hubbard, P.K. Jimack, Nonlinear multigrid methods for second order differential operators with nonlinear diffusion coefficient, Comput. Math. Appl. 68 (2014) 1619–1634. [11] S.C. Brenner, H. Li, L. Sung, Multigrid methods for saddle point problems: Oseen system, Comput. Math. Appl. 74 (2017) 2056–2067. [12] J. Zhao, T. Liu, G. Feng, A nonlinear multigrid method for inversion of two-dimensional acoustic wave equation, J. Inverse Ill-Posed Probl. 22 (2014) 429–448. [13] J. Zhao, T. Liu, An adaptive multigrid conjugate gradient method for permeability identifcation of nonlinear diffusion equation, J. Inverse Ill-Posed Probl. 24 (2016) 89–97. [14] T. Liu, An adaptive multigrid conjugate radient method for the inversion of a nonlinear convection–diffusion equation, J. Inverse Ill-Posed Probl. 26 (2018) 623–631. [15] U. Ascher, E. Haber, A multigrid method for distributed parameter estimation problems, ETNA 18 (2003) 1–18. [16] T. Liu, A nonlinear multigrid method for inverse problem in the multiphase porous media flow, Appl. Math. Comput. 320 (2018) 271–281. [17] S.F. McCormick, J.G. Wade, Multigrid solution of a linearized, regularized least-squares problem in electrical impedance tomography, Inverse Problems 9 (1993) 697–713. [18] L. Borcea, A nonlinear multigrid for imaging electrical conductivity and permittivity at low frequency, Inverse Problems 17 (2001) 329–359. [19] J.C. Ye, C.A. Bouman, K.J. Webb, R.P. Millane, Nonlinear multigrid algorithms for Bayesian optical diffusion tomography, IEEE Trans. Image Process. 10 (2001) 909–922. [20] S. Oh, A.B. Milstein, C.A. Bouman, K.J. Webb, A general framework for nonlinear multigrid inversion, IEEE Trans. Image Process. 14 (2005) 125–140. [21] S.J. Liao, Homotopy analysis method: a new analytic method for nonlinear problems, Appl. Math. Mech. 19 (1998) 957–962. [22] J.H. He, Homotopy perturbation technique, Comput. Methods Appl. Mech. Engrg. 178 (1999) 257–262. [23] J.H. He, New interpretation of homotopy perturbation method, Internat. J. Modern Phys. B 20 (2006) 1141–1199. [24] J. Biazar, F. Badpeima, F. Azimi, Application of the homotopy perturbation method to Zakharov-Kuznetsov equations, Comput. Math. Appl. 58 (2009) 2391–2394. [25] A. Ebaid, Remarks on the homotopy perturbation method for the peristaltic flow of Jeffrey fluid with nano-particles in an asymmetric channel, Comput. Math. Appl. 68 (2014) 77–85. [26] V. Marinca, N. Herişanu, Nonlinear dynamic analysis of an electrical machine rotor-bearing system by the optimal homotopy perturbation method, Comput. Math. Appl. 61 (2011) 2019–2024. [27] V. Marinca, N. Herişanu, Application of optimal homotopy asymptotic method for solving nonlinear equations arising in heat transfer, Int. Commun. Heat Mass Transfer 35 (2008) 710–715. [28] V. Marinca, N. Herişanu, Determination of periodic solutions for the motion of a particle on a rotating parabola by means of the optimal homotopy asymptotic method, J. Sound Vib. 329 (2010) 1450–1459. [29] N. Herişanu, V. Marinca, Accurate analytical solutions to oscillators with discontinuities and fractional-power restoring force by means of the optimal homotopy asymptotic method, Comput. Math. Appl. 60 (2010) 1607–1615. [30] V. Marinca, N. Herişanu, On the flow of a Walters-type B’ viscoelastic fluid in a vertical channel with porous wall, Int. J. Heat Mass Transfer 79 (2014) 146–165. [31] M.S. Hashmi, N. Khan, S. Iqbal, Numerical solutions of weakly singular Volterra integral equations using the optimal homotopy asymptotic method, Comput. Math. Appl. 64 (2012) 1567–1574. [32] H.B. Keller, D.J. Perozzi, Fast seismic ray tracing, SIAM J. Appl. Math. 43 (1983) 981–992. [33] D.W. Vasco, Singularity and branching: a path-following formalism for geophysical inverse problems, Geophys. J. Int. 119 (1994) 809–830. [34] M.E. Everett, Homotopy, polynomial equations, gross boundary data, and small Helmholtz systems, J. Comput. Phys. 124 (1996) 431–441. [35] M.D. Jegen, E.M. Everett, A. Schultz, Using homotopy to invert geophysical data, Geophysics 66 (2001) 1749–1760. [36] T. Liu, S.S. Liu, Identification of diffusion parameters in a non-linear convection–diffusion equation using adaptive homotopy perturbation method, Inverse Probl. Sci. Eng. 26 (2018) 464–478. [37] B. Han, H.S. Fu, Z. Li, A homotopy method for the inversion of a two-dimensional acoustic wave equation, Inverse Probl. Sci. Eng. 13 (2005) 411–431. [38] B. Han, H.S. Fu, H. Liu, A homotopy method for well-log constraint waveform inversion, Geophysics 72 (2007) R1–R7. [39] T.K. Nilssen, T. Mannseth, X.C. Tai, Permeability estimation with the augmented Lagrangian method for a nonlinear diffusion equation, Comput. Geosci. 7 (2003) 27–47. [40] L. Cao, B. Han, Convergence analysis of the homotopy perturbation method for solving nonlinear ill-posed operator equations, Comput. Math. Appl. 61 (2011) 2058–2061. [41] T. Liu, A wavelet multiscale method for the inverse problem of a nonlinear convection–diffusion equation, J. Comput. Appl. Math. 330 (2018) 165–176.
Please cite this article as: T. Liu, A multigrid–homotopy method for nonlinear inverse problems, Computers and Mathematics with Applications (2019), https://doi.org/10.1016/j.camwa.2019.09.023.
Contents lists available at ScienceDirect
Computers and Mathematics with Applications journal homepage: www.elsevier.com/locate/camwa
A multigrid–homotopy method for nonlinear inverse problems Tao Liu School of Mathematics and Statistics, Northeastern University at Qinhuangdao, Qinhuangdao 066004, PR China
article
info
a b s t r a c t In the present contribution, we develop a novel method combining the multigrid idea and the homotopy technique for nonlinear inverse problems, in which the forward problems are modeled by some forms of partial differential equations. The method first attempts to use the multigrid method to decompose the original inverse problem into a sequence of sub-inverse problems which depend on the grid variables and are solved in proper order according to the grid size from the coarsest to the finest, and then carries out the inversion on the coarsest grid by the homotopy method. The strategy may give a rapidly and globally convergent method. As a practical application, this method is used to solve the nonlinear inverse problem of a nonlinear convection–diffusion equation, which is the saturation equation within the two-phase porous media flow. We demonstrate the effectiveness and merits of the multigrid–homotopy method on two actual model problems. © 2019 Elsevier Ltd. All rights reserved.
Article history: Received 2 April 2019 Received in revised form 25 September 2019 Accepted 30 September 2019 Available online xxxx Keywords: Multigrid Homotopy Regularization Nonlinear inverse problems Convection–diffusion equation
1. Introduction Generally speaking, we can formulate the nonlinear inverse problem as the nonlinear operator equation Φ(q) = f,
(1.1)
or the nonlinear least-square problem min ∥Φ(q) − f∥2 ,
(1.2)
q
where q is the parameter to be inverted, f is the measurement data, and Φ is the forward model, which can be computed from the appropriate governing partial differential equation. By using a standard procedure of Newton-type iterative method to linearize Eq. (1.1), we can get Φ′ (q)(q − qi ) = − [Φ(qi ) − f] ,
(1.3) ′
which can be used to solve Eq. (1.1). But it is general that Φ is not invertible, so we must regularize it [1]. The most famous regularization approach is the Tikhonov’s regularization approach [2]. Many researchers [3,4] have studied its analysis for the nonlinear ill-posed problems, and Bakushinskii [5] developed an iteratively regularized Gauss–Newton method qi+1 = qi − Φ′ (qi )∗ Φ′ (qi ) + βi I
[
]−1 [
i = 0, 1, 2, . . .
Φ′ (qi )∗ (Φ(qi ) − f) + βi (qi − q0 ) ,
]
(1.4)
E-mail address: [email protected]. https://doi.org/10.1016/j.camwa.2019.09.023 0898-1221/© 2019 Elsevier Ltd. All rights reserved.
Please cite this article as: T. Liu, A multigrid–homotopy method for nonlinear inverse problems, Computers and Mathematics with Applications (2019), https://doi.org/10.1016/j.camwa.2019.09.023.
2
T. Liu / Computers and Mathematics with Applications xxx (xxxx) xxx
Here, superscript ∗ refers to the adjoint of an operator. Introducing Tikhonov’s regularization to Eq. (1.3), we can obtain the Levenberg Marquardt method [6] qi+1 = qi − Φ′ (qi )∗ Φ′ (qi ) + βi I
[
]−1
Φ′ (qi )∗ [Φ(qi ) − f] ,
i = 0, 1, 2, . . .
(1.5)
Furthermore, the Landweber method is in essence another interesting regularization method. For example, Hanke et al. [7] proposed a Landweber iterative scheme to solve Eq. (1.1) qi+1 = qi − Φ′ (qi )∗ [Φ(qi ) − f] ,
i = 0, 1, 2, . . .
(1.6)
To solve inverse problems, such Newton-type methods work by performing the computations using a fixed discrete grid. While tremendous advance has been made in aspect of reducing the computational complexity of the fixed grid methods, calculation cost still deserves great attention. Perhaps more importantly, all these methods are locally convergent. In another word, good initial estimates have to be provided when these methods are used to solve Eq. (1.1). In order to accelerate the convergence rate and widen the convergence region, this paper combines the multigrid idea with the homotopy technique, and thereby proposes a multigrid–homotopy method for nonlinear inverse problems, in which the forward problems are modeled by some forms of partial differential equations. The multigrid method is known to be one of the best methods to solve the forward problems of partial differential equations [8–11]. By computing corrections on a coarser grid level and then prolongating them back to the original problem, it can reduce the computation complexities and remove the smooth error components. In recent years, this method has been extended to the study of inverse problems. Successful studies of the multigrid method include the inverse problems of the acoustic wave equation [12], the nonlinear diffusion equation [13], the nonlinear convection– diffusion equation [14], the three-dimension elliptic equation [15], the multiphase porous media flow [16], the electrical impedance tomography [17,18] and the Bayesian optical diffusion tomography [19,20]. It is shown in these papers that the multigrid method is effective on inverse problems. The homotopy method, due to its global convergence property, has been extensively used in the field of nonlinear problems. In 1998, Liao [21] employed the basic ideas of the homotopy in topology to propose a homotopy analysis method for nonlinear problems. Different from the homotopy analysis method, the homotopy perturbation method proposed by He [22,23] is actually a new perturbation technique coupled with the homotopy technique, and has been successfully used to solve nonlinear differential equations [24,25]. Based on the He’s homotopy perturbation method, Marinca and Herişanu [26] proposed an analytical technique, called the optimal homotopy perturbation method, and employed it to study a nonlinear dynamic analysis of an electrical machine rotor-bearing nonlinear system. The optimal homotopy asymptotic method established by Marinca and Herişanu [27] combines the features of the homotopy concept with an efficient computational algorithm which provides a simple and rigorous procedure to control the convergence of the solution. A series of papers by Marinca and Herişanu [28–30], Hashmi et al. [31] have proved the effectiveness, generalization and reliability of this method and obtained solutions of currently important application in science and engineering. Recently, the homotopy method has been also used to solve a wide range of inverse problems. Keller and Perozzi [32] may be the first authors who applied homotopy to the inverse problems. Then, Vasco [33] presented the homotopy method for the inverse problem and demonstrated its use for the traveltime tomography. Everett [34] transformed a discrete inverse problem of the two-dimension Helmholtz equation to a polynomial system, and used the homotopy method to solve it. Jegen et al. [35] used the homotopy method to solve the geophysical inversion of electromagnetic data, which depended on the normal equations. Liu et al. [36] used the homotopy method to identify the diffusion parameters in a nonlinear convection–diffusion equation. Han and Fu developed the homotopy method for the inverse problem of the acoustic wave equation [37] and waveform inversion constrained by well-log data [38]. The rest of the paper is organized as follows. In Sections 2 and 3, the homotopy method and the multigrid method are described respectively. Next, in Section 4 the proposed multigrid–homotopy method is presented. In Section 5, we demonstrate the merits and effectiveness of the proposed method by the inversion of a nonlinear convection–diffusion equation, which is the saturation equation within the two-phase porous media flow. Finally, Section 6 summarizes the conclusions. 2. Homotopy method If we have a good prior estimate q˜ for the solution q† of Eq. (1.1), we can introduce the Tikhonov regularization to approximate Eqs. (1.1) or (1.2), that is min ∥Φ(q) − f∥2 + ϵ∥q − q˜ ∥2 ,
[
]
q
(2.1)
where ϵ is the regularization parameter. Then, the regularized Gauss–Newton method (1.4) can be used to solve Eq. (2.1). However, in practical situation, this assumption is not realistic. In general, a good prior estimate q˜ cannot be provided for Eq. (2.1), so min ∥Φ(q) − f∥2 + ϵ∥q∥2 ,
[
q
]
(2.2)
Please cite this article as: T. Liu, A multigrid–homotopy method for nonlinear inverse problems, Computers and Mathematics with Applications (2019), https://doi.org/10.1016/j.camwa.2019.09.023.
T. Liu / Computers and Mathematics with Applications xxx (xxxx) xxx
3
is used to replace Eq. (2.1). Then, the Levenberg Marquardt method (1.5) can be used, and many trials of initial estimates q0 are made to obtain a satisfactory approximation to q† . This is a serious problem especially for the reservoir history matching. In the field of reservoir history matching, one is interested in inferring as many parameters as there are grid cells in the simulator, and this number usually exceeds 105 for a field simulation. A single field-model reservoir simulation may take several hours to complete and many such simulator runs are required to obtain an inversion result with inversion methods such as Levenberg Marquardt method, so many trials of initial estimates will result in the very high computational cost [39]. To overcome this difficulty, we use the homotopy method to solve Eq. (2.2). First, consider the normal equation of the objective functional in Eq. (2.2) Φ′ (q)∗ (Φ(q) − f) + ϵ q = 0,
(2.3)
and then construct its fixed-point homotopy equation H(q, η) = η Φ′ (q)∗ (Φ(q) − f) + ϵ q + (1 − η)(q − q0 ) = 0,
[
]
(2.4)
where η ∈ [0, 1] is a homotopy parameter, q is an optional initial estimate. When η = 0, the solution of H(q, 0) = 0 is the known initial estimate q0 . When η = 1, H(q, 1) = 0 is just Eq. (2.3). In order to solve Eq. (2.3), we divide the interval [0, 1] into 0 = η0 < η1 < · · · < ηK = 1. For η = ηk (k = 1, 2, . . . , K ), use the iterative method to solve Eq. (2.4) in sequence from k = 1 to k = K . As the solution q0 of H(q, η0 ) = 0 is known, it can be used as the initial estimate to the next equation H(q, η1 ) = 0. Assuming that the approximation solution qk−1 of H(q, ηk−1 ) = 0 has already been obtained, a successive linearization method could be used to compute the solution qk of H(q, ηk ) = 0. Suppose that the ith approximation qki to qk is known. To compute the next approximation qki+1 , in equation H(q, ηk ) = 0, we use the linear function Φ′ (qki )(q − qki ) + Φ(qki ) to replace Φ(q), then have 0
[ ] ηk Φ′ (qki )∗ (Φ′ (qki )(q − qki ) + Φ(qki ) − f) + ϵ q + (1 − ηk )(q − q0 ) = 0. The solution of Eq. (2.5) is just the (i + 1)th approximation qki+1
to q
(2.5)
k
[ ]−1 = qki − ηk Φ′ (qki )∗ Φ′ (qki ) + (ϵ + 1 − ηk )I
[ ] · ηk Φ′ (qki )∗ (Φ(qki ) − f) + ϵ qki + (1 − ηk )(qki − q0 ) , qk0
qki+1
k−1
=q
,
k
q =
qkNk +1
,
i = 0, 1, . . . , Nk ,
(2.6)
k = 1, 2, . . . , K .
Eq. (2.6) has a fast convergence rate similar to the regularized Gauss–Newton method, therefore a good approximation to qk can be obtained by only one iteration when ηk −ηk−1 is small enough. In order to save unnecessary computational cost, we can let ηk = Kk (k = 0, 1, . . . , K ) by an isometric division, and Nk ≡ 0. The solution of Eq. (2.4) will not be affected if K is appropriately chosen. Then, Eq. (2.6) can be transformed into a simplified version
) ]−1 k ϵ+1− I K K [ ( ) ] k ′ k ∗ k · Φ (q ) (Φ(qk ) − f) + ϵ qk + 1 − (qk − q0 ) ,
qk+1 = qk −
[
k
K
Φ′ (qk )∗ Φ′ (qk ) +
(
K
(2.7) k = 0, 1, . . . , K − 1.
In fact, the final iteration result qK of Eq. (2.7) is not the solution to H(q, ηK ) = 0 (also Eq. (2.3)), but it is the one iteration result of Eq. (2.6) with the initial estimate qK −1 and ηK −1 = K K−1 . qK is such a good approximation to the solution of the normal Eq. (2.3) that it can be used as the initial estimate q˜ to the regularized Gauss–Newton method (1.4). Therefore, in order to find the solution q† , we need to perform a regularized Gauss–Newton method. Then, Eqs. (2.7) and (1.4) can form a globally convergent homotopy method for the nonlinear inverse problem (1.1). For the convergence analysis of the homotopy method for nonlinear inverse problems see [40]. 3. Multigrid method In most inverse problems, it is important to choose a discrete grid spacing to describe both the forward problem and its subsequent inverse operation. Usually, a fine grid is attractive because it can reduce modeling error and increase the resolution of the solution, but these are improved at the cost of a tremendous increase in calculation amount. To overcome the problem of calculation complexity, we can use the multigrid method to solve the nonlinear inverse problems. From the beginning we introduce not only one grid but a sequence of grids {Ωl } (l = 0, 1, . . . , L). Let Ω0 be the finest discretization grid, and the grid spacing of Ωl (l = 1, . . . , L) can be determined by multiplying the grid spacing of Ω0 by 2l . The nonlinear operator equation obtained by discretizing the inverse problem on grid Ωl is ΦΩl (qΩl ) = fΩl ,
(3.1)
and the corresponding least-square problem with regularization is min JΩl (qΩl ) = ∥ΦΩl (qΩl ) − fΩl ∥2 + ϵΩl ∥qΩl ∥2 . qΩ l
(3.2)
Please cite this article as: T. Liu, A multigrid–homotopy method for nonlinear inverse problems, Computers and Mathematics with Applications (2019), https://doi.org/10.1016/j.camwa.2019.09.023.
4
T. Liu / Computers and Mathematics with Applications xxx (xxxx) xxx
old Let qold Ωl be some given approximation on grid Ωl , and a coarser resolution representation qΩl+1 can be obtained from qold Ωl by the relation l+1 old qold Ωl+1 = Il qΩl ,
(3.3)
where Ill+1 is the restriction operator. We would like to improve the finer grid approximation qold Ωl by first performing the fixed grid iteration method on the coarser grid Ωl+1 , starting with the initial estimate qold Ωl+1 w old qne Ωl+1 ←− Fixed_Grid_Iteration(qΩl+1 , JΩl (·)),
(3.4)
w where qne Ωl+1 is the iteration result, and then using this result to correct the finer grid approximation old new l old w qne Ωl ←− qΩl + Il+1 (qΩl+1 − qΩl+1 ),
(3.5)
where Ill+1 is the prolongation operator. The fixed grid iteration method can be chosen as any one of the above three common iterative methods: regularized Gauss–Newton method (1.4), Levenberg Marquardt method (1.5), and Landweber method (1.6). old w It is ideal that the new approximation qne Ωl should be better than the old approximation qΩl , specifically, w old JΩl (qne Ωl ) ≤ JΩl (qΩl ).
However, this would not be the case unless the objective functionals are consistent. In reality, if we naively choose a group of objective functionals such as JΩl (·) in Eq. (3.2), the approximation after correction will easily be far away from the optimal solution. To eliminate the problem of inconsistent objective functionals, we use the conditions introduced in [20], which guarantee that the multigrid inversion method is monotonely convergent. Firstly, the adjusted objective functionals are defined as follows adj
JΩl (qΩl ) = JΩl (qΩl ) − rΩl qΩl
(3.6)
= ∥ΦΩl (qΩl ) − fΩl ∥2 + ϵΩl ∥qΩl ∥2 − rΩl qΩl ,
where rΩl is a row vector, whose function is to adjust the functional’s gradient. On the finest discretization grid Ω0 , all adj variables choose their fine grid values and rΩ0 = 0, such that JΩ0 (qΩ0 ) = JΩ0 (qΩ0 ). Secondly, the initial error between the forward model and measurement data on the coarse grid is equal to the one on the fine grid l+1
ΦΩl+1 (Il
qΩl ) − fΩl+1 = ΦΩl (qΩl ) − fΩl .
(3.7)
This leads to the recursive expression for fΩl+1 fΩl+1 ←− fΩl − ΦΩl (qΩl ) − ΦΩl+1 (Ill+1 qΩl ) .
]
[
(3.8)
Thirdly, the initial regularization term on the coarse grid is equal to the one on the fine grid
ϵΩl+1 ∥Ill+1 qΩl ∥2 = ϵΩl ∥qΩl ∥2 .
(3.9)
This leads to the recursive expression for ϵΩl+1
ϵΩl+1 ←−
∥qΩl ∥2 ∥Ill+1 qΩl ∥2
ϵΩl .
(3.10)
Finally, the gradients of the coarse and fine adjusted objective functionals are the same at the current values of Ill+1 qΩl and qΩl adj l+1 l ∇ Jadj Ωl+1 (Il qΩl ) = ∇ JΩl (qΩl )Il+1 .
(3.11)
This leads to the recursive expression for rΩl+1 rΩl+1 ←− ∇ JΩl+1 (Ill+1 qΩl ) − ∇ JΩl (qΩl ) − rΩl Ill+1 .
(
)
(3.12)
The multigrid-V method is obtained by simply replacing the fixed grid iteration method on Ωl+1 of this two-grid method with another two-grid method recursively. After initialization of fΩ0 ← f, ϵΩ0 ← ϵ , and rΩ0 ← 0, we can solve the nonlinear inverse problem (1.1) by the multigrid-V method, which moves from the finest grid to the coarsest grid, and then moves back to the finest grid, as described in Algorithm 3.1. Please cite this article as: T. Liu, A multigrid–homotopy method for nonlinear inverse problems, Computers and Mathematics with Applications (2019), https://doi.org/10.1016/j.camwa.2019.09.023.
T. Liu / Computers and Mathematics with Applications xxx (xxxx) xxx
5
Algorithm 3.1. The multigrid algorithm. (1) Initialization. Given q0Ω0 ← q0 , fΩ0 ← f, ϵΩ0 ← ϵ , rΩ0 ← 0. (2) For l = 1, 2, . . . , L, compute q0Ωl ←− Ill−1 q0Ωl−1 ,
]
[
fΩl ←− fΩl−1 − ΦΩl−1 (q0Ωl−1 ) − ΦΩl (q0Ωl ) ,
ϵΩl ←−
∥q0Ωl−1 ∥2 ∥q0Ωl ∥2
ϵΩl−1 ,
)
(
rΩl ←− ∇ JΩl (q0Ωl ) − ∇ JΩl−1 (q0Ωl−1 ) − rΩl−1 Ill−1 , and obtain the adjusted objective functional adj
JΩl (qΩl ) = ∥ΦΩl (qΩl ) − fΩl ∥2 + ϵΩl ∥qΩl ∥2 − rΩl qΩl . †
adj
(3) For l = L, L − 1, . . . , 1, use the common iterative method to compute the solution qΩl of min JΩl (qΩl ) starting from the initial estimate q0Ωl , and carry out the correction
(
†
)
q0Ωl−1 ←− q0Ωl−1 + Ill−1 qΩl − q0Ωl . (4)
Use the common iterative method to compute the solution q† of min JΩ0 (qΩ0 ) starting from the initial estimate q0Ω0 .
From Algorithm 3.1, it can be seen that the gradients of JΩl (l = 0, 1, . . . , L) are only calculated once at the beginning of the multigrid algorithm. The gradients on the coarser grids are less costly, and the most computationally costly part is the gradient on the finest grid. The common iterative method on the finest grid Ω0 also requires the calculation of the gradient of JΩ0 at every iteration step. In fact, the computational cost of this part of the multigrid algorithm, which is part (2) in Algorithm 3.1, is less than that of one iteration of the common iterative method on the finest grid, and it does not influence how many grids are chosen. 4. Multigrid–homotopy method The multigrid inversion method is one special kind of multiresolution techniques, which decomposes the original inverse problem into a sequence of sub-inverse problems which depend on the different grid variables, from the coarsest to the finest. Since the coarse grid problems involve the reduced number of variable dimension, the multigrid inversion method can efficiently compute the solution for a desired fine grid problem. Importantly, the multigrid inversion method can greatly accelerate convergence and reduce computational cost because both the forward and inverse problems are more coarsely discretized on the coarser grids. In [20], Oh et al. gave the following monotone convergence theorem of the multigrid inversion method, where condition (1) states that the monotone convergence of the multigrid inversion method depends on the choice of the fixed grid iteration method. When the fixed grid iteration method is chosen to be the locally convergent common iterative method, the multigrid inversion method will not achieve convergence if the initial estimate q0 is far enough away from the solution q† . To enlarge the convergence region, we make the combination of the multigrid method and the homotopy method, and thereby propose a novel inversion method, named as the multigrid–homotopy method. Theorem 4.1 (Multigrid Monotone Convergence). For 0 ≤ l ≤ L − 1, define the functional ξΩl+1 as
ξΩl+1 (qΩl+1 ) = JΩl+1 (qΩl+1 ) − JΩl (qΩl + Ill+1 (qΩl+1 − Ill+1 qΩl )), where JΩl (·) and JΩl+1 (·) are continuously differentiable. Assume that the following conditions are satisfied. (1) (2) (3) (4)
The fixed grid iteration method is monotone for 0 ≤ l ≤ L.
ξΩl+1 (·) is convex for 0 ≤ l < L.
The adjustment vector rΩl+1 is given by Eq. (3.12) for 0 ≤ l < L. The iteration number on grid Ωl is greater than or equal to 1 for 0 ≤ l ≤ L.
Then, the multigrid inversion method of Algorithm 3.1 is monotone for JΩ0 (·). In the inversion process, we first discretize the inverse problem on the different discretization grids Ωl (l = 0, 1, . . . , L) to obtain the objective functionals JΩl (qΩl ) (l = 0, 1, . . . , L), and then adjust these objective functionals at the values of q0Ω0 = q0 , q0Ωl = Ill−1 q0Ωl−1 (l = 1, . . . , L) in the manner the same as that described in Section 3, thereby leading to the adj
adjusted objective functionals JΩl (qΩl ) (l = 0, 1, . . . , L). Next, we use the homotopy method described in Section 2 to adj
†
minimize the adjusted objective functional JΩL (qΩL ) starting from the initial estimate q0ΩL , and use the obtained result qΩL Please cite this article as: T. Liu, A multigrid–homotopy method for nonlinear inverse problems, Computers and Mathematics with Applications (2019), https://doi.org/10.1016/j.camwa.2019.09.023.
6
T. Liu / Computers and Mathematics with Applications xxx (xxxx) xxx
to correct the finer grid value q0ΩL−1 . On the second coarsest grid ΩL−1 , the corrected result serves as the initial estimate, adj
and the common iterative method can be used to minimize JΩL−1 (qΩL−1 ). Successively fining upward in this way succeeds in finding the solution q† of the original fine grid inverse problem. The corresponding algorithm is summarized as follows: Algorithm 4.1. The multigrid–homotopy algorithm. (1) Initialization. Given q0Ω0 ← q0 , fΩ0 ← f, ϵΩ0 ← ϵ , rΩ0 ← 0. (2) For l = 1, 2, . . . , L, compute q0Ωl ←− Ill−1 q0Ωl−1 ,
]
[
fΩl ←− fΩl−1 − ΦΩl−1 (q0Ωl−1 ) − ΦΩl (q0Ωl ) ,
ϵΩl ←−
∥q0Ωl−1 ∥2 ∥q0Ωl ∥2
ϵΩl−1 ,
)
(
rΩl ←− ∇ JΩl (q0Ωl ) − ∇ JΩl−1 (q0Ωl−1 ) − rΩl−1 Ill−1 , and obtain the adjusted objective functional adj
JΩl (qΩl ) = ∥ΦΩl (qΩl ) − fΩl ∥2 + ϵΩl ∥qΩl ∥2 − rΩl qΩl . †
adj
(3) Use the homotopy method to compute the solution qΩL of min JΩL (qΩL ) starting from the initial estimate q0ΩL , and carry out the correction
(
†
)
q0ΩL−1 ←− q0ΩL−1 + ILL−1 qΩL − q0ΩL . †
adj
(4) For l = L − 1, L − 2, . . . , 1, use the common iterative method to compute the solution qΩl of min JΩl (qΩl ) starting from the initial estimate q0Ωl , and carry out the correction
(
†
)
q0Ωl−1 ←− q0Ωl−1 + Ill−1 qΩl − q0Ωl . (5)
Use the common iterative method to compute the solution q† of min JΩ0 (qΩ0 ) starting from the initial estimate q0Ω0 .
On the coarsest grid, there is no prior information about the initial estimate q0ΩL , so the convergence may not be achieved if directly applying the locally convergent common iterative method. To overcome this difficulty, the globally convergent homotopy method is used on the coarsest grid. On the other grids, the initial estimates are all corrected, so it is better to use the common iterative method. If the homotopy method is applied on all grids, unnecessary computational cost will be introduced. 5. Inversion of a nonlinear convection–diffusion equation The saturation equation within the two-phase porous media flow can be described by the following nonlinear convection–diffusion equation ut +
∂ ∂ ϕ (u) + ψ (u) − ∇ · (q(x, y)N(u)∇ u) = s(x, y, t), ∂x ∂y
t > 0,
(5.1)
where q(x, y) is the permeability at (x, y) in the medium, s(x, y, t) is the piecewise smooth source function, N(u) is the positive nonlinear diffusion function, and ϕ (u), ψ (u) are the nonlinear S-shaped flux functions of Buckley–Leverett type in the x and y directions respectively. For simplicity, we will study the problem in the domain Ω = [0, 1] × [0, 1], with the boundary condition u(x, y, t)|(x,y)∈∂ Ω = 0,
(5.2)
and the initial condition u(x, y, 0) = φ (x, y).
(5.3)
Eqs. (5.1)–(5.3) constitute the forward problem of the nonlinear convection–diffusion equation. But in practice, the permeability q(x, y) is not known. What we know is just some measurement data of the concentration field, e.g., u(xp , yp , t) = f (xp , yp , t),
p = 1, 2, . . . , P ,
0 ≤ t ≤ T.
(5.4)
We should reconstruct the permeability from Eqs. (5.1)–(5.3) and the measurement data f , then Eqs. (5.1)–(5.3) and (5.4) constitute the inverse problem of the nonlinear convection–diffusion equation. Please cite this article as: T. Liu, A multigrid–homotopy method for nonlinear inverse problems, Computers and Mathematics with Applications (2019), https://doi.org/10.1016/j.camwa.2019.09.023.
T. Liu / Computers and Mathematics with Applications xxx (xxxx) xxx
7
Fig. 1. The exact model and inversion results on different grids by MGHT in the first example.
After the finite-difference discretization, Eqs. (5.1)–(5.4) can be written as
⎧ m m−1 ui,j −ui,j −1 m−1 m m m ⎪ ⎪ + ∇ · (ϕ (um ⎪ ∆t i,j ), ψ (ui,j )) − ∇ · (qi,j Ni,j ∇ ui,j ) = si,j , ⎪ ⎪ ⎪ i = 1, 2, . . . , I − 1; j = 1, 2, . . . , J − 1, ⎪ ⎪ ⎨ m i = 0, 1, . . . , I , ui,0 = um i,1 = 0, m m ⎪ u = u = 0 , j = 0, 1, . . . , J , ⎪ 0 , j 1 , j ⎪ ⎪ ⎪ 0 ⎪ u = φ , i = 0 , 1, . . . , I ; j = 0, 1, . . . , J , i , j ⎪ i,j ⎪ ⎩ m uxp ,yp = fxmp ,yp , p = 1, 2, . . . , P ,
(5.5)
m m where um i,j = u(i∆x, j∆y, m∆t), qi,j = q(i∆x, j∆y), si,j = s(i∆x, j∆y, m∆t), φi,j = φ (i∆x, j∆y), fxp ,yp = f (xp , yp , m∆t), ∆t,
∆x, ∆y are respectively the time and spatial step lengths, and I =
1
∆x
,J =
1
∆y
−1 m−1 m m . ∇ · (ϕ (um i,j ), ψ (ui,j )) and ∇ · (qi,j Ni,j ∇ ui,j )
are the discretization expressions of the convection term and the diffusion term, respectively, which have been described in detail in our previous paper [41]. The difference Eq. (5.5) defines a vector-valued function Φ : q → f, where q and f are the vectors composed of qi,j and fxmp ,yp q⊤ = (q1,1 , q1,2 , . . . , q1,J , q2,1 , q2,2 , . . . , q2,J , . . . , qI ,1 , qI ,2 , . . . , qI ,J ), f⊤ = (fx11 ,y1 , fx12 ,y2 , . . . , fx1P ,yP , . . . , fxM1 ,y1 , fxM2 ,y2 , . . . , fxMP ,yP ), where M =
T
∆t
. Let f¯xmp ,yp denote the measurement data and form the vector f¯ in the same sequence as f
f¯⊤ = (f¯x11 ,y1 , f¯x12 ,y2 , . . . , f¯x1P ,yP , . . . , f¯xM1 ,y1 , f¯xM2 ,y2 , . . . , f¯xMP ,yP ), then the inversion of q is reduced to the nonlinear operator equation Φ(q) = f¯.
(5.6)
Please cite this article as: T. Liu, A multigrid–homotopy method for nonlinear inverse problems, Computers and Mathematics with Applications (2019), https://doi.org/10.1016/j.camwa.2019.09.023.
8
T. Liu / Computers and Mathematics with Applications xxx (xxxx) xxx
Fig. 2. The inversion results by HT and MG in the first example.
Fig. 3. ε with respect to the CPU times with different initial estimate values in the first example.
To examine the effectiveness of the proposed multigrid–homotopy method, we perform the numerical simulations for u2 (1−5(1−u2 )) two models on a 2.6 GHz PC with 4 GB RAM in MATLAB 2012b environment under Windows 7. Take ϕ (u) = u2 +(1−u)2 ,
ψ (u) =
u2 , u2 +(1−u)2
N(u) = u2 − u + 1, s(x, y, t) = 0, φ (x, y) = sin(π x) sin(π y), T = 0.1, ∆t = 0.005, and ϵ = 10−8 . By
trial-and-error, we set K = 10 in the homotopy method. The initial estimate q0 is set as four different constant values 1 1.51, 3.89, 6.59, and 9.05. The inverse problem is discretized on three grids Ω0 , Ωl , Ω2 , that is L = 2, and ∆x = ∆y = 24 on the finest grid Ω0 . The measurements points are located at each grid points of Ω2 , and the measurement data are taken from t = 0 to t = 0.1 with time span of ∆t = 0.005. To examine the ability of noise suppression of the methods, Please cite this article as: T. Liu, A multigrid–homotopy method for nonlinear inverse problems, Computers and Mathematics with Applications (2019), https://doi.org/10.1016/j.camwa.2019.09.023.
T. Liu / Computers and Mathematics with Applications xxx (xxxx) xxx
9
Table 1 The CPU run times (seconds) by four different methods. Example number
Initial estimate
MGHT
HT
MG
GN
1
1.51 3.89 6.59 9.05
783.347 874.550 783.936 867.263
1265.565 1409.426 1297.970 1336.795
852.307 No convergence No convergence No convergence
1112.221 No convergence No convergence No convergence
2
1.51 3.89 6.59 9.05
483.956 464.332 452.011 446.015
1144.559 1037.194 984.816 933.361
479.844 471.176 463.266 No convergence
915.021 837.895 757.788 No convergence
Table 2 σ by four different methods. Example number
Initial estimate
MGHT
HT
MG
GN
1
1.51 3.89 6.59 9.05
0.0835 0.0851 0.0863 0.0863
0.0836 0.0851 0.0863 0.0864
0.0851 No convergence No convergence No convergence
0.0852 No convergence No convergence No convergence
2
1.51 3.89 6.59 9.05
0.0338 0.0333 0.0334 0.0339
0.0340 0.0338 0.0341 0.0342
0.0346 0.0346 0.0341 No convergence
0.0342 0.0339 0.0341 No convergence
we add 1% Gaussian noise to the measurement data. The common iterative method in the multigrid–homotopy method and the multigrid method is chosen to be the regularized Gauss–Newton method. The error ε of the iterated inversion result and the exact permeability model is computed by
ε = ∥qiter − qexact ∥, where qiter and qexact are respectively the iterated inversion result and exact permeability model. The relative error σ of the inversion result is computed by
σ =
∥qresu − qexact ∥ , ∥qexact ∥
where qresu is the inversion result. In the first example, the exact permeability model is a stratified medium consisting of three vertical layers, as shown in Fig. 1(a), and we respectively use the multigrid–homotopy method (MGHT), the homotopy method (HT), the multigrid method (MG), and the regularized Gauss–Newton method (GN). The typical inversion results on grids Ω0 , Ω1 , and Ω2 by MGHT are respectively shown in Figs. 1(b), 1(c), and 1(d), and the typical inversion results by HT and MG when the convergence can be achieved are shown in Fig. 2. To compare the performance of MGHT, HT, MG, and GN, the errors ε of the iterated inversion results and the exact permeability model vs. the CPU times with different initial estimate values 1.51, 3.89, 6.59, and 9.05 are plotted in Fig. 3. If the initial estimate is poor (i.e. the cases of 3.89, 6.59, and 9.05), GN and MG fail to converge, while MGHT and HT can succeed to converge. Therefore, there are no black lines of GN and green lines of MG in Figs. 3(b), 3(c), and 3(d). The reason that MG achieves no convergence is that the condition (1) in Theorem 4.1 is not satisfied. It is shown that MGHT has a larger convergence region than MG and GN. It is also shown from Fig. 3 that MGHT requires less computational cost than HT and GN. In the second example, the exact permeability model has two anomalous bodies in a homogeneous medium, as shown in Fig. 4(a), and we respectively use MGHT, HT, MG, and GN. The typical inversion results on grids Ω0 , Ω1 , and Ω2 by MGHT are respectively shown in Figs. 4(b), 4(c), and 4(d), and the typical inversion results by HT and MG when the convergence can be achieved are shown in Fig. 5. In the second example, the errors ε of the iterated inversion results and the exact permeability model vs. the CPU times with different initial estimate values 1.51, 3.89, 6.59, and 9.05 are plotted in Fig. 6. In addition, the CPU run times and the relative errors σ of the inversion results by these four methods in above two examples are respectively listed in Tables 1 and 2. These results further demonstrate that (1) the multigrid–homotopy method is stable, fast, and globally convergent; (2) the multigrid–homotopy method has a larger convergence region than the multigrid method and regularized Gauss–Newton method; (3) the multigrid–homotopy method requires less computational cost than the homotopy method and regularized Gauss–Newton method. Please cite this article as: T. Liu, A multigrid–homotopy method for nonlinear inverse problems, Computers and Mathematics with Applications (2019), https://doi.org/10.1016/j.camwa.2019.09.023.
10
T. Liu / Computers and Mathematics with Applications xxx (xxxx) xxx Table 3 σ by MGHT for several values of ϵ with the initial estimate value 1.51. Example number
ϵ = 10−5
ϵ = 10−6
ϵ = 10−7
ϵ = 10−8
1 2
0.1584 0.0408
0.1367 0.0370
0.1111 0.0355
0.0835 0.0338
Fig. 4. The exact model and inversion results on different grids by MGHT in the second example.
Remark 5.1. From Eq. (3.8), it can be noted that fΩl+1 and fΩl are of the same length. To do this, the measurements points are located at the grid points of the coarsest grid. Therefore, the choice of L is related to the locations of the measurements points. Remark 5.2. The regularization parameter ϵ is determined by trial-and-error, i.e., we test several ϵ values and check which performs best. The relative errors σ of the inversion results by MGHT for several values of ϵ with the initial estimate value 1.51 are listed in Table 3. It can be seen from Table 3 that the regularization parameter ϵ = 10−8 performs best. 6. Conclusions For the general nonlinear inverse problems, this paper proposed a multigrid–homotopy method which combines the advantages of the multigrid method and homotopy method. This method first uses the multigrid idea to decompose the original inverse problem into a sequence of sub-inverse problems which depend on the grid variables and are solved in proper order according to the grid size from the coarsest to the finest, and then carries out the inversion on the coarsest grid by the homotopy technique. As a practical application, this new method has been used successfully to solve the inverse problem of the nonlinear convection–diffusion equation, which can be seen as the saturation equation within the two-phase porous media flow. The numerical results show the effectiveness of the multigrid–homotopy method in the aspects of stability, global convergence, and fast computation. Most important of all, the comparative numerical examples show that the multigrid–homotopy method not only has a larger convergence region than the multigrid Please cite this article as: T. Liu, A multigrid–homotopy method for nonlinear inverse problems, Computers and Mathematics with Applications (2019), https://doi.org/10.1016/j.camwa.2019.09.023.
T. Liu / Computers and Mathematics with Applications xxx (xxxx) xxx
11
Fig. 5. The inversion results by HT and MG in the second example.
Fig. 6. ε with respect to the CPU times with different initial estimate values in the second example.
method and regularized Gauss–Newton method, but also requires less computational cost than the homotopy method and regularized Gauss–Newton method. The proposed multigrid–homotopy method can be well suited for other nonlinear inverse problems. Acknowledgments The author would like to thank Editor-in-Chief Daniele Boffi and anonymous referees for their valuable help. This work was supported by the National Natural Science Foundation of China (11271102), the Natural Science Foundation of He-bei Province of China (A2017501001), the Fundamental Research Funds for the Central Universities, PR China (N182304023). Please cite this article as: T. Liu, A multigrid–homotopy method for nonlinear inverse problems, Computers and Mathematics with Applications (2019), https://doi.org/10.1016/j.camwa.2019.09.023.
12
T. Liu / Computers and Mathematics with Applications xxx (xxxx) xxx
References [1] H.W. Engl, M. Hanke, A. Neubauer, Regularization of Inverse Problems, Kluwer, Dordrecht, 1996. [2] A.N. Tikhonov, V.Y. Arsenin, Solution of Ill-posed Problems, John Wiley and Sons, New York, 1977. [3] H.W. Engl, M. Hanke, A. Neubauer, Convergence rates for tikhonov regularization of nonlinear ill-posed problem, Inverse Problems 5 (1989) 523–540. [4] B. Kaltenbacher, Some Newton-type methods for the regularization of nonlinear ill-posed problems, Inverse Problems 13 (1997) 729–753. [5] A.B. Bakushinskii, The problem of the convergence of the iteratively regularized Gauss-Newton method, Comput. Math. Phys. 32 (1992) 1353–1359. [6] M. Hanke, A regularization Levenberg–Marquardt scheme, with application to inverse groundwater filtration problems, Inverse Problems 13 (1997) 79–95. [7] M. Hanke, A. Neubaner, O. Scherzer, A convergence analysis of the Landweber iteration for nonlinear ill-posed problems, Numer. Math. 72 (1995) 21–27. [8] J. Sogn, S. Takacs, Robust multigrid solvers for the biharmonic problem in isogeometric analysis, Comput. Math. Appl. 77 (2019) 105–124. [9] R. Dai, P. Lin, J. Zhang, An EXCMG accelerated multiscale multigrid computation for 3D Poisson equation, Comput. Math. Appl. 77 (2019) 2051–2060. [10] K.J. Brabazon, M.E. Hubbard, P.K. Jimack, Nonlinear multigrid methods for second order differential operators with nonlinear diffusion coefficient, Comput. Math. Appl. 68 (2014) 1619–1634. [11] S.C. Brenner, H. Li, L. Sung, Multigrid methods for saddle point problems: Oseen system, Comput. Math. Appl. 74 (2017) 2056–2067. [12] J. Zhao, T. Liu, G. Feng, A nonlinear multigrid method for inversion of two-dimensional acoustic wave equation, J. Inverse Ill-Posed Probl. 22 (2014) 429–448. [13] J. Zhao, T. Liu, An adaptive multigrid conjugate gradient method for permeability identifcation of nonlinear diffusion equation, J. Inverse Ill-Posed Probl. 24 (2016) 89–97. [14] T. Liu, An adaptive multigrid conjugate radient method for the inversion of a nonlinear convection–diffusion equation, J. Inverse Ill-Posed Probl. 26 (2018) 623–631. [15] U. Ascher, E. Haber, A multigrid method for distributed parameter estimation problems, ETNA 18 (2003) 1–18. [16] T. Liu, A nonlinear multigrid method for inverse problem in the multiphase porous media flow, Appl. Math. Comput. 320 (2018) 271–281. [17] S.F. McCormick, J.G. Wade, Multigrid solution of a linearized, regularized least-squares problem in electrical impedance tomography, Inverse Problems 9 (1993) 697–713. [18] L. Borcea, A nonlinear multigrid for imaging electrical conductivity and permittivity at low frequency, Inverse Problems 17 (2001) 329–359. [19] J.C. Ye, C.A. Bouman, K.J. Webb, R.P. Millane, Nonlinear multigrid algorithms for Bayesian optical diffusion tomography, IEEE Trans. Image Process. 10 (2001) 909–922. [20] S. Oh, A.B. Milstein, C.A. Bouman, K.J. Webb, A general framework for nonlinear multigrid inversion, IEEE Trans. Image Process. 14 (2005) 125–140. [21] S.J. Liao, Homotopy analysis method: a new analytic method for nonlinear problems, Appl. Math. Mech. 19 (1998) 957–962. [22] J.H. He, Homotopy perturbation technique, Comput. Methods Appl. Mech. Engrg. 178 (1999) 257–262. [23] J.H. He, New interpretation of homotopy perturbation method, Internat. J. Modern Phys. B 20 (2006) 1141–1199. [24] J. Biazar, F. Badpeima, F. Azimi, Application of the homotopy perturbation method to Zakharov-Kuznetsov equations, Comput. Math. Appl. 58 (2009) 2391–2394. [25] A. Ebaid, Remarks on the homotopy perturbation method for the peristaltic flow of Jeffrey fluid with nano-particles in an asymmetric channel, Comput. Math. Appl. 68 (2014) 77–85. [26] V. Marinca, N. Herişanu, Nonlinear dynamic analysis of an electrical machine rotor-bearing system by the optimal homotopy perturbation method, Comput. Math. Appl. 61 (2011) 2019–2024. [27] V. Marinca, N. Herişanu, Application of optimal homotopy asymptotic method for solving nonlinear equations arising in heat transfer, Int. Commun. Heat Mass Transfer 35 (2008) 710–715. [28] V. Marinca, N. Herişanu, Determination of periodic solutions for the motion of a particle on a rotating parabola by means of the optimal homotopy asymptotic method, J. Sound Vib. 329 (2010) 1450–1459. [29] N. Herişanu, V. Marinca, Accurate analytical solutions to oscillators with discontinuities and fractional-power restoring force by means of the optimal homotopy asymptotic method, Comput. Math. Appl. 60 (2010) 1607–1615. [30] V. Marinca, N. Herişanu, On the flow of a Walters-type B’ viscoelastic fluid in a vertical channel with porous wall, Int. J. Heat Mass Transfer 79 (2014) 146–165. [31] M.S. Hashmi, N. Khan, S. Iqbal, Numerical solutions of weakly singular Volterra integral equations using the optimal homotopy asymptotic method, Comput. Math. Appl. 64 (2012) 1567–1574. [32] H.B. Keller, D.J. Perozzi, Fast seismic ray tracing, SIAM J. Appl. Math. 43 (1983) 981–992. [33] D.W. Vasco, Singularity and branching: a path-following formalism for geophysical inverse problems, Geophys. J. Int. 119 (1994) 809–830. [34] M.E. Everett, Homotopy, polynomial equations, gross boundary data, and small Helmholtz systems, J. Comput. Phys. 124 (1996) 431–441. [35] M.D. Jegen, E.M. Everett, A. Schultz, Using homotopy to invert geophysical data, Geophysics 66 (2001) 1749–1760. [36] T. Liu, S.S. Liu, Identification of diffusion parameters in a non-linear convection–diffusion equation using adaptive homotopy perturbation method, Inverse Probl. Sci. Eng. 26 (2018) 464–478. [37] B. Han, H.S. Fu, Z. Li, A homotopy method for the inversion of a two-dimensional acoustic wave equation, Inverse Probl. Sci. Eng. 13 (2005) 411–431. [38] B. Han, H.S. Fu, H. Liu, A homotopy method for well-log constraint waveform inversion, Geophysics 72 (2007) R1–R7. [39] T.K. Nilssen, T. Mannseth, X.C. Tai, Permeability estimation with the augmented Lagrangian method for a nonlinear diffusion equation, Comput. Geosci. 7 (2003) 27–47. [40] L. Cao, B. Han, Convergence analysis of the homotopy perturbation method for solving nonlinear ill-posed operator equations, Comput. Math. Appl. 61 (2011) 2058–2061. [41] T. Liu, A wavelet multiscale method for the inverse problem of a nonlinear convection–diffusion equation, J. Comput. Appl. Math. 330 (2018) 165–176.
Please cite this article as: T. Liu, A multigrid–homotopy method for nonlinear inverse problems, Computers and Mathematics with Applications (2019), https://doi.org/10.1016/j.camwa.2019.09.023.