This article is a brief survey of David Gottlieb's extraordinary research career. It is impossible to give a thorough presentation of all his research and the impact of his work, but we shall describe the main contributions and give examples of the results in some of his papers. David was for many years the dominating person in the development of spectral methods, and we devote much of the space in this article to this field.

Given a piecewise smooth function, it is possible to construct a global expansion in some complete orthogonal basis, such as the Fourier basis. However, the local discontinuities of the function will destroy the convergence of global approximations, even in regions for which the underlying function is analytic. The global expansions are contaminated by the presence of a local discontinuity, and the result is that the partial sums are oscillatory and feature non-uniform convergence. This characteristic behavior is called the Gibbs phenomenon. However, David Gottlieb and Chi-Wang Shu showed that these slowly and non-uniformly convergent global approximations retain within them high order information which can be recovered with suitable postprocessing. In this paper we review the history of the Gibbs phenomenon and the story of its resolution.

In many problems, one wishes to solve the Helmholtz equation with variable coefficients within the Laplacian-like term and use a high order accurate method (e.g., fourth order accurate) to alleviate the points-per-wavelength constraint by reducing the dispersion errors. The variation of coefficients in the equation may be due to an inhomogeneous medium and/or non-Cartesian coordinates. This renders existing fourth order finite difference methods inapplicable. We develop a new compact scheme that is provably fourth order accurate even for these problems. We present numerical results that corroborate the fourth order convergence rate for several model problems.

The stochastic collocation method using sparse grids has become a popular choice for performing stochastic computations in high dimensional (random) parameter space. In addition to providing highly accurate stochastic solutions, the sparse grid collocation results naturally contain sensitivity information with respect to the input random parameters. In this paper, we use the sparse grid interpolation and cubature methods of Smolyak together with combinatorial analysis to give a computationally efficient method for computing the global sensitivity values of Sobol'. This method allows for approximation of all main effect and total effect values from evaluation of f on a single set of sparse grids. We discuss convergence of this method, apply it to several test cases and compare to existing methods. As a result which may be of independent interest, we recover an explicit formula for evaluating a Lagrange basis interpolating polynomial associated with the Chebyshev extrema. This allows one to manipulate the sparse grid collocation results in a highly efficient manner.

The incompatibilities between the initial and boundary data will cause singularities at the time-space corners, which in turn adversely affect the accuracy of the numerical schemes used to compute the solutions. We study the corner singularity issue for nonlinear evolution equations in 1D, and propose two remedy procedures that effectively recover much of the accuracy of the numerical scheme in use. Applications of the remedy procedures to the 1D viscous Burgers equation, and to the 1D nonlinear reaction-diffusion equation are presented. The remedy procedures are applicable to other nonlinear diffusion equations as well.

In this article, we discuss a least-squares/fictitious domain method for the solution of linear elliptic boundary value problems with Robin boundary conditions. Let Ω and ω be two bounded domains of R^d such that ω⊂Ω. For a linear elliptic problem in Ω\ω with Robin boundary condition on the boundary γ of ω, our goal here is to develop a fictitious domain method where one solves a variant of the original problem on the full Ω, followed by a well-chosen correction over ω. This method is of the virtual control type and relies on a least-squares formulation making the problem solvable by a conjugate gradient algorithm operating in a well chosen control space. Numerical results obtained when applying our method to the solution of two-dimensional elliptic and parabolic problems are given; they suggest optimal order of convergence.

We investigate the effectiveness of two extrapolation-based methods aiming to approximate the initial state required by an iterative solver in simulations of unsteady flow problems. The methods lead to about a ten-fold reduction in the iteration count while requiring only negligible computational overhead. They are particularly suitable for parallel computing since they are based almost exclusively on data stored locally on each processor. Performance has been evaluated in simulations of turbulent flow in a stenosed carotid artery and also in laminar flow in a very large domain containing the human intracranial arterial tree.

A recent work of Li et al. [Numer. Math. Theor. Meth. Appl., 1(2008), pp. 92-112] proposed a finite volume solver to solve 2D steady Euler equations. Although the Venkatakrishnan limiter is used to prevent the non-physical oscillations nearby the shock region, the overshoot or undershoot phenomenon can still be observed. Moreover, the numerical accuracy is degraded by using Venkatakrishnan limiter. To fix the problems, in this paper the WENO type reconstruction is employed to gain both the accurate approximations in smooth region and non-oscillatory sharp profiles near the shock discontinuity. The numerical experiments will demonstrate the efficiency and robustness of the proposed numerical strategy.

We present a fast Poisson solver on spherical shells. With a special change of variable, the radial part of the Laplacian transforms to a constant coefficient differential operator. As a result, the Fast Fourier Transform can be applied to solve the Poisson equation with O(N³logN) operations. Numerical examples have confirmed the accuracy and robustness of the new scheme.

A novel Eulerian Gaussian beam method was developed in [8] to compute the Schrödinger equation efficiently in the semiclassical regime. In this paper, we introduce an efficient semi-Eulerian implementation of this method. The new algorithm inherits the essence of the Eulerian Gaussian beam method where the Hessian is computed through the derivatives of the complexified level set functions instead of solving the dynamic ray tracing equation. The difference lies in that, we solve the ray tracing equations to determine the centers of the beams and then compute quantities of interests only around these centers. This yields effectively a local level set implementation, and the beam summation can be carried out on the initial physical space instead of the phase plane. As a consequence, it reduces the computational cost and also avoids the delicate issue of beam summation around the caustics in the Eulerian Gaussian beam method. Moreover, the semi-Eulerian Gaussian beam method can be easily generalized to higher order Gaussian beam methods, which is the topic of the second part of this paper. Several numerical examples are provided to verify the accuracy and efficiency of both the first order and higher order semi-Eulerian methods.

We consider constraint preserving multidimensional evolution equations. A prototypical example is provided by the magnetic induction equation of plasma physics. The constraint of interest is the divergence of the magnetic field. We design finite volume schemes which approximate these equations in a stable manner and preserve a discrete version of the constraint. The schemes are based on reformulating standard edge centered finite volume fluxes in terms of vertex centered potentials. The potential-based approach provides a general framework for faithful discretizations of constraint transport and we apply it to both divergence preserving as well as curl preserving equations. We present benchmark numerical tests which confirm that our potential-based schemes achieve high resolution, while being constraint preserving.

A new adaptive cell average spectral element method (SEM) is proposed to solve the time-dependent Wigner equation for transport in quantum devices. The proposed cell average SEM allows adaptive non-uniform meshes in phase spaces to reduce the high-dimensional computational cost of Wigner functions while preserving exactly the mass conservation for the numerical solutions. The key feature of the proposed method is an analytical relation between the cell averages of the Wigner function in the k-space (local electron density for finite range velocity) and the point values of the distribution, resulting in fast transforms between the local electron density and local fluxes of the discretized Wigner equation via the fast sine and cosine transforms. Numerical results with the proposed method are provided to demonstrate its high accuracy, conservation, convergence and a reduction of the cost using adaptive meshes.

The rotational incremental pressure-correction (RIPC) scheme, described in Timmermans et al. [Int. J. Numer. Methods. Fluids., 22 (1996)] and Shen et al. [Math. Comput., 73 (2003)] for non-rotational Navier-Stokes equations, is extended to rotating incompressible flows. The method is implemented in the context of a pseudo Fourier-spectral code and applied to several rotating laminar and turbulent flows. The performance of the scheme and the computational results are compared to the so-called diagonalization method (DM) developed by Morinishi et al. [Int. J. Heat. Fluid. Flow., 22 (2001)]. The RIPC predictions are in excellent agreement with the DM predictions, while being simpler to implement and computationally more efficient. The RIPC scheme is not in anyway limited to implementation in a pseudo-spectral code or periodic boundary conditions, and can be used in complex geometries and with other suitable boundary conditions.

In this exploratory study, we present a new method of approximating a large system of ODEs by one with fewer equations, while attempting to preserve the essential dynamics of a reduced set of variables of interest. The method has the following key elements: (i) put a (simple, ad-hoc) probability distribution on the phase space of the ODE; (ii) assert that a small set of replacement variables are to be unknown linear combinations of the not-of-interest variables, and let the variables of the reduced system consist of the variables-of-interest together with the replacement variables; (iii) find the linear combinations that minimize the difference between the dynamics of the original system and the reduced system. We describe this approach in detail for linear systems of ODEs. Numerical techniques and issues for carrying out the required minimization are presented. Examples of systems of linear ODEs and variable-coefficient linear PDEs are used to demonstrate the method. We show that the resulting approximate reduced system of ODEs gives good approximations to the original system. Finally, some directions for further work are outlined.

We construct a well-conditioned hierarchical basis for triangular H(curl)- conforming elements with selected orthogonality. The basis functions are grouped into edge and interior functions, and the later is further grouped into normal and bubble functions. In our construction, the trace of the edge shape functions are orthonormal on the associated edge. The interior normal functions, which are perpendicular to an edge, and the bubble functions are both orthonormal among themselves over the reference element. The construction is made possible with classic orthogonal polynomials, viz., Legendre and Jacobi polynomials. For both the mass matrix and the quasi-stiffness matrix, better conditioning of the new basis is shown by a comparison with the basis previously proposed by Ainsworth and Coyle [Comput. Methods. Appl. Mech. Engrg., 190 (2001), 6709-6733].

In this paper we consider two commonly used classes of finite volume weighted essentially non-oscillatory (WENO) schemes in two dimensional Cartesian meshes. We compare them in terms of accuracy, performance for smooth and shocked solutions, and efficiency in CPU timing. For linear systems both schemes are high order accurate, however for nonlinear systems, analysis and numerical simulation results verify that one of them (Class A) is only second order accurate, while the other (Class B) is high order accurate. The WENO scheme in Class A is easier to implement and costs less than that in Class B. Numerical experiments indicate that the resolution for shocked problems is often comparable for schemes in both classes for the same building blocks and meshes, despite of the difference in their formal order of accuracy. The results in this paper may give some guidance in the application of high order finite volume schemes for simulating shocked flows.

A quasi non-overlapping hybrid scheme that combines the finite-difference time-domain (FDTD) method and the finite-element time-domain (FETD) method with nonconforming meshes is developed for time-domain solutions of Maxwell's equations. The FETD method uses mixed-order basis functions for electric and magnetic fields, while the FDTD method uses the traditional Yee's grid; the two methods are joined by a buffer zone with the FETD method and the discontinuous Galerkin method is used for the domain decomposition in the FETD subdomains. The main features of this technique is that it allows non-conforming meshes and an arbitrary numbers of FETD and FDTD subdomains. The hybrid method is completely stable for the time steps up to the stability limit for the FDTD method and FETD method. Numerical results demonstrate the validity of this technique.