Two fractional regularization methods for identifying the radiogenic source of the Helium production-diffusion equation

1.   Introduction

The study of helium diffusion dynamics has attracted the attention of many scholars in recent years. Compared with other isotope systems, helium has a high yield. In addition, high-precision and high-sensitivity helium analysis is relatively easy, and the (U-Th)/He isotopes dating has reached a relatively low-temperature condition, which is very important for the study of thermochronology. In ^[1], the author gives the helium diffusion and low-temperature thermochronometry of apatite. The effects of long alpha-stopping distance on (U-Th)/He ages are studied in literature ^[2]. For the application of helium isotope as thermochronometer in terrestrial and extrater restrial materials, refer to the paper ^[3]. For more important applications of helium isotopes in physics, please refer to the literature^[4,5,6].

This article will consider the prediction of helium concentration as a function of the spatial variable source term, which is closely related to helium isotope dating and is a low-temperature thermochronometry method. The helium production-diffusion model is as follows:

where $0 < a_{0}\leq a(t)\leq a_{1}(a_{0}$ and $a_{1}$ are constants), $r$ is the dimensional radial variable, $R$ is the radius of the spherical diffusion domain, given the final observation of the helium concentration as follows

The object is to reconstruct the source term $f(r)$.

Inverse source problem is a typical ill-posed problem, which has been studied by many scholars. For the discussion of the existence, uniqueness and stability of solutions, please refer to reference ^{[7,8,9,10,11]}. In reference ^[12,13], Isakov has given some theoretical studies on the inverse source problem, and in ^[14], Isakov explains the general inverse problem of partial differential equations. In ^[15], Bao et al. has used the Tikhonov method and the truncation method with a-priori parameter choice rule to study the problem of the inverse radiogenic source identification problem. In ^[16], Zhang and Yan applied a-posteriori truncation method to the radiogenic source identification for the helium production-diffusion equation, which is closely connected to the helium isotopes dating as one method of the low-temperature thermochronometry. For more literature on numerical methods, see ^{[17,18,19,20,21,22,23,24,25]}.

Regarding the research on the regularization theory of inverse radiogenic source problem, most of the results are discussed in the context of a-priori parameter selection, the numerical results of posterior parameters are more less. The a-priori method is based on the smoothness conditions of the solution, which is convenient for theoretical analysis, but it is difficult to verify. Therefore, in actual calculations, the regularization method of the posterior parameter selection rule is more widely used. In order to better identify the radiation source problem of helium production-diffusion Eq (1.1), we will give two fractional regularization methods, namely weighted fractional Tikhonov regularization method (WTRM) and fractional Landweber regularization method (FLRM). For related research on fractional regularization methods, please refer to the literatures ^{[26,27,28,29,30,31,32]}. Also, for the application of more fractional regularization methods please refer to ^{[33,34,35,36]}.

The outline of the paper is as follows. In Section 2, we give some preliminary results; In order to overcome the difficulty, a novel a-priori bound is introduced. In Section 3, the ill-posedness of inverse radiogenic source promblem $(1.1)$ is also given. In Section 4, we constructed the regularization solution by a WTRM. Moreover, we have performed a detailed convergence analysis and given the a-priori and a-posterior regularization parameters choice rule. In Section 5, we use the FLRM to give the regularization solution of the inverse source problem, and give the regularization parameter choice rule of the provided method and the corresponding error estimate. In Section 6, Numerical examples are given, and the numerical results show that the proposed method is accurate and effective. This article summarizes some general conclusions in Section 7.

2.   Preliminaries

The following lemmas will be used.

Lemma 2.1. Given $0 < \alpha\leq1$, for constants $x > 0$, $\underline{c} = \frac{R^2}{2a_1 \pi^2}$ and $\overline{c} = \frac{R^2}{a_0 \pi^2}$, $a_0$ and $a_1$ are the diffusivity ^[15], we hold the following inequality

where $c_1 = \frac{\underline{c}^\alpha \alpha^{\frac{\alpha}{\alpha+1}}}{\overline{c}^\alpha(\alpha+1)} > 0$ are independent of $\alpha, \overline{c}, \underline{c}$.

Proof. For $0 < \alpha\leq1$, we define the following function:

where $x > 0$, $h_1(x)$ has a unique point $x_0 = \Big(\frac{\overline{c}^{\alpha+1}}{\alpha \mu}\Big)^{\frac{1}{\alpha+1}}\geq 0$ such that $h_1'(x_0) = 0$.

Clearly we have

Lemma 2.2. For constants $x > 0$, $0 < \alpha\leq1$ and $p > 0$, we have

where $c_2 = \overline{c}^{-\frac{p}{2}}\Big (1-\frac{p}{2\alpha+2}\Big)\Big(\frac{p}{2\alpha+2-p}\Big)^{\frac{p}{2\alpha+2}}, c_3 = \underline{c}^{-(\alpha+1)}$.

Proof. Given the following function:

where $x > 0$.

If $0 < p < 2\alpha+2$, then $\lim\limits_{x\rightarrow 0}h_2(x) = \lim\limits_{x\rightarrow \infty}h_2(x) = 0$. Thus there exists a $x_0 = \overline{c}^{\frac{1}{2}}\Big(\frac{2\alpha+2-p}{\mu p}\Big)^{\frac{1}{2\alpha+2}}\geq 0$ which is a global maximizer such that $h_2'(x_0) = 0$, we have

Thus, we have

If $p\geq2\alpha+2$, for $x\geq1$ then we have

Lemma 2.3. For constants $x > 0$ and $0 < \alpha\leq1$, we have

where $c_4 = \frac{\overline{c}(p+2)^{2+p}(2\alpha-p)^{\frac{2\alpha-p}{2\alpha+2}}}{\underline{c}^{1+\frac{p}{2}}(2\alpha+2)}, c_5 = \frac{\overline{c}}{\underline{c}^{\alpha+1}} > 0$.

Proof. Define the following function:

If $0 < p < 2\alpha$, then $\lim\limits_{x\rightarrow 0}h_3(x) = \lim\limits_{x\rightarrow \infty}h_3(x) = 0$. Thus there exists a unique $x_0 = \underline{c}^{\frac{1}{2}}\bigg(\frac{2\alpha-p}{\mu (p+2)}\bigg)^{\frac{1}{2\alpha+2}}\geq 0$ which is a global maximizer such that $h_3'(x_0) = 0$, we have

Therefore, we have

If $p\geq2\alpha$, $x\geq1$, we have

Lemma 2.4.^[37,38] For $0 < \lambda < 1$, $\upsilon > 0, \, m\in \mathbb{N}$, let $r_{m}(\lambda): = (1-\lambda)^{m}$, there holds:

where

Lemma 2.5. For $k_n > 0$ and $\frac{1}{2} < \alpha < 1$, $0 < \beta k_{n}^{2} < 1, m\geq 1$, we have

Proof. We define two functions with $\tau^{2}: = \beta k_{n}^{2}$:

and

Obviously $\psi(\tau) = \beta\phi(\tau)$. These two functions are continuous in $\tau\in (0, 1)$.

For $\frac{1}{2} < \alpha < 1$ and $\tau\in (0, 1)$, using the Lemma 3.3 in ^[28], we have

Therefore,

Lemma 2.6. For $m\geq1, k_{n} > 0, 0 < \beta k_{n}^{2} < 1$, we have

where the constant $c_6 = \Big(\frac{p}{4\beta}\Big)^{\frac{p}{4}}$.

Proof. We introduce a new variable $x: = k_{n}^{2}, x < 1/\beta$, and let

It is easy to see that there exists a unique $x_0 = \frac{z}{\beta(z+m)}$ with $z = \frac{p}{4}$ such that $h_4'(x_0) = 0$. We find that

3.   Problem formulation and Ill-posedness

In this section, we derive an analytical solution for the inverse radiogenic source problem based on the eigenfunction expansion, and give an analysis on the ill-posedness of the inverse source problem $(1.1)$.

Throughout this paper, the Hilbert space of square integrable functions on $[0, R]$ is denoted as $L^{2}([0, R])$. $\langle \cdot, \cdot\rangle$ and $\|\cdot\|$ are the inner product and norm on $L^{2}([0, R])$ respectively, introduced as follows

In order to solve the problem (1.1), we introduce a new function $\omega(r, t)$ under the substitution

It follows from the inverse radiogenic source problem (1.1), that $\omega$ satisfies:

The corresponding final observation of the helium concentration becomes

Applying the method of separation of variables, consider the solution of problem (3.1) of the form

substitute it into the $(3.1)$, we obtain the following Sturm-Liouville problem

where $\lambda$ is an unknown constant.

Through calculation, we can get the eigenvalues of (3.4) as follows

and the corresponding eigenfunctions are

From the orthogonality and completeness of the eigenfunction system $\{sin(\frac{n\pi r}{R})\}_{n = 1}^{\infty}$ in $L^2([0, R])$, we get the solution $\omega(r, t)$ and the source term $r^{1}f(r)$ can be represented as

where

Substituting (3.7) and (3.8) into (3.1), we have

Solving the above initial-value problem yields the solution

Therefore, the solution of (3.1) can be written as the infinite series

Evaluating (3.10) at t = T on both sides and using the final helium concentration (3.2) give

Define

It is easy to verify that the eigenfunctions $\{\varphi_n(r)\}_{n = 1}^{\infty}$ form an orthonormal basis in $L^{2}([0, R])$. Using the eigenfunctions as a basis, formula (3.8) can be rewritten as:

For convenience, we denote

To get $f_{r}$, define an operator $K:f_r\rightarrow g_r$, then the inverse source problem can be represented by the following linear operator equation:

Using (3.13), it holds

with

Therefore, the analytical solution of the inverse source problem is:

Because measurement errors exist in the data function $g_r$, the solution has to be reconstructed from noisy data $g_r^{\delta}$ which is assumed to satisfy

Here $\delta > 0$ represents the noise level, and both $g_r(r)$ and $g_r^{\delta}(r)$ are assumed to be functions in $L^{2}([0, R])$.

To study the ill-posed nature of the inverse problem, it is sufficient to investigate the decay property of the eigenvalues. From ^[15], we can see that the upper and lower bounds of $k_n$ are as follows

where the constant $\underline{c} = \frac{R^2}{2a_1 \pi^2}$ and $\overline{c} = \frac{R^2}{a_0 \pi^2}$.

From (3.20), we note that when $n\rightarrow \infty$, the eigenvalue $k_n\rightarrow 0$[15,p7]. Fixed size data error can be amplified arbitrarily much by the factors $k_n^{-1}$. Therefore, the problem of identifying $f(r)$ is ill-posed. In the following, we will use two fractional regularization methods to solve the inverse radiogenic source promblem.

To obtain the error estimates, it is necessary to assume certain regularity of the exact source function. Here we assume that there exists an a priori estimate for the source function $f_r(r)$, i.e.,

where $E > 0$ is a constant. where the norm is defined in terms of the eigenfunctions

It is easy to check $\|f_{r}\|_{0} = \|f_{r}\|$.

4.   Regularization solution of WTRM and convergence analysis

In this section, we propose a WTRM to solve the ill-posed problem $(1.1)$ and give convergence estimate under the a-priori regularization parameter choice rule.

Then we consider WTRM to solve the ill-posed problem, the regularization solution is

where $\mu > 0$ plays the role of regularization parameter, we define $\alpha$ as the fractional parameter. When $\alpha = 1$, it expresses the classic Tikhonov method. However, for $0 < \alpha < 1$ we can see it prevent the effect of oversmoothing and obtain a more accurate numerical results for the discontinuity of solution.

4.1.   Convergence analysis of WTRM under an a-priori parameter choice rule

Theorem 4.1. Suppose the a priori condition (3.21) and the noise assumption (3.19) hold,

(1) If $0 < p < 2\alpha+2$ and choice $\mu = \Big(\frac{\delta}{E}\Big)^{\frac{2\alpha+2}{p+2}}$, we have a convergence estimate

(2) If $p\geq2\alpha+2$ and choice $\mu = \Big(\frac{\delta}{E}\Big)^{\frac{\alpha+1}{\alpha+2}}$, we have a convergence estimate

Proof. By the triangle inequality, we know

We first give the estimate of $I_1$, with Lemma 2.1 and (3.19), we have

Now we estimate $I_2$, by Lemma 2.2 and a priori bound condition (3.21), we can deduce that

Combining the above two inequality and choose the regularization parameter $\mu$, we obtain

4.2.   Convergence analysis of WTRM under an a-posteriori parameter selection rule

In this section, we give the regularization parameter choice rule of the posterior fractional regularization method. We can also obtain a convergence rate for the regularized solution (4.1) under this parameter choice rule. The most general a posteriori rule is the Morozov's discrepancy principle ^[39].

We use the discrepancy principle in the following form:

where $0 < \alpha\leq1, \tau > 1$ is a constant, $\mu > 0$ is regularization parameter, $K$ is defined by the operator Eq (3.15). According to the following lemma, we know there exists a unique solution for Eq (4.5) if $0 < \tau\delta < \|g_r^\delta\|$.

Lemma 4.2. Let $d(\mu) = \|Kf^{\delta, \mu}_{r}(r)-g_{r}^{\delta}(r)\|$, then we have the following conclusions: (1) $d(\mu)$ is a continuous function; (2) $\lim\limits_{\mu\rightarrow 0}d(\mu) = \|g_r^\delta(r)\|$; (3) $\lim\limits_{\mu\rightarrow \infty}d(\mu) = 0$; (4) $d(\mu)$ is a strictly decreasing function over $(0, \infty)$.

Proof. The proofs are straightforward results by virtue of

Lemma 4.3. If $\mu$ is the solution of Eq (4.5), we also obtain the following inequality:

Proof. From (4.5), and according to Lemma $4.2$, we obtain

According to Lemma $2.3$, we have

This yields

Theorem 4.4. Suppose the a priori condition (3.21) and the noise assumption (3.19) hold, and the regularization parameter $\mu$ is chosen by discrepancy principle (4.5), then,

(1) If $0 < p < 2\alpha$, we have a convergence estimate

(2) If $p\geq2\alpha$, we have a convergence estimate

Proof. By the triangle inequality, we know

(1) For $0 < p < 2\alpha$. We first give the estimate of $I_1$, with Lemma 4.3 we have

Now we estimate $I_2$, by (3.19) and (3.21), we can deduce that

Combining (4.9), (4.10) and (4.11), we obtain

(2) For $p\geq2\alpha$. From (4.9), we first give the estimate of $I_1$, with Lemma 4.2 we have

Then we estimate $I_2$, by Lemma 2.3 and (3.21), we known

Combining (4.9), (4.12) and (4.13), we obtain

This completes the proof.

5.   Regularization solution of FLRM and convergence analysis

In this section, we propose a FLRM to solve the ill-posed problem $(1.1)$ and give convergence estimate under the a-priori regularization parameter choice rule.

We denote regularization solution of FLRM with the noisy data as follows:

where $m\ge 1$ plays the role of regularization parameter, $0 < \beta < \frac{2}{k_n^{2}}$, $\alpha$ is called the fractional parameter. When $\alpha = 1$, it is the classic Landweber iterative method.

5.1.   Convergence analysis of FLRM under an a-priori parameter selection rule

Now we give the main result of this section.

Theorem 5.1. Suppose the a priori condition (3.21) and the noise assumption (3.19) hold, let $m = \lfloor(\frac{E}{\delta})^{\frac{4}{p+2}}\rfloor$, we have the convergence estimate

where $\lfloor s\rfloor$ denotes the largest integer smaller than or equal to $s$, $c_5$ are the positive constants depending on $\beta, p, \alpha,$ and $\underline{c}$.

Proof. By the triangle inequality, we know

As in the estimate of $I_1$, by Lemma 2.5 and (3.19), we have

Now we estimate $J_2$, by Lemma 2.6 and the a-priori bound condition (3.21), we can deduce that

Combining the above two inequalities, we obtain

Choose the regularization parameter $m$ by

then we have the following result

5.2.   Convergence analysis of FLRM under an a-posteriori parameter selection rule

Due to the semi-convergence property of iterative methods for ill-posed problems, we need a reliable stopping rule for detecting the moment from convergence to divergence. In this section, we give the a-posteriori parameter choice rule for the FLRM. We can obtain a convergence rate for the regularized solution (4.1) under this parameter choice rule. The most general a-posteriori rule is the Morozov's discrepancy principle ^[39].

We use the discrepancy principle in the following form:

where $\tau > 1$ is a user-supplied constant independent on $\delta$, $m > 0$ is regularization parameter which makes (5.4) hold at the first time, $K$ is defined by the operator Eq (3.15).

Lemma 5.2. Let $d(m) = \|Kf^{\delta, m}_{r}(r)-g_r^\delta(r)\Big\|$, then we have the following conclusions: (1) $d(m)$ is a continuous function; (2) $\lim\limits_{m\rightarrow 0}d(m) = \|g_r^\delta(r)\|$; (3) $\lim\limits_{m\rightarrow \infty}d(m) = 0$; (4) $d(m)$ is a strictly decreasing function over $(0, \infty)$.

Proof. By (5.1) and (5.4), we have

Obviously, $\lim\limits_{m\to 0}d(m) = \Big(\sum_{n = 1}^\infty (g_r^\delta)^{2}\Big)^{\frac{1}{2}} = \|g_r^\delta(r)\|.$ Therefore the conclusions (1)–(4) are obvious.

Remark 5.3. We assume that the noisy data $\|g_r^\delta\|$ is large enough such that $0 < \tau\delta < \|g_r^\delta\|$, thus according to Lemma 5.2, there exists a unique minimum solution for inequality (5.4).

Lemma 5.4. If $m$ is the solution of Eq (5.4), we can obtain the following inequality:

where

Proof. From (5.4), and according to Lemma $2.4$, we obtain

then

This yields

Theorem 5.5. Suppose the a priori condition (3.21) and the noise assumption (3.19) hold, then we have the convergence estimate

Proof. By the triangle inequality, we know

On the estimate of $J_1$, by Lemma 5.4, we have

Now we estimate $J_2$, by the a priori bound condition (3.21), we can deduce that

where we have used the Hölder inequality. Therefore, by the triangle inequality, we obtain

Combining (5.7), (5.8) and (5.10), we can get the conclusion.

6.   Numerical examples

In this section, three simple numerical examples are presented to show the validity of the two fractional regularization methods. the simulated data aregenerated as:

where $g$ is the solution of the forward problem, and randn is the white noise, the magnitude $\delta$ indicates the noise level of the measured data. In the numerical experiment, we use the spatial discretization number $M = 401$, and we fix $a(t) = 1, T = 1, R = 1$. $f(r)$ is the exact solution, $f^{\delta, \mu}(r)$ and $f^{\delta, m}(r)$ are regularized solutions of the WTRM and FLRM, respectively. The relative error calculations of the WTRM and FLRM are as follows

where $\|\cdot\|$ is the $L^{2}([0, R])$ norm.

In order to obtain the artificial data $g^\delta$, we need to solve the following forward problem:

Since the diffusivity $a(t)$ is taken to be a constant, the specific representation of the eigenvalues, eigenfunctions and the regularization solution could be calculated as follows

and the eigenfunctions

Substituting the eigenfunctions gives a formula to compute the coefficients

Note that the right hand side is essentially the sine transform of the function $rg(r)$, which can be efficiently implemented by using a version of the fast Fourier transform for real functions ^[40]. Based on the explicit expressions for the eigensystems, the regularization solution of WTRM could be calculated as follows

the regularization solution of FLRM could be calculated as follows

In practical computation, we take the first $n = 40$ terms of the sum $(6.1)$ and (6.2).

Example 6.1. Consider a smooth exact solution:

in interval $[0, 1]$.

Figures 1–3 show the numerical results of Example 6.1. Figure 1 shows the solution of forward problem and the unperturbed data function $g(r)$. Figure 2 shows the comparison for the regularization parameter chosen by both the a-priori WTRM and the a-posteriori WTRM with respect to different noise level. Figure 3 shows the comparison for the regularization parameter chosen by both the a-priori FLRM and the a-posteriori FLRM with respect to different noise level. Table 1(a) and 1(b) show the relative error between the exact solution and the approximate solution calculated by the WTRM and the FLRM, respectively.

From the numerical results, the approximate solution calculated by the a-posterior method is better than the a-prior method, and the relative error of the a-posterior calculation is smaller than that of the a-priori parameter choice. As the noise level increases, the numerical performance deteriorates.

Example 6.2. Consider the oscillating source function

in interval $[0, 1]$.

The numerical experiments of Example 6.2 is similar to those of Example 6.1. The result of Example 6.2 is shown in Figures 4–6 and Table 2.

Example 6.3. Consider the nonsmooth source function:

In this example, we compare the numerical results of two fractional regularization methods and classical regularization methods under the same parameter choice. Figure 7 shows the solution of forward problem and the unperturbed data function $g(r)$ of Example 6.3. Figure 8 and Table 3 show the comparison of the numerical results of the weighted fractional Tikhonov regularization method and the classical Tikhonov regularization method. The result shows that the weighted fractional Tikhonov method outperforms the classical Tikhonov method under the same parameters, and the classic Tikhonov regularized solution oversmooths.

Figure 9 and Table 4 show the comparison of the numerical results of the fractional Landweber regularization method and the classical Landweber regularization method. We can see that the fractional Landweber method provides better numerical result than the classical Landweber method under the same iterative steps.

7.   Conclusions

As we have seen, the fractional methods include the classical method by introducing a new fractional parameter. We have proved the error estimates for the fractional regularization methods under the the a-priori parameter choice rule and the Morozov's parameter choice rule. The error estimates are order-optimal. This shows that in theory the fractional regularization methods are not inferior to the classical regularization method. Numerical results are also displays that the fractional methods can overcome the disadvantage of over-smoothing of the classical methods. The future research will be generalize the fractional regularization methods to some fractional differential equations with important application background.

Acknowledgments

This research is partially supported by the Fundamental Research Funds for the Central Universities, No. lzujbky-2020-12. The author is very grateful to the reviewers for their constructive comments and valuable suggestions, which greatly improved the quality of our papers.

Conflict of interest

All authors declare no conflict of interest in this paper.