Abstract: Publication date: Available online 26 August 2016
Source:Applied Numerical Mathematics
Author(s): Federica Porta, Anastasia Cornelio, Valeria Ruggiero
It is well known that there is a strong connection between time integration and convex optimization. In this work, inspired by the equivalence between the forward Euler scheme and the gradient descent method, we broaden our analysis to the family of Runge-Kutta methods and show that they enjoy a natural interpretation as first-order optimization algorithms. The strategies intrinsically suggested by Runge-Kutta methods are exploited in order to detail novel proposal for either scaling or preconditioning gradient-like approaches, whose convergence is ensured by the stability condition for Runge-Kutta schemes. The theoretical analysis is supported by the numerical experiments carried out on some test problems arising from suitable applications where the proposed techniques can be efficiently employed.

Abstract: Publication date: Available online 26 August 2016
Source:Applied Numerical Mathematics
Author(s): T. Diogo, P.M. Lima, A. Pedas, G. Vainikko
This work is concerned with the construction and analysis of high order numerical methods for solving initial value problems for linear Volterra integro-differential equations with different types of singularities. Using an integral reformulation of the initial value problem, a smoothing transformation is applied so that the exact solution of the resulting equation does not contain any singularities in its derivatives up to a certain order. After that, the regularized equation is solved by a piecewise polynomial collocation method on a uniform or mildly graded grid. Finally, the obtained spline approximations can be used to define (typically non-polynomial) approximations for the initial value problem. The theoretical results are tested by some numerical examples.

Abstract: Publication date: Available online 26 August 2016
Source:Applied Numerical Mathematics
Author(s): Haiyun Dong, Maojun Li
In this paper, we present a class of high order reconstructed central discontinuous Galerkin-finite element methods for the fully nonlinear weakly dispersive Green-Naghdi model, which describes a large spectrum of shallow water waves. In the proposed methods, we first reformulate the Green-Naghdi model into conservation laws coupled with an elliptic equation, and then discretize the conservation laws with reconstructed central discontinuous Galerkin methods and the elliptic equation with continuous FE methods. The reconstructed central discontinuous Galerkin methods can be viewed as a class of fast central discontinuous Galerkin methods, in which we replace the standard formula for the numerical solution defined on the dual mesh in the central discontinuous Galerkin method with a projection equation in the L 2 sense. The proposed methods reduce the computational cost of the traditional methods by nearly half but still maintain the formal high order accuracy. We study the L 2 stability and an L 2 a priori error estimate for smooth solutions of the reconstructed central discontinuous Galerkin method for linear hyperbolic equation. Numerical tests are presented to illustrate the accuracy and computational efficiency of the proposed method.

Abstract: Publication date: Available online 26 August 2016
Source:Applied Numerical Mathematics
Author(s): Alexandru Mihai Bica, Mircea Curila, Sorin Curila
A new iterative numerical method to solve two-point boundary value problems associated to functional differential equations of even order is proposed. The method uses a cubic spline interpolation procedure activated at each iterative step. The convergence of the method is proved and it is tested on some numerical experiments. The notion of numerical stability with respect to the choice of the first iteration is introduced proving that the proposed method is numerically stable in this sense.

Abstract: Publication date: Available online 27 August 2016
Source:Applied Numerical Mathematics
Author(s): Paulo B. Vasconcelos
Data-sparse representation techniques are emerging on computing approximate solutions for large scale problems involving matrices with low numerical rank. This representation provides both low memory requirements and cheap computational costs. In this work we consider the numerical solution of a large dimensional problem resulting from a finite rank discretization of an integral radiative transfer equation in stellar atmospheres. The integral operator, defined through the first exponential-integral function, is of convolution type and weakly singular. Hierarchically semiseparable representation of the matrix operator with low-rank blocks is built and data-sparse matrix computations can be performed with almost linear complexity. This representation of the original fully populated matrix is an algebraic multilevel structure built from a specific hierarchy of partitions of the matrix indices. Numerical tests illustrate the benefits of this matrix technique compared to standard storage schemes, dense and sparse, in terms of computational cost as well as memory requirements. This approach is particularly useful when a fine discretization of the integral equation is required and the resulting linear system of equations is of large dimension and numerically difficult to solve.

Abstract: Publication date: Available online 20 August 2016
Source:Applied Numerical Mathematics
Author(s): Sebastian Franz
It is well known that continuous Galerkin methods lack stability for singularly perturbed convection-diffusion problems. One approach to overcome this behaviour is to use discontinuous Galerkin methods instead. Unfortunately, this increases the number of degrees of freedom and thus the computational costs. We analyse discontinuous Galerkin methods of anisotropic polynomial order and discrete discontinuous spaces. By enforcing continuity in the vertices of a mesh, the number of unknowns can be reduced while the convergence order in the dG-norm is still sustained. Numerical experiments for several polynomial elements and finite element spaces support our theoretical results.

Abstract: Publication date: November 2016
Source:Applied Numerical Mathematics, Volume 109
Author(s): M.I.M. Copetti, M. Aouadi
The problem of thermoviscoelastic quasi-static contact between a rod and a rigid obstacle, when the diffusion effect is taken into account, is modeled and analyzed. The contact is modeled by the Signorini's condition and the stress–strain constitutive equation is of the Kelvin–Voigt type. In the quasi-static case, the governing equations correspond to the coupling of an elliptic and two parabolic equations. It poses some new mathematical difficulties due to the nonlinear boundary conditions. The existence of solutions is proved as the limit of solutions to a penalized problem. Moreover, we show that the weak solution converges to zero exponentially as time goes to infinity. Finally, we give some computational results where the influence of diffusion and viscosity are illustrated in contact.

Abstract: Publication date: November 2016
Source:Applied Numerical Mathematics, Volume 109
Author(s): Philipp Bader, David I. McLaren, G.R.W. Quispel, Marcus Webb
It is a classical theorem of Liouville that Hamiltonian systems preserve volume in phase space. Any symplectic Runge–Kutta method will respect this property for such systems, but it has been shown by Iserles, Quispel and Tse and independently by Chartier and Murua that no B-Series method can be volume preserving for all volume preserving vector fields. In this paper, we show that despite this result, symplectic Runge–Kutta methods can be volume preserving for a much larger class of vector fields than Hamiltonian systems, and discuss how some Runge–Kutta methods can preserve a modified measure exactly.

Abstract: Publication date: November 2016
Source:Applied Numerical Mathematics, Volume 109
Author(s): Wei-ping Shen, Chong Li, Xiao-qing Jin, Jen-chih Yao
We consider the convergence problem of some Newton-type methods for solving the inverse singular value problem with multiple and positive singular values. Under the nonsingularity assumption of the relative generalized Jacobian matrices at the solution c ⁎ , a convergence analysis for the multiple and positive case is provided and the superlinear or quadratical convergence properties are proved. Moreover, numerical experiments are given in the last section and comparisons are made.

Abstract: Publication date: November 2016
Source:Applied Numerical Mathematics, Volume 109
Author(s): J.B. Collins, P.A. Gremaud
In this paper we analyze the convergence of the domain decomposition method applied to transport problems on networks. In particular, we derive estimates for the number of required iterations for linear problems. These estimates can be used to determine when the implementation of domain decomposition methods would be beneficial for this type of problems.

Abstract: Publication date: November 2016
Source:Applied Numerical Mathematics, Volume 109
Author(s): Thomas Wolf, Heiko K.F. Panzer, Boris Lohmann
A new version of the alternating directions implicit (ADI) iteration for the solution of large-scale Lyapunov equations is introduced. It generalizes the hitherto existing iteration, by incorporating tangential directions in the way they are already available for rational Krylov subspaces. Additionally, first strategies to adaptively select shifts and tangential directions in each iteration are presented. Numerical examples emphasize the potential of the new results.

Abstract: Publication date: November 2016
Source:Applied Numerical Mathematics, Volume 109
Author(s): Jingjun Zhao, Rui Zhan, Alexander Ostermann
The aim of this paper is to analyze stability properties of explicit exponential integrators for three kinds of delay differential equations. First, linear autonomous delay differential equations are studied, and sufficient conditions for P- and GP-contractivity of explicit exponential Runge–Kutta methods are given. Second, for linear nonautonomous test equations, PN- and GPN-stability of a particular Magnus integrator is investigated. It is shown that the Magnus integrator is GPN-stable and convergent of order two. Finally, for semilinear delay differential equations, RN- and GRN-stability of explicit exponential Runge–Kutta methods is studied and sufficient conditions for GRN-stability are derived. Some examples of P-, GP-, RN- and GRN-stable exponential integrators are given, and numerical experiments that illustrate the theoretical results are included.

Abstract: Publication date: November 2016
Source:Applied Numerical Mathematics, Volume 109
Author(s): Matthias Schlottbom
We use a diffuse interface method for solving Poisson's equation with a Dirichlet condition on an embedded curved interface. The resulting diffuse interface problem is identified as a standard Dirichlet problem on approximating regular domains. We estimate the errors introduced by these domain perturbations, and prove convergence and convergence rates in the H 1 -norm, the L 2 -norm and the L ∞ -norm in terms of the width of the diffuse layer. For an efficient numerical solution we consider the finite element method for which another domain perturbation is introduced. These perturbed domains are polygonal and non-convex in general. We prove convergence and convergences rates in the H 1 -norm and the L 2 -norm in terms of the layer width and the mesh size. In particular, for the L 2 -norm estimates we present a problem adapted duality technique, which crucially makes use of the error estimates derived for the regularly perturbed domains. Our results are illustrated by numerical experiments, which also show that the derived estimates are sharp.

Abstract: Publication date: November 2016
Source:Applied Numerical Mathematics, Volume 109
Author(s): Y. Mesri, M. Khalloufi, E. Hachem
In this paper we derive a multi-dimensional mesh adaptation method which produces optimal meshes for quadratic functions, positive semi-definite. The method generates anisotropic adaptive meshes as quasi-uniform ones in some metric space, with the tensor metric being computed based on interpolation error estimates. It does not depend, a priori, on the PDEs at hand in contrast to residual methods. The estimated error is then used to steer local modifications of the mesh in order to reach a prescribed level of error in L p -norm or a prescribed number of elements. The L p -norm of the estimated error is then minimized in order to get an optimal mesh. Numerical examples in 2D and 3D for analytic challenging problems and an application to a Computational Fluid Dynamics problem are presented and discussed in order to show how the proposed method recovers optimal convergence rates as well as to demonstrate its computational performance.

Abstract: Publication date: November 2016
Source:Applied Numerical Mathematics, Volume 109
Author(s): María G. Armentano, Verónica Moreno
The goal of this work is to introduce a local and a global interpolator in Jacobi-weighted spaces, with optimal order of approximation in the context of the p-version of finite element methods. Then, an a posteriori error indicator of the residual type is proposed for a model problem in two dimensions and, in the mathematical framework of the Jacobi-weighted spaces, the equivalence between the estimator and the error is obtained on appropriate weighted norm.

Abstract: Publication date: November 2016
Source:Applied Numerical Mathematics, Volume 109
Author(s): Mehdi Dehghan, Mostafa Abbaszadeh
The element free Galerkin technique is a meshless method based on the variational weak form in which the test and trial functions are the shape functions of moving least squares approximation. Since the shape functions of moving least squares approximation do not have the Kronecker property thus the Dirichlet boundary condition can not be applied directly and also in this case obtaining an error estimate is not simple. The main aim of the current paper is to propose an error estimate for the extracted numerical scheme from the element free Galerkin method. To this end, we select the fractional cable equation with Dirichlet boundary condition. Firstly, we obtain a time-discrete scheme based on a finite difference formula with convergence order O ( τ 1 + min { α , β } ) , then we use the meshless element free Galerkin method, to discrete the space direction and obtain a full-discrete scheme. Also, for calculating the appeared integrals over the boundary and the domain of problem the Gauss–Legendre quadrature rule has been used. In the next, we change the main problem with Dirichlet boundary condition to a new problem with Robin boundary condition. Then, we show that the new technique is unconditionally stable and convergent using the energy method. We show convergence orders of the time discrete scheme and the full discrete scheme are O ( τ 1 + min { α , β } ) and O ( r p + 1 + τ 1 + min { α , β } ) , respectively. So, we can say that the main aim of this paper is as follows, (1) Transferring the main problem with Dirichlet boundary condition (old problem) to a problem with Robin boundary condition (new problem), (2) Showing that with special condition (when σ → + ∞ ) the solution of the new problem is convergent to the solution of the old problem, (3) Obtaining an error estimate for the new problem. Numerical examples confirm the efficiency and accuracy of the proposed scheme.

Abstract: Publication date: December 2016
Source:Applied Numerical Mathematics, Volume 110
Author(s): Franco Dassi, Luca Formaggia, Stefano Zonca
Standard 3D mesh generation algorithms may produce a low quality tetrahedral mesh, i.e., a mesh where the tetrahedra have very small dihedral angles. In this paper, we propose a series of operations to recover these badly-shaped tetrahedra. In particular, we will focus on the shape of these undesired mesh elements by proposing a novel method to distinguish and classify them. For each of these configurations, we apply a suitable sequence of operations to get a higher mesh quality. Finally, we employ a random algorithm to avoid locks and loops in the procedure. The reliability of the proposed mesh optimization algorithm is numerically proved with several examples.

Abstract: Publication date: November 2016
Source:Applied Numerical Mathematics, Volume 109
Author(s): Yali Gao, Liquan Mei
In this paper, implicit–explicit multistep Galerkin methods are studied for two-dimensional nonlinear Schrödinger equations and coupled nonlinear Schrödinger equations. The spatial discretization is based on Galerkin method using linear and quadratic basis functions on triangular and rectangular finite elements. And the implicit–explicit multistep method is used for temporal discretization. Linear and nonlinear numerical tests are presented to verify the validity and efficiency of the numerical methods. The numerical results record that the optimal order of the error in L 2 and L ∞ norm can be reached.

Abstract: Publication date: Available online 3 August 2016
Source:Applied Numerical Mathematics
Author(s): Arvet Pedas, Enn Tamme, Mikk Vikerpuur
We consider a class of boundary value problems for linear fractional weakly singular integro-differential equations which involve Caputo-type derivatives. Using an integral equation reformulation of the boundary value problem, we first study the regularity of the exact solution. Based on the obtained regularity properties and spline collocation techniques, the numerical solution of the boundary value problem by suitable non-polynomial approximations is discussed. Optimal global convergence estimates are derived and a super-convergence result for a special choice of grid and collocation parameters is given. A numerical illustration is also presented.

Abstract: Publication date: Available online 3 August 2016
Source:Applied Numerical Mathematics
Author(s): Changbum Chun, Beny Neta
In this paper we analyze Murakami's family of fifth order methods for the solution of nonlinear equations. We show how to find the best performer by using a measure of closeness of the extraneous fixed points to the imaginary axis. We demonstrate the performance of these members as compared to the two members originally suggested by Murakami. We found several members for which the extraneous fixed points are on the imaginary axis, only one of these has 6 such points (compared to 8 for the other members). We show that this member is the best performer.

Abstract: Publication date: Available online 3 August 2016
Source:Applied Numerical Mathematics
Author(s): Francesco Fambri, Michael Dumbser
In this paper two new families of arbitrary high order accurate spectral discontinuous Galerkin (DG) finite element methods are derived on staggered Cartesian grids for the solution of the incompressible Navier–Stokes (NS) equations in two and three space dimensions. The discrete solutions of pressure and velocity are expressed in the form of piecewise polynomials along different meshes. While the pressure is defined on the control volumes of the main grid, the velocity components are defined on edge-based dual control volumes, leading to a spatially staggered mesh. Thanks to the use of a nodal basis on a tensor-product domain, all discrete operators can be written efficiently as a combination of simple one-dimensional operators in a dimension-by-dimension fashion. In the first family, high order of accuracy is achieved only in space, while a simple semi-implicit time discretization is derived by introducing an implicitness factor θ ∈ [ 0.5 , 1 ] for the pressure gradient in the momentum equation. The real advantages of the staggering arise after substituting the discrete momentum equation into the weak form of the continuity equation. In fact, the resulting linear system for the pressure is symmetric and positive definite and either block penta-diagonal (in 2D) or block hepta-diagonal (in 3D). As a consequence, the pressure system can be solved very efficiently by means of a classical matrix-free conjugate gradient method. From our numerical experiments we find that the pressure system appears to be reasonably well-conditioned, since in all test cases shown in this paper the use of a preconditioner was not necessary. This is a rather unique feature among existing implicit DG schemes for the Navier–Stokes equations. In order to avoid a stability restriction due to the viscous terms, the latter are discretized implicitly using again a staggered mesh approach, where the viscous stress tensor is also defined on the dual mesh. The second family of staggered DG schemes proposed in this paper achieves high order of accuracy also in time by expressing the numerical solution in terms of piecewise space–time polynomials. In order to circumvent the low order of accuracy of the adopted fractional stepping, a simple iterative Picard procedure is introduced, which leads to a space–time pressure-correction algorithm. In this manner, the symmetry and positive definiteness of the pressure system are not compromised. The resulting algorithm is stable, computationally very efficient, and at the same time arbitrary high order accurate in both space and time. These features are typically not easy to obtain all at the same time for a numerical method applied to the incompressible Navier–Stokes equations. The new numerical method has been thoroughly validated for approximation polynomials of degree up to N = 11 , using a large set of non-trivial test problems in two and three space dimensions, for which either analytical, numerical or experimental reference solutions exist.

Abstract: Publication date: Available online 12 August 2016
Source:Applied Numerical Mathematics
Author(s): John D. Towers
This paper presents a finite volume scheme for a scalar one-dimensional fluid-particle interaction model. When devising a finite volume scheme for this model, one difficulty that arises is how to deal with the moving source term in the PDE while maintaining a fixed grid. The fixed grid requirement comes from the ultimate goal of accommodating two or more particles. The finite volume scheme that we propose addresses the moving source term in a novel way. We use a modified computational stencil, with the lower part of the stencil shifted during those time steps when the particle crosses a mesh point. We then employ an altered convective flux to compensate the stencil shifts. The resulting scheme uses a fixed grid, preserves total momentum, and enforces several stability properties in the single-particle case. The single-particle scheme is easily extended to multiple particles by a splitting method.

Abstract: Publication date: Available online 12 August 2016
Source:Applied Numerical Mathematics
Author(s): V. Nijimbere, L.J. Campbell
This paper examines the development of a time-dependent nonreflecting boundary condition (or radiation condition) for use in simulations of the propagation of internal gravity waves in a two-dimensional geophysical fluid flow configuration. First, a linear radiation condition, originally derived by Campbell and Maslowe, is implemented in some linear test cases. It involves the computation of a Laplace convolution integral which is nonlocal in time and thus requires values of the dependent variable at all previous time levels. An approximation for the integral is implemented here to reduce the expense of the computation and the results obtained are shown to be more accurate than those obtained using steady boundary conditions. For larger amplitude waves, nonlinear equations are required and the application of the linear radiation condition gives rise to instabilities. A new nonlinear time-dependent nonreflecting boundary condition is introduced which takes into account wave mean flow interactions in the vicinity of the outflow boundary by including a component corresponding to the vertical divergence of the horizontal momentum flux. This prevents the development of numerical instabilities and gives more accurate results in a nonlinear test problem than the results obtained using the linear radiation condition.

Abstract: Publication date: October 2016
Source:Applied Numerical Mathematics, Volume 108
Author(s): D.S. Kamenetskiy
The paper traces the relation between the Embedded Discontinuous Galerkin method and the Streamline Upwind Petrov–Galerkin method and its generalizations, such as the Discontinuous Residual-free Bubble method, for the Euler and Navier–Stokes equations. Popular choices of the stabilization terms for both discretizations are related on the analytical level. The conservation property of the Embedded Discontinuous Galerkin discretization on dual volumes of the computational grid is established.

Abstract: Publication date: Available online 3 June 2016
Source:Applied Numerical Mathematics
Author(s): Leonardo Di G. Sigalotti, Jaime Klapp, Otto Rendón, Carlos A. Vargas, Franklin Peña-Polo
The problem of consistency of smoothed particle hydrodynamics (SPH) has demanded considerable attention in the past few years due to the ever increasing number of applications of the method in many areas of science and engineering. A loss of consistency leads to an inevitable loss of approximation accuracy. In this paper, we revisit the issue of SPH kernel and particle consistency and demonstrate that SPH has a limiting second-order convergence rate. Numerical experiments with suitably chosen test functions validate this conclusion. In particular, we find that when using the root mean square error as a model evaluation statistics, well-known corrective SPH schemes, which were thought to converge to second, or even higher order, are actually first-order accurate, or at best close to second order. We also find that observing the joint limit when N → ∞ , h → 0 , and n → ∞ , as was recently proposed by Zhu et al., where N is the total number of particles, h is the smoothing length, and n is the number of neighbor particles, standard SPH restores full C 0 particle consistency for both the estimates of the function and its derivatives and becomes insensitive to particle disorder.

Abstract: Publication date: Available online 7 June 2016
Source:Applied Numerical Mathematics
Author(s): Po-Hsien Lin, S.-T. John Yu
A set of model equations are proposed to simulate waves generated by unsteady, low-speed, nearly incompressible air and water flows. The equations include the continuity and momentum equations with pressure and velocity as the unknowns. Compressibility effect associated with waves motion is directly tracked by time-accurate calculation of pressure fluctuations. The corresponding density changes are modeled by using the bulk modulus of the medium. The three-dimensional equations are shown to be hyperbolic by analyzing eigenvalues and eigenvectors of the composite Jacobian matrix of the equations. Specifically, the matrix is shown to be diagonalizable and have a real spectrum. Moreover, an analytical form of the Riemann invariants of the one-dimensional equations are derived. To validate the model equations, the space-time Conservation Element and Solution Element (CESE) method and the SOLVCON code are employed to solve the two-dimensional equations. Aeolian tones generated by air and water flows passing a cylinder and over an open cavity are simulated. Numerical results compare well with previously reported data.

Abstract: Publication date: Available online 9 June 2016
Source:Applied Numerical Mathematics
Author(s): Veronika Schleper
We present a new HLL-type approximate Riemann solver for a compressible two-phase flow model with phase transition and surface forces such as surface tension or electric forces. The solver is obtained following the main ideas of the HLL-approach. Due to the nonlinearity in the kinetic relation driving the phase transition, this solver involves the solution of one single nonlinear equation, contrary to a single-phase HLL-type solver, where no nonlinear equations have to be solved. We present some illustrative numerical examples to show the performance and accuracy of the new solver, comparing it to the relaxation Riemann solver from the literature as well as to the two-shock approximation, which is the most accurate interface Riemann solver currently available for the present situation, but suffers from high computational costs.

Abstract: Publication date: Available online 11 June 2016
Source:Applied Numerical Mathematics
Author(s): M. Kordy, E. Cherkaev, P. Wannamaker
This work develops a model order reduction method for a numerical solution of an inverse multi-frequency eddy current problem using a rational interpolation of the transfer function in the complex plane. We use the Pade interpolation in the complex frequency plane; this allows us to speed up the calculation of the frequency-dependent Jacobian in the inversion procedure without loosing accuracy. Interpolating frequencies are chosen adaptively to reduce the maximal approximation error. We use the error indicator that is equivalent to a seminorm of the residual. The efficiency of the developed approach is demonstrated by applying it to the inverse magnetotelluric problem, which is a geophysical electromagnetic remote sensing method used in mineral, geothermal, and groundwater exploration. In this application, the transfer function values are needed for shifts in a purely imaginary interval. Thus we consider the interpolating shifts in the same interval as well as in a purely real interval, containing the spectrum of the operator. Numerical tests show an excellent performance of the proposed methods characterized by a significant reduction of computational time without loss of accuracy of the calculated Jacobian.

Abstract: Publication date: September 2016
Source:Applied Numerical Mathematics, Volume 107
Author(s): S. Khodayari-Samghabadi, S.H. Momeni-Masuleh, A. Malek
In this paper, we present a stabilized explicit-extended penalty Galerkin method based on the implicit pressure and explicit saturation method to find the global solution for the two-phase flow in porous media at each time step. The bubble functions are employed as basis of the spatial dimensions for the extended penalty Galerkin method. The forward Euler method is applied to the temporal discretization. Since the accuracy of numerical simulations flow through porous media depends on the modeling of the injection and production well, we propose a new well model for the presented method. The details of the stability analysis for the proposed method are provided and suitable values of the penalty term and time steps are calculated. The efficiency of the method is illustrated by simulations of a waterflood in a heterogeneous oil reservoir. Comparisons are made with available literature which show the efficiency and accuracy of the proposed method.

Abstract: Publication date: Available online 11 June 2016
Source:Applied Numerical Mathematics
Author(s): C.F. Bracciali, A. Sri Ranga, A. Swaminathan
When a nontrivial measure μ on the unit circle satisfies the symmetry d μ ( e i ( 2 π − θ ) ) = − d μ ( e i θ ) then the associated orthogonal polynomials on the unit circle, say Φ n , are all real. In this case, in 1986, Delsarte and Genin have shown that the two sequences of para-orthogonal polynomials { z Φ n ( z ) + Φ n ⁎ ( z ) } and { z Φ n ( z ) − Φ n ⁎ ( z ) } , where Φ n ⁎ ( z ) = z n Φ n ( 1 / z ‾ ) ‾ , satisfy three term recurrence formulas and have also explored some further consequences of these sequences of polynomials such as their connections to sequences of orthogonal polynomials on the interval [ − 1 , 1 ] . The same authors, in 1988, have also provided a means to extend these results to cover any nontrivial measure on the unit circle. However, only recently the extension associated with the para-orthogonal polynomials z Φ n ( z ) − Φ n ⁎ ( z ) was thoroughly explored, especially from the point of view of three term recurrence and chain sequences. The main objective of the present article is to provide the theory surrounding the extension associated with the para-orthogonal polynomials z Φ n ( z ) + Φ n ⁎ ( z ) for any nontrivial measure on the unit circle. As an important application of the theory, a characterization for the existence of the integral ∫ 0 2 π e i θ − w − 2 d μ ( e i θ ) , where w is such that w = 1 , is given in terms of the coefficients α n − 1 = − Φ n ( 0 ) ‾ ... PubDate: 2016-06-16T18:07:53Z

Abstract: Publication date: October 2016
Source:Applied Numerical Mathematics, Volume 108
Author(s): Sanjib Kumar Acharya, Ajit Patel
In this article, a class of second order parabolic initial-boundary value problems in the framework of primal hybrid principle is discussed. The interelement continuity requirement for standard finite element method has been alleviated by using primal hybrid method. Finite elements are constructed and used in spatial direction, and backward Euler scheme is used in temporal direction for solving fully discrete scheme. Optimal order estimates for both the semidiscrete and fully discrete method are derived with the help of modified projection operator. Numerical results are obtained in order to verify the theoretical analysis.

Abstract: Publication date: September 2016
Source:Applied Numerical Mathematics, Volume 107
Author(s): Tomoyuki Miyaji, Paweł Pilarczyk, Marcio Gameiro, Hiroshi Kokubu, Konstantin Mischaikow
We study the usefulness of two most prominent publicly available rigorous ODE integrators: one provided by the CAPD group (capd.ii.uj.edu.pl), the other based on the COSY Infinity project (cosyinfinity.org). Both integrators are capable of handling entire sets of initial conditions and provide tight rigorous outer enclosures of the images under a time-T map. We conduct extensive benchmark computations using the well-known Lorenz system, and compare the computation time against the final accuracy achieved. We also discuss the effect of a few technical parameters, such as the order of the numerical integration method, the value of T, and the phase space resolution. We conclude that COSY may provide more precise results due to its ability of avoiding the variable dependency problem. However, the overall cost of computations conducted using CAPD is typically lower, especially when intervals of parameters are involved. Moreover, access to COSY is limited (registration required) and the rigorous ODE integrators are not publicly available, while CAPD is an open source free software project. Therefore, we recommend the latter integrator for this kind of computations. Nevertheless, proper choice of the various integration parameters turns out to be of even greater importance than the choice of the integrator itself.

Abstract: Publication date: September 2016
Source:Applied Numerical Mathematics, Volume 107
Author(s): Jochen Schütz, Klaus Kaiser
In this publication, we consider IMEX methods applied to singularly perturbed ordinary differential equations. We introduce a new splitting into stiff and non-stiff parts that has a direct extension to systems of conservation laws and investigate its performance analytically and numerically. We show that this splitting can in some cases improve the order of convergence, demonstrating that the phenomenon of order reduction is not only a consequence of the method but also of the splitting.

Abstract: Publication date: October 2016
Source:Applied Numerical Mathematics, Volume 108
Author(s): Haibing Wang, Jijun Liu
Consider the scattering of long ocean tidal waves by an island taking into account the influence of daily rotation of the Earth, which is modeled by an exterior boundary value problem for the two-dimensional Helmholtz equation with generalized oblique derivative boundary condition. In this paper, we are concerned with a corresponding inverse scattering problem which is to reconstruct the unknown obstacle (island) from the far-field data. After proving the unique solvability of the direct scattering problem in a suitable function space required for our inverse scattering problem, we establish the linear sampling method (LSM) for reconstructing the boundary of the obstacle from the far-field data. To clarify the validity of such a sampling-type method which essentially depends on the solvability of an interior boundary value problem, we show that, except a discrete set of wave numbers, such an interior problem has a unique solution. Finally, some numerical examples are presented to demonstrate the efficiency of the reconstruction scheme.

Abstract: Publication date: October 2016
Source:Applied Numerical Mathematics, Volume 108
Author(s): Yifen Ke, Changfeng Ma
In this note, a technical error is pointed out in the proof of Theorem 3 in the paper M. Benzi and X.-P. Guo [1]. A correct proof of this theorem is given.

Abstract: Publication date: October 2016
Source:Applied Numerical Mathematics, Volume 108
Author(s): Igor Boglaev
The paper deals with numerical solution of coupled systems of nonlinear parabolic equations based on a nonlinear ADI scheme. The convergence of the nonlinear ADI scheme to the continuous solution is proved. A monotone iterative ADI method is constructed. The existence and uniqueness of a solution of the nonlinear ADI scheme are established. An analysis of convergence of the monotone iterative ADI method to the solution of the nonlinear ADI scheme on the whole time interval is given. Numerical experiments are presented.

Abstract: Publication date: October 2016
Source:Applied Numerical Mathematics, Volume 108
Author(s): John T. Betts, Stephen L. Campbell, Karmethia C. Thompson
The numerical treatment of optimal control problems with state and control delays is important in a wide variety of scientific and technical applications. Solutions to these types of problems are difficult to obtain via analytic techniques since the system may be nonlinear and subjected to complicated inputs and constraints. There are several numerical methods available to compute the solutions of optimal control problems without delays. One such popular method is direct transcription. Although the numerical solutions of optimal control delay problems are important, less literature and software exists in this area. A general purpose industrial grade direct transcription code that can handle optimal control problems with both state and control constraints and delays is under development. Control delays pose a special challenge. A new technique for treating control delays when using a direct transcription approach is investigated in this paper.

Abstract: Publication date: October 2016
Source:Applied Numerical Mathematics, Volume 108
Author(s): Cong Shi, Chen Wang, Ting Wei
In this paper, we consider a class of severely ill-posed backward problems for linear parabolic equations. We use a convolution regularization method to obtain a stable approximate initial data from the noisy final data. The convergence rates are obtained under an a priori and an a posteriori regularization parameter choice rule in which the a posteriori parameter choice is a new generalized discrepancy principle based on a modified version of Morozov's discrepancy principle. The log-type convergence order under the a priori regularization parameter choice rule and log log -type order under the a posteriori regularization parameter choice rule are obtained. Two numerical examples are tested to support our theoretical results.

Abstract: Publication date: October 2016
Source:Applied Numerical Mathematics, Volume 108
Author(s): Raul Borsche
In this paper we propose a procedure to extend classical numerical schemes for hyperbolic conservation laws to networks of hyperbolic conservation laws. At the junctions of the network we solve the given coupling conditions and minimize the contributions of the outgoing numerical waves. This flexible procedure allows us to also use central schemes at the junctions. Several numerical examples are considered to investigate the performance of this new approach compared to the common Godunov solver and exact solutions.

Abstract: Publication date: October 2016
Source:Applied Numerical Mathematics, Volume 108
Author(s): Jing An
In this paper we present an efficient spectral method based on the Legendre–Galerkin approximation for the transmission eigenvalue problem. A rigorous error analysis is presented by using the minmax principle for the generalized eigenvalue problems associated to a transmission eigenvalue problem. However, this formulation can only compute real eigenvalues. Thus, we also present another formulation based on second order equations and construct an appropriate set of basis functions such that the matrices in the discrete variational form are sparse. For the case of constant medium, we derive the matrix formulations based on the tensor-product for the discrete variational form in two and three-dimensional cases, respectively. In addition, we also establish an optimization scheme based on the Legendre–Galerkin approximation. With this scheme we can estimate the index of refraction of an inhomogeneous medium. We also present ample numerical results to show that our method is very effective and high accurate.

Abstract: Publication date: October 2016
Source:Applied Numerical Mathematics, Volume 108
Author(s): Changfeng Ma, Na Huang
By reformulating a class of weakly nonlinear complementarity problems as implicit fixed-point equations based on splitting of the system matrix, a modified modulus-based matrix splitting algorithm is presented. The convergence analysis of proposed algorithm is established for the case that the splitting of the system matrix is an H-splitting. Numerical experiments on two model problems are given to illustrate the theoretical results and examine the numerical effectiveness.

Abstract: Publication date: October 2016
Source:Applied Numerical Mathematics, Volume 108
Author(s): Peter Benner, Matthias Heinkenschloss, Jens Saak, Heiko K. Weichelt
This paper improves the inexact Kleinman–Newton method for solving algebraic Riccati equations by incorporating a line search and by systematically integrating the low-rank structure resulting from ADI methods for the approximate solution of the Lyapunov equation that needs to be solved to compute the Kleinman–Newton step. A convergence result is presented that tailors the convergence proof for general inexact Newton methods to the structure of Riccati equations and avoids positive semi-definiteness assumptions on the Lyapunov equation residual, which in general do not hold for low-rank approaches. In the convergence proof of this paper, the line search is needed to ensure that the Riccati residuals decrease monotonically in norm. In the numerical experiments, the line search can lead to substantial reduction in the overall number of ADI iterations and, therefore, overall computational cost.

Abstract: Publication date: Available online 16 May 2016
Source:Applied Numerical Mathematics
Author(s): John W. Pearson
In this manuscript we consider the development of fast iterative solvers for Stokes control problems, an important class of PDE-constrained optimization problems. In particular we wish to develop effective preconditioners for the matrix systems arising from finite element discretizations of time-dependent variants of such problems. To do this we consider a suitable rearrangement of the matrix systems, and exploit the saddle point structure of many of the relevant sub-matrices involved – we may then use this to construct representations of these sub-matrices based on good approximations of their ( 1 , 1 ) -block and Schur complement. We test our recommended iterative methods on a distributed control problem with Dirichlet boundary conditions, and on a time-periodic problem.

Abstract: Publication date: Available online 11 May 2016
Source:Applied Numerical Mathematics
Author(s): Joachim Rang
It is well-known that one-step methods have order reduction if they are applied on stiff ODEs such as the example of Prothero–Robinson. In this paper we analyse the local error of Runge–Kutta and Rosenbrock–Wanner methods. We derive new order conditions and define with them B P R -consistency. We show that for strongly A-stable methods B P R -consistency implies B P R -convergence. Finally we analyse methods from literature, derive new B P R -consistent methods and present numerical examples. The numerical and analytical results show the influence of different properties of the methods and of different order conditions on the numerical error and on the numerical convergence order.

Abstract: Publication date: Available online 9 May 2016
Source:Applied Numerical Mathematics
Author(s): Antoine Tambue, Jean Medard T. Ngnotchouye
We consider a finite element approximation of a general semi-linear stochastic partial differential equation (SPDE) driven by space-time multiplicative and additive noise. We examine the full weak convergence rate of the exponential Euler scheme when the linear operator is self adjoint and also provide the full weak convergence rate for non-self-adjoint linear operator with additive noise. Key part of the proof does not rely on Malliavin calculus. For non-self-adjoint operators, we analyse the optimal strong error for spatially semi discrete approximations for both multiplicative and additive noise with truncated and non-truncated noise. Depending on the regularity of the noise and the initial solution, we found that in some cases the rate of weak convergence is twice the rate of the strong convergence. Our convergence rate is in agreement with some numerical results in two dimensions.

Abstract: Publication date: Available online 10 May 2016
Source:Applied Numerical Mathematics
Author(s): YuFeng Shi, Yan Guo
In this paper, we apply a maximum-principle-satisfying finite volume compact weighted scheme to numerical modelling traffic flow problems on networks. Road networks can be numerically model as a graph, whose edges are a finite number of roads that join at junctions. The evolution on each road is described by a scalar hyperbolic conservation law, and traffic distribution matrices are used to formulate coupling conditions at the network junctions. In order to achieve maximum-principle of the traffic density on each road, the maximum-principle-satisfying polynomial rescaling limiter is adopted. Numerical results for road networks with rich solution structures are presented in this work and indicate that the finite volume compact weighted scheme produces essentially non-oscillatory, maximum principle preserving and high resolution solutions.

Abstract: Publication date: Available online 10 May 2016
Source:Applied Numerical Mathematics
Author(s): Sabine Le Borne, Lusine Shahmuradyan
In several production processes, the distribution of particles dispersed in an environmental phase may be mathematically described by the solution of population balance equations. We are concerned with the development of efficient numerical techniques for the aggregation process: It invokes an integral term that is usually numerically expensive to evaluate and often dominates the total simulation cost. We describe an approach on locally refined nested grids to evaluate both the source and the sink terms in almost linear complexity (instead of quadratic complexity resulting from a direct approach). The key is to switch from a nodal to a wavelet basis representation of the density function. We illustrate the numerical performance of this approach, both in comparison to a discretization of piecewise constant functions on a uniform grid as well as to the fixed pivot method on a geometric grid.

Abstract: Publication date: Available online 2 May 2016
Source:Applied Numerical Mathematics
Author(s): G. Kreiss, B. Krank, G. Efraimsson
A zone of increasingly stretched grid is a robust and easy-to-use way to avoid unwanted reflections at artificial boundaries in wave propagating simulations. In such a buffer zone there are two main damping mechanisms, dissipation and under-resolution that turns a traveling wave into an evanescent wave. We present analysis in one and two space dimensions showing that evanescent decay through under-resolution is a very efficient way to damp waves. The analysis is supported by numerical computations.