Abstract: Publication date: Available online 28 April 2016
Source:Applied Numerical Mathematics
Author(s): Giuseppe Izzo, Zdzislaw Jackiewicz
In this paper we use the theoretical framework of General Linear Methods (GLMs) to analyze and generalize the class of Cash's Modified Extended Backward Differentiation Formulae (MEBDF). Keeping the structure of MEBDF and their computational cost we propose a more general class of methods that can be viewed as a composition of modified linear multistep methods. These new methods are characterized by smaller error constants and possibly larger angles of A ( α ) -stability. Numerical experiments which confirm the good performance of these methods on a set of stiff problems are also reported.

Abstract: Publication date: Available online 29 April 2016
Source:Applied Numerical Mathematics
Author(s): Zhiqiang Li, Yubin Yan, Neville J. Ford
In this paper, we first introduce an alternative proof of the error estimates of the numerical methods for solving linear fractional differential equations proposed in Diethelm [6] where a first-degree compound quadrature formula was used to approximate the Hadamard finite-part integral and the convergence order of the proposed numerical method is O ( Δ t 2 − α ) , 0 < α < 1 , where α is the order of the fractional derivative and Δt is the step size. We then use a similar idea to prove the error estimates of the high order numerical method for solving linear fractional differential equations proposed in Yan et al. [37], where a second-degree compound quadrature formula was used to approximate the Hadamard finite-part integral and we show that the convergence order of the numerical method is O ( Δ t 3 − α ) , 0 < α < 1 . Nnumerical examples are given to show that the numerical results are consistent with the theoretical results.

Abstract: Publication date: July 2016
Source:Applied Numerical Mathematics, Volume 105
Author(s): Hai Bi, Hao Li, Yidu Yang
This paper proposes and analyzes an a posteriori error estimator for the finite element multi-scale discretization approximation of the Steklov eigenvalue problem. Based on the a posteriori error estimates, an adaptive algorithm of shifted inverse iteration type is designed. Finally, numerical experiments comparing the performances of three kinds of different adaptive algorithms are provided, which illustrate the efficiency of the adaptive algorithm proposed here.

Abstract: Publication date: July 2016
Source:Applied Numerical Mathematics, Volume 105
Author(s): Julia Leibinger, Michael Dumbser, Uwe Iben, Isabell Wayand
Flexible tubes are widely used in modern industrial hydraulic systems as connections between different components like valves, pumps and actuators. For the design and the analysis of the temporal behavior of a hydraulic system, one therefore needs an accurate mathematical model that describes the fluid flow in a compliant duct. Hence, in this paper we want to model the fluid–structure-interaction (FSI) problem given by the axially symmetric flow of a compressible barotropic fluid that flows through flexible tubes made of vulcanized rubber. The material of the tube can be described by using a visco-elastic rheology, which takes into account the strain relaxation of the material. The resulting mathematical model consists in a one-dimensional system of nonlinear hyperbolic partial differential equations (PDE) with non-conservative products and algebraic source terms. To solve this system numerically, we apply the DOT method, which is a generalized path-conservative Osher-type Riemann solver for conservative and non-conservative hyperbolic PDE recently proposed in [23] and [22]. We provide numerical evidence that the proposed DOT Riemann solver is well-balanced for the governing PDE system under consideration. The method is compared to available quasi-exact solutions of the Riemann problem in the case of an elastic wall described by the Laplace law. It is also compared to available experimental data and exact solutions obtained in the frequency domain for a linear visco-elastic wall behavior. In all cases under investigation the proposed path-conservative finite volume scheme based on the DOT Riemann solver is able to produce very accurate results.

Abstract: Publication date: Available online 21 April 2016
Source:Applied Numerical Mathematics
Author(s): Zhiping Mao, Sheng Chen, Jie Shen
We consider numerical approximation of the Riesz Fractional Differential Equations (FDEs), and construct a new set of generalized Jacobi functions, J n − α , − α ( x ) , which are tailored to the Riesz fractional PDEs. We develop optimal approximation results in non-uniformly weighted Sobolev spaces, and construct spectral Petrov–Galerkin algorithms to solve the Riesz FDEs with two kinds of boundary conditions (BCs): (i) homogeneous Dirichlet boundary conditions, and (ii) Integral BCs. We provide rigorous error analysis for our spectral Petrov–Galerkin methods, which show that the errors decay exponentially fast as long as the data (right-hand side function) is smooth, despite that fact that the solution has singularities at the endpoints. We also present some numerical results to validate our error analysis.

Abstract: Publication date: Available online 22 April 2016
Source:Applied Numerical Mathematics
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 16 April 2016
Source:Applied Numerical Mathematics
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: Available online 16 April 2016
Source:Applied Numerical Mathematics
Author(s): Dajana Conte, Beatrice Paternoster
The purpose of this paper is to employ graphics processing units for the numerical solution of large systems of weakly singular Volterra Integral Equations (VIEs), by means of Waveform Relaxation (WR) methods. A CUDA solver based on different kinds of WR iterations is developed. Numerical results on large systems of VIEs arising from the semi-discretization in space of fractional diffusion-wave equations are presented, showing the obtained speed-up.

Abstract: Publication date: Available online 13 April 2016
Source:Applied Numerical Mathematics
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: Available online 8 April 2016
Source:Applied Numerical Mathematics
Author(s): Peiqi Huang, Mingchao Cai, Feng Wang
In this paper, we propose a two-grid finite element method for solving the mixed Navier-Stokes/Darcy model with the Beavers-Joseph-Saffman interface condition. After solving a coupled nonlinear problems on a coarse grid, we sequentially solve decoupled and linearized subproblems on a fine grid and then correct the solution on the same grid. Compared with the existing work on the two-grid methods for the coupled model, it is shown that our two-grid method allows a much higher order scaling between the coarse grid size H and the fine grid size h. Specifically, if a k-th order discretization is applied, by using h = H 2 k + 1 k for k = 1 , 2 and h = H k + 1 k − 1 for k ≥ 3 , the final step two-grid solution errors in the energy norm are still optimal. Numerical experiments are also given to confirm the theoretical analysis.

Abstract: Publication date: Available online 7 April 2016
Source:Applied Numerical Mathematics
Author(s): K. Jbilou, A. Messaoudi
In the present paper we introduce new block extrapolation methods as generalizations of the well known vector extrapolation methods. We give expressions of the obtained approximations via the Schur complement and also propose an efficient implementation of these methods. Applications to linearly generated sequences are given and extensions to nonlinear problems are also given. Applications of the proposed block extrapolation methods to some nonlinear matrix equations are considered and some numerical examples are given.

Abstract: Publication date: Available online 4 April 2016
Source:Applied Numerical Mathematics
Author(s): F. Yılmaz, A. Çıbık
In this work, we apply a variational multiscale stabilization (VMS) to the optimal control problems of Navier-Stokes equations. We first obtain the optimality conditions by using Lagrange approach. After stating the continuous optimality system, we formulate the discrete problem by a projection-based stabilized finite element scheme. We give the stability estimates for both state and adjoint state variables. Then, we present the a priori error analysis. The main issue of this paper is to show the efficiency of the variational multiscale stabilization for this optimal control problem. We verify our results by numerical examples.

Abstract: Publication date: Available online 5 April 2016
Source:Applied Numerical Mathematics
Author(s): Mahboub Baccouch
We develop and analyze a new residual-based a posteriori error estimator for the discontinuous Galerkin (DG) method for nonlinear ordinary differential equations (ODEs). The a posteriori DG error estimator under investigation is computationally simple, efficient, and asymptotically exact. It is obtained by solving a local residual problem with no boundary condition on each element. We first prove that the DG solution exhibits an optimal O ( h p + 1 ) convergence rate in the L 2 -norm when p-degree piecewise polynomials with p ≥ 1 are used. We further prove that the DG solution is O ( h 2 p + 1 ) superconvergent at the downwind points. We use these results to prove that the p-degree DG solution is O ( h p + 2 ) super close to a particular projection of the exact solution. This superconvergence result allows us to show that the true error can be divided into a significant part and a less significant part. The significant part of the discretization error for the DG solution is proportional to the ( p + 1 ) -degree right Radau polynomial and the less significant part converges at O ( h p + 2 ) rate in the L 2 -norm. Numerical experiments demonstrate that the theoretical rates are optimal. Based on the global superconvergent approximations, we construct asymptotically exact a posteriori error estimates and prove that they converge to the true errors in the L 2 -norm under mesh refinement. The order of convergence is proved to be p + 2 . Finally, we prove that the global effectivity index in the L 2 -norm converges to unity at O ( h ) rate. Several numerical examples are provided to illustrate the global superconvergence results and the convergence of the proposed estimator under mesh refinement. A local adaptive procedure that makes use of our local a posteriori error estimate is also presented.

Abstract: Publication date: Available online 1 April 2016
Source:Applied Numerical Mathematics
Author(s): Norikazu Saito, Yoshiki Sugitani, Guanyu Zhou
We consider the stationary Stokes equations under a unilateral boundary condition of Signorini's type, which is one of artificial boundary conditions in flow problems. Well-posedness is discussed through its variational inequality formulation. We also consider the finite element approximation for a regularized penalty problem. The well-posedness, stability and error estimates of optimal order are established. The lack of a coupled Babuška and Brezzi's condition makes analysis difficult. We offer a new method of analysis. Particularly, our device to treat the pressure is novel and of some interest. Numerical examples are presented to validate our theoretical results.

Abstract: Publication date: Available online 30 March 2016
Source:Applied Numerical Mathematics
Author(s): Zhendong Gu, Xiaojing Guo, Daochun Sun
We propose series expansion method to solve VIEs (Volterra integral equations) with smooth given functions, including weakly singular VIEs possessing unsmooth solution. The key step in proposed method is to approximate given functions by their own Chebyshev Gauss-Lobatto interpolation polynomials. Lubich's results play an important role in the proposed method for weakly singular VIEs. Convergence analysis is provided for proposed method. Numerical experiments are carried out to confirm theoretical results.

Abstract: Publication date: Available online 30 March 2016
Source:Applied Numerical Mathematics
Author(s): I. Arun, Murugesan Venkatapathi
Sommerfeld integrals relate a spherical wave from a point source to a convolution set of plane and cylindrical waves. This relation does not have analytical solutions but it submits to a solution by numerical integration. Among others, it is significant for theoretical studies of many optical and radiation phenomena involving surfaces. This approach is preferred over discretized computational models of the surface because of the many orders of increased computations involved in the latter. One of the most widely used and accurate methods to compute these solutions is the numerical integration of the Sommerfeld integrand over a complex contour. We have analyzed the numerical advantages offered by this method, and have justified the optimality of the preferred contour of integration and the choice of two eigenfunctions used. In addition to this, we have also analyzed four other approximate methods to compute the Sommerfeld integral and have identified their regions of validity, and numerical advantages, if any. These include the high relative permittivity approximation, the short distance approximation, the exact image theory and Fourier expansion of the reflection coefficient. We also finally compare these five methods in terms of their computational cost.

Abstract: Publication date: June 2016
Source:Applied Numerical Mathematics, Volume 104
Author(s): Georgios D. Akrivis, Vassilios A. Dougalis, Efstratios Gallopoulos, Apostolos Hadjidimos, Ilias S. Kotsireas, Dimitrios Noutsos, Yiannis G. Saridakis, Michael N. Vrahatis

Abstract: Publication date: June 2016
Source:Applied Numerical Mathematics, Volume 104
Author(s): Martin P. Arciga Alejandre, Francisco J. Ariza Hernandez, Jorge Sanchez Ortiz
In this work, we consider an initial boundary-value problem for a stochastic evolution equation with Riesz-fractional spatial derivative and white noise on the half-line, { u t ( x , t ) = D x α u ( x , t ) + N u ( x , t ) + B ˙ ( x , t ) , x > 0 , t ∈ [ 0 , T ] , u ( x , 0 ) = u 0 ( x ) , x > 0 , u x ( 0 , t ) = g 1 ( t ) , t ∈ [ 0 , T ] , where D x α is the Riesz-fractional derivative, α ∈ ( 2 , 3 ) , N is a Lipschitzian operator and B ˙ ( x , t ) is the white noise. To construct the integral representation of solutions we use the Fokas method and Picard scheme to prove existence and uniqueness. Moreover, Monte Carlo methods are implemented to approximate solutions.

Abstract: Publication date: Available online 25 March 2016
Source:Applied Numerical Mathematics
Author(s): G. Landi, E. Loli Piccolomini, I. Tomba
We present a discrepancy-like stopping criterium for iterative regularization methods for the solution of linear discrete ill-posed problems. The presented criterium terminates the iterations of the iterative method when the residual norm of the computed solution becomes less or equal to the residual norm of a regularized Truncated Singular Value Decomposition (TSVD) solution. We present two algorithms for the automatic computation of the TSVD residual norm using the Discrete Picard Condition. The first algorithm uses the SVD coefficients while the second one uses the Fourier coefficients. In this work, we mainly focus on the Conjugate Gradient Least Squares method, but the proposed criterium can be used for terminating the iterations of any iterative regularization method. Many numerical tests on some selected one dimensional and image deblurring problems are presented and the results are compared with those obtained by state-of-the-art parameter selection rules. The numerical results show the efficiency and robustness of the proposed criterium.

Abstract: Publication date: June 2016
Source:Applied Numerical Mathematics, Volume 104
Author(s): Robert Nürnberg, Andrea Sacconi
Microelectronic circuits usually contain small voids or cracks, and if those defects are large enough to sever the line, they cause an open circuit. A fully practical finite element method for the temporal analysis of the migration of voids in the presence of surface diffusion, electric loading and elastic stress is presented. We simulate a bulk–interface coupled system, with a moving interface governed by a fourth-order geometric evolution equation and a bulk where the electric potential and the displacement field are computed. The method presented here follows a fitted approach, since the interface grid is part of the boundary of the bulk grid. A detailed analysis, in terms of experimental order of convergence (when the exact solution to the free boundary problem is known) and coupling operations (e.g., smoothing/remeshing of the grids, intersection between elements of the two grids), is carried out. A comparison with a previously introduced unfitted approach (where the two grids are totally independent) is also performed, along with several numerical simulations in order to test the accuracy of the methods.

Abstract: Publication date: June 2016
Source:Applied Numerical Mathematics, Volume 104
Author(s): S. González-Pinto, D. Hernández-Abreu
A family of splitting methods for the time integration of evolutionary Advection Diffusion Reaction Partial Differential Equations (PDEs) semi-discretized in space by Finite Differences is obtained. The splitting is performed in the Jacobian matrix by using the Approximate Matrix Factorization (AMF) and by considering up to three inexact Newton Iterations applied to the two-stage Radau IIA method along with a very simple predictor. The overall process allows to reduce the storage and the algebraic costs involved in the numerical solution of the multidimensional linear systems to the level of 1D-dimensional linear systems with small bandwidths. Some specific AMF-Radau methods are constructed after studying the expression for the local error in semi-linear equations, and their linear stability properties are widely studied. The wedge of stability of the methods depends on the number of splittings used for the Jacobian matrix of the spatial semidiscretized ODEs, J h = ∑ j = 1 d J h , j , where h stands for the spatial grid resolution. A-stability is proven for the cases d = 1 , 2 , and A ( 0 ) -stability for any d ≥ 1 . Numerical experiments on a 3D semi-linear advection diffusion reaction test problem and a 2D-combustion model are presented. The experiments show that the methods compare well with standard classical methods in parabolic problems and can also be successfully used for advection dominated problems when some diffusion or stiff reactions are present. In the latter case the stability imposes restrictions on the number of splitting terms (d).

Abstract: Publication date: June 2016
Source:Applied Numerical Mathematics, Volume 104
Author(s): Helen Christodoulidi, Tassos Bountis, Lambros Drossos
We study numerically classical 1-dimensional Hamiltonian lattices involving inter-particle long range interactions that decay with distance like 1 / r α , for α ≥ 0 . We demonstrate that although such systems are generally characterized by strong chaos, they exhibit an unexpectedly organized behavior when the exponent α < 1 . This is shown by computing dynamical quantities such as the maximal Lyapunov exponent, which decreases as the number of degrees of freedom increases. We also discuss our numerical methods of symplectic integration implemented for the solution of the equations of motion together with their associated variational equations. The validity of our numerical simulations is estimated by showing that the total energy of the system is conserved within an accuracy of 4 digits (with integration step τ = 0.02 ), even for as many as N = 8000 particles and integration times as long as 106 units.

Abstract: Publication date: June 2016
Source:Applied Numerical Mathematics, Volume 104
Author(s): Wanyok Atisattapong, Pasin Maruphanton
In this work we improve the accuracy and the convergence of the 1 / t algorithm for multidimensional numerical integration. The proposed strategy is to introduce a new approximation method which obviates the bin width effect of the conventional 1 / t algorithm by using the average of y values, which varies as the number of Monte Carlo trials changes, instead of the fixed value of y. The non-convergence of the 1 / t algorithm and the convergence of the new method are proved by theoretical analysis. The potential of the method is illustrated by the evaluation of one-, two- and multi- dimensional integrals up to six dimensions. Our results show that the numerical estimates from our method converge to their exact values without either error saturation or the bin with effect, in contrast with the conventional 1 / t algorithm.

Abstract: Publication date: June 2016
Source:Applied Numerical Mathematics, Volume 104
Author(s): A.C.L. Ashton, K.M. Crooks
In this paper we examine a numerical implementation of Fokas' unified method for elliptic boundary value problems on convex polygons. Within this setting the unified method provides a reconstruction of the unknown boundary values from the known boundary data for an important class of elliptic boundary value problems. We demonstrate that this approach can be used to study the classical Dirichlet eigenvalue problem for the Laplacian on convex polygons. We show that this problem is equivalent to a nonlinear eigenvalue problem for a semi-Fredholm operator which has holomorphic dependence on the eigenvalue itself. We also study the classical Dirichlet problem for the Helmholtz operator. We provide a new Galerkin scheme for the underlying Dirichlet–Neumann map, and prove that for sufficiently regular Dirichlet data this scheme converges with spectral accuracy.

Abstract: Publication date: June 2016
Source:Applied Numerical Mathematics, Volume 104
Author(s): Chris G. Antonopoulos, Tassos Bountis, Lambros Drossos
We investigate dynamically and statistically diffusive motion in a chain of linearly coupled 2-dimensional symplectic McMillan maps and find evidence of subdiffusion in weakly and strongly chaotic regimes when all maps of the chain possess a saddle point at the origin and the central map is initially excited. In the case of weak coupling, there is either absence of diffusion or subdiffusion with q > 1 -Gaussian probability distributions, characterizing weak chaos. However, for large enough coupling and already moderate number of maps, the system exhibits strongly chaotic ( q ≈ 1 ) subdiffusive behavior, reminiscent of the subdiffusive energy spreading observed in a disordered Klein–Gordon Hamiltonian. Our results provide evidence that coupled symplectic maps can exhibit physical properties similar to those of disordered Hamiltonian systems, even though the local dynamics in the two cases is significantly different.

Abstract: Publication date: June 2016
Source:Applied Numerical Mathematics, Volume 104
Author(s): G.V. Kozyrakis, A.I. Delis, G. Alexandrakis, N.A. Kampanis
A bed-load sediment transport model is used to describe realistic cases of the morphodynamics in coastal areas. The hydrodynamic equations are based on the well-known, two-dimensional depth-averaged non-linear shallow water equations, with bathymetry forces and friction, which are subsequently coupled to the Exner equation to describe the morphological evolution. Different forms of the bed-load transport flux are considered in the Exner equation and certain relations between them are established. The numerical model is expressed in a fully-coupled form where a single system of equations is solved by a high-resolution two-dimensional finite volume scheme of the relaxation type. The relaxation is performed by classical models where neither approximate Riemann solvers nor characteristic decompositions are needed. The overall numerical scheme is validated in benchmark problems, and for a realistic application of a coastal area on the northern side of the island of Crete. The validity of these results is established by comparisons made with the well-known MIKE Software by DHI Group.

Abstract: Publication date: June 2016
Source:Applied Numerical Mathematics, Volume 104
Author(s): D. Mantzavinos, M.G. Papadomanolaki, Y.G. Saridakis, A.G. Sifalakis
Gliomas are among the most aggressive forms of brain tumors. Over the last years mathematical models have been well developed to study gliomas growth. We consider a simple and well established mathematical model focused on proliferation and diffusion. Due to the heterogeneity of the brain tissue (white and grey matter) the diffusion coefficient is considered to be discontinuous. Fokas transform approach for the solution of linear PDE problems, apart from the fact that it avoids solving intermediate ODE problems, yields novel integral representations of the solution in the complex plane that decay exponentially fast and converge uniformly at the boundaries. To take advantage of these properties for the solution of the model problem at hand, we have successfully implemented Fokas transform method in the multi-domain environment induced by the interface discontinuities of our problem's domain. The fact that the integral representation of the solution at any time–space point of our problem's domain is independent on any other points of the domain, except of course on initial data, coupled with a simple composite trapezoidal rule, implemented on appropriately chosen integration contours, yields a fast and efficient analytical–numerical technique capable of producing directly high-order approximations of the solution at any point of the domain requiring no prior knowledge of the solution at any other time instances or space information.

Abstract: Publication date: June 2016
Source:Applied Numerical Mathematics, Volume 104
Author(s): Panagiotis D. Michailidis, Konstantinos G. Margaritis
Numerical linear algebra is one of the most important forms of scientific computation. The basic computations in numerical linear algebra are matrix computations and linear systems solution. These computations are used as kernels in many computational problems. This study demonstrates the parallelisation of these scientific computations using multi-core programming frameworks. Specifically, the frameworks examined here are Pthreads, OpenMP, Intel Cilk Plus, Intel TBB, SWARM, and FastFlow. A unified and exploratory performance evaluation and a qualitative study of these frameworks are also presented for parallel scientific computations with several parameters. The OpenMP and SWARM models produce good results running in parallel with compiler optimisation when implementing matrix operations at large and medium scales, whereas the remaining models do not perform as well for some matrix operations. The qualitative results show that the OpenMP, Cilk Plus, TBB, and SWARM frameworks require minimal programming effort, whereas the other models require advanced programming skills and experience. Finally, based on an extended study, general conclusions regarding the programming models and matrix operations for some parameters were obtained.

Abstract: Publication date: June 2016
Source:Applied Numerical Mathematics, Volume 104
Author(s): Marco Donatelli, Matteo Molteni, Vincenzo Pennati, Stefano Serra-Capizzano
In this paper we propose a scheme based on cubic splines for the solution of the second order two point boundary value problems. The solution of the algebraic system is computed by using optimized multigrid methods. In particular the transformation of the stiffness matrix essentially in a block Toeplitz matrix and its spectral analysis allow to choose smoothers able to reduce error components related to the various frequencies and to obtain an optimal method. The main advantages of our strategy can be listed as follows: (i) a fourth order of accuracy combined with a quadratic conditioning matrix, (ii) a resulting matrix structure whose eigenvalues can be compactly described by a symbol (this information is the key for designing an optimal multigrid method). Finally, some numerics that confirm the predicted behavior of the method are presented and a discussion on the two dimensional case is given, together with few 2D numerical experiments.

Abstract: Publication date: June 2016
Source:Applied Numerical Mathematics, Volume 104
Author(s): Andrea Cangiani, Emmanuil H. Georgoulis, Max Jensen
A discontinuous Galerkin (dG) method for the numerical solution of initial/boundary value multi-compartment partial differential equation (PDE) models, interconnected with interface conditions, is analysed. The study of interface problems is motivated by models of mass transfer of solutes through semi-permeable membranes. The case of fast reactions is also included. More specifically, a model problem consisting of a system of semilinear parabolic advection–diffusion–reaction partial differential equations in each compartment with only local Lipschitz conditions on the nonlinear reaction terms, equipped with respective initial and boundary conditions, is considered. General nonlinear interface conditions modelling selective permeability, congestion and partial reflection are applied to the compartment interfaces. The interior penalty dG method for this problem, presented recently, is analysed both in the space-discrete and in fully discrete settings for the case of, possibly, fast reactions. The a priori analysis shows that the method yields optimal a priori bounds, provided the exact solution is sufficiently smooth. Numerical experiments indicate agreement with the theoretical bounds.

Abstract: Publication date: Available online 21 March 2016
Source:Applied Numerical Mathematics
Author(s): M. Hubenthal, D. Onofrei
In previous works we considered the Helmholtz equation with fixed frequency k outside a discrete set of resonant frequencies, where it is implied that, given a source region D a ⊂ R d ( d = 2 , 3 ‾ ) and u 0 , a solution of the homogeneous scalar Helmholtz equation in a set containing the control region D c ⊂ R d , there exists an infinite class of boundary data on ∂ D a so that the radiating solution to the corresponding exterior scalar Helmholtz problem in R d ∖ D a will closely approximate u 0 in D c . Moreover, it will have vanishingly small values beyond a certain large enough “far-field” radius R. In this paper we study the minimal energy solution of the above problem (e.g. the solution obtained by using Tikhonov regularization with the Morozov discrepancy principle) and perform a detailed sensitivity analysis. In this regard we discuss the stability of the minimal energy solution with respect to measurement errors as well as the feasibility of the active scheme (power budget and accuracy) depending on: the mutual distances between the antenna, control region and far field radius R; value of the regularization parameter; frequency; location of the source.

Abstract: Publication date: Available online 22 March 2016
Source:Applied Numerical Mathematics
Author(s): Mansur I. Ismailov, Ibrahim Tekin
In this paper, the direct and inverse initial boundary value problems for a first order system of two hyperbolic equations are considered. The method of characteristics and the finite difference method are applied to the theoretical and numerical solutions of the direct problem, respectively. Moreover the suitability of the method of characteristics for the inverse problem of finding solely space-dependent coefficients and the finite difference method for solely time-dependent coefficients of the first order hyperbolic system are shown. The stability of the numerical method is supported by the examples.

Abstract: Publication date: Available online 21 March 2016
Source:Applied Numerical Mathematics
Author(s): Alexandra Koulouri, Ville Rimpiläinen, Mike Brookes, Jari P. Kaipio
In the inverse source problem of the Poisson equation, measurements on the domain boundaries are used to reconstruct sources inside the domain. The problem is an ill-posed inverse problem and it is sensitive to modelling errors of the domain. These errors can be boundary, structure and material property errors, for example. In this paper, we investigate whether the recently proposed Bayesian approximation error (BAE) approach could be used to alleviate the source estimation errors when an approximate model for the domain is employed. The BAE is based on postulating a probabilistic model for the uncertainties, in this case the geometry and structure of the domain, and to carry out approximate marginalization over these nuisance parameters. We particularly consider electroencephalography (EEG) source imaging as an application. EEG is a diagnostic brain imaging modality, and it can be used to reconstruct neural sources in the brain from electric potential measurements along the scalp. In the feasibility study, we assess to which degree one can recover from the modelling errors that are induced by the use of the three concentric circle head model instead of an anatomically accurate head model. The studied domain modelling errors include errors in the geometry of the exterior boundary and the structure of the interior. We show that, in particular with superficial dipole sources, the BAE yields estimates that can in some cases be considered adequately accurate. This would avoid the need for the extraction of the accurate head features which is conventionally carried out via expensive and time consuming auxiliary imaging modalities such as magnetic resonance imaging.

Abstract: Publication date: Available online 24 March 2016
Source:Applied Numerical Mathematics
Author(s): V.S. Sizikov, D.N. Sidorov
We propose the generalized quadrature methods for numerical solution of singular integral equation of Abel type. We overcome the singularity using the analytic computation of the singular integral. The problem of solution of singular integral equation is reduced to nonsingular system of linear algebraic equations without shift meshes techniques employment. We also propose generalized quadrature method for solution of Abel equation using the singular integral. Relaxed errors bounds are derived. In order to improve the accuracy we use Tikhonov regularization method. We demonstrate the efficiency of proposed techniques on infrared tomography problem. Numerical experiments show that it make sense to apply regularization in case of highly noisy (about 10%) sources only. That is due to the fact that singular integral equations enjoy selfregularization property.

Abstract: Publication date: Available online 22 March 2016
Source:Applied Numerical Mathematics
Author(s): M.A.V. Pinto, C. Rodrigo, F.J. Gaspar, C.W. Oosterlee
In this work, incomplete factorization techniques are used as smoothers within a geometric multigrid algorithm on triangular grids. A local Fourier analysis is proposed to study the smoothing properties of these methods, as well as the asymptotic convergence of the whole multigrid procedure. With this purpose, two- and three-grid local Fourier analysis are performed. Several two-dimensional diffusion problems, including different kinds of anisotropy are considered to demonstrate the robustness of this type of methods.

Abstract: Publication date: Available online 2 March 2016
Source:Applied Numerical Mathematics
Author(s): Haihua Qin, Xiaodong Liu
We are concerned with the reconstruction of both the penetrable inhomogeneous medium and the buried impenetrable obstacle. Firstly, the classical linear sampling method is used to recover the support of the inhomogeneous medium, and then a modification of the linear sampling method is proposed for objects buried in a known layered medium. The main feature of our method is that it avoids using knowledge of the Green's function for the background media. Finally, some numerical experiments are presented to demonstrate the feasibility and effectiveness of our method.

Abstract: Publication date: Available online 3 March 2016
Source:Applied Numerical Mathematics
Author(s): C. Rodrigo, F.J. Gaspar, F.J. Lisbona
A general local Fourier analysis for overlapping block smoothers on triangular grids is presented. This analysis is explained in a general form for its application to problems with different discretizations. This tool is demonstrated for two different problems: a stabilized linear finite element discretization of Stokes equations and an edge-based discretization of the curl-curl operator by lowest-order Nédélec finite element method. In this latter, special Fourier modes have to be considered in order to perform the analysis. Numerical results comparing two- and three-grid convergence factors predicted by the local Fourier analysis to real asymptotic convergence factors are presented to confirm the predictions of the analysis and show their usefulness.

Abstract: Publication date: Available online 15 February 2016
Source:Applied Numerical Mathematics
Author(s): Yiming Bu, Bruno Carpentieri, Zhaoli Shen, Tingzhu Huang
In this paper we introduce an algebraic recursive multilevel incomplete factorization preconditioner, based on a distributed Schur complement formulation, for solving general linear systems. The novelty of the proposed method is to combine factorization techniques of both implicit and explicit type, recursive combinatorial algorithms, multilevel mechanisms and overlapping strategies to maximize sparsity in the inverse factors and consequently reduce the factorization costs. Numerical experiments demonstrate the good potential of the proposed solver to precondition effectively general linear systems, also against other state-of-the-art iterative solvers of both implicit and explicit form.

Abstract: Publication date: July 2016
Source:Applied Numerical Mathematics, Volume 105
Author(s): Chang-tao Sheng, Zhong-qing Wang, Ben-yu Guo
In this paper, we propose a multistep Legendre–Gauss spectral collocation method for the nonlinear Volterra functional integro-differential equations (VFIDEs) with vanishing delays. This method is easy to implement and possesses the high-order accuracy. We also provide a rigorous convergence analysis of the hp-version of the multistep spectral collocation method under H 1 -norm. Numerical results confirm the theoretical analysis.

Abstract: Publication date: Available online 10 February 2016
Source:Applied Numerical Mathematics
Author(s): Pandelitsa Panaseti, Antri Zouvani, Niall Madden, Christos Xenophontos
We consider a fourth order singularly perturbed boundary value problem (BVP) in one-dimension and the approximation of its solution by the hp version of the Finite Element Method (FEM). The given problem's boundary conditions are not suitable for writing the BVP as a second order system, hence the approximation will be sought from a finite dimensional subspace of the Sobolev space H 2 . We construct suitable C 1 hierarchical basis functions for the approximation and we show that the hp FEM on the Spectral Boundary Layer Mesh yields a robust approximation that converges exponentially in the energy norm, as the number of degrees of freedom is increased. Numerical examples that validate (and extend) the theory are also presented.

Abstract: Publication date: Available online 11 February 2016
Source:Applied Numerical Mathematics
Author(s): A.R. Manapova, F.V. Lubyshev
In this work we consider optimization problems for processes described by semi-linear partial differential equations of elliptic type with discontinuous coefficients and solutions (with imperfect contact matching conditions), with controls involved in the coefficients. Finite difference approximations of optimization problems are constructed. For the numerical implementation of finite optimization problems differentiability and Lipshitz-continuity of the grid functional of the approximating grid problems are proved. An iterative method for solving boundary value problems of contact for PDEs of elliptic type with discontinuous coefficients and solutions is developed and validated. The convergence of the iterative process is investigated. And the convergence rate of iterations (with calculated constants) is estimated.

Abstract: Publication date: July 2016
Source:Applied Numerical Mathematics, Volume 105
Author(s): Franziska Nestler
We present an efficient method to compute the electrostatic fields, torques and forces in dipolar systems, which is based on the fast Fourier transform for nonequispaced data (NFFT). We consider 3d-periodic, 2d-periodic, 1d-periodic as well as 0d-periodic (open) boundary conditions. The method is based on the corresponding Ewald formulas, which immediately lead to an efficient algorithm only in the 3d-periodic case. In the other cases we apply the NFFT based fast summation in order to approximate the contributions of the nonperiodic dimensions in Fourier space. This is done by regularizing or periodizing the involved functions, which depend on the distances of the particles regarding the nonperiodic dimensions. The final algorithm enables a unified treatment of all types of periodic boundary conditions, for which only the precomputation step has to be adjusted.

Abstract: Publication date: Available online 27 January 2016
Source:Applied Numerical Mathematics
Author(s): Kasra Mohaghegh, Roland Pulch, Jan ter Maten
Nowadays electronic circuits comprise about a hundred million components on slightly more than one square centimeter. The model order reduction (MOR) techniques are among the most powerful tools to conquer this complexity and scale, although the nonlinear MOR is still an open field of research. On the one hand, the MOR techniques are well developed for linear ordinary differential equations (ODEs). On the other hand, we deal with differential algebraic equations (DAEs), which result from models based on network approaches. There are the direct and the indirect strategy to convert a DAE into an ODE. We apply the direct approach, where an artificial parameter is introduced in the linear system of DAEs. This results in a singular perturbed problem. On compact domains, uniform convergence of the transfer function of the regularized system towards the transfer function of the system of DAEs is proved in the general linear case. This convergence is for the transfer functions of the full model. We apply and investigate two different ways of MOR techniques in this context. We have two test examples, which are both TL models.

Abstract: Publication date: Available online 21 January 2016
Source:Applied Numerical Mathematics
Author(s): Farshid Dabaghi, Adrien Petrov, Jérôme Pousin, Yves Renard
This paper deals with a one-dimensional elastodynamic contact problem and aims to highlight some new numerical results. A new proof of existence and uniqueness results is proposed. More precisely, the problem is reformulated as a differential inclusion problem, the existence result follows from some a priori estimates obtained for the regularized problem while the uniqueness result comes from a monotonicity argument. An approximation of this evolutionary problem combining the finite element method as well as the mass redistribution method which consists on a redistribution of the body mass such that there is no inertia at the contact node, is introduced. Then two benchmark problems, one being new with convenient regularity properties, together with their analytical solutions are presented and some possible discretizations using different time-integration schemes are described. Finally, numerical experiments are reported and analyzed.

Abstract: Publication date: Available online 15 January 2016
Source:Applied Numerical Mathematics
Author(s): M.J. Ruijter, C.W. Oosterlee
We develop a Fourier method to solve quite general backward stochastic differential equations (BSDEs) with second-order accuracy. The underlying forward stochastic differential equation (FSDE) is approximated by different Taylor schemes, such as the Euler, Milstein, and Order 2.0 weak Taylor schemes, or by exact simulation. A θ-time-discretization of the time-integrands leads to an induction scheme with conditional expectations. The computation of the conditional expectations appearing relies on the availability of the characteristic function for these schemes. We will use the characteristic function of the discrete forward process. The expected values are approximated by Fourier cosine series expansions. Numerical experiments show rapid convergence of our efficient probabilistic numerical method. Second-order accuracy is observed and also proved. We apply the method to, among others, option pricing problems under the Constant Elasticity of Variance and Cox-Ingersoll-Ross processes.

Abstract: Publication date: Available online 8 January 2016
Source:Applied Numerical Mathematics
Author(s): C. Reisinger, P.A. Forsyth
An advantageous feature of piecewise constant policy timestepping for Hamilton-Jacobi–Bellman (HJB) equations is that different linear approximation schemes, and indeed different meshes, can be used for the resulting linear equations for different control parameters. Standard convergence analysis suggests that monotone (i.e., linear) interpolation must be used to transfer data between meshes. Using the equivalence to a switching system and an adaptation of the usual arguments based on consistency, stability and monotonicity, we show that if limited, potentially higher order interpolation is used for the mesh transfer, convergence is guaranteed. We provide numerical tests for the mean-variance optimal investment problem and the uncertain volatility option pricing model, and compare the results to published test cases.