Abstract: Publication date: August 2016
Source:Applied Numerical Mathematics, Volume 106
Author(s): Fikriye Yılmaz, Aytekin Çı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 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 17 May 2016
Source:Applied Numerical Mathematics
Author(s): A. Cardone, R. D'Ambrosio, B. Paternoster
The present paper illustrates the construction of direct quadrature methods of arbitrary order for Volterra integral equations with periodic solution. The coefficients of these methods depend on the parameters of the problem, following the exponential fitting theory. The convergence of these methods is analyzed, and some numerical experiments are illustrated to confirm theoretical expectations and for comparison with other existing methods.

Abstract: Publication date: August 2016
Source:Applied Numerical Mathematics, Volume 106
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 problem 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, 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: August 2016
Source:Applied Numerical Mathematics, Volume 106
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: August 2016
Source:Applied Numerical Mathematics, Volume 106
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 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.

Abstract: Publication date: Available online 4 May 2016
Source:Applied Numerical Mathematics
Author(s): Somayeh Gh. Bardeji, Isabel N. Figueiredo, Ercília Sousa
An optical flow variational model is proposed for a sequence of images defined on a domain in R 2 . We introduce a regularization term given by the L 1 norm of a fractional differential operator. To solve the minimization problem we apply the split Bregman method. Extensive experimental results, with performance evaluation, are presented to demonstrate the effectiveness of the new model and method and to show that our algorithm performs favorably in comparison to another existing method. We also discuss the influence of the order α of the fractional operator in the estimation of the optical flow, for 0 ≤ α ≤ 2 . We observe that the values of α for which the method performs better depend on the geometry and texture complexity of the image. Some extensions of our algorithm are also discussed.

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 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 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.