A gradient weighted moving finite element method (GWMFE) based on piecewise polynomial of any degree is developed to solve time-dependent problems in two space dimensions. Numerical experiments are employed to test the accuracy and effciency of the proposed method with nonlinear Burger equation.

1. Introduction

Many problems in science and engineering are formulated in terms of time-dependent partial differential equations (PDEs). It is well known that due to the moving steep fronts present in the solution, these problems present serious numerical difficulties. We present an approach where the mesh moves dynamical to capture the sharp front with a small number of space nodes.

Moving finite-element method (MFE) is a discretization technique on continuously deforming spatial grids introduced by K. Miller and R. N. Miller [1, 2] to deal with time-dependent partial differential equations involving fine scale phenomena such as moving fronts, pulses and shocks [37]. Much significant work on MFE has been done by Baines and Wathen and their collaborators [1, 4, 5, 8, 9]. In particular, we mention the Baines book [8] and its excellent bibliography. For more information about moving finite-element method and its aspects and applications, see [918].

In all of these works, the method is based on a minimization of the PDE residual that is obtained by approximating the solution with piecewise linear elements. In [19], Coimbra et. al. introduced the MFE in two dimensions in which the degree of approximation polynomial was greater than 1. In [6, 7], Carlson and Miller introduced the GWMFE in two dimension in which the approximation polynomial was linear. In this paper, we present a formulation of GWMFE with approximation of higher degree to solve two dimensional time-dependent partial differential equations. This formulation of the GWMFE, which builds on the original approach of Miller [1, 2], uses piecewise polynomial approximations of any degree of the 2D spatial domain. Numerical investigations are presented to show the accuracy and effectiveness of our method.

2. The GWMFE in 2D

Our formulation of GWMFE has been designed to solve a PDE of the type where is a first or second order differential operator on the 2D spatial interval , if and be zero or not, respectively, and under Dirichlet, Neumann, or Robin boundary conditions and initial conditions satisfying , .

2.1. Description of the Method

The GWMFE is a numerical procedure which allows the local gradient adaptation of the finite-element approximation space with time. For the space discretization, we consider a hexagonally connected triangularization of , where is the number of elements on . In each triangle , the solution approximated by a polynomial of degree greater than 1. We define the polynomial approximation to as where is the th Lagrange basis function on the th element, is the number of interpolations points in an element, and is the value of at the th interpolation point of the th element. Because of minimizing the algebra requirements in the formulation of GWMFE for computing , the physical coordinates on triangular element are introduced. So let be the physical coordinates of . The Lagrangian interpolation functions are given by where , and with and that denote the degree of approximation and is th physical coordinate of the th interpolation point.

For example, for node 1 of triangle element with 6 nodes, we have then For the node 2, we have

then For other nodes, we have A weighted form of the variational formulation is often recommended, in particular, when the method is overly sensitive to specific features of the physical problem such as steep fronts. Such weighting replaces the inner products by inner products with respect to a given, positive, weight function . Then, in one dimension, we will have: In two dimension, a standard form of weighted function is the positive function , in meaning of the Galerkin approach, given rise to minimization where

It is also possible to use a solution-dependent weight function , which depends on or those known function, say , or its first derivatives, say or .

In GWMFE, this weight function is taken to be and we will have where is the gradient with respect to the physical variable and .

The argument for the use of this weight function is that it de-emphasizes those parts of the integral where and are large and therefore reduce the effect of minimization in steep parts of the solution. These moving node methods are especially suited to problems which develop sharp moving fronts, especially problems where one needs to resolve the fine-scale structure of the fronts.

We add the penalization term to the objective function (2.14) in order to prevent singularities depending on the area of each element , . So we will have Here is the time derivative of . The internodal viscosity function and the internodal spring function associated to an element are chosen in a closed form to the original proposed by Miller [1]. We consider, respectively, . The penalty constants , and are small constants supplied by the user. Penalty functions do not interfere on the solution, but exclusively on the movement of the nodes in order to prevent singularities. Their disadvantage is that it is not possible to set up a priori a relation between them and the problem we are solving.

The discretization of space-variables transforms each PDE in a system of ODE. To accomplish the discretization of problem (2.1) some overwriting may be necessary in order to apply the appropriate boundary conditions. The full discretization of (2.1) is obtained solving the ODE system by a suitable ODE solver such as LSODI [20] in FORTRAN software or ODE15S [16, 21] in MATLAB software.

2.2. Time Derivative of

The approximation to is dependent on the nodal amplitudes and on the nodal position . So we can write In order to define the system of ODEs generated by space discretization, it is necessary to evaluate the derivatives for all , and and for all .

Let us consider a global node , which is either a vertex of a triangle or an interpolation point that belongs to an edge or is placed inside the triangle. Let be the number of these global nodes. The support of a global node is the union of triangles to whom belongs, say , be the number of elements surrounding node . Suppose that is the global node associated with , the th interpolation point of th element of . Defining the global basis function , as After some simple computations, from (2.3) and (2.18), we have where is the global node associated to the th interpolation point of th element.

The computation of and for are similar. Consider After some computations [22], we have where is the th physical coordinates of .

2.3. The GWMFE Equations

Our GWMFE discretization leads to a large ODE system in which , and . Now consider the first equations in (2.22). From (2.15) and (2.17), we have for , . Considering the support of the global node , say , associated with the th interpolation point of the th element we have

for .

For the second equation in (2.22), if we suppose that is the th node space, then for .

Similarly, the third equation is for .

2.4. Second-Order Terms

Second order terms such as the Laplacian need to be interpreted in GWMFE in the sense of smoothing. That is, we imagine the corners of our GWMFE to be ever so slightly smoothed.

Based on the idea of smoothing, there are basically three techniques for dealing with this problem. The -mollification of Miller [1], the application of Greens theorem to reduce the order of the differentiation [14], and the idea of recovery [9] which requires constructing a function from the piecewise polynomial approximation with sufficient continuity to ensure that all the integrals involving second-order derivatives exist and may be evaluated. In order to define and evaluate the integrals involving second-order derivatives over the support, , of a global node, , we consider a -neighbourhood of an edge. In each of -neighbourhood of edge adjacent to the node , we need to calculate where is the -neighbourhood of an edge adjacent to the node and is the recovered function of . We take the Hermite cubic recovered function that has continuity [8, 23]. With changing the coordinates, we can describe in terms of and , the length of the edge, . Then we have where is the rotation angle. Let , in which and , the recovered function, defined by where is the Hermite cubic polynomial defined by the values of and at and .

So the integral (2.28) may be approximate without difficulties in the usual way by some quadrature rule, for example, the mid-point rule. Thus when , with some computation [19], we get where is the quadrature weight associated to the quadrature point and where and denote the right and left triangle, respectively (Figure 1).

2.5. General Equations of the GWMFE

Consider the global node , , and its support . Assume that is the th interpolation point in the element , . Denoting by the length of the th edge adjacent to , the GWMFE equation (2.24) associated to the nodal amplitude is for .

If , be a spatial node defining the th vertex of the element surrounding the node . For , (2.25) leads to

and (2.26) leads to

for , where

However (2.31), (2.32), and (2.33) define the general form of the system of ODE generated by the GWMFE. According to the boundary conditions, some of the equations may have to be rewritten.

This system of ODE can be written as follows: where .

This system of ODE has a stiff mass matrix and appropriate methods are thus required. In the present work, we use the ODE15S package [16, 21] under All cab software for integrating in time.

3. Local Time Step Refinement

Let time steps of the problem have the form where .

Now, we apply the refinement scheme at each time step [17], for example, on the first time step . Set , , Thus the time integration from to involves sub step such that

solves ODE systems similar (2.35) to approximate , and with ODE15S(@Function,,,Option).

In each ODE system, we need the initial conditions which are obtained by solving the previous ODE system. In other words, the initial condition of the th ODE system is the approximation of the th ODE system at .

Generally, suppose that we are at time level and want to move toward , similarly, consider such that The values of , and are obtained by solving ODE system This process continues until .

3.1. Algorithm

The local time step refinement (LTSR) method may be derived as follows

Step 1. Set , and , and obtain with initial value , , and set .

Step 2. Set , and .

Step 3. Solve ODE system (2.35) at time level as follows: = ODE15S(@Function,,,Option),the initial value of which in this step is obtained by solving ODE system (2.35) at .

Step 4. Set , , , , , and .

Step 5. If then go to Step 3, or else .

Step 6. If then go to Step 2, or else break.

So our solution, will be obtained in desired .

4. Numerical results

We present a numerical example to illustrate the performance of our GWMFE. The integrals that appear in the system of ODE, say (2.31), (2.32), and (2.33), are evaluated by 1D and 2D Gauss quadrature formulaes with 7 interior points [24]. For integration in time, we have used the integrator ODE15S. We select the standard choices for both absolute and relative time tolerances for the ODE15S integrator, . In order to define the penalty functions the user must supply the penalty constants , , and . The value of corresponds to the minimum allowed for an element area and in all computations we consider .

4.1. Burgers Equation in 2D

Some of the more difficult and interesting real life problems in which adaptive algorithms are needed arise in transport phenomena in which steep fronts propagate through the domain. The special case of the nonlinear Burger equation is often used to test numerical methods so we consider the nonlinear evolution equation We assume that initial and Dirichlet boundary condition are chosen to correspond to the analytic solution This problem is solved from to with with quadratic polynomial as approximation function. Solution of this problem can be obtained with similar computational effort for smaller values of . All numerical results shown here are obtained on a Pentium IV processor at 3.00 GHz.

Figure 2 presents the adaptivity and nodes movement for to seconds with quadratic elements at some cases. Figures 36 present nodes movement and their solutions with quadratic elements for second, second, second and seconds, respectively.

Let us consider to the average error, where is number of spacial nodes in , is a space node, is the GWMFE solution at and is the exact solution at .

Table 1 present CPU time, number of function evaluation (NFE) and average error () for GWMFE at seconds seconds with quadratic element and , and as penalty constants in some meshes. As Table 1 shows, when number of nodes in each direction are increased, then NFE in ODE15S and therefore CPU time has been increased and our average error has been decreased.

5. Conclusions

In this paper, we presented a gradient weighted moving finite-element method based on polynomial approximations of high degree for the solution of time-dependent PDEs on two-dimensional space domains. We used a solution-dependent weight function for original MFE formulation to have better performance and mesh adaptivity. These moving nodes method are especially suited to problems which develop sharp moving fronts, especially problems where one needs to resolve the fine-scale structure of the fronts.

A careful treatment of the general second order terms is carried out. Moreover, by using numerical evaluations of all integrals, we can solve a large class of problems without extra calculations. The GWMFE is applied to the Burger test equation for transport process with quadratic polynomial as interpolation function. One can solve this problem with other nonlinear approximation function as well as other penalty constants. Numerical results are given to illustrate the good behavior of the GWMFE when using some cases of penalty constants.