Numerical investigations of traveling singular sources problems via moving mesh method

This paper studies the numerical solution of traveling singular sources problems. In such problems, a big challenge is the sources move with different speeds, which are described by some ordinary differential equations. A predictor-corrector algorithm is presented to simulate the position of singular sources. Then a moving mesh method in conjunction with domain decomposition is derived for the underlying PDE. According to the positions of the sources, the whole domain is splitted into several subdomains, where moving mesh equations are solved respectively. On the resulting mesh, the computation of jump $[\dot{u}]$ is avoided and the discretization of the underlying PDE is reduced into only two cases. In addition, the new method has a desired second-order of the spatial convergence. Numerical examples are presented to illustrate the convergence rates and the efficiency of the method. Blow-up phenomenon is also investigated for various motions of the sources.

as the model problem in this paper. Here q > 0 is the number of singular sources. The initial value u 0 (x) is taken to be continuous and compatible with the boundary conditions, i.e. u 0 (x) → 0 as |x| → ∞. The local source functions F i (t, x, u) (i = 0, 1, . . . , q − 1) might be given a priori or can be determined from some additional constraints on the solution. The traveling sources are located at with the mesh composed of the local mesh on each subdomain. Taking the advantages of domain decomposition [21], MMPDEs could be solved efficiently by parallel computing. Moreover, It can be found that the computation of [u] is avoided, thus the discretization scheme for the underlying PDE becomes very simple. In addition, our method has an expected second-order convergence in space. The organization of the paper is as follows. In section 2, we introduce the moving mesh method in conjunction with domain decomposition for the model problem. In section 3, the discretization schemes for the physical problem will be derived in detail. In section 4, several numerical examples are given to demonstrate the numerical efficiency and accuracy of our method. The conclusions are presented in the last section.

Moving mesh method in conjunction with domain decomposition
In the last decade, moving mesh method in conjunction with a Schwarz domain decomposition has been developed by Haynes and his co-workers (see i.e. [22,23] and references therein). And in this section, we will introduce a slight different method, that is, moving mesh method in conjunction with a non-overlapping domain decomposition.
Denote the observed domain by [x l , x r ] and assume it containing all sources, that is, x l < α 0 (t) < α 1 (t) < · · · < α q−1 (t) < x r . Here x l , x r are either constants or variables of t. Then, the observed domain is divided into q + 1 subdomains [α i−1 , α i ] (i = 0, 1, . . . , q) with α −1 = x l , α q = x r , by the q sources respectively. Obviously, the sizes of the subdomains are variables of t too.
Let x and ξ denote physical and computational coordinates, respectively. Without loss of generality we assume the computational domain is [0, 1]. Then an one-to-one coordinate transformation between the observed domain [x l , x r ] and the computational domain [0, 1] is defined by with For a given uniform mesh, ξ j = j N , j = 0, 1, . . . , N , on the computational domain, the corresponding mesh on the observed domain [x l , x r ] is In our method, the coordinate transformation (6) is determined as a piecewise smooth function. On each subdomain [α i−1 , α i ], i = 0, 1, . . . , q, it is the solution of an MMPDE which is derived from the equidistribution principle. In the literature, the following MMPDEs which known as MMPDE4, MMPDE5 and MMPDE6, respectively, are popularly used after they were originally established and analyzed in [14]. Here M = M (x, t) is the monitor function giving some measure of the solution error on the physical domain and τ > 0 is a parameter representing a timescale for adjusting the mesh toward equidistribution. In the asymptotic case t → ∞, the solution of MMPDE4, MMPDE5 and MMPDE6 would satisfy the equidistribution principle, which is stated that [14] ∂ ∂ξ M ∂x ∂ξ = 0.
For more details about MMPDE, one can refer to [14] or the recent book [16]. In this paper, MMPDE6 (9) with the boundary condition is employed as an example to describe our moving mesh strategy in conjunction with domain decomposition. Here j s i is some fixed index satisfying 0 < j s i < N . The resulting mesh, used to solve the model problem on [x l , x r ], satisfies the property that a fixed mesh point is located on each source during the time in consideration, i.e. x j s i ≡ α i (t).
Given the old mesh x n j on the observed domain [x l , x r ] and the corresponding solution on the mesh.   Figure 1 shows the moving mesh strategy in conjunction with domain decomposition. Here the computation of the monitor function will be presented in section 4. And MMPDE6 (9) is solved by the following finite difference scheme in our numerical examples, where ∆t n = t n+1 − t n and M j+ 1 2 = (M j+1 + M j )/2. The new mesh could be obtained very efficiently by parallel computing based on domain decomposition methods [21]. And it is best in the sense of equidistribution on each subdomain. On the other hand, we will found in the next section that the computation of the jump [u] is avoided, hence the discretization scheme for the physical PDE (1) becomes very simple.

Model discretization and final algorithm
In this section, we derive the discretization schemes for the physical model problem (1)−(3) and (4) on the observed domain [x l , x r ] with appropriate boundary conditions. Then present a full algorithm of moving mesh method for the model problem.

Discretization schemes
For an arbitrary function Through the coordinate transformation (6), we can rewrite equation (1) on the computational coordinates asu Since the right-hand side of (13) vanishes when x = α i (t), that is, we conduct the discretization scheme for (13) on the above equation as [19], with each term on the left-hand side of (14) containing the information of jumps when they cross the sources. Physically, the value u of ith source changes smoothly as time evolution, which means the jump of the directional derivative of u(x, t) along the vector (α ′ i (t), 1) is zero [4,19], i.e., Recalling that By using the above equation, we can deduce from (14) that Then we obtain immediately by taking (4) and (5) into (17). Similarly to [19], the discretization scheme for (13) are divided into two cases due to x j s i ≡ α i (t) during time integration. For j = j s i , i = 0, 1, . . . , q −1, the mesh point not located at the source, (13) is discretized by standard center difference for spatial variable and backward difference for temporal variable, that is, where h n j = x n j − x n j−1 . Here x n j , u n j are the mesh and the solution on it at time step t n , respectively. For j = j s i , the mesh point just located at the source, the jump informations should be incorporated into the discretization scheme. For this case, the discretization scheme for (13) reads where ψ n+1 In the above schemes, we need the source position α i (t n+1 ) at time step t n+1 . For the general movement (4), it is computed by the following Crank-Nicolson scheme as in [20]. If ψ i (t, α i (t), u), i = 0, 1, . . . , q − 1, are independent of u, the source position α n+1 i and the speed ψ n+1 i can be calculated in advance before solving the discretization schemes for MMPDE6 (9) and physical PDE (13). Otherwise, the resulting system would be too complicated to be solved. In this case, we decouple the discretization system by a predictor-corrector algorithm. For the predictor step, assume ψ n+1 i = ψ n i and solve (21) to get an approximate variable α * i of α n+1 i . Then substituting ψ n+1 i and α * i into the discretization schemes for (9) and (13) to obtain an approximate solution u * of u n+1 . For the corrector step, compute ψ n+1 i = ψ i (t n+1 , α * i , u * ) and solve (21) to get α n+1 i , then obtain the solution u n+1 at time step t n+1 by the discretization schemes for (9) and (13). To complete the discretization schemes, we require an appropriate condition for u on the boundary of the observed domain [x l , x r ]. For the observed domain is small enough, we employ a third-order local absorbing boundary condition (LABC) proposed in [24] for (1) as in [19,20]. Here s 0 is an user-defined parameter, the plus sign in "±" corresponds to the LABC at the right boundary x r , and the minus sign corresponds to the one at the left boundary x l . Under the map (6), we get the LABC for (13) as followṡ where the plus sign in "±" corresponds to the right boundary, and the minus sign corresponds to the left boundary. According to [19,20], a finite difference scheme for (23) is on the left boundary, and on the right boundary. Here two ghost points x −1 and x N +1 are used. On the other hand, if the observed domain is big enough or else, Dirichlet boundary conditions are employed.

Full algorithm
We close this section with a full algorithm in Figure 2 for the model problem (1)− (3) and (4). Here the choice of the time step ∆t n will be determined in the following concrete examples, and T ol > 0 is set to be 10 −16 .

Numerical examples
In this section, we present some numerical examples to verify the convergence rate and illustrate efficiency of the full algorithm in Figure 2.
Example 1. We consider a nonlinear moving interface problem with the following exact solution for some choice of ω 1 and ω 2 . The interface α 0 (t) is determined by solving the scalar equation so that u(x, t) is continuous across the interface.
The equation (27) has a unique solution on [0, 1] if we take, for example, π < ω 1 , ω 2 < 2π. Then we have the ordinary differential equation for the motion of the interface Prepare the initial values x 0 , u 0 and the terminate time T . Let n = 0, t n = 0.
Moving mesh strategy in conjunction with domain decomposition (see Figure 1).
Solve the discretization schemes for the physical PDE (13) to calculate u * .

Moving mesh strategy in conjunction with
domain decomposition (see Figure 1).
Solve the discretization schemes for the physical PDE (13) to calculate u n+1 . Based on the jump conditions, the source function F 0 (t, x, u) is Same as in [2], we take ω 1 = 5π/4, ω 2 = 7π/4. The observed domain is set by [0, 1], where the initial position of the interface is α 0 (0) = 0.58333. Since we have the exact solution, Dirichlet boundary conditions are employed. In this example, we simply use the uniform time step, i.e. ∆t n ≡ const, and the total number of the time meshes is L. The monitor function for MMPDE6 (9) takes the form where 0 < θ < 1, 0 < ε ≪ 1. This is consistent with the choice in [19,20]. In practice, smoothing the monitor function can improve the accuracy of the numerical solution, and we utilize the smoothing technique proposed in [25]. Here the parameters in MMPDE6 (9) and the monitor function (30) are given by τ = 10 −3 , θ = 0.5, and ε = 10 3 /N 4 . Since backward Euler scheme is used to solve the physical PDE in this paper, the truncation error for time discretization is only first-order. To verify our algorithm has a second-order convergence rate for space, the number of L should be fourfold when N is double in the convergence test. Computational results with different number of N and L at the time T = 0.1 are listed in Table 1, where the errors are defined as Here, u e is the true solution, α * 0 used as the exact interface is the solution of a zero-finding MATLAB function fzero for (27). The numerical solution U is obtained by the algorithm where α 0 (t) and α ′ 0 (t) are exactly calculated. AndŨ ,α 0 represent respectively the numerical solution and interface, solved with the full predictor-corrector algorithm. The ratios in Table 1 are E 2N,4L /E N,L ,Ẽ 2N,4L /Ẽ N,L and E α 2N,4L /E α N,L , respectively. It is shown that our algorithm solves the solution and the interface very well, and has a second-order convergence rate for space, i.e. O(1/N 2 ). Additionally, compared the corresponding results in [3], our algorithm is better than the method proposed in [3].    Figure 3 presents the profiles of the solution in physical variable and computational variable and the evolving mesh from t = 0 to t = 0.1. The number of the mesh is N = 24, with half mesh points on each side of the interface. We can see that we get excellent resolution of the example even with a grid as coarse as N = 24.
The rest examples are from traveling heat sources problems with the solution may be blowup [4,17,19,20]. If not specifically pointed out, the initial value is given by the observed domain is set by [−10, 10] with u(−10, t) = u(10, t) = 0, and the source functions F i (t, x, u) are simply specified by The resulting nonlinear system is solved by Newton iteration with the tolerance tol = 10 −8 . The monitor function for MMPDE6 (9) takes the form where the parameters 0 < θ i < 1, q+1 i=0 θ i = 1, 0 < ε ≪ 1, p > 0 will be determined later. For non-blowup case, the following graded time steps [4,17] t n = n T L 2 , n = 0, 1, . . . , L, are used with [0, T ] the time integration interval and L the number of time meshes. While for blow-up case, the time step ∆t n = t n+1 − t n is chosen to be [4,26] ∆t n = min where ε is same in the monitor function, µ is a small positive constant with µ = 10 −3 in the test.
The profiles of the computed solution in physical variable and computational variable and the evolving mesh are presented in Figure 4 for q = 1 and in Figure 5 for q = 2. For simplicity, each subdomain has 50 mesh points, i.e. N = 100 for q = 1 and N = 150 for q = 2. The numerical results are coincide with that in [4,17], and the blow-up time is 2.039708648680643 at the first source x = 4.079417297361286, corresponding to the maximum value of u max = 3.16 × 10 6 .   Example 3 (Sin-type moving sources). We now consider two sources case, in which the sources move periodically with the same speed while separated by a constant distance d 1 = 2.5, that is, , α 0 (0) = 0. The Blow-up phenomenon is studied in [19] for different amplitudes A. Here we only give the numerical results for A = π (see Figure 6), since all results are similar to those in [19]. All parameters are chosen the same as those in last example. The blow-up occurs at t = 1.689611393639939 on the second source with the maximum value of u max = 3.16 × 10 6 .   Example 4 (Symmetric periodic moving sources). We consider the case for two sources, which move periodically and symmetrically. The motion are described by α ′ 0 (t) = A cos(πt), α 0 (0) = −2.0, and α 1 (t) = −α 0 (t), with e.g. A = π.
To our best knowledge, there has no theoretical results for multi-sources with different speeds and this is the first time numerically investigating the phenomenon for this case. It is shown in Figure  7 that blow-up occurs on both sources at t = 2.496881990359248, corresponding to the maximum value of u max = 3.16 × 10 6 .
If local absorbing boundary conditions (23) are used, the observed domain can be chosen more smaller while the results do not be influenced almost. See Figure 8 as an example, where the observed domain is set by [α 0 (t) − 4.0, α 1 (t) + 4.0], changed as time evolution. Now blow-up occurs on both sources at t = 2.496370241342059 with the maximum value of u max = 3.16 × 10 6 .

Conclusions
In this paper, our work focus on the problem of traveling singular sources with different speeds. A new moving mesh method in conjunction with a non-overlapping domain decomposition is proposed  for solving this problems. The whole domain is splitted into q + 1 subdomains by the q sources, whose positions are gotten by a predictor-corrector algorithm. Taking the advantages of the domain decomposition, the computation of jump [u] is avoided and there are only two different cases discussed in the discretization of the physical PDE. Thus, it is easy for the implementation to solve the problems with two traveling sources or more. Moreover, the moving mesh method of MMPDEs can be applied into each sub-domain respectively. The second-order of the spatial convergence can be proved for the new method under a special time marching implementation. The good performance of the new method for the blow-up phenomenon is demonstrated through a number of examples with two sources. Furthermore, using the new method, we successfully simulate the solutions of two sources with different speeds. To our best knowledge, this is the first time investigation for this case. The case of three sources or more can be implemented similarly.