Two fractional regularization methods for identifying the radiogenic source of the Helium production-di ﬀ usion equation

: The predication of the helium di ﬀ usion concentration as a function of a source term in di ﬀ usion equation is an ill-posed problem. This is called inverse radiogenic source problem. Although some classical regularization methods have been considered for this problem, we propose two new fractional regularization methods for the purpose of reducing the over-smoothing of the classical regularized solution. The corresponding error estimates are proved under the a-priori and the a-posteriori regularization parameter choice rules. Some numerical examples are shown to display the necessarity of the methods.


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 highsensitivity 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: ∂r ] + f (r), t ∈ (0, T ), 0 < r < R, u(r, 0) = 0, 0 < r < R, u(R, t) = 0, t ∈ (0, T ), lim r→0 u(r, t) bounded, t ∈ (0, T ), (1.1) where 0 < a 0 ≤ a(t) ≤ 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 u(r, T ) = g(r), 0 < r < R. (1. 2) 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.
Proof. We introduce a new variable x := k 2 n , x < 1/β, and let It is easy to see that there exists a unique x 0 = z β(z+m) with z = p 4 such that h 4 (x 0 ) = 0. We find that 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]). ·, · and · are the inner product and norm on L 2 ([0, R]) respectively, introduced as follows f, g = R 0 f (r)g(r)dr, and u = (  In order to solve the problem (1.1), we introduce a new function ω(r, t) under the substitution ω(r, t) = ru(r, t).
It follows from the inverse radiogenic source problem (1.1), that ω 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 X (r) + λX(r) = 0. 0 < r < R, X(R) = 0, lim r→0 r −1 X(r) bounded, (3.4) where λ is an unknown constant.
Through calculation, we can get the eigenvalues of (3.4) as follows 5) and the corresponding eigenfunctions are From the orthogonality and completeness of the eigenfunction system {sin( nπr , we get the solution ω(r, t) and the source term r 1 f (r) can be represented as 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 It is easy to verify that the eigenfunctions {ϕ n (r)} ∞ n=1 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 f r (r) = r f (r) and g r (r) = rg(r). (3.14) To get f r , define an operator K : f r → g r , then the inverse source problem can be represented by the following linear operator equation: 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 which is assumed to satisfy Here δ > 0 represents the noise level, and both g r (r) and g δ r (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 c n 2 ≤ k n ≤ c n 2 , as n → ∞, (3.20) where the constant c = R 2 2a 1 π 2 and c = R 2 a 0 π 2 . From (3.20), we note that when n → ∞, the eigenvalue k n → 0 [15, p7]. Fixed size data error can be amplified arbitrarily much by the factors k −1 n . 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., f r p ≤ E, for p > 0, (3.21) 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 .

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 µ > 0 plays the role of regularization parameter, we define α as the fractional parameter. When α = 1, it expresses the classic Tikhonov method. However, for 0 < α < 1 we can see it prevent the effect of oversmoothing and obtain a more accurate numerical results for the discontinuity of solution.

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, 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 µ, we obtain

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 Morozovs discrepancy principle [39].
We use the discrepancy principle in the following form: where 0 < α ≤ 1, τ > 1 is a constant, µ > 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 < τδ < g δ r . Lemma 4.3. If µ 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 µ is chosen by discrepancy principle (4.5), then, (1) If 0 < p < 2α, we have a convergence estimate (4.7) (2) If p ≥ 2α, we have a convergence estimate Proof. By the triangle inequality, we know (4.9) (1) For 0 < p < 2α. We first give the estimate of I 1 , with Lemma 4.3 we have (4.10) Now we estimate I 2 , by (3.19) and (3.21), we can deduce that (4.11) Combining (4.9), (4.10) and (4.11), we obtain (2) For p ≥ 2α. From (4.9), we first give the estimate of I 1 , with Lemma 4.2 we have (4.12) Then we estimate I 2 , by Lemma 2.3 and (3.21), we known (4.13) Combining (4.9), (4.12) and (4.13), we obtain This completes the proof.

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 ≥ 1 plays the role of regularization parameter, 0 < β < 2 k 2 n , α is called the fractional parameter. When α = 1, it is the classic Landweber iterative method.

Convergence analysis of FLRM under an a-priori parameter selection rule
Now we give the main result of this section. where s denotes the largest integer smaller than or equal to s, c 5 are the positive constants depending on β, p, α, and 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

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 τ > 1 is a user-supplied constant independent on δ, m > 0 is regularization parameter which makes (5.4) hold at the first time, K is defined by the operator Eq (3.15).
Obviously, lim m→0 d(m) = ∞ n=1 (g δ r ) 2 1 2 = g δ r (r) . Therefore the conclusions (1)-(4) are obvious. Remark 5.3. We assume that the noisy data g δ r is large enough such that 0 < τδ < g δ r , 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: Proof. From (5.4), and according to Lemma 2.4, we obtain 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 (5.10) Combining (5.7), (5.8) and (5.10), we can get the conclusion.

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: g δ = g(1 + δ · randn(size(g))).
where g is the solution of the forward problem, and randn is the white noise, the magnitude δ 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 δ,µ (r) and f δ,m (r) are regularized solutions of the WTRM and FLRM, respectively. The relative error calculations of the WTRM and FLRM are as follows where · is the L 2 ([0, R]) norm.
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 k n = T 0 e −λ n T τ a(t)dt dτ = 1 − e −aλ n T aλ n = 1 − e −a(nπ) 2 a(nπ) 2 , and the eigenfunctions ϕ n (r) = √ 2 sin(nπr).
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).  Relative error of FLRM (α = 0.8). 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 apriori 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.     (b) Relative error of classic Tikhonov method (α = 1).
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.   (b) Relative error of classic Landweber method (α = 1). 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.

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.