Numerical Analysis Seminar

Abstract: In many practical applications, for instance, in computational electromagnetics, the excitation is time-harmonic. Switching from the time domain to the frequency domain allows us to avoid expensive time-stepping schemes by solving a simple elliptic equation for the amplitude. This is possible for linear problems, but not for nonlinear problems. However, due to the periodicity of the solution, we can expand the solution in a Fourier series. Truncating this Fourier series and approximating the Fourier coefficients by finite elements, we arrive at a large-scale coupled nonlinear system for determining the finite element approximation to the Fourier coefficients. The construction of fast solvers for such systems is very crucial for the efficiency of this multiharmonic approach. In this talk, we construct and analyze nearly optimal solvers for the Jacobi systems arising from the Newton linearization of the large-scale coupled nonlinear system. Numerical experiments demonstrate the performance of the solver.

We have calculated a six-dimensional ab initio potential energy surface for the OC-HF dimer at the CCSD(T)/aug-cc-pVTZ level of theory. A least-squares fitting method was used to fit the angular part of the potential, and the six-dimensional surface was obtained by interpolating the angular potential on a grid of R, rCO, and rHF points, using a three-dimensional reproducing kernel Hilbert spaces. In order to reduce the dimensionality of the system from six to four dimensions, the C-O and H-F stretching motions were adiabatically separated from the bending and stretching motions of the complex, using the vibrational self-consistent field method. The Hamiltonian operator has been split into six different terms, where the kinetic energy operators are diagonal in the spectral representation of the wave function. On the other hand, the potential energy is diagonal in the grid representation of the wave function. The spectral and grid representation are related through a unitary transformation. The total representation of the Hamiltonian operator in the spectral basis is obtained by transforming the vector in the grid representation to a spectral representation. Then the rovibrational energy levels were calculated using a Lanczos based iterative diagonalization scheme. This scheme requires repetitively acting the Hamiltonian operator on a vector, while avoiding the problem of constructing the full matrix.

Finite element method (FEM) is a powerful way to solve partial differential equations numerically. On the other hand, multigrid method gives fast solution to the matrix system resulting from FEM. Hybridized FEM has many advantages over tradition FEM, but its unique features make it difficult to discover corresponding multigrid algorithm. In this talk, I will give a general introduction to hybridized FEM (Raviart-Thmoas), and a newly-developed efficient multigrid method.

Iteratively least squares and gradient iterations intertwined with thresholding operations have been recently investigated for addressing inverse problems whose solutions are characterized by a few significant degrees of freedom.
We retrace some of the history of these algorithms and known results, and also address a variety of improved methods.

While the convergence of these algorithms is quite clarified, convergence rates and complexity are known only in special situations. In this talk we would like to focus on the complexity of compressive algorithms when addressing certain infinite dimensional problems. It is known that they may perform "arbitrarily bad" when applied for the regularized inversion of compact operators. Indeed for such operators the (infinite) matrix representation with respect to a "good basis", in the sense that it quasi-diagonalizes the operator, turns out to be diagonal dominant with fast decaying diagonal entries. The rate of convergence of the algorithms is related to the "local conditioning" of such a matrix, i.e., how well-conditioned is any relatively small group of columns. This is the case, for instance, when we deal with potential operators, such as in magnetic tomography, and matrix representations with respect to multiscale bases or wavelets. We discuss how to precondition these problems in order to obtain a uniform condition number of the resulting matrices over any small group of columns. In particular, we will show how block-diagonal preconditioning will produce infinite matrices with a "Restricted Isometry Property (RIP)", as the one introduced for finite dimensional situations in compressed sensing problems. We will use this property in order to show how adaptive numerical iterations can be performed guaranteeing a controlled linear convergence of these algorithms.

We discuss a conservative level set model for the simulation of two-phase flow based on a finite element discretization. The method consists of an incompressible Navier-Stokes solver to determine the fluid flow, coupled to an advection equation for the level set function that describes the position of the interface separating the two fluids. In our model, the level set function is set to be a smoothed indicator function, which provides a natural framework for discrete approximations of forces localized on the interface. Moreover, the (re-)initialization of the profile can be implemented such that the volume fractions of the individual fluids are conserved.

Numerical inaccuracies in our level set implementation can be traced back to errors in the approximation of normal vectors and curvature. The finite element spaces for approximating pressure and level set function need to match, in order not to introduce an additional imbalance between interface forces and pressure gradients. We analyze these inaccuracies on a drop in equilibrium, where numerical errors give rise to unphysical velocities.

A shortcoming of the level set model is its deficiency to model physics at contact lines, i.e., locations where the interface meets solids and aligns at material-dependent contact angles. We propose to tackle this problem by extending the level set model by a so-called phase field model close to boundaries, which does include a mechanism for setting oblique contact angles.