A Full Predictor-Corrector Finite Element Method for the One-Dimensional Heat Equation with Time-Dependent Singularities ()
1. Introduction
The heat equation is a partial differential equation known to model problems in thermodynamics, Brownian motions and particle diffusion. In the last seven decades, the increasing availability of powerful computers has been a major propulsion to the implementation of numerical methods for partial differential equations. Concerning the numerical solution of the heat equation, the usual methods are the finite difference methods (see [1] [2] [3] ), the methods of lines (see [4] [5] [6] ) and more recently the space-time finite element methods (FEMs) (see [7] - [12] ).
The singular behavior of solutions of a partial differential equation is known to affect the convergence rate of the numerical approximation. The singularity can be geometric (due to the existence of reentrant corners or edges, the change of boundary conditions), or can be directly related to a singularity of a component of the equation (sign-changing coefficient, component singular at a point...). The goal of this paper is to use classical FEM to solve the one-dimensional heat problem: Find u such that:
(1)
where
,
,
and the source term
is some function singular at a point
.
The space-time FEM is known to give an approximation of the solution of Problem (1) with optimal convergence rates, but the method either uses discontinuous Galerkin methods (see [12] ) that increase the number of degrees of freedom, or a particular interpolation error estimate in order to guarantee the optimal convergence (see [7] ). Those techniques are efficient but involve particular adaptations’ techniques.
Singular complement methods and predictor-corrector FEM (PCFEM) have been introduced to solve the problem of geometric singularities of elliptic boundary value problems (see [13] - [23] ) and more recently on geometric singularities of hyperbolic problems (see [24] ). The treatment of time-dependent singularities by means of PCFEM is also introduced in [25] using the method of lines. The main advantage of PCFEM is its application on any standard type of mesh (structured or unstructured) and on usual interpolation operators. This paper presents a PCFEM for the one-dimensional heat Problem (1). The method relies on a Fourier decomposition of the solution where computable formulas of the singular functions are derived and the solution is recomputed with a corrector algorithm in order to cancel the effect of the singularity caused by the source term f.
This paper is organized as follows: Section 2 presents the initial boundary value problem with the function spaces, the existence and uniqueness of the solution, a stability result and the Fourier decomposition of the solution with an extraction formula of the coefficients of the singularity. Section 3 presents the FEM, the predictor-corrector algorithm and some error estimates. Section 4 gives some numerical results that show the efficiency of the method and Section 5 is the conclusion of the paper.
2. Initial Boundary Value Problem
2.1. Function Spaces and Stability Results
In order to study Problem (1), it is crucial to introduce the following spaces:
1) Given a domain
,
, the space
of classes of Lebesgue-measurable and square-integrable functions over Ω with the norm:
where
can be a vector or a scalar depending on n.
2) Given an interval
and a Hilbert space V,
will denote the space of bounded k-continuously differentiable functions of the form
with the norm:
where
is a multi-index,
and
.
3) Given an interval
, a Hilbert space V and
, define the space
be the space of k-continuously differentiable functions of the form
such that
,
, is α-Hölder. The space
is equipped with the norm:
4) For a given domain Ω, the space
,
,
denotes the space of Lebesgue-measurable classes of functions u such that
and
for
, equipped with the norm:
The space
where
represents the boundary of Ω. The space
is equipped with the norm of
.
The space
is the dual space of
in the sense of distributions.
5) For an interval
and a domain Ω, the space
is the completion of
with respect to the norm:
The space
with the norm:
where
The space
with the norm:
6) The space
7) The space
.
The variational formulation of Problem (1) is: Find
such that:
(2)
where
.
Theorem 1. ( [26] p. 382 (Theorem 11.3))
Assume that
. Then, Problem (1) has a unique solution
.
Remark. From the Sobolev embedding theorem (see [27] , p. 98]), the space
, so the solution u of Problem (1) is at least
-Hölder continuous in time.
Define the linear continuous map named as the Newton potential (see [7] , p. 5])
by
such that
and
for any
, if one uses the Cauchy-Schwartz inequality and the Hölder inequality,
Hence,
is bounded and
(3)
so if
, one has:
and using the definition of
, the Cauchy-Schwartz inequality and (3),
(4)
Hence, from (3) and (4),
and from (3),
Then,
(5)
and one has the stability condition:
(6)
2.2. Fourier Decomposition of the Solution
Assuming that the solution u of (1) is an element of
,
, and has a line singularity at a point
, based on the boundary conditions one can set
and obtain use separation of variables to obtain the following result.
Theorem 2. The solution u of Problem (1) can be written as
where
, and setting
,
1) If
,
(7)
and
.
2) If
,
(8)
and
.
3. Finite Element Discretization
3.1. Finite ELement Method
Let
be a triangulation of the domain
with quasi-uniform meshes of diameter h and define:
with the norm of
.
The Galerkin approximation consists to find
such that:
(9)
The variational Problem (9) has a unique solution according to Theorem 1 and since
, for any
, from the inequality (5),
(10)
where
is a constant depending on
.
Furthermore, one has the Galerkin orthogonality:
(11)
The following theorem is a lemma equivalent to the Céa’s Lemma known for elliptic boundary value problems. This theorem has been proved in [7] in a more general case.
Theorem 3. Let u be the solution of the variational Problem (2) and
be the solution of the discretized variational Problem (9). Then, there exists
depending on
such that:
Proof. From the discrete stability condition (10) and the Galerkin orthogonality (11), for any
,
Then,
for any
where
is an arbitrary constant depending on
.
Hence,
.£
The following theorem proves the energy error estimate when the solution is smooth. This theorem has been proved in [7] in a more general case.
Theorem 4. Let u be the unique solution of Problem (1) and
be the unique solution of the equivalent discretized variational Problem (9). Assume
, then
Proof.
where
is the unique solution of the variational problem,
Note that for any
such that
and
for
or
,
Then,
£
3.2. Predictor-Corrector Algorithm
The predictor-corrector method presented is based on the decomposition of the approximate solution
where
is the solution of (9),
and
are defined as in Theorem 2.
Predictor-Corrector Algorithm
1) Solve the problem: Find
such that:
2) For
fixed, do:
a) For
, compute
and
using Formula (7) or (8).
b) If
solve the problem: Find
such that:
c) Compute
.
3.3. Error Estimates
From [28] , if
,
, then for any
,
,
(12)
where
and
are defined in Theorem 2 and
is an arbitrary constant.
The convergence result in the energy norm is given by the following theorem.
Theorem 5. Let
,
be the solution of Problem (1) and let
,
,
be the approximate solution of u using the Predictor-Corrector Algorithm 3.2 on a quasi-uniform and triangulation
of the domain
. Assume further that
, then there is
such that:
Proof. Suppose
is the decomposition of the solution u as in Theorem 2. Set
be the approximation of u using the Predictor-Corrector Algorithm 3.2.
From Theorem 4, (12) and the assumption
, one has:
£
4. Numerical Experiments
The main purpose of this section is to present some computational results carried out with the Predictor-Corrector Algorithm 3.2. Even though the examples considered may not represent any real world situation, the difficulties due to singularities in the solution are well reflected. We consider two examples with known solutions. All computations are carried out using a code developped in python. The resulting linear systems are solved iteratively with a generalized minimal residual method (GMRES) and all the integrals are computed using a 15-point Legendre Gauss formula on the stifness, mass matrices and load vectors, and a 150-points Legendre Gauss formula for integrals used to compute the Fourier coefficients.
In all the computations, starting from a triangulation
, the refinement are perfomed by dividing the concerned lentgh by two (
). The order or rate of convergence
is computed using the formula:
In all the following examples,
is the expected order of convergence. Figure 1 presents the type of mesh refinements used in the examples.
4.1. Example 1
In this example, we consider Problem (9) with the exact solution
and
,
,
.
The function has a line singularity at
and belongs to
. Table 1 presents the error estimates in the energy norm
and the required number of Fourier coefficients N in order to obtain an optimal convergence
Figure 1. The uniform meshes: Refinements 1, 2 and 3
Table 1. Error estimates in the energy norm.
rate. Note that we need
at each refinement
. A convergence of the classical FEM is observed and the convergence rate that is greater than 0.5, which is better than the expected convergence rate 0.25. However the convergence rate of the classical FEM is still not optimal. From the results obtained with the PCFEM, the optimal convergence rate is recovered.
Figures 2-4 present the plots of the exact solution, the approximate solution before adaptation and the aproximate solution using the Predictor-Corrector Algorithm 3.2 respectively at the fixed times
,
,
and
. Figures 5-7 present the 3-dimensional line countour plots of the exact solution, the approximate solution before adaptation and the approximate solution after using the predictor-corrector adaptation (PCFEM) with 75 lines. It is observed that the predictor-corrector solution has a better approximation with a least number of nodes in the refinement.
4.2. Example 2
In this example, we consider Problem (9) with the exact solution
and
,
,
.
The solution has a line singularity at
and is an element of
so the number of Fourier coefficients at each refinement on a triangulation of diameter h should satisfy
. Figure 8 shows the
Figure 2. The exact solution, the FEM solution and the PCFEM solution at Refinement 3 (
).
Figure 3. The exact solution, the FEM solution and the PCFEM solution at Refinement 4 (
).
Figure 4. The exact solution, the FEM solution and the PCFEM solution at Refinement 5 (
).
Figure 5. 3D plots of the exact solution, the FEM solution and the PCFEM solution at Refinement 3.
Figure 6. 3D plots of the exact solution, the FEM solution and the PCFEM solution at Refinement 4.
Figure 7. 3D plots of the exact solution, the FEM solution and the PCFEM solution at Refinement 5.
Figure 8. The graph of
against
,
,
.
error estimates of the FEM and the PCFEM. The convergence rate of the classical FEM is not optimal while the one of the PCFEM is optimal.
5. Conclusions
This paper presents an adaptive predictor-corrector finite element method (PCFEM) for the one-dimensional heat equation with zero initial and boundary conditions. The method is based on a Fourier decomposition with a correction of the solution by a predictor-corrector algorithm. The PCFEM is shown to be efficient on quasi-uniform meshes and then can be used without a particular mesh treatment or a special interpolation error in the energy norm. In the predictor-corrector algorithm, the stiffness matrices are the same, only the load vector is modified to correct the solution, and then the complexity of the PCFEM algorithm and the classical FEM algorithm are of the same order.
The PCFEM developed can be applied to the heat equation with non-zero initial and boundary conditions, and the proof of the stability conditions and the Homogenization methods can be adapted using the proofs developed in [7] provided that the initial and boundary conditions are regular enough. It is worth saying that the method developed is another motivation to continue using standard FEM to solve partial differential equations.
List of Notations
The following are some variables used in the paper.