Mathematics - Numerical Analysis Publications (50)


Mathematics - Numerical Analysis Publications

We propose an alternative method to generate samples of a spatially correlated random field with applications to large-scale problems for forward propagation of uncertainty. A classical approach for generating these samples is the Karhunen-Lo\`{e}ve (KL) decomposition. However, the KL expansion requires solving a dense eigenvalue problem and is therefore computationally infeasible for large-scale problems. Read More

In this work, we propose a parallel-in-time solver for linear and nonlinear ordinary differential equations. The approach is based on an efficient multilevel solver of the Schur complement related to a multilevel time partition. For linear problems, the scheme leads to a fast direct method. Read More

We consider a quasilinear degenerate diffusion-reaction system that describes biofilm formation. The model exhibits two non-linear effects: a power law degeneracy as one of the dependent variables vanishes and a super diffusion singularity as it approaches unity. Biologically relevant solutions are characterized by a moving interface and gradient blow-up there. Read More

While reduced-order models (ROMs) have been popular for efficiently solving large systems of differential equations, the stability of reduced models over long-time integration is of present challenges. We present a greedy approach for ROM generation of parametric Hamiltonian systems that captures the symplectic structure of Hamiltonian systems to ensure stability of the reduced model. Through the greedy selection of basis vectors, two new vectors are added at each iteration to the linear vector space to increase the accuracy of the reduced basis. Read More

A new numerical method is proposed for a 1-D inverse medium scattering problem with multi-frequency data. This method is based on the construction of a weighted cost functional. The weight is a Carleman Weight Function (CWF). Read More

We discuss various universality aspects of numerical computations using standard algorithms. These aspects include empirical observations and rigorous results. We also make various speculations about computation in a broader sense. Read More

In this paper we present two strategies to enable "parallelization across the method" for spectral deferred corrections (SDC). Using standard low-order time-stepping methods in an iterative fashion, SDC can be seen as preconditioned Picard iteration for the collocation problem. Typically, a serial Gauss-Seidel-like preconditioner is used, computing updates for each collocation node one by one. Read More

This work reformulates the complete electrode model of electrical impedance tomography in order to enable its more efficient numerical solution. The model traditionally assumes constant contact conductances on all electrodes, which leads to a discontinuous Robin boundary condition since the gaps between the electrodes can be described by vanishing conductance. As a consequence, the regularity of the electromagnetic potential is limited to less than two square-integrable weak derivatives, which negatively affects the convergence of, e. Read More

In this paper we demonstrate that the framework of nonlinear spectral decompositions based on total variation (TV) regularization is very well suited for image fusion as well as more general image manipulation tasks. The well-localized and edge-preserving spectral TV decomposition allows to select frequencies of a certain image to transfer particular features, such as wrinkles in a face, from one image to another. We illustrate the effectiveness of the proposed approach in several numerical experiments, including a comparison to the competing techniques of Poisson image editing, linear osmosis, wavelet fusion and Laplacian pyramid fusion. Read More

Local adaptivity and mesh refinement are key to the efficient simulation of wave phenomena in heterogeneous media or complex geometry. Locally refined meshes, however, dictate a small time-step everywhere with a crippling effect on any explicit time-marching method. In [18] a leap-frog (LF) based explicit local time-stepping (LTS) method was proposed, which overcomes the severe bottleneck due to a few small elements by taking small time-steps in the locally refined region and larger steps elsewhere. Read More

A new approach is introduced for deriving a mixed variational formulation for Kirchhoff plate bending problems with mixed boundary conditions involving clamped, simply supported, and free boundary parts. Based on a regular decomposition of an appropriate nonstandard Sobolev space for the bending moments, the fourth-order problem can be equivalently written as a system of three (consecutively to solve) second-order problems in standard Sobolev spaces. This leads to new discretization methods, which are flexible in the sense, that any existing and well-working discretization method and solution strategy for standard second-order problems can be used as a modular building block of the new method. Read More

Newton method is one of the most powerful methods for finding solution of nonlinear equations. In its classical form it is applied for systems of $n$ equations with $n$ variables. However it can be modified for underdetermined equations (with $mRead More

We consider the statistical inverse problem to recover $f$ from noisy measurements $Y = Tf + \sigma \xi$ where $\xi$ is Gaussian white noise and $T$ a compact operator between Hilbert spaces. Considering general reconstruction methods of the form $\hat f_\alpha = q_\alpha \left(T^*T\right)T^*Y$ with an ordered filter $q_\alpha$, we investigate the choice of the regularization parameter $\alpha$ by minimizing an unbiased estimate of the predictive risk $\mathbb E\left[\Vert Tf - T\hat f_\alpha\Vert^2\right]$. The corresponding parameter $\alpha_{\mathrm{pred}}$ and its usage are well-known in the literature, but oracle inequalities and optimality results in this general setting are unknown. Read More

Assuming that the absence of perturbations guarantees weak or strong convergence to a common fixed point, we study the behavior of perturbed products of an infinite family of nonexpansive operators. Our main result indicates that the convergence rate of unperturbed products is essentially preserved in the presence of perturbations. This, in particular, applies to the linear convergence rate of dynamic string averaging projection methods, which we establish here as well. Read More

We discuss the construction of robust preconditioners for finite element approximations of Biot's consolidation model in poroelasticity. More precisely, we study finite element methods based on generalizations of the Hellinger-Reissner principle of linear elasticity, where the stress tensor is one of the unknowns. The Biot model has a number of applications in science, medicine, and engineering. Read More

We are concerned with the efficient implementation of symplectic implicit Runge-Kutta (IRK) methods applied to systems of (non-necessarily Hamiltonian) ordinary differential equations by means of Newton-like iterations. We pay particular attention to symmetric symplectic IRK schemes (such as collocation methods with Gaussian nodes). For a $s$-stage IRK scheme used to integrate a $d$-dimensional system of ordinary differential equations, the application of simplified versions of Newton iterations requires solving at each step several linear systems (one per iteration) with the same $sd \times sd$ real coefficient matrix. Read More

We introduce a new finite element (FE) discretization framework applicable for covariant split equations. The introduction of additional differential forms (DF) that form pairs with the original ones permits the splitting of the equations into topological momentum and continuity equations and metric-dependent closure equations that apply the Hodge-star operator. Our discretization framework conserves this geometrical structure and provides for all DFs proper FE spaces such that the differential operators hold in strong form. Read More

We propose a method to construct Riesz MRA on local fields of positive characteristic and corresponding scaling step functions that generate it. Read More

In this paper, a two-grid method is proposed to linearize and symmetrize the steady-state Poisson-Nernst-Planck equations. The computational system is decoupled to linearize and symmetrize equations by using this method, which can improve the computational efficiency compared with the finite element method. Error estimates are derived for the proposed method. Read More

We collect examples of boundary-value problems of Dirichlet and Dirichlet-Neumann type which we found instructive when designing and analysing numerical methods for fully nonlinear elliptic partial differential equations. In particular, our model problem is the Monge-Amp\`ere equation, which is treated through its equivalent reformulation as a Hamilton-Jacobi-Bellman equation. Our examples illustrate how the different notions of boundary conditions appearing in the literature may admit different sets of viscosity sub- and supersolutions. Read More

We propose general computational procedures based on descriptor state-space realizations to compute coprime factorizations of rational matrices with minimum degree denominators. Enhanced recursive pole dislocation techniques are developed, which allow to successively place all poles of the factors into a given "good" domain of the complex plane. The resulting McMillan degree of the denominator factor is equal to the number of poles lying in the complementary "bad" region and therefore is minimal. Read More

We propose and analyze a heterogenous multiscale method for the efficient integration of constant-delay differential equations subject to fast periodic forcing. The stroboscopic averaging method (SAM) suggested here may provide approximations with \(\mathcal{O}(H^2+1/\Omega^2)\) errors with a computational effort that grows like \(H^{-1}\) (the inverse of the stepsize), uniformly in the forcing frequency \(\Omega\). Read More

The paper discusses the construction of high dimensional spatial discretizations for arbitrary multivariate trigonometric polynomials, where the frequency support of the trigonometric polynomial is known. We suggest a construction based on the union of several rank-1 lattices as sampling scheme. We call such schemes multiple rank-1 lattices. Read More

The Green's function of a transformer is essential for prediction of its vibration. As the Green's function cannot be measured directly and completely, the finite element analysis (FEA) is typically used for its estimation. However, because of the complexity of the transformer structure, the calculations involved in FEA are time consuming. Read More

For time-dependent partial differential equations, parallel-in-time integration using the "parallel full approximation scheme in space and time" (PFASST) is a promising way to accelerate existing space-parallel approaches beyond their scaling limits. Inspired by the classical Parareal method and multigrid ideas, PFASST allows to integrate multiple time-steps simultaneously using a space-time hierarchy of spectral deferred correction sweeps. While many use cases and benchmarks exist, a solid and reliable mathematical foundation is still missing. Read More

We develop a finite element method for the Laplace--Beltrami operator on a surface described by a set of patchwise parametrizations. The patches provide a partition of the surface and each patch is the image by a diffeomorphism of a subdomain of the unit square which is bounded by a number of smooth trim curves. A patchwise tensor product mesh is constructed by using a structured mesh in the reference domain. Read More

Given a tetrahedral mesh and objective functionals measuring the mesh quality which take into account the shape, size, and orientation of the mesh elements, our aim is to improve the mesh quality as much as possible. In this paper, we combine the moving mesh smoothing, based on the integration of an ordinary differential equation coming from a given functional, with the lazy flip technique, a reversible edge removal algorithm to modify the mesh connectivity. Moreover, we utilize radial basis function (RBF) surface reconstruction to improve tetrahedral meshes with curved boundary surfaces. Read More

In recent years, the use of sparse recovery techniques in the approximation of high-dimensional functions has garnered increasing interest. In this work, we present a survey of recent progress in this emerging topic. Our main focus is on the computation of polynomial approximations of high-dimensional functions on d-dimensional hypercubes. Read More

In this work, we outline the entropy viscosity method and discuss how the choice of scaling influences the size of viscosity for a simple shock problem. We present examples to illustrate the performance of the entropy viscosity method under two distinct scalings. Read More

A new class of Hermite methods for solving nonlinear conservation laws is presented. While preserving the high order spatial accuracy for smooth solutions in the existing Hermite methods, the new methods come with better stability properties. Artificial viscosity in the form of the entropy viscosity method is added to capture shocks. Read More

Estimating the norm of the solution of the linear difference equation $u(\theta)-u(\theta+\omega)=v(\theta)$ plays a fundamental role in KAM theory. Optimal (in certain sense) estimates for the solution of this equation were provided by R\"ussmann in the mid 70's. The aim of this paper is to compare the sharpness of these classic estimates with more specific estimates obtained with the help of the computer. Read More

In this paper, we present an efficient energy stable scheme to solve a phase field model incorporating contact line condition. Instead of the usually used Cahn-Hilliard type phase equation, we adopt the Allen-Cahn type phase field model with the static contact line boundary condition that coupled with incompressible Navier-Stokes equations with Navier boundary condition. The projection method is used to deal with the Navier-Stokes equa- tions and an auxiliary function is introduced for the non-convex Ginzburg-Landau bulk potential. Read More

This paper proposes an adaptive timestep construction for an Euler-Maruyama approximation of the ergodic SDEs with a drift which is not globally Lipschitz over an infinite time interval. If the timestep is bounded appropriately, we show not only the stability of the numerical solution and the standard strong convergence order, but also that the bound for moments and strong error of the numerical solution are uniform in T, which allow us to introduce the adaptive multilevel Monte Carlo. Numerical experiments support our analysis. Read More

This paper is concerned with minimization of a fourth-order linearized Canham--Helfrich energy subject to Dirichlet boundary conditions on curves inside the domain. Such problems arise in the modeling of the mechanical interaction of biomembranes with embedded particles. There, the curve conditions result from the imposed particle--membrane coupling. Read More

For a symmetric positive semidefinite linear system of equations $\mathcal{Q} {\bf x} = {\bf b}$, where ${\bf x} = (x_1,\ldots,x_s)$ is partitioned into $s$ blocks, with $s \geq 2$, we show that each cycle of the classical block symmetric Gauss-Seidel (block sGS) method exactly solves the associated quadratic programming (QP) problem but added with an extra proximal term of the form $\frac{1}{2} \| {\bf x}-{\bf x}^k \|_{\mathcal T}^2$, where ${\mathcal T}$ is a symmetric positive semidefinite matrix related to the sGS decomposition and ${\bf x}^k$ is the previous iterate. By leveraging on such a connection to optimization, we are able to extend the result (which we name as the block sGS decomposition theorem) for solving a convex composite QP (CCQP) with an additional possibly nonsmooth term in $x_1$, i.e. Read More

In this paper we describe two fully mass conservative, energy stable, finite difference methods on a staggered grid for the quasi-incompressible Navier-Stokes-Cahn-Hilliard (q-NSCH) system governing a binary incompressible fluid flow with variable density and viscosity. Both methods, namely the primitive method (finite difference method in the primitive variable formulation) and the projection method (finite difference method in a projection-type formulation), are so designed that the mass of the binary fluid is preserved, and the energy of the system equations is always non-increasing in time at the fully discrete level. We also present an efficient, practical nonlinear multigrid method - comprised of a standard FAS method for the Cahn-Hilliard equation, and a method based on the Vanka-type smoothing strategy for the Navier-Stokes equation - for solving these equations. Read More

In this paper, we present a variational integrator that is based on an approximation of the Euler--Lagrange boundary-value problem via Taylor's method. This can viewed as a special case of the shooting-based variational integrator. The Taylor variational integrator exploits the structure of the Taylor method, which results in a shooting method that is one order higher compared to other shooting methods based on a one-step method of the same order. Read More

This paper investigates gradient recovery schemes for data defined on discretized manifolds. The proposed method, parametric polynomial preserving recovery (PPPR), does not ask for the tangent spaces of the exact manifolds which have been assumed for some significant gradient recovery methods in the literature. Another advantage of the proposed method is that it removes the symmetric requirement from the existing methods for the superconvergence. Read More

In this paper, we deal with the acceleration of the convergence of infinite series $\sum^\infty_{n=1}a_n$, when the terms $a_n$ are in general complex and have asymptotic expansions that can be expressed in the form $$ a_n\sim[\Gamma(n)]^{s/m}\exp\left[Q(n)\right]\sum^\infty_{i=0}w_i n^{\gamma-i/m}\quad\text{as $n\to\infty$},$$ where $\Gamma(z)$ is the gamma function, $m\geq1$ is an arbitrary integer, $Q(n)=\sum^{m}_{i=0}q_in^{i/m}$ is a polynomial of degree at most $m$ in $n^{1/m}$, $s$ is an arbitrary integer, and $\gamma$ is an arbitrary complex number. This can be achieved effectively by applying the $\tilde{d}^{(m)}$ transformation of the author to the sequence $\{A_n\}$ of the partial sums $A_n=\sum^n_{k=1}a_k$, $n=1,2,\dots\ .$ We give a detailed review of the properties of such series and of the $\tilde{d}^{(m)}$ transformation and the recursive W-algorithm that implements it. Read More

We study the effect of adaptive mesh refinement on a parallel domain decomposition solver of a linear system of algebraic equations. These concepts need to be combined within a parallel adaptive finite element software. A prototype implementation is presented for this purpose. Read More

We show that the concept of $H^2$-gradient flow for the Willmore energy and other functionals that depend at most quadratically on the second fundamental form is well-defined in the space of immersions of Sobolev class $W^{2,p}$ from a compact, $n$-dimensional manifold into Euclidean space, provided that $p \geq 2$ and $p>n$. We also discuss why this is not true for Sobolev class $H^2=W^{2,2}$. In the case of equality constraints, we provide sufficient conditions for the existence of the projected $H^2$-gradient flow and demonstrate its usability for optimization with several numerical examples. Read More

The conditioning of implicit Runge-Kutta (RK) integration for linear finite element approximation of diffusion equations on general anisotropic meshes is investigated. Bounds are established for the condition number of the resulting linear system with and without diagonal preconditioning for the implicit Euler (the simplest implicit RK method) and general implicit RK methods. It is shown that the conditioning of an implicit RK method behaves like that of the implicit Euler method. Read More

We develop a Monte-Carlo based numerical method for solving discrete-time stochastic optimal control problems with inventory. These are optimal control problems in which the control affects only a deterministically evolving inventory process on a compact state space while the random underlying process manifests itself through the objective functional. We propose a Regress Later modification of the traditional Regression Monte Carlo which allows to decouple inventory levels in two successive time steps and to include in the basis functions of the regression the dependence on the inventory levels. Read More

Kernel quadratures and other kernel-based approximation methods typically suffer from prohibitive cubic time and quadratic space complexity in the number of function evaluations. The problem arises because a system of linear equations needs to be solved. In this article we show that the weights of a kernel quadrature rule can be computed efficiently and exactly for up to tens of millions of nodes if the kernel, integration domain, and measure are fully symmetric and the node set is a union of fully symmetric sets. Read More

A new highly efficient method is developed for computation of traveling periodic waves (Stokes waves) on the free surface of deep water. A convergence of numerical approximation is determined by the complex singularites above the free surface for the analytical continuation of the travelling wave into the complex plane. An auxiliary conformal mapping is introduced which moves singularities away from the free surface thus dramatically speeding up numerical convergence by adapting the numerical grid for resolving singularities while being consistent with the fluid dynamics. Read More

The (fast) component-by-component (CBC) algorithm is an efficient tool for the construction of generating vectors for quasi-Monte Carlo rank-1 lattice rules in weighted reproducing kernel Hilbert spaces. We consider product weights, which assigns a weight to each dimension. These weights encode the effect a certain variable (or a group of variables by the product of the individual weights) has. Read More

We define a generalized finite element method for the discretization of elliptic partial differential equations in heterogeneous media. An adaptive local finite element basis (AL basis) on a coarse mesh which does not resolve the matrix of the media is constructed by solving finite-dimensional localized problems. The method requires $O(log(1/H)^{d+1})$ basis functions per mesh point. Read More

Unfitted finite element methods have a great potential for large scale simulations, since avoid the generation of body-fitted meshes and the use of graph partitioning techniques, two main bottlenecks for problems with non-trivial geometries. However, the linear systems that arise from these discretizations can be much more ill-conditioned, due to the so-called small cut cell problem. The state-of-the-art approach is to rely on sparse direct methods, which have quadratic complexity and thus, are not well-suited for large scale simulations. Read More

We consider Trace finite element methods for the linear membrane problem on second order tetrahedral elements. To accomplish this, zero-level set reconstruction methods for second order tetrahedra are considered. For the higher order membrane model a corresponding stabilization is proposed and numerically evaluated. Read More

The dense matrix resulting from an integral equation (IE) based solution of Maxwell's equations can be compactly represented by an ${\cal H}^2$-matrix. Given a general dense ${\cal H}^2$-matrix, prevailing fast direct solutions involve approximations whose accuracy can only be indirectly controlled. In this work, we propose new accuracy-controlled direct solution algorithms, including both factorization and inversion, for solving general ${\cal H}^2$-matrices, which does not exist prior to this work. Read More