Meshless Local Petrov–Galerkin Formulation of Inverse Stefan Problem via Moving Least Squares Approximation

Abstract: In this paper, we study the meshless local Petrov–Galerkin (MLPG) method based on the moving least squares (MLS) approximation for finding a numerical solution to the Stefan free boundary problem. Approximation of this problem, due to the moving boundary, is difficult. To overcome this difficulty, the problem is converted to a fixed boundary problem in which it consists of an inverse and nonlinear problem. In other words, the aim is to determine the temperature distribution and free boundary. The MLPG method using the MLS approximation is formulated to produce the shape functions. The MLS approximation plays an important role in the convergence and stability of the method. Heaviside step function is used as the test function in each local quadrature. For the interior nodes, a meshless Galerkin weak form is used while the meshless collocation method is applied to the the boundary nodes. Since MLPG is a truly meshless method, it does not require any background integration cells. In fact, all integrations are performed locally over small sub-domains (local quadrature domains) of regular shapes, such as intervals in one dimension, circles or squares in two dimensions and spheres or cubes in three dimensions. A two-step time discretization method is used to deal with the time derivatives. It is shown that the proposed method is accurate and stable even under a large measurement noise through several numerical experiments.


Introduction
A free boundary problem (FBP) is a partial differential equation with initial and boundary conditions, in which a part of the boundary of domain, called a free boundary, is unknown at the outset of the problem. FBPs have many applications in science and engineering. FBPs usually happen in phase separation problems, which can be either stationary or moving free boundaries. The Stefan problem is one kind of the free boundary problems which describes the process of melting and solidification. In this paper, the numerical solution of these problems are considered. The determination of the temperature function and the position of the free boundary are desired [1][2][3][4][5]. The existence and uniqueness of the solution to these problems are investigated in References [2,3,5]. In recent years, several methods have been employed for solving the Stefan problems numerically, such as the homotopy analysis method [6,7], Lie-group shooting method [8], finite difference and finite element methods [9] and the variational iteration method [10]. Grzymkowski and Slota [11,12] applied the Adomian decomposition method (ADM) to solve one-phase Stefan problems and Slota [13] used the homotopy perturbation method for one-phase inverse Stefan problems. In Reference [14] the method of fundamental solutions was applied to the one-dimensional Stefan problem.
For many years, the finite element method (FEM) has been considered a standard and effective technique for numerically solving many applied problems in science and engineering [15,16]. Due to several limitations, these techniques cannot solve some of the complex problems of today's world. For this reason, the development and formulation of new and effective numerical techniques in recent years has been an interesting field for some engineers and mathematicians. In recent years, meshless methods have gained considerable attention, in engineering and applied mathematics. Flexibility and simplicity are the advantages of these methods. Meshless methods overcome the shortcomings of the mesh-based techniques [17]. In these methods, a system of algebraic equations is created using a set of scattered nodes-called field nodes-within the domain and its boundary for representation (but not for discretization) the whole domain of the problem and its boundary, therefore it is not necessary to use a predefined mesh for the domain discretization.
Meshless methods are generally divided into three categories. The first category includes methods that use integration and are based on weak forms of PDEs, such as the element free Galerkin method [18][19][20][21][22][23][24]. The second methods are based on the strong forms of PDEs and use integration, for example, the meshless collocation method based on radial basis functions (RBFs) [25][26][27][28] are in this category. The third category is a set of methods based on the combination of weak forms and strong forms. The meshless methods based on strong form are truly meshless methods and the implementation of these methods are usually simple. They are also computationally efficient. In spite of several advantages, they also have some shortcomings, such as numerical instability and less accuracy.
The meshless weak form methods are those that use the global and local weak form. The stability and accuracy of these methods make them more attractive. In these methods, the global weak forms are used and numerical integrations are carried out on the global background cells in solving the algebraic equations. In meshless local weak form methods, it does not require any background integration cells for field nodes. The meshless local Petrov-Galerkin (MLPG) method [29][30][31][32][33][34][35][36][37][38][39][40] is based on the local weak form of PDEs. In the MLPG method, the numerical integrations are performed over a local small sub-domain defined for each node. The local sub-domains usually have a regular shape, such as interval, circle, square, sphere, cube, and so forth.
The moving least square (MLS) approximations have an important role in using the MLPG method. By considering a local sub-domain for each field node, the MLS approximates the solving function at each field node. In this paper, a kind of MLPG method using the MLS approximation is applied for numerically solving the Stefan problem.
The layout of the paper is as follows. In the next section, we give the formulation of the inverse Stefan problem. Section 3 briefly describes the MLS approximation. In Section 4 we present the time discretization of the problem. In Section 5 the local weak form formulation of the discretized problem is presented. We present the MLPG discretization of the problem in Section 6. Numerical examples are given and solved to observe the performance of the proposed method in Section 7. At last, we give a conclusion in Section 8.

Statement of the Problem
Consider the following heat conduction equation subject to the initial condition: v(y, 0) = f (y), 0 < y < s 0 , s 0 = s(0), and boundary condition v(0, t) = g(t), where α denotes the thermal diffusivity and v(y, t), t and x denote the temperature, time and spatial location, respectively. On the free boundary where α is the thermal diffusivity, β is the thermal conductivity, κ is the latent heat of fusion per unit volume, h(t) is the temperature of the phase change, v(y, t) is temperature and t and x refer to time and spatial location, respectively. Equation (4) represents the continuity of temperature and Equation (5) is the Stefan condition.
In this problem, we try to find v(y, t), the temperature distribution in given domain and s(t), the free boundary.This is a nonlinear problem due to Stefan condition [4].
By using the change of variable, we the free boundary problem is transformed into a fixed boundary problem. Let Then Equations (1)-(5) are changed to In this paper, an approach based on the MLPG method and MLS approximations is applied to the Equation (7), which is subjected to the initial condition (8) and over-specified boundary conditions (9)-(11).

The MLS Approximation Technique
In this section, the formulation of MLS approximation is explained. The trial functions at each node is represented by the MLS approximation. Consider the sub-domain Ω s , with the boundary ∂Ω s , of problem global domain Ω around point x. In fact, Ω s is the domain of definition (or support) of the MLS approximation for the trial function at x. Let q t (x) = [q 1 (x), q 2 (x), ..., q m (x)] be a complete monomial basis in the space coordinate x. For example, the linear basis for one-dimensional is and the quadratic basis function is For all x belong to Ω s the MLS approximation u h (x) of u in Ω s , over a set of random nodes x i (i = 1, 2, ..., n) located in Ω s , is given as where ] is a vector of coefficients. In order to determine the unknown coefficient vector λ(x), we define a function I(λ(x)) as follows where the matrices Q and W(x) in Equation (15) are defined as In the above relations, w i (x), i = 1, 2, ..., n, is the weight function corresponding to the node x i , so that for each x in the support of w i (x) we have w i (x) > 0, n is the number of nodes in Ω s for which the weight functions w i (x) > 0 andû t = [û 1 ,û 2 , ...,û n ] is the vector of fictitious nodal values. It is necessary to mention thatû i , i = 1, 2, ..., n, are not equal to nodal values u i , i = 1, 2, ..., n, of the unknown trial function u h (x) in general ( Figure 1). The stationarity of I(λ(x)) in Equation (15) where F(x) and G(x) are matrices defined as follows The MLS approximation is well-defined only when the matrix F in Equation (16) is non-singular, that is, if and only if the rank of Q equals m. A necessary condition to have a well-defined MLS approximation is that at least m weight functions are non-zero (i.e., n > m) for each sample point x ∈ Ω. Computing λ(x) from Equation (16) and substituting it into Equation (14), gives where or The function ψ i is usually called the shape function of the MLS approximation corresponding to nodal point x i .
The partial derivative of ψ i (x) with respect to x is defined as x denotes the derivative with respect to x. In this paper, the following Gaussian weight function is used

The Time Discretization of the Problem
We use the following finite difference approximations for the time derivative operators where ., M and t = T M . Also by using the Crank-Nicolson technique, we have the following approximations: Using the above approximations, Equations (7) and (11) can be respectively written as: Suppose that λ = 2 t , then we have

The Local Weak Form Formulation
Let Ω i q be a sub-domain associated with the nodal point x i , i = 1, 2, ..., N, (called local quadrature cell) in the global domain Ω. Ω i q i = 1, 2, ..., N, overlap each other and union of them cover the whole global domain Ω. In this paper Ω i q are intervals centered at x i of radius r q . By applying the MLPG method, the local weak form is obtained over local quadrature cells Ω i q . For each node x i ∈ Ω i q the local weak form of Equation (29) is represented as follows where the Heaviside step function [41,42] ν is used as the test function. Applying the integration by parts in Equation (31) the following local weak form equation is obtained: where ∂Ω i q is the boundary of Ω i q .

MLPG Discretization
In this section, we obtain a system of algebraic equations from discretization of the Equation (33), by employing MLS approximation. Equation (33) is discretized for this purpose. Consider the N regularly points x i , i = 1, 2, ...N, in the domain of the problem and its boundary such that x i+1 − x i = h. Suppose that u(x i , t k ), is determined and u(x i , t k+1 ), is unknown for i = 1, 2, ...N. In order to determine the N unknown quantities u(x i , t k+1 ) we need to have N equations. For interior nodes x i of the domain Ω, by replacing MLS approximation Formula (19) in the Equation (33), the following discrete equations are obtained For boundary nodes x = 0 and x = 1 we set together with Equation (30) for which its discrete form is written as The matrix form of Equations (34)- (36) for all N nodal points can be represented as follows where Assuming that and U = (u i ) N×1 , Equations (37) and (38) yield the following system of equations According to the boundary conditions (35) and (36) we have for each step At the first step, when k = 0, due to the initial conditions, we have the following assumptions

Numerical Experiments
In this section, we test the described meshless method with two examples. In numerical computations, the input Cauchy data are considered with noisẽ where δ denotes the level of noise, and R(i) are random numbers in [−1, 1]. In these examples, the domain integrals are approximated using the 4 points Gaussian quadrature rule. In order to investigate the accuracy of computed approximations and the efficiency of the presented method, the following root mean square (RMS) error and absolute error formulas are applied Absolute error s = |s exact (t j ) − s approx (t j )|.
In implementing the meshless local weak form, each local quadrature domain Ω i q is taken as interval centered at x i of radius r q = 0.7h where h = x i+1 − x i , i = 0, 1, 2, ...N − 1. Also the radius of support domain Ω s is r s = 4r q and the quadratic basis functions (13) is used in Equation (14).
The results of using the proposed method are obtained with t = 0.01, h = 0.01. Figure 2 presents the RMS error for u(x, t) and Absolute error for s(t) versus the shape parameter c at t = 1. For other values of t the results are almost the same. In this example, the interval (0.0097, 0.01155) is suggested for choosing c. It is necessary to note that, ill conditioning occurs by increasing c. Hereafter, we fix it at c = 1.1h = 0.011. In Figure 3 RMS versus N is plotted at t = 1. It can be seen that in this figure, the error values decrease by increasing N. Values of RMS error for u(x, t) and Absolute error for s(t) with δ = 0 and δ = 0.1 are presented in Table 1. It shows that the numerical results are more accurate when there exists no noise on the input data. Under a noise level δ = 0.1, the numerical result obtained by the MLPG method is also acceptable. Exact solution, numerical solution and Absolute errors for u(x, t) and s(t) are plotted in Figures 4-6.

Example 2.
In this example, we consider the problem (1)-(5) with α = 1, β = −1, κ = 1, T = 1 and The analytical solutions of problem are given by v(y, t) = exp(1 − The approximation results of using the proposed method are obtained with t = 0.01, h = 0.01. Figure 7 presents the RMS error for u(x, t) and absolute error for s(t) versus the shape parameter c at t = 1. For other values of t the results are similar. In this example, the interval (0.0097, 0.01155) is also suggested for choosing c. It is necessary to note that ill conditioning occurs by increasing c. Hereafter we fix it at c = 1.1h = 0.011. In Figure 8 RMS versus N is plotted at t = 1. We observe that in this figure, the error values decrease by increasing N. Values of RMS error for u(x, t) and absolute error for s(t) with δ = 0 and δ = 0.1 are presented in Table 2. We see that the numerical results are more accurate when there exist no noise on the input data. Under a noise level δ = 0.1, the numerical result obtained by the MLPG method is also acceptable. The exact solution, numerical solution and absolute errors for u(x, t) and s(t) are plotted in Figures 9-11, respectively.

Example 3.
With an example, we compare the average values of the absolute errors between the present method and the Adomian decomposition method and the fourth-order Runge-Kutta method obtained in Reference [43]. We consider the problem (1)-(5) with α = 0.1, κ = 10, β = −1, T = 0.5 and The analytical solutions of problem are v(y, t) = exp(1 − y + 0.1t), s(t) = 0.1t + 1, and u(x, t) = v(xs(t), t) = exp(1 − x(0.1t + 1) + t), s(t) = 0.1t + 1. Tables 3 and 4 show the results of calculations related to the reconstruction of the moving boundary and the temperature distribution using the present method, the Adomian decomposition method and the fourth-order Runge-Kutta method. It can be seen that the presented procedure in this paper is useful and efficient in finding solutions of the considered problem.

Conclusions
In this paper, a kind of MLPG method using the MLS approximation to represent the trial function at each field node, is applied for numerically solving a nonlinear one-phase Stefan problem. Nonlinearity of this problem is due to the Stefan condition. The free boundary problem is transformed into a fixed boundary by the change of variable. In the presented method, all integrations are performed over small local quadrature domains so it does not require using any background integration cells. In the proposed method, the shape functions are produced by the MLS approximation technique. A two-step time discretization method is used to approximate the time derivatives operators. The Heaviside step function was used as the test function in the local weak form method in MLPG. Numerical results show that the proposed method is accurate and stable, although under a large measurement noise.
Author Contributions: All authors contributed equally to this work.