Some Research Of Numerical Solution For Inverse Problem In Parabolic Equation

Inverse problem in parabolic equation is mainly discussed in heat conduction mod-eling and contaminant transport modeling, while the homotopy method is a global con-vergent method with a lot of application. In this article, we applied the homotopy method to the inverse problem in parabolic equation(s), focused on the problem of de-tecting the source term and permeability coefficients. The algorithm and numerical experiments are supplied to show the benefit of this method.1. Inverse problem for source term in a parabolic equationThe inverse heat problem, which we will discuss, IP, can be stated as following:let φ{x) and θ(x) be given functions defined in (0,1) with θ{x) satisfying the homogeneous Neumann boundary condition. Determine the functions u(x, t) and q(x) satisfying andThe L1-norm estimate of q(x) in terms of the initial and final conditions is pro-posed. We assume that φ(x), q(x) and θ(x) satisfy the following conditions(Al) φ(x) ∈C1,φ{x)≥0 and φx(0)=φx(1)=0,x∈[0,1],(A2) q(x) ∈C0, q(x)≥0,x∈[0,1],(A3) θ(x)∈C0,x∈[0,1].The extremum principle for parabolic equation was used to obtain the following resultsLemma 1. Assume that φ(x) and q(x) are given and satisfy, respectively, conditions (Al) and (A2). Then the direct problem of finding u(x,t) solving equation (1) has a unique solution satisfying consequentlyNext, we give the stability estimate of q(x) in L1-norm.Lemma 2. Ifφ(x) and q(x) satisfy, respectively, the conditions (Al) and (A2). More-over if then there exists a constant γ0>0 such thatTheorem 1. Let the assumptions in Lemma (2) hold, then we have the estimate of q(x) as followsWe discretize the IP using the finite-difference method, and then present the model in a discrete setting.First, the equation (1) is approximated by the following finite-difference scheme: where r and h are, respectively, the time step size and the space step size. Let xn= nh, tm= mr, qn= q(xn), N=1/h, M=T/τ, and denote u(xn, tm) by unm. The equation (2) is discretized as following For simplicity, we define a vector in RN+1 as follows: q＝(q0,q,,…,qN)T, the inversion of q can be reduced to solve a nonlinear operator equation: F(q)= 0. (9) As the the estimation of q(x), we will give an estimate of q for the discrete problem (7)-(8).Lemma 3. Let q and φ be nonnegative vectors, and satisfy Then the solution um of the equation (8) is nonnegative and satisfies ‖um‖∞≤‖φ‖∞,m=1,2,...,MTheorem 2. Let the assumptions in Lemma (3) hold, and there exists a positive con-stant γ1 such that Then the following estimate about q holdswe construct a Newton homotopy map H:RN+1×[0,1] ×RN+1�'RN+1 by H(q,λ,q0)= F(q)-(1-λ)F(q0), (14) where A ∈ [0,1], and q0 is the initial value of the process.The theoretical foundation of the probability of one of the globally convergent ho-motopy methods was given in the Parameterized Sard’s theorem in previous reports. The main result of this theorem is that the zero set Hq0-1 (0) consists of smooth, nonin-tersecting curves in U x [0,1). These curves either are closed loops, or have endpoints in U x 0 or U x 1, or go to infinity. Another important consequence is that these curves have a finite arc length in any compact subset of U×[0,1).In the end, we prove the existence of the homotopy curve for the equation (9). We introduce some definitions first: where E0, E are positive constants.Theorem 3. Let E0 be chosen suitably, that Jacobian matrix DF(q0) has full rank at q0 E ∈QE0 And the assumptions in Theorem (2) hold. Then we have a zero curve (q(s),λ(s)) as Hq0{q, λ)= 0 for almost all q0∈QE0.Along this curve, the Jacobian matrix has full rank, and we could follow this zero curve from (q0,0) to (q(smax), λ(smax)). λ(smax) has the two following situations:1.λ(smax)<1,{q(s):s ∈[0,smax]} runs across SE and out of QE,2. λ(smaa;)= 1,{q{s):s ∈ [0, smax]} C QE and F(q{Smax))= 0.In this section, an adaptive regularized homotopy method is designed successfully to solve the radiative coefficient identification problem of the heat equation. We focused on the theoretical analysis of the homotopy method applying to IP. Evidence for the existence of q is given by the existence of a zero curve of the homotopy map, and then the numerical inversion of q could be made by tracing the zero curve of the homotopy map. The theoretical details for this is also discussed. The numerical results indicate that the homotopy method is globally convergent, efficient and stability. In addition, it was effective for the discontinuous q(x), even q(x) is similar to a 5-function.2. Inverse problem for permeability coefficient in parabolic equa-tions systemIn this section we discuss the inverse problem for a(x) in the following PDE sys-tem: with initial and boundary conditions subject to the observation data θ at the final time where Γ1 and Γ2 are on the boundary of Ω,and not intersected.For simplicity,Ω is assumed to be the unit square.We have a estimate of α(x)as followsTheorem 4.Suppose that (A1-A3) re satisfied,thereexists a positive constant γ0>0, s.t. S is the proportion of Ω.we havePDE System(15)can be approximated by the following finite-differenee equa-tions. where ui,jn=u(iΔx,jΔy,nΔt),vi,jn=v(iΔx,jΔy,nΔt),qij=1(iΔx,jΔy),N= 1/Δx=1Δy,M=T/τ,τ is the time step size and Δx,△y are the step sizes of the square grid in the x-and y-direction,then V.(αi,j(?)uijn)is the discretization of the diffusion term.Let represent the discrete derivatives in the x.direction for the difference approximations, respectively.A corresponding notation is used for y-direction derivatives.The mean values with the discretized permeability a are defined with where Then the discretization of the diffusion term is expressed as Initial and boundary conditions are discretized as Equation (20) can be discretized asThe estimation of a is as followsTheorem 5. Suppose there exists a positive constant γ1, s.t. S is the proportion of Ω, here S= 1, and we haveThen the inversion, for a is reduced to a nonlinear operator equation As is well known, the nonlinear operator Equation (30) is an ill-posed problem, we introduce the Tikhonov regularization, we use the following optimal problem instead of solving (30) where‖·‖ is the L2 norm, α1,α2 are the regularization parameters, and M1, M2 are the second-order smoothing matrices in the x-and y-direction. M1 is as follows:The corresponding Euler equation of (31), which has the same solution, where "T" denotes the matrix transposition. We use a successive method to construct a basic iterative method.Assume that the k -th approximation ak of Equation (32) has been found, to get the (k+1)-th approximation ak+1, according to Gauss-Newton’s method Avoiding the impact of the second derivative, we have the iterative method below which can be called as the iteratively regularized Gauss-Newton’s method. It has a fast convergence speed and good stability. But it is only of local convergence. To overcome the weakness of the local convergence of Gauss-Newton’s method, we use the homotopy method to solve (31). Construct the homotopy equation where λ∈ [0,1] is the homotopy parameter,α0 is the given initial value of the whole computation. We also need to regularize (35) because of its ill-posednessLet us choose an isometric division λω=ω/W, ω= 0,1,..., W. The corresponding ω-th Euler equation is Assume that the solution aω-1 of the (ω-1)-th equation has already been found, it could be used as the initial guess of the ω-th equation. The iterative formula can then be constructed as After αw is obtained, we also make correction by using (34). The method constituted by (38) and (34) can be called as a simplified regularized homotopy method.The homotopy parameters λw in (38) play a very important role. If the value of W is too large, there will be large burden of calculation, if a too small value is selected, the solution may no converge. Thus we properly modify the simplified regularized homotopy method by choosing the homotopy parameters adaptively.Algorithm The adaptive regularized homotopy method.(1) Choose an initial value a0, small threshold ∈, initial homotopy step size δ=1/W, regularizationparameter α1,α2, and set λ1= δ, λ2= λ1+δ, ω= 2. Calculate α1 by(2) Calculate αw from(3) Calculate and check if κ<ε, and then set δ=2δ; if κ≥ε, then set δ= 0.5δ. Let ω= ω+1 and λω=λω+δ. If λω< 1, go to step (2).(4) Make a correction by using Gauss-Newton on the initial estimate αω until the global minimum is obtained.From the results of numerical simulations, it seems clear that the proposed method is globally convergent, computationally efficient and stable with respect to noise. In addition, it is necessary in contaminant transport modeling because of the greatly large jumps in the permeability.3. Laplace transform for inverse problem of permeability coefficient in parabolic systemIn this section we discuss the inverse problem for a(x) in the following PDE sys-tem: with initial and boundary conditions subject to the observation data h(t) asThis is a seriously ill-posed inverse problem, we use the Laplace transform to make this problem easier to solve. Let u(x, s) be the Laplace transform of u(x, t) with t,The equation system (39) becomes where Define map G(s; a, k) as thenThe inverse problem (39) could become the next nonlinear operator equationWe solve (49) with homotopy method. From the results of numerical simulations, it seems clear that the proposed method is globally convergent, computationally efficient and stable with respect to noise. In addition, it is necessary to be improved because of bad result in greatly large jumps in the permeability.