The present paper conjectures a topological condition which classifies a T-spline into standard, semi-standard and non-standard. We also provide the basic framework to prove the conjecture on the classification of semi-standard T-splines and give the proof for the semi-standard of bi-degree (1, d) and (d, 1) T-splines.

In this paper, we present a unified finite volume method preserving discrete maximum principle (DMP) for the conjugate heat transfer problems with general interface conditions. We prove the existence of the numerical solution and the DMP-preserving property. Numerical experiments show that the nonlinear iteration numbers of the scheme in [24] increase rapidly when the interfacial coefficients decrease to zero. In contrast, the nonlinear iteration numbers of the unified scheme do not increase when the interfacial coefficients decrease to zero, which reveals that the unified scheme is more robust than the scheme in [24]. The accuracy and DMP-preserving property of the scheme are also verified in the numerical experiments.

In this article, we consider the primal-dual path-following method and the trust-region updating strategy for the standard linear programming problem. For the rank-deficient problem with the small noisy data, we also give the preprocessing method based on the QR decomposition with column pivoting. Then, we prove the global convergence of the new method when the initial point is strictly primal-dual feasible. Finally, for some rankdeficient problems with or without the small noisy data from the NETLIB collection, we compare it with other two popular interior-point methods, i.e. the subroutine pathfollow.m and the built-in subroutine linprog.m of the MATLAB environment. Numerical results show that the new method is more robust than the other two methods for the rank-deficient problem with the small noise data.

In this paper, we consider the knot placement problem in B-spline curve approximation. A novel two-stage framework is proposed for addressing this problem. In the first step, the l_∞,1-norm model is introduced for the sparse selection of candidate knots from an initial knot vector. By this step, the knot number is determined. In the second step, knot positions are formulated into a nonlinear optimization problem and optimized by a global optimization algorithm — the differential evolution algorithm (DE). The candidate knots selected in the first step are served for initial values of the DE algorithm. Since the candidate knots provide a good guess of knot positions, the DE algorithm can quickly converge. One advantage of the proposed algorithm is that the knot number and knot positions are determined automatically. Compared with the current existing algorithms, the proposed algorithm finds approximations with smaller fitting error when the knot number is fixed in advance. Furthermore, the proposed algorithm is robust to noisy data and can handle with few data points. We illustrate with some examples and applications.

In this paper, ETD3-Padé and ETD4-Padé Galerkin finite element methods are proposed and analyzed for nonlinear delayed convection-diffusion-reaction equations with Dirichlet boundary conditions. An ETD-based RK is used for time integration of the corresponding equation. To overcome a well-known difficulty of numerical instability associated with the computation of the exponential operator, the Padé approach is used for such an exponential operator approximation, which in turn leads to the corresponding ETD-Padé schemes. An unconditional L² numerical stability is proved for the proposed numerical schemes, under a global Lipshitz continuity assumption. In addition, optimal rate error estimates are provided, which gives the convergence order of ^O(k³ + h^r) (ETD3- Padé) or O(k⁴ + h^r) (ETD4-Padé) in the L² norm, respectively. Numerical experiments are presented to demonstrate the robustness of the proposed numerical schemes.

In this paper, we present a two-grid discretization scheme for semilinear parabolic integro-differential equations by H¹-Galerkin mixed finite element methods. We use the lowest order Raviart-Thomas mixed finite elements and continuous linear finite element for spatial discretization, and backward Euler scheme for temporal discretization. Firstly, a priori error estimates and some superclose properties are derived. Secondly, a two-grid scheme is presented and its convergence is discussed. In the proposed two-grid scheme, the solution of the nonlinear system on a fine grid is reduced to the solution of the nonlinear system on a much coarser grid and the solution of two symmetric and positive definite linear algebraic equations on the fine grid and the resulting solution still maintains optimal accuracy. Finally, a numerical experiment is implemented to verify theoretical results of the proposed scheme. The theoretical and numerical results show that the two-grid method achieves the same convergence property as the one-grid method with the choice h = H ².

A combined scheme of the improved two-grid technique with the block-centered finite difference method is constructed and analyzed to solve the nonlinear time-fractional parabolic equation. This method is considered where the nonlinear problem is solved only on a coarse grid of size H and two linear problems based on the coarse-grid solutions and one Newton iteration is considered on a fine grid of size h. We provide the rigorous error estimate, which demonstrates that our scheme converges with order O(Δt^2-α+h²+H⁴) on non-uniform rectangular grid. This result indicates that the improved two-grid method can obtain asymptotically optimal approximation as long as the mesh sizes satisfy h=O(H²). Finally, numerical tests confirm the theoretical results of the presented method.

This is one of our series works on numerical methods for mean-field forward backward stochastic differential equations (MFBSDEs). In this work, we propose an explicit multistep scheme for MFBSDEs which is easy to implement, and is of high order rate of convergence. Rigorous error estimates of the proposed multistep scheme are presented. Numerical experiments are carried out to show the efficiency and accuracy of the proposed scheme.

This note introduces a method for sampling Ising models with mixed boundary conditions. As an application of annealed importance sampling and the Swendsen-Wang algorithm, the method adopts a sequence of intermediate distributions that keeps the temperature fixed but turns on the boundary condition gradually. The numerical results show that the variance of the sample weights is relatively small.

A common way to handle the Tikhonov regularization method for the first kind Fredholm integral equations, is first to discretize and then to work with the final linear system. This unavoidably inflicts discretization errors which may lead to disastrous results, especially when a quadrature rule is used. We propose to regularize directly the integral equation resulting in a continuous Tikhonov problem. The Tikhonov problem is reduced to a simple least squares problem by applying the Golub-Kahan bidiagonalization (GKB) directly to the integral operator. The regularization parameter and the iteration index are determined by the discrepancy principle approach. Moreover, we study the discrete version of the proposed method resulted from numerical evaluating the needed integrals. Focusing on the nodal values of the solution results in a weighted version of GKB-Tikhonov method for linear systems arisen from the Nyström discretization. Finally, we use numerical experiments on a few test problems to illustrate the performance of our algorithms.

This paper concerns the reconstruction of a scalar coefficient of a second-order elliptic equation in divergence form posed on a bounded domain from internal data. This problem finds applications in multi-wave imaging, greedy methods to approximate parameterdependent elliptic problems, and image treatment with partial differential equations. We first show that the inverse problem for smooth coefficients can be rewritten as a linear transport equation. Assuming that the coefficient is known near the boundary, we study the well-posedness of associated transport equation as well as its numerical resolution using discontinuous Galerkin method. We propose a regularized transport equation that allow us to derive rigorous convergence rates of the numerical method in terms of the order of the polynomial approximation as well as the regularization parameter. We finally provide numerical examples for the inversion assuming a lower regularity of the coefficient, and using synthetic data.

Flexible GMRES (FGMRES) is a variant of preconditioned GMRES, which changes preconditioners at every Arnoldi step. GMRES often has to be restarted in order to save storage and reduce orthogonalization cost in the Arnoldi process. Like restarted GMRES, FGMRES may also have to be restarted for the same reason. A major disadvantage of restarting is the loss of convergence speed. In this paper, we present a heavy ball flexible GMRES method, aiming to recoup some of the loss in convergence speed in the restarted flexible GMRES while keep the benefit of limiting memory usage and controlling orthogonalization cost. Numerical tests often demonstrate superior performance of the proposed heavy ball FGMRES to the restarted FGMRES.

By employing EQ₁^rot nonconforming finite element, the numerical approximation is presented for multi-term time-fractional mixed sub-diffusion and diffusion-wave equation on anisotropic meshes. Comparing with the multi-term time-fractional sub-diffusion equation or diffusion-wave equation, the mixed case contains a special time-space coupled derivative, which leads to many difficulties in numerical analysis. Firstly, a fully discrete scheme is established by using nonconforming finite element method (FEM) in spatial direction and L1 approximation coupled with Crank-Nicolson (L1-CN) scheme in temporal direction. Furthermore, the fully discrete scheme is proved to be unconditional stable. Besides, convergence and superclose results are derived by using the properties of EQ₁^rot nonconforming finite element. What's more, the global superconvergence is obtained via the interpolation postprocessing technique. Finally, several numerical results are provided to demonstrate the theoretical analysis on anisotropic meshes.

In this paper, we systematically construct two classes of structure-preserving schemes with arbitrary order of accuracy for canonical Hamiltonian systems. The one class is the symplectic scheme, which contains two new families of parameterized symplectic schemes that are derived by basing on the generating function method and the symmetric composition method, respectively. Each member in these schemes is symplectic for any fixed parameter. A more general form of generating functions is introduced, which generalizes the three classical generating functions that are widely used to construct symplectic algorithms. The other class is a novel family of energy and quadratic invariants preserving schemes, which is devised by adjusting the parameter in parameterized symplectic schemes to guarantee energy conservation at each time step. The existence of the solutions of these schemes is verified. Numerical experiments demonstrate the theoretical analysis and conservation of the proposed schemes.

We consider a nonlinear stochastic Volterra integral equation with time-dependent delay and the corresponding Euler-Maruyama method in this paper. Strong convergence rate (at fixed point) of the corresponding Euler-Maruyama method is obtained when coefficients f and g both satisfy local Lipschitz and linear growth conditions. An example is provided to interpret our conclusions. Our result generalizes and improves the conclusion in[J. Gao, H. Liang, S. Ma, Strong convergence of the semi-implicit Euler method for nonlinear stochastic Volterra integral equations with constant delay, Appl. Math. Comput., 348 (2019), 385-398.]

Numerical integration poses greater challenges in Galerkin meshless methods than finite element methods owing to the non-polynomial feature of meshless shape functions. The reproducing kernel gradient smoothing integration (RKGSI) is one of the optimal numerical integration techniques in Galerkin meshless methods with minimum integration points. In this paper, properties, quadrature rules and the effect of the RKGSI on meshless methods are analyzed. The existence, uniqueness and error estimates of the solution of Galerkin meshless methods under numerical integration with the RKGSI are established. A procedure on how to choose quadrature rules to recover the optimal convergence rate is presented.

In this paper, we apply the Wasserstein-Fisher-Rao (WFR) metric from the unbalanced optimal transport theory to the earthquake location problem. Compared with the quadratic Wasserstein (W₂) metric from the classical optimal transport theory, the advantage of this method is that it retains the important amplitude information as a new constraint, which avoids the problem of the degeneration of the optimization objective function near the real earthquake hypocenter and origin time. As a result, the deviation of the global minimum of the optimization objective function based on the WFR metric from the true solution can be much smaller than the results based on the W₂ metric when there exists strong data noise. Thus, we develop an accurate earthquake location method under strong data noise. Many numerical experiments verify our conclusions.

In this paper, a new hybridized mixed formulation of weak Galerkin method is studied for a second order elliptic problem. This method is designed by approximate some operators with discontinuous piecewise polynomials in a shape regular finite element partition. Some discrete inequalities are presented on discontinuous spaces and optimal order error estimations are established. Some numerical results are reported to show super convergence and confirm the theory of the mixed weak Galerkin method.

In this paper, we consider the generalized Nash equilibrium with shared constraints in the stochastic environment, and we call it the stochastic generalized Nash equilibrium. The stochastic variational inequalities are employed to solve this kind of problems, and the expected residual minimization model and the conditional value-at-risk formulations defined by the residual function for the stochastic variational inequalities are discussed. We show the risk for different kinds of solutions for the stochastic generalized Nash equilibrium by the conditional value-at-risk formulations. The properties of the stochastic quadratic generalized Nash equilibrium are shown. The smoothing approximations for the expected residual minimization formulation and the conditional value-at-risk formulation are employed. Moreover, we establish the gradient consistency for the measurable smoothing functions and the integrable functions under some suitable conditions, and we also analyze the properties of the formulations. Numerical results for the applications arising from the electricity market model illustrate that the solutions for the stochastic generalized Nash equilibrium given by the ERM model have good properties, such as robustness, low risk and so on.

In this paper, we present the backward stochastic Taylor expansions for a Ito process, including backward Ito-Taylor expansions and backward Stratonovich-Taylor expansions. We construct the general full implicit strong Taylor approximations (including Ito-Taylor and Stratonovich-Taylor schemes) with implicitness in both the deterministic and the stochastic terms for the stiff stochastic differential equations (SSDE) by employing truncations of backward stochastic Taylor expansions. We demonstrate that these schemes will converge strongly with corresponding order 1, 2, 3, . . . . Mean-square stability has been investigated for full implicit strong Stratonovich-Taylor scheme with order 2, and it has larger meansquare stability region than the explicit and the semi-implicit strong Stratonovich-Taylor schemes with order 2. We can improve the stability of simulations considerably without too much additional computational effort by using our full implicit schemes. The full implicit strong Taylor schemes allow a larger range of time step sizes than other schemes and are suitable for SSDE with stiffness on both the drift and the diffusion terms. Our numerical experiment show these points.

A nonlinear parabolic system is derived to describe compressible miscible displacement in a porous medium. The concentration equation is treated by a mixed finite element method with characteristics (CMFEM) and the pressure equation is treated by a parabolic mixed finite element method (PMFEM). Two-grid algorithm is considered to linearize nonlinear coupled system of two parabolic partial differential equations. Moreover, the L^q error estimates are conducted for the pressure, Darcy velocity and concentration variables in the two-grid solutions. Both theoretical analysis and numerical experiments are presented to show that the two-grid algorithm is very effective.

To reduce the computational cost, we propose a regularizing modified LevenbergMarquardt scheme via multiscale Galerkin method for solving nonlinear ill-posed problems. Convergence results for the regularizing modified Levenberg-Marquardt scheme for the solution of nonlinear ill-posed problems have been proved. Based on these results, we propose a modified heuristic parameter choice rule to terminate the regularizing modified Levenberg-Marquardt scheme. By imposing certain conditions on the noise, we derive optimal convergence rates on the approximate solution under special source conditions. Numerical results are presented to illustrate the performance of the regularizing modified Levenberg-Marquardt scheme under the modified heuristic parameter choice.

We study from a numerical point of view a multidimensional problem involving a viscoelastic body with two porous structures. The mechanical problem leads to a linear system of three coupled hyperbolic partial differential equations. Its corresponding variational formulation gives rise to three coupled parabolic linear equations. An existence and uniqueness result, and an energy decay property, are recalled. Then, fully discrete approximations are introduced using the finite element method and the implicit Euler scheme. A discrete stability property and a priori error estimates are proved, from which the linear convergence of the algorithm is derived under suitable additional regularity conditions. Finally, some numerical simulations are performed in one and two dimensions to show the accuracy of the approximation and the behaviour of the solution.

In this paper, based on discrete gradient, a dissipation-preserving integrator for weakly dissipative perturbations of oscillatory Hamiltonian system is established. The solution of this system is a damped nonlinear oscillator. Basically, lots of nonlinear oscillatory mechanical systems including frictional forces lend themselves to this approach. The new integrator gives a discrete analogue of the dissipation property of the original system. Meanwhile, since the integrator is based on the variation-of-constants formula for oscillatory systems, it preserves the oscillatory structure of the system. Some properties of the new integrator are derived. The convergence is analyzed for the implicit iterations based on the discrete gradient integrator, and it turns out that the convergence of the implicit iterations based on the new integrator is independent of ||M||, where M governs the main oscillation of the system and usually ||M|| ≫ 1. This significant property shows that a larger stepsize can be chosen for the new schemes than that for the traditional discrete gradient integrators when applied to the oscillatory Hamiltonian system. Numerical experiments are carried out to show the effectiveness and efficiency of the new integrator in comparison with the traditional discrete gradient methods in the scientific literature.

In this paper, we consider the Euler-Maruyama method for a class of stochastic Volterra integral equations (SVIEs). It is known that the strong convergence order of the Euler-Maruyama method is $\frac{1}{2}$. However, the strong superconvergence order 1 can be obtained for a class of SVIEs if the kernels σi(t, t) = 0 for i = 1 and 2; otherwise, the strong convergence order is $\frac{1}{2}$ . Moreover, the theoretical results are illustrated by some numerical examples.

In this paper, we derive and analyse waveform relaxation (WR) methods for solving differential equations evolving on a Lie-group. We present both continuous-time and discrete-time WR methods and study their convergence properties. In the discrete-time case, the novel methods are constructed by combining WR methods with Runge-KuttaMunthe-Kaas (RK-MK) methods. The obtained methods have both advantages of WR methods and RK-MK methods, which simplify the computation by decoupling strategy and preserve the numerical solution of Lie-group equations on a manifold. Three numerical experiments are given to illustrate the feasibility of the new WR methods.

In this paper, we develop an active set identification technique. By means of the active set technique, we present an active set adaptive monotone projected Barzilai-Borwein method (ASAMPBB) for solving nonnegative matrix factorization (NMF) based on the alternating nonnegative least squares framework, in which the Barzilai-Borwein (BB) step sizes can be adaptively picked to get meaningful convergence rate improvements. To get optimal step size, we take into account of the curvature information. In addition, the larger step size technique is exploited to accelerate convergence of the proposed method. The global convergence of the proposed method is analysed under mild assumption. Finally, the results of the numerical experiments on both synthetic and real-world datasets show that the proposed method is effective.

The main objective of this paper is to present an efficient structure-preserving scheme, which is based on the idea of the scalar auxiliary variable approach, for solving the twodimensional space-fractional nonlinear Schrödinger equation. First, we reformulate the equation as an canonical Hamiltonian system, and obtain a new equivalent system via introducing a scalar variable. Then, we construct a semi-discrete energy-preserving scheme by using the Fourier pseudo-spectral method to discretize the equivalent system in space direction. After that, applying the Crank-Nicolson method on the temporal direction gives a linearly-implicit scheme in the fully-discrete version. As expected, the proposed scheme can preserve the energy exactly and more efficient in the sense that only decoupled equations with constant coefficients need to be solved at each time step. Finally, numerical experiments are provided to demonstrate the efficiency and conservation of the scheme.

This note introduces the double flip move to accelerate the Swendsen-Wang algorithm for Ising models with mixed boundary conditions below the critical temperature. The double flip move consists of a geometric flip of the spin lattice followed by a spin value flip. Both symmetric and approximately symmetric models are considered. We prove the detailed balance of the double flip move and demonstrate its empirical efficiency in mixing.

In this paper, we study the well-posedness and solution regularity of a multi-term variable-order time-fractional diffusion equation, and then develop an optimal Galerkin finite element scheme without any regularity assumption on its true solution. We show that the solution regularity of the considered problem can be affected by the maximum value of variable-order at initial time t = 0. More precisely, we prove that the solution to the multi-term variable-order time-fractional diffusion equation belongs to C²([0,T ]) in time provided that the maximum value has an integer limit near the initial time and the data has sufficient smoothness, otherwise the solution exhibits the same singular behavior like its constant-order counterpart. Based on these regularity results, we prove optimalorder convergence rate of the Galerkin finite element scheme. Furthermore, we develop an efficient parallel-in-time algorithm to reduce the computational costs of the evaluation of multi-term variable-order fractional derivatives. Numerical experiments are put forward to verify the theoretical findings and to demonstrate the efficiency of the proposed scheme.

This paper revisits the classical problem "Can we hear the density of a string?", which can be formulated as an inverse spectral problem for a Sturm-Liouville operator. Instead of inverting the map from density to spectral data directly, we propose a novel method to reconstruct the density based on inverting a sequence of trace formulas which bridge the density and its spectral data clearly in terms of a series of nonlinear integral equations. Numerical experiments are presented to verify the validity and effectiveness of the proposed numerical algorithm. The impact of different parameters involved in the algorithm is also discussed.

Direct and inverse problems for the scattering of cracks with mixed oblique derivative boundary conditions from the incident plane wave are considered, which describe the scattering phenomenons such as the scattering of tidal waves by spits or reefs. The solvability of the direct scattering problem is proven by using the boundary integral equation method. In order to show the equivalent boundary integral system is Fredholm of index zero, some relationships concerning the tangential potential operator is used. Due to the mixed oblique derivative boundary conditions, we cannot employ the factorization method in a usual manner to reconstruct the cracks. An alternative technique is used in the theoretical analysis such that the far field operator can be factorized in an appropriate form and fulfills the range identity theorem. Finally, we present some numerical examples to demonstrate the feasibility and effectiveness of the factorization method.

This paper concerns the convex optimal control problem governed by multiscale elliptic equations with arbitrarily rough L^∞ coefficients, which has not only complex coupling between nonseparable scales and nonlinearity, but also important applications in composite materials and geophysics. We use one of the recently developed numerical homogenization techniques, the so-called Rough Polyharmonic Splines (RPS) and its generalization (GRPS) for the efficient resolution of the elliptic operator on the coarse scale. Those methods have optimal convergence rate which do not rely on the regularity of the coefficients nor the concepts of scale-separation or periodicity. As the iterative solution of the nonlinearly coupled OCP-OPT formulation for the optimal control problem requires solving the corresponding (state and co-state) multiscale elliptic equations many times with different right hand sides, numerical homogenization approach only requires one-time pre-computation on the fine scale and the following iterations can be done with computational cost proportional to coarse degrees of freedom. Numerical experiments are presented to validate the theoretical analysis.

A simple criterion is studied for the first time for identifying the discrete energy dissipation of the Crank-Nicolson scheme for Maxwell's equations in a Cole-Cole dispersive medium. Several numerical formulas that approximate the time fractional derivatives are investigated based on this criterion, including the L1 formula, the fractional BDF-2, and the shifted fractional trapezoidal rule (SFTR). Detailed error analysis is provided within the framework of time domain mixed finite element methods for smooth solutions. The convergence results and discrete energy dissipation law are confirmed by numerical tests. For nonsmooth solutions, the method SFTR can still maintain the optimal convergence order at a positive time on uniform meshes. Authors believe this is the first appearance that a second-order time-stepping method can restore the optimal convergence rate for Maxwell's equations in a Cole-Cole dispersive medium regardless of the initial singularity of the solution.

Two nonconforming penalty methods for the two-dimensional stationary Navier-Stokes equations are studied in this paper. These methods are based on the weakly continuous P¹ vector fields and the locally divergence-free (LDF) finite elements, which respectively penalize local divergence and are discontinuous across edges. These methods have no penalty factors and avoid solving the saddle-point problems. The existence and uniqueness of the velocity solution are proved, and the optimal error estimates of the energy norms and L²-norms are obtained. Moreover, we propose unified pressure recovery algorithms and prove the optimal error estimates of L²-norm for pressure. We design a unified iterative method for numerical experiments to verify the correctness of the theoretical analysis.

In this paper, we present a method for generating Bézier surfaces from the boundary information based on a general second order functional and a third order functional associated with the triharmonic equation. By solving simple linear equations, the internal control points of the resulting Bézier surface can be obtained as linear combinations of the given boundary control points. This is a generalization of previous works on Plateau-Bézier problem, harmonic, biharmonic and quasi-harmonic Bézier surfaces. Some representative examples show the effectiveness of the presented method.

A low order nonconforming mixed finite element method (FEM) is established for the fully coupled non-stationary incompressible magnetohydrodynamics (MHD) problem in a bounded domain in 3D. The lowest order finite elements on tetrahedra or hexahedra are chosen to approximate the pressure, the velocity field and the magnetic field, in which the hydrodynamic unknowns are approximated by inf-sup stable finite element pairs and the magnetic field by H¹(?)-conforming finite elements, respectively. The existence and uniqueness of the approximate solutions are shown. Optimal order error estimates of L²(H¹)-norm for the velocity field, L²(L²)-norm for the pressure and the broken L²(H¹)-norm for the magnetic field are derived.

Two-phase image segmentation is a fundamental task to partition an image into foreground and background. In this paper, two types of nonconvex and nonsmooth regularization models are proposed for basic two-phase segmentation. They extend the convex regularization on the characteristic function on the image domain to the nonconvex case, which are able to better obtain piecewise constant regions with neat boundaries. By analyzing the proposed non-Lipschitz model, we combine the proximal alternating minimization framework with support shrinkage and linearization strategies to design our algorithm. This leads to two alternating strongly convex subproblems which can be easily solved. Similarly, we present an algorithm without support shrinkage operation for the nonconvex Lipschitz case. Using the Kurdyka-Lojasiewicz property of the objective function, we prove that the limit point of the generated sequence is a critical point of the original nonconvex nonsmooth problem. Numerical experiments and comparisons illustrate the effectiveness of our method in two-phase image segmentation.

In this paper, we apply local discontinuous Galerkin (LDG) methods for pattern formation dynamical model in polymerizing actin flocks. There are two main difficulties in designing effective numerical solvers. First of all, the density function is non-negative, and zero is an unstable equilibrium solution. Therefore, negative density values may yield blow-up solutions. To obtain positive numerical approximations, we apply the positivitypreserving (PP) techniques. Secondly, the model may contain stiff source. The most commonly used time integration for the PP technique is the strong-stability-preserving Runge-Kutta method. However, for problems with stiff source, such time discretizations may require strictly limited time step sizes, leading to large computational cost. Moreover, the stiff source any trigger spurious filament polarization, leading to wrong numerical approximations on coarse meshes. In this paper, we combine the PP LDG methods with the semi-implicit Runge-Kutta methods. Numerical experiments demonstrate that the proposed method can yield accurate numerical approximations with relatively large time steps.

We propose an accurate and energy-stable parametric finite element method for solving the sharp-interface continuum model of solid-state dewetting in three-dimensional space. The model describes the motion of the film/vapor interface with contact line migration and is governed by the surface diffusion equation with proper boundary conditions at the contact line. We present a weak formulation for the problem, in which the contact angle condition is weakly enforced. By using piecewise linear elements in space and backward Euler method in time, we then discretize the formulation to obtain a parametric finite element approximation, where the interface and its contact line are evolved simultaneously. The resulting numerical method is shown to be well-posed and unconditionally energystable. Furthermore, the numerical method is generalized to the case of anisotropic surface energies in the Riemannian metric form. Numerical results are reported to show the convergence and efficiency of the proposed numerical method as well as the anisotropic effects on the morphological evolution of thin films in solid-state dewetting.