A Gradient Weighted Moving Finite-Element Method with Polynomial Approximation of Any Degree

Abstract

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 [3–7]. 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 [9–18].

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 𝜕𝑢𝜕𝑡=𝐿(𝑢), (2.1) where 𝐿(𝑢)=𝐹(𝑥,𝑦,𝑡,𝑢,𝜕𝑢,𝜕𝑥𝜕𝑢)𝜕𝜕𝑦²𝑢𝜕𝑥²+𝐺(𝑥,𝑦,𝑡,𝑢,𝜕𝑢,𝜕𝑥𝜕𝑢)𝜕𝜕𝑦²𝑢𝜕𝑦²+𝐻(𝑥,𝑦,𝑡,𝑢,𝜕𝑢,𝜕𝑥𝜕𝑢𝜕𝑦) (2.2) 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 𝑢(𝑥,𝑦,0)=𝑢₀(𝑥,𝑦), (𝑥,𝑦)∈Ω.

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 Ω=⋃^𝑛𝑒_𝑗=1Ω_𝑗, 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 𝑈_𝑗(𝑥,𝑦,𝑡)=^𝑛𝑝_𝑘=1𝑏_𝑗,𝑘(𝑥,𝑦)𝑈^𝑘_𝑗(𝑡), (2.3) 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 𝑏_𝑗,𝑘(𝑥,𝑦)=𝐵^(𝑘)(𝐿₁)𝐵^(𝑘)(𝐿₂)𝐵^(𝑘)(𝐿₃), (2.4) where 1≤𝑗≤𝑛_𝑒, 1≤𝑘≤𝑛_𝑝, and 𝐵^(𝑘)(𝐿_𝑁)=^{𝑝𝐿𝑁(𝑘)}_𝑠=11𝑠(𝑝𝐿_𝑁−𝑠+1)𝑝𝐿_𝑁^(𝑘)≥1,1𝑝𝐿_𝑁^(𝑘)=0, (2.5) with 𝑁=1,2,3 and 𝑘=1,…,𝑛_𝑝 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 𝐵⁽¹⁾(𝐿₁)=^2×1_𝑠=11𝑠(2𝐿₁−𝑠+1)=𝐿₁(2𝐿₁𝐵−1),₍₁₎(𝐿₂)=1,𝐵⁽¹⁾(𝐿₃)=1, (2.6) then 𝑏_𝑗,1=𝐿₁(2𝐿₁−1). (2.7) For the node 2, we have

𝐵⁽²⁾(𝐿₁)=^2×(1/2)_𝑠=11𝑠(2𝐿₁−𝑠+1)=2𝐿₁,𝐵⁽²⁾(𝐿₂)=^2×(1/2)_𝑠=11𝑠(2𝐿₂−𝑠+1)=2𝐿₂,𝐵⁽²⁾(𝐿₃)=1, (2.8) then 𝑏_𝑗,2=4𝐿₁𝐿₂. (2.9) For other nodes, we have 𝑏_𝑗,3=𝐿₂(2𝐿₂−1),𝑏_𝑗,4=4𝐿₂𝐿₃,𝑏_𝑗,5=𝐿₃(2𝐿₃−1),𝑏_𝑗,6=4𝐿₁𝐿₃. (2.10) 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: (𝑓,𝑔)_𝑤=_Ω𝑓(𝑥)𝑔(𝑥)𝑤(𝑥)𝑑Ω. (2.11) In two dimension, a standard form of weighted function is the positive function 𝑤(𝑥,𝑦), in meaning of the Galerkin approach, given rise to minimization min_{̇𝑈𝑖,̇𝑋𝑖,̇𝑌𝑖Ω}𝑅²𝑤(𝑥,𝑦)𝑑Ω, (2.12) where

𝑅=𝜕𝑈𝜕𝑡−𝐹(𝑥,𝑦,𝑡,𝑈,𝜕𝑈,𝜕𝑥𝜕𝑈)𝜕𝜕𝑦²𝑈𝜕𝑥²−𝐺(𝑥,𝑦,𝑡,𝑈,𝜕𝑈,𝜕𝑥𝜕𝑈)𝜕𝜕𝑦²𝑈𝜕𝑦²−𝐻(𝑥,𝑦,𝑡,𝑈,𝜕𝑈,𝜕𝑥𝜕𝑈𝜕𝑦). (2.13) 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 1/1+‖∇𝑈‖²₂ and we will have min_{̇𝑈𝑖,̇𝑋𝑖,̇𝑌𝑖}𝜙=_Ω𝑅²1+‖∇𝑈‖²₂𝑑Ω, (2.14) 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 ∑^𝑛𝑒_𝑗=1(𝜀_𝑗Δ̇_𝑗−𝑆_𝑗) to the objective function (2.14) in order to prevent singularities depending on the area of each element Δ_𝑗, 𝑗=1,…,𝑛_𝑒. So we will have min_{̇𝑈𝑖,̇𝑋𝑖,̇𝑌𝑖}𝜙=_Ω𝑅²1+‖∇𝑈‖²₂𝑑Ω+^𝑛𝑒_𝑗=1(𝜀_𝑗Δ̇_𝑗−𝑆_𝑗). (2.15) 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, 𝑆_𝑗=𝐶₁Δ_𝑗−𝐶₃𝐶(1+³Δ_𝑗−𝐶₃)²,𝜀_𝑗𝐶=(¹Δ_𝑗−𝐶₃+𝐶₂𝐶)(1+³Δ_𝑗−𝐶₃)², (2.16)𝑗=1,…,𝑛_𝑒. 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 𝜕𝑈=𝜕𝑡^𝑛𝑒𝑛_𝑗=1^𝑝_𝑘=1𝜕𝑈𝜕𝑈^𝑘_𝑗𝑈̇^𝑘_𝑗+^𝑛𝑠_ℓ=1𝜕𝑈𝜕𝑋_ℓ𝑋̇_ℓ+^𝑛𝑠_ℓ=1𝜕𝑈𝜕𝑌_ℓ𝑌̇_ℓ. (2.17) In order to define the system of ODEs generated by space discretization, it is necessary to evaluate the derivatives 𝜕𝑈/𝜕𝑈^𝑘_𝑗 for all 1≤𝑗≤𝑛_𝑒, 1≤𝑘≤𝑛_𝑝, and 𝜕𝑈/𝜕𝑋_ℓ and 𝜕𝑈/𝜕𝑌_ℓ for all 1≤ℓ≤𝑛_𝑠.

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 Ω^𝑖=⋃^𝐼(𝑖)_𝑠=1Ω_𝑖𝑠, 𝐼(𝑖) 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 𝜑_𝑖, 1≤𝑖≤𝑛_𝑔 as 𝜑_𝑖=0(𝑥,𝑦)∉Ω^𝑖,𝑏_𝑗,𝑘(𝑥,𝑦)∈Ω_𝑗⊂Ω^𝑖. (2.18) After some simple computations, from (2.3) and (2.18), we have 𝜕𝑈𝜕𝑈^𝑘_𝑗=𝜑_𝑖, (2.19) where 𝜑_𝑖 is the global node associated to the 𝑘th interpolation point of 𝑗th element.

The computation of 𝜕𝑈/𝜕𝑋_ℓ and 𝜕𝑈/𝜕𝑋_ℓ for 1≤ℓ≤𝑛_𝑠 are similar. Consider 𝜕𝑈𝜕𝑋_ℓ=^𝑛𝑝_𝑘=1𝜕𝑏_𝑘𝜕𝑋_ℓ(𝑥,𝑦)𝑈^𝑘_𝑗,𝜕𝑈𝜕𝑌_ℓ=^𝑛𝑝_𝑘=1𝜕𝑏_𝑘𝜕𝑌_ℓ(𝑥,𝑦)𝑈^𝑘_𝑗. (2.20) After some computations [22], we have 𝜕𝑈𝜕𝑋_ℓ=−𝐿_𝑘1𝜕𝑈,𝜕𝑥𝜕𝑈𝜕𝑌_ℓ=−𝐿_𝑘1𝜕𝑈𝜕𝑦, (2.21) where 𝐿_𝑘1 is the 𝐾₁th physical coordinates of (𝑥,𝑦)∈Ω_𝑗.

2.3. The GWMFE Equations

Our GWMFE discretization leads to a large ODE system 𝜕𝜙𝜕𝑈̇^𝑘_𝑗=0,𝜕𝜙𝜕𝑋̇_ℓ=0,𝜕𝜙𝜕𝑌̇_ℓ=0, (2.22) in which 1≤𝑗≤𝑛_𝑒, 1≤𝑘≤𝑛_𝑝, and 1≤ℓ≤𝑛_𝑠. Now consider the first equations in (2.22). From (2.15) and (2.17), we have _Ω𝑅²1+‖∇𝑈‖²₂𝜕𝑈𝜕𝑈^𝑘_𝑗𝑑Ω=0, (2.23) for 1≤𝑗≤𝑛_𝑒, 1≤𝑘≤𝑛_𝑝. Considering the support of the global node 𝑣_𝑖, say Ω^𝑖, associated with the 𝑘th interpolation point of the 𝑗th element we have

_Ω𝑖11+‖∇𝑈_𝑖‖²₂𝜑_𝑖𝜕𝑈𝜕𝑡𝑑Ω=_Ω𝑖11+‖∇𝑈_𝑖‖²₂𝜑_𝑖𝐹(𝑥,𝑦,𝑡,𝑈,𝜕𝑈,𝜕𝑥𝜕𝑈)𝜕𝜕𝑦²𝑈𝜕𝑥²+𝑑Ω_Ω𝑖11+‖∇𝑈_𝑖‖²₂𝜑_𝑖𝐺(𝑥,𝑦,𝑡,𝑈,𝜕𝑈,𝜕𝑥𝜕𝑈)𝜕𝜕𝑦²𝑈𝜕𝑦²+𝑑Ω_Ω𝑖11+‖∇𝑈_𝑖‖²₂𝜑_𝑖𝐻(𝑥,𝑦,𝑡,𝑈,𝜕𝑈,𝜕𝑥𝜕𝑈𝜕𝑦)𝑑Ω, (2.24) for 𝑖=1,…,𝑛_𝑔.

For the second equation in (2.22), if we suppose that 𝑣_𝑖 is the 𝑙th node space, then _Ω𝑖11+‖∇𝑈_𝑖‖²₂𝜕𝑈𝜕𝑋_ℓ𝜕𝑈𝜕𝑡𝑑Ω=_Ω𝑖11+‖∇𝑈_𝑖‖²₂𝜕𝑈𝜕𝑋_ℓ𝐹(𝑥,𝑦,𝑡,𝑈,𝜕𝑈,𝜕𝑥𝜕𝑈)𝜕𝜕𝑦²𝑈𝜕𝑥²+𝑑Ω_Ω𝑖11+‖∇𝑈_𝑖‖²₂𝜕𝑈𝜕𝑋_ℓ𝐺(𝑥,𝑦,𝑡,𝑈,𝜕𝑈,𝜕𝑥𝜕𝑈)𝜕𝜕𝑦²𝑈𝜕𝑦²+𝑑Ω_Ω𝑖11+‖∇𝑈_𝑖‖²₂𝜕𝑈𝜕𝑋_ℓ𝐻(𝑥,𝑦,𝑡,𝑈,𝜕𝑈,𝜕𝑥𝜕𝑈−𝜕𝑦)𝑑Ω_{𝐼(𝑖)𝑠=1}12𝜀_𝑖𝑠Δ𝑦_𝑘1(𝜀_𝑖𝑠Δ̇_𝑖𝑠−𝑆_𝑖𝑠), (2.25) for ℓ=1,…,𝑛_𝑠.

Similarly, the third equation is _Ω𝑖11+‖∇𝑈_𝑖‖²₂𝜕𝑈𝜕𝑌_ℓ𝜕𝑈𝜕𝑡𝑑Ω=_Ω𝑖11+‖∇𝑈_𝑖‖²₂𝜕𝑈𝜕𝑌_ℓ𝐹(𝑥,𝑦,𝑡,𝑈,𝜕𝑈,𝜕𝑥𝜕𝑈)𝜕𝜕𝑦²𝑈𝜕𝑥²+𝑑Ω_Ω𝑖11+‖∇𝑈_𝑖‖²₂𝜕𝑈𝜕𝑌_ℓ𝐺(𝑥,𝑦,𝑡,𝑈,𝜕𝑈,𝜕𝑥𝜕𝑈)𝜕𝜕𝑦²𝑈𝜕𝑦²+𝑑Ω_Ω𝑖11+‖∇𝑈_𝑖‖²₂𝜕𝑈𝜕𝑌_ℓ𝐻(𝑥,𝑦,𝑡,𝑈,𝜕𝑈,𝜕𝑥𝜕𝑈−𝜕𝑦)𝑑Ω_{𝐼(𝑖)𝑠=1}12𝜀_𝑖𝑠Δ𝑥_𝑘1(𝜀_𝑖𝑠Δ̇_𝑖𝑠−𝑆_𝑖𝑠), (2.26) for ℓ=1,…,𝑛_𝑠.

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 𝐼_𝛼=_𝑅𝜕𝑈{𝐹(𝑥,𝑦,𝑡,𝑈,,𝜕𝑈𝜕𝑥)𝜕𝜕𝑦²𝑈𝜕𝑥²𝜕𝑈+𝐺(𝑥,𝑦,𝑡,𝑈,,𝜕𝑈𝜕𝑥)𝜕𝜕𝑦²𝑈𝜕𝑦²}𝑑𝑥𝑑𝑦, (2.27) where 𝑅 is the 𝛿-neighbourhood of an edge adjacent to the node 𝑣_𝑖 and 𝑈 is the recovered function of 𝑈. We take the 1𝐷 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, 𝑅=[−𝛿,𝛿]×[0,𝛼]. Then we have 𝐼_𝛼=^𝛼₀^𝛿_−𝛿𝜕𝑈{𝐹(𝜉,𝜂,𝑡,𝑈,,𝜕𝑈𝜕𝜉𝜕𝜂)cos²𝜕𝑈(𝜃)+𝐺(𝜉,𝜂,𝑡,𝑈,,𝜕𝑈𝜕𝜉𝜕𝜂)sin²𝜕(𝜃)}²𝑈𝜕𝜉²𝑑𝜉𝑑𝜂, (2.28) where 𝜃 is the rotation angle. Let 𝑝=(0,𝜂_𝛼), in which 1≤𝜂_𝛼≤𝛼 and 𝑈, the recovered function, defined by 𝑈(𝜉,𝜂_𝛼)=𝑈𝜉∉[−𝛿,𝛿],𝐻𝜉∈[−𝛿,𝛿], (2.29) where 𝐻 is the 1𝐷 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 𝛿→0, with some computation [19], we get 𝐼_𝛼=𝛼_𝑝𝜕𝑈𝑤(𝑝)(|𝜕𝜉_{𝑝∈Ω𝑅}−𝜕𝑈|𝜕𝜉_{𝑝∈Ω𝐿})𝜕𝑈×{𝐹(𝜉,𝜂,𝑡,𝑈,,𝜕𝑈𝜕𝜉𝜕𝜂)cos²𝜕𝑈(𝜃)+𝐺(𝜉,𝜂,𝑡,𝑈,,𝜕𝑈𝜕𝜉𝜕𝜂)sin²(𝜃)}|_𝑝, (2.30) where 𝑤(𝑝) is the 1𝐷 quadrature weight associated to the quadrature point 𝑝 and where Ω_𝑅 and Ω_𝐿 denote the right and left triangle, respectively (Figure 1).

Figure 1 

𝛿-neighbourhood of an edge adjacent to the node 𝑣𝑖.

2.5. General Equations of the GWMFE

Consider the global node 𝑣_𝑖, 𝑖=1,…,𝑛_𝑔, and its support Ω^𝑖. Assume that 𝑣_𝑖 is the 𝑘_𝑖𝑠th interpolation point in the element Ω_𝑖𝑠, 𝑠=1,…,𝐼(𝑖). Denoting by 𝑎_𝑖𝑠 the length of the 𝑠th edge adjacent to 𝑣_𝑖, the GWMFE equation (2.24) associated to the nodal amplitude is ^{𝐼(𝑖)𝑛}_𝑠=1^𝑝_{𝑘=1Ω𝑖𝑠}11+‖∇𝑈_𝑖𝑠‖²₂𝑏_𝑖𝑠𝑏_𝑘𝑑Ω𝑈̇^𝑘_𝑗+^𝐼(𝑖)3_{𝑠=1𝑁=1Ω𝑖𝑠}−11+‖∇𝑈_𝑖𝑠‖²₂𝐿_𝑘𝑁𝑏_𝑖𝑠𝜕𝑈̇𝑋𝜕𝑥𝑑Ω_ℓ𝑁+^𝐼(𝑖)3_{𝑠=1𝑁=1Ω𝑖𝑠}−11+‖∇𝑈_𝑖𝑠‖²₂𝐿_𝑘𝑁𝑏_𝑖𝑠𝜕𝑈̇𝑌𝜕𝑦𝑑Ω_ℓ𝑁=^𝐼(𝑖)_{𝑠=1Ω𝑖𝑠}11+‖∇𝑈_𝑖𝑠‖²₂𝑏_𝑖𝑠𝐹(𝑥,𝑦,𝑡,𝑈,𝜕𝑈,𝜕𝑥𝜕𝑈)𝜕𝜕𝑦²𝑈𝜕𝑥²+𝑑Ω_{𝐼(𝑖)𝑠=1Ω𝑖𝑠}11+‖∇𝑈_𝑖𝑠‖²₂𝑏_𝑖𝑠𝐺(𝑥,𝑦,𝑡,𝑈,𝜕𝑈,𝜕𝑥𝜕𝑈)𝜕𝜕𝑦²𝑈𝜕𝑦²+𝑑Ω_{𝐼(𝑖)𝑠=1Ω𝑖𝑠}11+‖∇𝑈_𝑖𝑠‖²₂𝑏_𝑖𝑠𝐻(𝑥,𝑦,𝑡,𝑈,𝜕𝑈,𝜕𝑥𝜕𝑈+𝜕𝑦)𝑑Ω_{𝐼(𝑖)𝑠=1𝑎𝑖𝑠0}𝑏_𝑘𝑖𝑠(𝜕𝑈|𝜕𝜉_{𝑝∈Ω𝑅}−𝜕𝑈|𝜕𝜉_{𝑝∈Ω𝐿})𝜕𝑈×{𝐹(𝜉,𝜂,𝑡,𝑈,,𝜕𝑈𝜕𝜉𝜕𝜂)cos²𝜕𝑈(𝜃)+𝐺(𝜉,𝜂,𝑡,𝑈,,𝜕𝑈𝜕𝜉𝜕𝜂)sin²(𝜃)}𝑑𝜂, (2.31) for 𝑖=1,…,𝑛_𝑔.

If 𝜉_ℓ=(𝑋_ℓ,𝑌_ℓ), ℓ=1,…,𝑛_𝑠 be a spatial node defining the 𝑘₁th vertex of the element surrounding the node 𝑖. For ℓ=1,…,𝑛_𝑠, (2.25) leads to

^{𝐼(𝑖)𝑛}_𝑠=1^𝑝_{𝑘=1Ω𝑖𝑠}−11+‖∇𝑈_𝑖𝑠‖²₂𝐿_𝑘1𝑏_𝑘𝜕𝑈̇𝑈𝜕𝑥𝑑Ω^𝑘_𝑗+^𝐼(𝑖)3_{𝑠=1𝑁=1Ω𝑖𝑠}11+‖∇𝑈_𝑖𝑠‖²₂𝐿_𝑘𝑁𝐿_𝑘1(𝜕𝑈)𝜕𝑥²𝑑Ω𝑋̇_ℓ𝑁+^𝐼(𝑖)3_{𝑠=1𝑁=1Ω𝑖𝑠}11+‖∇𝑈_𝑖𝑠‖²₂𝐿_𝑘𝑁𝐿_𝑘1𝜕𝑈𝜕𝑥𝜕𝑈̇𝑌𝜕𝑦𝑑Ω_ℓ𝑁+^𝐼(𝑖)_𝑠=112𝜀²_𝑖𝑠Δ𝑦_𝑘1Δ̇_𝑖𝑠=^𝐼(𝑖)_{𝑠=1Ω𝑖𝑠}−11+‖∇𝑈_𝑖𝑠‖²₂𝐿_𝑘1𝜕𝑈𝜕𝑥𝐹(𝑥,𝑦,𝑡,𝑈,𝜕𝑈,𝜕𝑥𝜕𝑈)𝜕𝜕𝑦²𝑈𝜕𝑥²+𝑑Ω_{𝐼(𝑖)𝑠=1Ω𝑖𝑠}−11+‖∇𝑈_𝑖𝑠‖²₂𝐿_𝑘1𝜕𝑈𝜕𝑥𝐺(𝑥,𝑦,𝑡,𝑈,𝜕𝑈,𝜕𝑥𝜕𝑈)𝜕𝜕𝑦²𝑈𝜕𝑦²+𝑑Ω_{𝐼(𝑖)𝑠=1Ω𝑖𝑠}−11+‖∇𝑈_𝑖𝑠‖²₂𝐿_𝑘1𝜕𝑈𝜕𝑥𝐻(𝑥,𝑦,𝑡,𝑈,𝜕𝑈,𝜕𝑥𝜕𝑈𝜕𝑦)𝑑Ω+^𝐼(𝑖)_𝑠=112𝜀_𝑖𝑠𝑆_𝑖𝑠Δ𝑦_𝑘1+^𝐼(𝑖)_𝑠=1^𝑎𝑖𝑠₀−𝐿_𝑘11{(2(𝜕𝑈|𝜕𝜉_{𝑝∈Ω𝑅}+𝜕𝑈|𝜕𝜉_{𝑝∈Ω𝐿}𝜕𝑈))cos(𝜃)−𝜕𝜂sin(𝜃)}𝑄𝑑𝜂, (2.32) and (2.26) leads to

^{𝐼(𝑖)𝑛}_𝑠=1^𝑝_{𝑘=1Ω𝑖𝑠}−11+‖∇𝑈_𝑖𝑠‖²₂𝐿_𝑘1𝑏_𝑘𝜕𝑈̇𝑈𝜕𝑦𝑑Ω^𝑘_𝑗+^𝐼(𝑖)3_{𝑠=1𝑁=1Ω𝑖𝑠}11+‖∇𝑈_𝑖𝑠‖²₂𝐿_𝑘𝑁𝐿_𝑘1𝜕𝑈𝜕𝑥𝜕𝑈̇𝑋𝜕𝑦𝑑Ω_ℓ𝑁+^𝐼(𝑖)3_{𝑠=1𝑁=1Ω𝑖𝑠}11+‖∇𝑈_𝑖𝑠‖²₂𝐿_𝑘𝑁𝐿_𝑘1(𝜕𝑈)𝜕𝑦²𝑑Ω𝑌̇_ℓ𝑁+^𝐼(𝑖)_𝑠=112𝜀²_𝑖𝑠Δ𝑥_𝑘1Δ̇_𝑖𝑠=^𝐼(𝑖)_{𝑠=1Ω𝑖𝑠}−11+‖∇𝑈_𝑖𝑠‖²₂𝐿_𝑘1𝜕𝑈𝜕𝑦𝐹(𝑥,𝑦,𝑡,𝑈,𝜕𝑈,𝜕𝑥𝜕𝑈)𝜕𝜕𝑦²𝑈𝜕𝑥²+𝑑Ω_{𝐼(𝑖)𝑠=1Ω𝑖𝑠}−11+‖∇𝑈_𝑖𝑠‖²₂𝐿_𝑘1𝜕𝑈𝜕𝑦𝐺(𝑥,𝑦,𝑡,𝑈,𝜕𝑈,𝜕𝑥𝜕𝑈)𝜕𝜕𝑦²𝑈𝜕𝑦²+𝑑Ω_{𝐼(𝑖)𝑠=1Ω𝑖𝑠}−11+‖∇𝑈_𝑖𝑠‖²₂𝐿_𝑘1𝜕𝑈𝜕𝑦𝐻(𝑥,𝑦,𝑡,𝑈,𝜕𝑈,𝜕𝑥𝜕𝑈𝜕𝑦)𝑑Ω+^𝐼(𝑖)_𝑠=112𝜀_𝑖𝑠𝑆_𝑖𝑠Δ𝑥_𝑘1+^𝐼(𝑖)_𝑠=1^𝑎𝑖𝑠₀−𝐿_𝑘11{(2(𝜕𝑈|𝜕𝜉_{𝑝∈Ω𝑅}+𝜕𝑈|𝜕𝜉_{𝑝∈Ω𝐿}𝜕𝑈))cos(𝜃)−𝜕𝜂sin(𝜃)}𝑄𝑑𝜂, (2.33) for ℓ=1,…,𝑛_𝑠, where 𝑄=^𝐼(𝑖)_𝑠=1^𝑎𝑖𝑠₀𝑏_𝑘𝑖𝑠(𝜕𝑈|𝜕𝜉_{𝑝∈Ω𝑅}−𝜕𝑈|𝜕𝜉_{𝑝∈Ω𝐿})𝜕𝑈×{𝐹(𝜉,𝜂,𝑡,𝑈,,𝜕𝑈𝜕𝜉𝜕𝜂)cos²𝜕𝑈(𝜃)+𝐺(𝜉,𝜂,𝑡,𝑈,,𝜕𝑈𝜕𝜉𝜕𝜂)sin²(𝜃)}𝑑𝜂. (2.34)

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: 𝑀(𝑡,𝑌)𝑌̇=𝐹(𝑡,𝑌), (2.35) 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 {𝑡₀,𝑡₁,…,𝑡_𝑘}, (3.1) where 𝑡_𝑘≤𝑇𝑓.

Now, we apply the refinement scheme at each time step [17], for example, on the first time step [𝑡₀,𝑡₁]. Set Δ𝑡=((𝑡₁−𝑡₀)/𝑙₀), 𝑙₀∈𝑁, 𝑡^{0+(𝑙0−𝑖)Δ𝑡}=𝑡₀+(𝑙₀−𝑖)Δ𝑡,𝑖=𝑙₀,𝑙₀−1,…,0. (3.2) Thus the time integration from 𝑡₀ to 𝑡₁ involves 𝑙₀ sub step such that

𝑡=𝑡₀→𝑡=𝑡₀+Δ𝑡→𝑡=𝑡₀+2Δ𝑡→,…,𝑡=𝑡₀+𝑙₀Δ𝑡=𝑡₁ (3.3) 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 (𝑘−1)th ODE system at 𝑡_𝑘.

Generally, suppose that we are at time level 𝑡=𝑡_𝑘 and want to move toward 𝑡=𝑡_𝑘+1, similarly, consider Δ𝑡=(𝑡_𝑘+1−𝑡_𝑘)/𝑙_𝑘 such that 𝑡^{𝑘+(𝑙𝑘−𝑖)Δ𝑡}=𝑡_𝑘+(𝑙_𝑘−𝑖)Δ𝑡,𝑖=𝑙_𝑘,𝑙_𝑘−1,…,0. (3.4) The values of 𝑈(𝑥,𝑦,𝑡_𝑘+1), 𝑥(𝑡_𝑘+1) and 𝑦(𝑡_𝑘+1) are obtained by solving 𝑙_𝑘 ODE system 𝑡=𝑡_𝑘→𝑡=𝑡^𝑘+Δ𝑡→𝑡=𝑡^𝑘+2Δ𝑡→,…,𝑡=𝑡^{𝑘+𝑙𝑘Δ𝑡}=𝑡_𝑘+1. (3.5) This process continues until 𝑡_𝑘+1≤𝑇𝑓.

3.1. Algorithm

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

Step 1. Set 𝑥_𝑖=𝑖/𝑛, 𝑖=1,2,…,𝑛, and 𝑦_𝑗=𝑗/𝑚, 𝑗=1,2,…,𝑚 and obtain 𝑈(𝑥,𝑦,0) with initial value 𝑥_𝑖, 𝑦_𝑗, 𝑌⁰=[𝑈⁰,𝑥⁰,𝑦⁰] and set 𝑘=0.

Step 2. Set Δ𝑡=(𝑡_𝑘+1−𝑡_𝑘)/𝑙_𝑘, 𝑙_𝑘∈𝑁 and 𝑠=1.

Step 3. Solve ODE system (2.35) at time level 𝑡=𝑡^{𝑘+𝑠Δ𝑡} as follows: [𝑡^{𝑘+𝑠Δ𝑡},𝑌^{𝑘+𝑠Δ𝑡}] = ODE15S(@Function,[𝑡^{𝑘+(𝑠−1)Δ𝑡},𝑡^{𝑘+𝑠Δ𝑡}],𝑌^{𝑘+(𝑠−1)Δ𝑡},Option),the initial value of which in this step is obtained by solving ODE system (2.35) at 𝑡=𝑡^{𝑘+(𝑠−1)Δ𝑡}.

Step 4. Set 𝑈^{𝑘+(𝑠−1)Δ𝑡}=𝑈^{𝑘+𝑠Δ𝑡}, 𝑥^{𝑘+(𝑠−1)Δ𝑡}=𝑥^{𝑘+𝑠Δ𝑡}, 𝑦^{𝑘+(𝑠−1)Δ𝑡}=𝑦^{𝑘+𝑠Δ𝑡}, Δ^{𝑘+(𝑠−1)Δ𝑡}=Δ^{𝑘+𝑠Δ𝑡}, 𝜀^{𝑘+(𝑠−1)Δ𝑡}=𝜀^{𝑘+𝑠Δ𝑡}, 𝑆^{𝑘+(𝑠−1)Δ𝑡}=𝑆^{𝑘+𝑠Δ𝑡}, and 𝑠=𝑠+1.

Step 5. If 𝑠≤𝑙_𝑘, then go to Step 3, or else 𝑘=𝑘+1.

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, Atol=Rtol=10⁻⁸. 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 𝐶₃=10⁻⁵.

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 𝜕𝑢𝜕𝜕𝑡=𝑝(²𝑢𝜕𝑥²+𝜕²𝑢𝜕𝑦²)−𝑢(𝜕𝑢+𝜕𝑥𝜕𝑢𝜕𝑦),0≤𝑥,𝑦≤1,𝑡≥0. (4.1) We assume that initial and Dirichlet boundary condition are chosen to correspond to the analytic solution 1𝑢(𝑥,𝑦,𝑡)=1+𝑒^{(𝑦+𝑥−𝑡)/2𝑝}. (4.2) This problem is solved from 𝑡=0 to 𝑡=2 with 𝑝=0.025 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 𝑡=0 to 𝑡=2 seconds with quadratic elements at some cases. Figures 3–6 present nodes movement and their solutions with quadratic elements for 𝑡=0 second, 𝑡=0.5 second, 𝑡=1 second and 𝑡=2 seconds, respectively.

(a)

(b)

(c)

(d)

(e)

(f)

Figure 2 

Nodes movement for Burger equation from 𝑡=0 to 𝑡=2 at some seconds.

(a)

(b)

Figure 3 

Nodes distribution and solution of Burgers equation with quadratic elements for uniform initial mesh.

(a)

(b)

Figure 4 

Nodes movement and related solution of Burgers equation with quadratic elements at 𝑡=0.5 second, 𝐶1=10−3, 𝐶2=10−4, 𝐶3=10−5.

(a)

(b)

Figure 5 

Nodes movement and related solution of Burgers equation with quadratic elements at 𝑡=1 second, 𝐶1=10−3, 𝐶2=10−4, 𝐶3=10−5.

(a)

(b)

Figure 6 

Nodes movement and related solution of Burgers equation with quadratic elements at 𝑡=2 seconds, 𝐶1=10−3, 𝐶2=10−4, 𝐶3=10−5.

Let us consider to the average error, 𝐸_ave=∑^𝑛𝑠_ℓ=1|𝑈(𝜉_ℓ)−𝑢(𝜉_ℓ)|𝑛_𝑠, (4.3) 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 (𝐸_ave) for GWMFE at 𝑡=2 seconds seconds with quadratic element and 𝐶₁=10⁻³, 𝐶₂=10⁻⁴, and 𝐶₃=10⁻⁵ 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.

Table 1 

CPU time, number of function evaluation (NFE), and average error (𝐸_𝐚𝐯𝐞) for GWMFE at t = 2 seconds with 𝐶₁=10⁻³, 𝐶₂=10⁻⁴, and 𝐶₃=10⁻⁵ in some meshes.

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.