![MathJax Logo](/templates/jsp/_style2/_tandf/pb2/images/math-jax.gif)
ABSTRACT
In this paper, we consider an inverse source problem for the helium production-diffusion equation. That is to determine a space-dependent source term in the helium production-diffusion equation from a noisy final data. Since the problem is ill-posed, we propose an effective iterative scheme to deal with the problem. More specifically, we use an acceleration technique to speed up a simple iterative scheme. We also give two kinds of convergence rates by using an a priori and a posteriori regularization parameter choice rule, respectively. Numerical examples are presented to illustrate the validity and effectiveness of the proposed method.
1. Introduction
Inverse source problems play an important role in many branches of engineering applications, such as migration of groundwater, identification and control of pollution source and environmental protection. Due to its application background, in recent years, the inverse source problem has received much attention. Studies on the existence, uniqueness and stability can be found in [Citation1–3]. Numerical methods may be found in [Citation4–9] and related references therein.
In this paper, we consider an inverse radiogenic source problem for the helium production-diffusion equation
(1)
(1)
where
(
and
are constants) is the time-dependent diffusion coefficient.
corresponds to the spatial variable-dependent radiogenic source production. The object of this inverse radiogenic source problem is to reconstruct the source term
by using the given data
(2)
(2)
which is the final observation of the helium concentration. Since there may be some measurement errors in the data function
, the source function
has to be reconstructed from noisy data
. This model is closely connected to the (U-Th)/He isotopes dating, which is one method of low-temperature thermochronometry. For a detailed physical background, we refer to the paper[Citation10–13]. Our model concerns the local helium concentration gradients resulting from the ejection of high-energy alpha particles from grain surfaces. It is assumed that the grain is of a spherical diffusion geometry, which is consistent with laboratory measurements of helium diffusion from apatite [Citation13].
The inverse radiogenic source problem is ill-posed [Citation11]. Therefore, in order to obtain a stable solution, some proper regularization methods are necessary. In [Citation11], two regularization methods, the Tikhonov regularization and the spectral cutoff regularization are considered to obtain the regularized solution, where the regularization parameters are selected by a priori parameter choice rule. Several numerical examples show the efficiency and accuracy of the proposed methods. One criticism for the a priori parameter choice rules is that the regularization parameter depends on the priori information about the exact solution which is unavailable in practical problems, hence the parameter can not be explicitly computed. Therefore, researchers are more interested in the a posterior parameter choice rule. In [Citation14], Zhang et al. formulate the spectral cutoff regularization method with the discrepancy principle in a uniform framework and adapt to solve the inverse radiogenic source problem. We note that both the Tikhonov and spectral cutoff regularized approximation solution in [Citation11,Citation14] are deduced in terms of the eigenfunction expansion, which may be not easy to calculate accurately if the diffusion coefficient is time dependent. Hence both the methods are not suitable for dealing with the inverse source problem (Equation1(1)
(1) ) – (Equation2
(2)
(2) ) with the diffusion coefficient depending on time variable.
Recently, a simple iterative scheme [Citation15] was used to solve the inverse source problem (Equation1(1)
(1) )–(Equation2
(2)
(2) ),
(3)
(3)
where the operator K is defined by
according to Equation (Equation1
(1)
(1) ), β is a relaxation factor satisfying
. This iterative process can be easily implemented, in each step, one only needs to solve a forward problem. It should be pointed out that the iterative scheme has been used to deal with sideways parabolic equation [Citation16], backward diffusion problem [Citation17] and inverse source problem [Citation18] of time-fractional diffusion equation. As reported in [Citation17,Citation18], the iterative scheme (Equation3
(3)
(3) ) is a very slow method, i.e. it needs far too many iteration steps to achieve the desired accuracy. Motivated by [Citation19,Citation20], in this paper, we construct the following iterative scheme,
(4)
(4)
with two initial elements
. The sequence
is chosen as
(5)
(5)
Compared with the iterative scheme (Equation3
(3)
(3) ), the iterative scheme (Equation4
(4)
(4) ) does not increase computation cost in each step except for one addition, which can be seen as an accelerated version of the iterative scheme (Equation3
(3)
(3) ). We will show that the iterative scheme (Equation4
(4)
(4) ) is a regularization method under appropriate stopping rules. Our theoretical analysis shows that, compared with the iterative scheme (Equation3
(3)
(3) ), the iterative scheme (Equation4
(4)
(4) ) can achieve the same convergence rate under both the a priori and a posteriori regularization parameter choice rules provided the parameter η in (Equation5
(5)
(5) ) is chosen large enough, but the number of iteration steps is much less than that of the iterative scheme (Equation3
(3)
(3) ). Therefore, it is a more efficient method for solving the inverse source problem (Equation1
(1)
(1) )–(Equation2
(2)
(2) ).
The rest of the paper is organized as follows. In Section 2, the iterative regularization method is proposed to obtain stable solution for the inverse source problem. Under a smoothness source condition on the exact source term, convergence rate estimates between the regularized solution and the exact solution are provided under both the a priori and a posteriori regularization parameter choice rules. In Section 3, three numerical examples are presented to illustrate the performance of the proposed method. Section 4 ends this paper with a conclusion and final remarks.
2. The iterative method and the related convergence rates
In this section, we will analyse the iterative regularization method and derive convergence rate estimates under a priori and a posteriori parameter choice rules. Throughout this paper, we denote by the Hilbert space of Lebesgue measurable function
with weight ρ on
. Denote by
and
the inner product and norm on
, respectively, given by
In order to derive error estimates and convergence rate results for the iterative method (Equation4
(4)
(4) ), we introduce a substitution
, then
satisfies
(6)
(6)
and the final observation of the helium concentration becomes
(7)
(7)
From the above discussion, we can reformulate the original inverse source problem (Equation1
(1)
(1) )-(Equation2
(2)
(2) ): to reconstruct the source
from the given data
, where
. If we define an operator
according to the system (Equation6
(6)
(6) ), then the inverse source problem (Equation6
(6)
(6) )–(Equation7
(7)
(7) ) can be formulated as the following operator equation:
(8)
(8)
Applying the method of separation of variable for (Equation6
(6)
(6) ) (see [Citation11] for details), there holds
with
(9)
(9)
and
where
and
are Bessel functions of fraction orders [Citation21].
Consequently, A is a linear self-adjoint compact operator with eigenvalues and eigenfunctions
. It is easy to check that the eigenfunctions
form an orthonormal basis in
. Due to the fact that the eigenvalues
tends to zero, the inverse problem is ill-posed [Citation22].
Based on the operator Equation (Equation8(8)
(8) ), we introduce the following iterative method,
(10)
(10)
with initial elements
, where β is a relaxation factor satisfying
. According to the iterative scheme (Equation10
(10)
(10) ), we have
where I is the identity operator. By induction, we can conclude that
belongs to the Krylov subspace
Hence for the noisy data
satisfying
(11)
(11)
by using the eigenfunction expansion, the regularized solution can be given by
(12)
(12)
where
is a polynomial of degree k−1. The residual satisfy
(13)
(13)
If we define
with
, based on the iterative scheme (Equation10
(10)
(10) ), straightforward computation shows that the residual polynomials
satisfy a three-term recurrence relation
(14)
(14)
with
. The expression for the residual polynomials
are proved in [Citation19].
Lemma 2.1
[Citation19]
Let satisfy the recurrence relation (Equation14
(14)
(14) ) with
as in (Equation5
(5)
(5) ), then
has the following form
(15)
(15)
with the Gegenbauer polynomials
.
The following two lemmas are important for us to prove the convergence rates.
Lemma 2.2
For polynomials with 0<s<1, then the following estimate hold for all
,
where
depend on
.
Proof.
The proof of this lemma can be found in the proof of Theorem 4 in [Citation19].
Lemma 2.3
For polynomials satisfying
with 0<s<1, then for all
,
Proof.
The proof of this lemma can be found in the proof of Proposition 1 in [Citation19].
Before deriving convergence rates under a priori and a posteriori parameter choice rules, we give the following elementary estimates for eigenvalues of operator A.
Lemma 2.4
For the eigenvalues defined in (Equation9(9)
(9) ), it hold
where
and
.
Proof.
Recall and
are the lower bound and the upper bound of
, simple calculation yields the conclusion.
In order to guarantee a certain convergence rate for , the set of the exact source function of problem (Equation6
(6)
(6) ) has to be restricted to certain source sets. Throughout this paper, we assume that
fulfils the following source condition:
(16)
(16)
where the norm is defined in terms of the eigenfunctions
It is easy to check that
.
2.1. Convergence rate under an a priori parameter choice rule
Before deriving convergence rate under an a priori parameter choice rule, we first give some lemmas about different errors.
Lemma 2.5
Let be the exact source function and fulfil the source condition (Equation16
(16)
(16) ),
is the solution of the iterative scheme (Equation10
(10)
(10) ) and
is the regularized solution given by (Equation12
(12)
(12) ). Then the approximation error and the data error satisfy
(17)
(17)
and
(18)
(18)
Proof.
We first estimate the approximation error. It follows the eigenfunction expansion of and the exact solution
that
(19)
(19)
Applying the source condition (Equation16
(16)
(16) ) and Lemma 2.4, we can obtain
An application of Lemma 2.2 can get the final estimate.
For the data error, according to the definition of the regularized solution (Equation12(12)
(12) ) and Lemma 2.3, we have
which completes the proof.
Remark 2.1
Following Lemma 2.5, we can understand the behaviour of the iteration. While the approximation error goes to 0 as , the data error may diverge to infinity. For a small value of k, the bound of the approximation error and the data error is of the order of δ only, whence the iteration seems to converge in the beginning, after a while the data error dominates the total error, the iteration may diverge. It is, therefore, necessary to terminate the iteration according to a stopping rule.
Theorem 2.6
Let be the exact radiogenic source,
is the regularized approximation given by (Equation12
(12)
(12) ). Assume the noisy data satisfies (Equation11
(11)
(11) ) and the exact solution satisfies the a priori condition (Equation16
(16)
(16) ).
If
and choosing
, then we have a convergence rate estimate
(20)
(20) with
.
If
and choosing
, then we have convergence rate estimate
Proof.
(1) For , by using the triangle inequality and Lemma 2.5, we can obtain
(2) For
, by using the triangle inequality and Lemma 2.5, we can obtain
which completes the proof.
Remark 2.2
In [Citation15], the authors obtained a convergence rate for all
with iterative scheme (Equation3
(3)
(3) ), where the regularization parameter k satisfies
. Theorem 2.6 shows that our method yields the same convergence rate provided that the parameter η is sufficiently large
, but only needs the square root of the number of iterations for the iterative scheme (Equation3
(3)
(3) ). For
, since
, the number of iteration steps of our method is still smaller than that for the iterative scheme (Equation3
(3)
(3) ).
2.2. Convergence rate under an a posteriori parameter choice rule
Since the a priori regularization parameter choice rule depends on the unavailable information about the exact source function, it can not be used in practical problems. Thus, it is of much interest to develop a posteriori parameter choice rule. In this subsection, we use the so-called discrepancy principle [Citation22] to choose the regularization parameter , i.e. choosing
such that
(22)
(22)
for some fixed number
.
The next lemma provides an upper bound for the regularization parameter following the discrepancy principle (Equation22
(22)
(22) ).
Lemma 2.7
Let be the exact source.
and
are the solutions of iterative scheme (Equation10
(10)
(10) ) with data
and
, respectively. Under the a priori condition (Equation16
(16)
(16) ) and noise assumption (Equation11
(11)
(11) ). If
is the smallest iteration number such that (Equation22
(22)
(22) ) holds, then we have the following upper bound.
If
, then there holds
If
, then there holds
Proof.
By the triangle inequality, we have
(23)
(23)
For the first term of the right-hand side of (Equation23
(23)
(23) ), according to the definition of the regularized solution (Equation12
(12)
(12) ) and Lemma 2.2, we can obtain
Now we estimate the second term of the right-hand side of (Equation23
(23)
(23) ). Applying the source condition (Equation16
(16)
(16) ) and Lemma 2.4, we have
Combining the above estimates, we can get
Noting that
, this yields
(24)
(24)
Applying Lemma 2.2, we have
(25)
(25)
Combining (Equation24
(24)
(24) ) and (Equation25
(25)
(25) ), we can obtain the conclusion.
Theorem 2.8
Let be the exact radiogenic source satisfying the a priori condition (Equation16
(16)
(16) ) and
be the regularized approximation given by the iterative scheme (Equation10
(10)
(10) ) with noisy data
satisfying the assumption (Equation11
(11)
(11) ). If the regularization parameter
is chosen according to the discrepancy principle (Equation22
(22)
(22) ), then we can obtain the following convergence rates:
If
, then there holds
where
If
, then there holds
where
Proof.
In virtue of the triangle inequality
(26)
(26)
we are led to estimate the two terms on the right-hand side. For the second term in (Equation26
(26)
(26) ), we derive from the definition of
,
Using the H
lder inequality, Lemma 2.2, Lemma 2.4 and the source condition (Equation16
(16)
(16) ), we can obtain
(27)
(27)
Since
(28)
(28)
combining (Equation27
(27)
(27) ) and (Equation28
(28)
(28) ) yields
Noting that
we have
(29)
(29)
Next, we estimate the data error
. Replacing k by
in the inequality (Equation18
(18)
(18) ) and using Lemma 2.5, we have
(30)
(30)
Using the upper bound for
provided in Lemma 2.7, we can get
(31)
(31)
By combining (Equation29
(29)
(29) ), (Equation31
(31)
(31) ), we can obtain the conclusion.
Remark 2.3
In [Citation15], the authors obtained a convergence rate for all
under the a posteriori parameter choice rule, the upper bound for the regularization parameter is
. Theorem 2.8 shows that our method yields the same convergence rate provide η is chosen large enough (
). For
, the convergence order is
, which is lower than that of iterative scheme (Equation3
(3)
(3) ). However, in both cases, the upper bound of the regularization parameter for our method is smaller than that for the iterative scheme (Equation3
(3)
(3) ), which means that, in general, our method requires fewer iteration steps.
3. Numerical results
In this section, we present some numerical simulations to illustrate the performance of the iterative method (Equation10(10)
(10) ) for the inverse radiogenic source problem. In our tests, unless otherwise specified, we shall take T = 1, R = 1,
. For comparison, we also give the numerical results of the iterative scheme (Equation3
(3)
(3) ). For the noisy data
, the regularized solutions are given by the iterative schemes
(32)
(32)
and
(33)
(33)
Since the a priori regularization parameter choice rule depends on the unknown information about the exact solution, we only give the numerical results under the a posteriori regularization parameter choice rule. According to the substitution
, once we have
, the regularized radiogenic source term
can be obtained by multiplying
. However, due to the singularity of the function
, this may cause some instability at
. Fortunately, this can be avoided in practice. Multiplying both sides of (Equation32
(32)
(32) ) and (Equation33
(33)
(33) ) by
in each step, we can find it is equivalent to solving
(34)
(34)
with
and
(35)
(35)
with two initial guess
, where
is defined by the following system
(36)
(36)
In the process of implementing iterative schemes (Equation34
(34)
(34) ) and (Equation35
(35)
(35) ) (we call the original method and the accelerated method, respectively), we need to solve the forward problem (Equation36
(36)
(36) ) repeatedly. Noticing that Equation (Equation36
(36)
(36) ) is singular at the origin
, in order to overcome the singularity, we will adopt the technique of shifting half mesh away from the pole [Citation23] to discrete the Equation (Equation36
(36)
(36) ). Based on this discretization on the spatial variable, Equation (Equation36
(36)
(36) ) can be solved by the method of lines [Citation24]. In order to avoid inverse crimes, the synthetic data
used in our experiments are generated on a much finer grid (N = 201) than the one used for the inversion (N = 51), where N = 51 is the total number of the spatial discrete points. The noisy data
is generated by adding a random perturbation to the ‘exact’ data, that is
where ϵ indicates the noise level of the measured data. To show the accuracy of numerical solution, we compute the relative error at the kth step given by
where f and
are the exact solution and the regularized approximation solution at the kth step, respectively.
In order to observe the behaviour of iterations, we also compute the residual which is defined by
where
is the final data with the source function
. All the computations are performed using MATLAB 2018a on an Intel-i5 desktop computer.
Example 3.1
We first consider the smooth source function
over the interval
.
We first investigate the influence of the parameter η which appears in (Equation5(5)
(5) ).For fixed noise level
, the relation between relative error
and η is plotted in Figure , which indicate that, in general, the value of η should not be selected too small, this is consistent with our theoretical results. In practical application, the parameter η is always selected as small as possible. Therefore in the following numerical experiments, we always choose
.
In Figure , we present the behaviour of residuals and relative errors with respect to the number of iteration steps using the noisy data with noise level
. We can see that the residuals corresponding to both the two methods decrease as k increase and the relative errors decrease first and then increase after several iteration steps, which indicate that choosing the appropriate number of stop iteration steps is very important for both methods. In the following, we take the regularization parameters based on (Equation22
(22)
(22) ) with
. The comparisons between the exact source function and the reconstructed solution obtained by both the original method (Equation34
(34)
(34) ) and the accelerated method (Equation35
(35)
(35) ) are presented in Figure . From Figure , it can be seen that, both the original method (Equation34
(34)
(34) ) and the accelerated method (Equation35
(35)
(35) ) can give satisfactory reconstruction results for different noise levels, but the required number of iteration steps and the computation time for the accelerated method are much less than that of the original method (see Table ).
Figure 3. Example 1: (a)–(c) The reconstructed source function obtained by the original method; (d)–(f) The reconstructed source function obtained by the accelerated method.
![Figure 3. Example 1: (a)–(c) The reconstructed source function obtained by the original method; (d)–(f) The reconstructed source function obtained by the accelerated method.](/cms/asset/7abb35dc-1178-4f15-89e7-4d4a90c021fa/gipe_a_2097224_f0003_oc.jpg)
Table 1. Numerical comparison on the original method and the accelerated method for Example 3.1.
Example 3.2
In this example, we consider the oscillating source function
over the interval
. The obtained numerical results are shown in Figure and Table . It is observed that, for this oscillating function, more iteration steps are needed to satisfy the discrepancy principle (Equation22
(22)
(22) ). Table further confirms that both the methods perform well although the number of iteration steps and computation time required by the accelerated method (Equation35
(35)
(35) ) are far less than that of the original method (Equation34
(34)
(34) ).
Figure 4. Example 2: (a)–(c) The reconstructed source function obtained by the original method; (d)–(f) the reconstructed source function obtained by the accelerated method.
![Figure 4. Example 2: (a)–(c) The reconstructed source function obtained by the original method; (d)–(f) the reconstructed source function obtained by the accelerated method.](/cms/asset/6150bd19-6bc5-40c7-aaf2-7b23465a91d6/gipe_a_2097224_f0004_oc.jpg)
Table 2. Numerical comparison on the original method and the accelerated method for Example 3.2.
Example 3.3
In this example, we consider the source function
This is a more difficult example since the function is not smooth at the point
. Figure shows the comparisons between the exact source function and the reconstructed solution obtained by both the original method (Equation34
(34)
(34) ) and the accelerated method (Equation35
(35)
(35) ). Further numerical results are shown in Table . These numerical results are very similar to that in the Example 3.1 and 3.2. Particularly, from Figure and Table , we can see that even for this nonsmooth function the accelerated method (Equation35
(35)
(35) ) still works well with fewer iteration steps.
Figure 5. Example 3: (a)–(c) The reconstructed source function obtained by the original method; (d)–(f) the reconstructed source function obtained by the accelerated method.
![Figure 5. Example 3: (a)–(c) The reconstructed source function obtained by the original method; (d)–(f) the reconstructed source function obtained by the accelerated method.](/cms/asset/04f6941c-8407-4f29-908c-c1dc2909e6e8/gipe_a_2097224_f0005_oc.jpg)
Table 3. Numerical comparison on the original method and the accelerated method for Example 3.3.
4. Conclusion
In this work, we consider an inverse radiogenic source problem associated to the helium production-diffusion equation. A fast iterative method has been proposed to solve the inverse source problem. Both the a priori and a posteriori regularization parameter choice rules were investigated for the proposed method. Furthermore, the convergence rates were provided for the regularized solutions generated by the iterative regularization method. Numerical results indicate that the proposed method is efficient, accurate and stable to compute the approximation of the unknown radiogenic source.
We conclude the paper by some remarks. In this paper, although our theoretical analysis is based on the eigenfunction expansion, our numerical experiments do not depend on the eigenvalues and eigenfunctions at all. In each step, we only need to numerically solve a forward problem. So our method can be easily extended to the problem for general geometry, which will be considered in the future.
Disclosure statement
No potential conflict of interest was reported by the author(s).
Additional information
Funding
References
- Choulli M, Yamamoto M. Conditional stability in determining a heat source. J Inverse Ill-posed Probl. 2004;12(3):233–243.
- Hettlich F, Rundell W. Identification of a discontinuous source in the heat equation. Inverse Probl. 2001;17(5):1465.
- Sakamoto K, Yamamoto M. Inverse heat source problem from time distributing overdetermination. Appl Anal. 2009;88(5):735–748.
- Dou FF, Fu CL, Yang F. Identifying an unknown source term in a heat equation. Inverse Probl Sci Eng. 2009;17(7):901–913.
- Hussein MS, Lesnic D, Johansson BT, et al. Identification of a multi-dimensional space-dependent heat source from boundary data. Appl Math Model. 2018;54:202–220.
- Johansson T, Lesnic D. Determination of a spacewise dependent heat source. J Comput Appl Math. 2007;209(1):66–80.
- Xie J, Zou J. Numerical reconstruction of heat fluxes. SIAM J Numer Anal. 2005;43(4):1504–1535.
- Yan L, Fu CL, Dou FF. A computational method for identifying a spacewise-dependent heat source. Int J Numer Method Biomed Eng. 2010;26(5):597–608.
- Yang F, Fu CL. The revised generalized Tikhonov regularization for the inverse time-dependent heat source problem. J Appl Math Comput. 2013;41(1–2):81–98.
- Bao G, Dou Y, Ehlers T, et al. Quantifying tectonic and geomorphic interpretations of thermochronometer data with inverse problem theory. Commun Comput Phys. 2011;9(1):129–146.
- Bao G, Ehlers TA, Li P. Radiogenic source identification for the helium production-diffusion equation. Commun Comput Phys. 2013;14(1):1–20.
- Shuster DL, Farley KA. 4he/3he thermochronometry: theory, practice, and potential complications. Rev Mineral Geochem. 2005;58(1):181–203.
- Wolf RA, Farley KA, Kass DM. Modeling of the temperature sensitivity of the apatite (ucth)/he thermochronometer. Chem Geol. 1998;148(1):105–114.
- Zhang YX, Yan L. The general a posteriori truncation method and its application to radiogenic source identification for the helium production–diffusion equation. Appl Math Model. 2017;43:126–138.
- Geng X, Cheng H, Liu M. Inverse source problem of heat conduction equation with time-dependent diffusivity on a spherical symmetric domain. Inverse Probl Sci Eng. 2021;29(11):1–16.
- Deng Y, Liu Z. Iteration methods on sideways parabolic equations. Inverse Probl. 2009;25(9):Article ID 095004, 14pp.
- Wang J-G, Wei T. An iterative method for backward time-fractional diffusion problem. Numer Methods Partial Differ Equ. 2014;30(6):2029–2041.
- Wang JG, Ran YH. An iterative method for an inverse source problem of time-fractional diffusion equation. Inverse Probl Sci Eng. 2018;26(10):1–13.
- Kindermann S. Optimal-order convergence of Nesterov acceleration for linear ill-posed problems. Inverse Probl. 2021;37(6):Article ID 065002, 21pp.
- Neubauer A. On Nesterov acceleration for Landweber iteration of linear ill-posed problems. J Inverse Ill-posed Probl. 2017;25(3):381–390.
- Watson GN. A treatise on the theory of Bessel functions. Cambridge: Cambridge University Press; 1995.
- Engl HW, Hanke M, Neubauer A. Regularization of inverse problems. Dordrecht: Kluwer Academic Publishers Group; 1996.
- Lai MC, Wang WC. Fast direct solvers for Poisson equation on 2d polar and spherical geometries. Numer Methods Partial Differ Equ: An Int J. 2002;18(1):56–68.
- Schiesser WE, Griffiths GW. A compendium of partial differential equation models: method of lines analysis with matlab. Cambridge: Cambridge University Press; 2009.