A tutorial on pseudospectral methods for computational optimal control

[English] This paper is a tutorial introduction to pseudospectral optimal control. With pseudospectral methods, a function is approximated as a linear combination of smooth basis functions, which are often chosen to be Legendre or Chebyshev polynomials. Collocation of the differential-algebraic equations is performed at orthogonal collocation points, which are selected to yield interpolation of high accuracy. Pseudospectral methods directly discretize the original optimal control problem to recast it into a nonlinear programming format. A numerical optimizer is then employed to find approximate local optimal solutions. The paper also briefly describes the functionality and implementation of PSOPT, an open source software package written in C++ that employs pseudospectral discretization methods to solve multi-phase optimal control problems. The software implements the Legendre and Chebyshev pseudospectral methods, and it has useful features such as automatic differentiation, sparsity detection, and automatic scaling. The use of pseudospectral methods is illustrated in two problems taken from the literature on computational optimal control. [Portuguese] Este artigo e um tutorial introdutorio sobre controle otimo pseudo-espectral. Em metodos pseudo-espectrais, uma funcao e aproximada como uma combinacao linear de funcoes de base suaves, tipicamente escolhidas como polinomios de Legendre ou Chebyshev. A colocacao de equacoes algebrico-diferenciais e realizada em pontos de colocacao ortogonal, que sao selecionados de modo a minimizar o erro de interpolacao. Metodos pseudoespectrais discretizam o problema de controle otimo original de modo a converte-lo em um problema de programa cao nao-linear. Um otimizador numerico e entao empregado para obter solucoes localmente otimas. Este artigo tambem descreve sucintamente a funcionalidade e a implementacao de um pacote computacional de codigo aberto escrito em C++ chamado PSOPT. Tal pacote emprega metodos de discretizacao pseudo-spectrais para resolver problemas de controle otimo com multiplas fase. O PSOPT permite a utilizacao de metodos de Legendre ou Chebyshev, e possui caractersticas uteis tais como diferenciacao automatica, deteccao de esparsidade e escalonamento automatico. O uso de metodos pseudo-espectrais e ilustrado em dois problemas retirados da literatura de controle otimo computacional.

Numerical computation of an analytic singular value decomposition of a matrix valued function

This paper extends the singular value decomposition to a path of matricesE(t). An analytic singular value decomposition of a path of matricesE(t) is an analytic path of factorizationsE(t)=X(t)S(t)Y(t) T whereX(t) andY(t) are orthogonal andS(t) is diagonal. To maintain differentiability the diagonal entries ofS(t) are allowed to be either positive or negative and to appear in any order. This paper investigates existence and uniqueness of analytic SVD's and develops an algorithm for computing them. We show that a real analytic pathE(t) always admits a real analytic SVD, a full-rank, smooth pathE(t) with distinct singular values admits a smooth SVD. We derive a differential equation for the left factor, develop Euler-like and extrapolated Euler-like numerical methods for approximating an analytic SVD and prove that the Euler-like method converges.

Plane wave discontinuous Galerkin methods for the 2D Helmholtz equation: analysis of the $p$-version

Plane wave discontinuous Galerkin (PWDG) methods are a class of Trefftz-type methods for the spatial discretization of boundary value problems for the Helmholtz operator $-\Delta-\omega^2$, $\omega>0$. They include the so-called ultra weak variational formulation from [O. Cessenat and B. Després, SIAM J. Numer. Anal., 35 (1998), pp. 255–299]. This paper is concerned with the a priori convergence analysis of PWDG in the case of $p$-refinement, that is, the study of the asymptotic behavior of relevant error norms as the number of plane wave directions in the local trial spaces is increased. For convex domains in two space dimensions, we derive convergence rates, employing mesh skeleton-based norms, duality techniques from [P. Monk and D. Wang, Comput. Methods Appl. Mech. Engrg., 175 (1999), pp. 121–136], and plane wave approximation theory.

Correlations between two sets of angular relation equations

It is shown here that the angular relation equations between direct and reciprocal vectors are very similar to the angular relation equations in Euler's theorem. These two sets of equations are usually treated separately as unrelated equations in different fields. In this careful study, the connection between the two sets of angular equations is revealed by considering the cosine rule for the spherical triangle. It is found that understanding of the correlation is hindered by the facts that the same variables are defined differently and different symbols are used to represent them in the two fields. Understanding the connection between different concepts is not only stimulating and beneficial, but also a fundamental tool in innovation and research, and has historical significance. The background of the work presented here contains elements of many scientific disciplines. This work illustrates the common ground of two theories usually considered separately and is therefore of benefit not only for its own sake but also to illustrate a general principle that a theory relevant to one discipline can often be used in another. The paper works with chemistry related concepts using mathematical methodologies unfamiliar to the usual audience of mainstream experimental and theoretical chemists.

Heat kernel smoothing via Laplace-Beltrami eigenfunctions and its application to subcortical structure modeling

We present a new subcortical structure shape modeling framework using heat kernel smoothing constructed with the Laplace-Beltrami eigenfunctions. The cotan discretization is used to numerically obtain the eigenfunctions of the Laplace-Beltrami operator along the surface of subcortical structures of the brain. The eigenfunctions are then used to construct the heat kernel and used in smoothing out measurements noise along the surface. The proposed framework is applied in investigating the influence of age (38-79 years) and gender on amygdala and hippocampus shape. We detected a significant age effect on hippocampus in accordance with the previous studies. In addition, we also detected a significant gender effect on amygdala. Since we did not find any such differences in the traditional volumetric methods, our results demonstrate the benefit of the current framework over traditional volumetric methods.

Introductory remarks

The very first numerical models which were developed more than 20 years ago were drastic simplifications of the real atmosphere and they were mostly restricted to describe adiabatic processes. For prediction of a day or two of the mid tropospheric flow these models often gave reasonable results but the result deteriorated quickly when the prediction was extended further in time. The prediction of the surface flow was unsatisfactory even for short predictions. It was evident that both the energy generating processes as well as the dissipative processes have to be included in numerical models in order to predict the weather patterns in the lower part of the atmosphere and to predict the atmosphere in general beyond a day or two. Present-day computers make it possible to attack the weather forecasting problem in a more comprehensive and complete way and substantial efforts have been made during the last decade in particular to incorporate the non-adiabatic processes in numerical prediction models. The physics of radiational transfer, condensation of moisture, turbulent transfer of heat, momentum and moisture and the dissipation of kinetic energy are the most important processes associated with the formation of energy sources and sinks in the atmosphere and these have to be incorporated in numerical prediction models extended over more than a few days. The mechanisms of these processes are mainly related to small scale disturbances in space and time or even molecular processes. It is therefore one of the basic characteristics of numerical models that these small scale disturbances cannot be included in an explicit way. The reason for this is the discretization of the model's atmosphere by a finite difference grid or the use of a Galerkin or spectral function representation. The second reason why we cannot explicitly introduce these processes into a numerical model is due to the fact that some physical processes necessary to describe them (such as the local buoyance) are a priori eliminated by the constraints of hydrostatic adjustment. Even if this physical constraint can be relaxed by making the models non-hydrostatic the scale problem is virtually impossible to solve and for the foreseeable future we have to try to incorporate the ensemble or gross effect of these physical processes on the large scale synoptic flow. The formulation of the ensemble effect in terms of grid-scale variables (the parameters of the large-scale flow) is called 'parameterization'. For short range prediction of the synoptic flow at middle and high latitudes, very simple parameterization has proven to be rather successful.

Symmetric stability of compressible zonal flows on a generalized equatorial β plane

Sufficient conditions are derived for the linear stability with respect to zonally symmetric perturbations of a steady zonal solution to the nonhydrostatic compressible Euler equations on an equatorial � plane, including a leading order representation of the Coriolis force terms due to the poleward component of the planetary rotation vector. A version of the energy–Casimir method of stability proof is applied: an invariant functional of the Euler equations linearized about the equilibrium zonal flow is found, and positive definiteness of the functional is shown to imply linear stability of the equilibrium. It is shown that an equilibrium is stable if the potential vorticity has the same sign as latitude and the Rayleigh centrifugal stability condition that absolute angular momentum increase toward the equator on surfaces of constant pressure is satisfied. The result generalizes earlier results for hydrostatic and incompressible systems and for systems that do not account for the nontraditional Coriolis force terms. The stability of particular equilibrium zonal velocity, entropy, and density fields is assessed. A notable case in which the effect of the nontraditional Coriolis force is decisive is the instability of an angular momentum profile that decreases away from the equator but is flatter than quadratic in latitude, despite its satisfying both the centrifugal and convective stability conditions.

Free gravity waves and balanced dynamics

It is shown how a renormalization technique, which is a variant of classical Krylov–Bogolyubov–Mitropol’skii averaging, can be used to obtain slow evolution equations for the vortical and inertia–gravity wave components of the dynamics in a rotating flow. The evolution equations for each component are obtained to second order in the Rossby number, and the nature of the coupling between the two is analyzed carefully. It is also shown how classical balance models such as quasigeostrophic dynamics and its second-order extension appear naturally as a special case of this renormalized system, thereby providing a rigorous basis for the slaving approach where only the fast variables are expanded. It is well known that these balance models correspond to a hypothetical slow manifold of the parent system; the method herein allows the determination of the dynamics in the neighborhood of such solutions. As a concrete illustration, a simple weak-wave model is used, although the method readily applies to more complex rotating fluid models such as the shallow-water, Boussinesq, primitive, and 3D Euler equations.

Coupled Kelvin-wave and mirage-wave instabilities in semigeostrophic dynamics

A weak instability mode, associated with phase-locked counterpropagating coastal Kelvin waves in horizontal anticyclonic shear, is found in the semigeostrophic (SG) equations for stratified flow in a channel. This SG instability mode approximates a similar mode found in the Euler equations in the limit in which particle-trajectory slopes are much smaller than f/N, where f is the Coriolis frequency and N > f the buoyancy frequency. Though weak under normal parameter conditions, this instability mode is of theoretical interest because its existence accounts for the failure of an Arnol’d-type stability theorem for the SG equations. In the opposite limit, in which the particle motion is purely vertical, the Euler equations allow only buoyancy oscillations with no horizontal coupling. The SG equations, on the other hand, allow a physically spurious coastal “mirage wave,” so called because its velocity field vanishes despite a nonvanishing disturbance pressure field. Counterpropagating pairs of these waves can phase-lock to form a spurious “mirage-wave instability.” Closer examination shows that the mirage wave arises from failure of the SG approximations to be self-consistent for trajectory slopes f/N.

Nonlinear stability of multilayer quasi-geostrophic flow

New nonlinear stability theorems are derived for disturbances to steady basic flows in the context of the multilayer quasi-geostrophic equations. These theorems are analogues of Arnol’d's second stability theorem, the latter applying to the two-dimensional Euler equations. Explicit upper bounds are obtained on both the disturbance energy and disturbance potential enstrophy in terms of the initial disturbance fields. An important feature of the present analysis is that the disturbances are allowed to have non-zero circulation. While Arnol’d's stability method relies on the energy–Casimir invariant being sign-definite, the new criteria can be applied to cases where it is sign-indefinite because of the disturbance circulations. A version of Andrews’ theorem is established for this problem, and uniform potential vorticity flow is shown to be nonlinearly stable. The special case of two-layer flow is treated in detail, with particular attention paid to the Phillips model of baroclinic instability. It is found that the short-wave portion of the marginal stability curve found in linear theory is precisely captured by the new nonlinear stability criteria.

On Arnol'd's second nonlinear stability theorem for two-dimensional quasi-geostrophic flow

Arnol'd's second hydrodynamical stability theorem, proven originally for the two-dimensional Euler equations, can establish nonlinear stability of steady flows that are maxima of a suitably chosen energy-Casimir invariant. The usual derivations of this theorem require an assumption of zero disturbance circulation. In the present work an analogue of Arnol'd's second theorem is developed in the more general case of two-dimensional quasi-geostrophic flow, with the important feature that the disturbances are allowed to have non-zero circulation. New nonlinear stability criteria are derived, and explicit bounds are obtained on both the disturbance energy and potential enstrophy which are expressed in terms of the initial disturbance fields. While Arnol'd's stability method relies on the second variation of the energy-Casimir invariant being sign-definite, the new criteria can be applied to cases where the second variation is sign-indefinite because of the disturbance circulations. A version of Andrews' theorem is also established for this problem.

A general method for finding extremal states of Hamiltonian dynamical systems, with applications to perfect fluids

In addition to the Hamiltonian functional itself, non-canonical Hamiltonian dynamical systems generally possess integral invariants known as ‘Casimir functionals’. In the case of the Euler equations for a perfect fluid, the Casimir functionals correspond to the vortex topology, whose invariance derives from the particle-relabelling symmetry of the underlying Lagrangian equations of motion. In a recent paper, Vallis, Carnevale & Young (1989) have presented algorithms for finding steady states of the Euler equations that represent extrema of energy subject to given vortex topology, and are therefore stable. The purpose of this note is to point out a very general method for modifying any Hamiltonian dynamical system into an algorithm that is analogous to those of Vallis etal. in that it will systematically increase or decrease the energy of the system while preserving all of the Casimir invariants. By incorporating momentum into the extremization procedure, the algorithm is able to find steadily-translating as well as steady stable states. The method is applied to a variety of perfect-fluid systems, including Euler flow as well as compressible and incompressible stratified flow.

First order k-th moment finite element analysis of nonlinear operator equations with stochastic data

We develop and analyze a class of efficient Galerkin approximation methods for uncertainty quantification of nonlinear operator equations. The algorithms are based on sparse Galerkin discretizations of tensorized linearizations at nominal parameters. Specifically, we consider abstract, nonlinear, parametric operator equations J(\alpha ,u)=0 for random input \alpha (\omega ) with almost sure realizations in a neighborhood of a nominal input parameter \alpha _0. Under some structural assumptions on the parameter dependence, we prove existence and uniqueness of a random solution, u(\omega ) = S(\alpha (\omega )). We derive a multilinear, tensorized operator equation for the deterministic computation of k-th order statistical moments of the random solution's fluctuations u(\omega ) - S(\alpha _0). We introduce and analyse sparse tensor Galerkin discretization schemes for the efficient, deterministic computation of the k-th statistical moment equation. We prove a shift theorem for the k-point correlation equation in anisotropic smoothness scales and deduce that sparse tensor Galerkin discretizations of this equation converge in accuracy vs. complexity which equals, up to logarithmic terms, that of the Galerkin discretization of a single instance of the mean field problem. We illustrate the abstract theory for nonstationary diffusion problems in random domains.

A finite element method for nonlinear elliptic problems

We present a Galerkin method with piecewise polynomial continuous elements for fully nonlinear elliptic equations. A key tool is the discretization proposed in Lakkis and Pryer, 2011, allowing us to work directly on the strong form of a linear PDE. An added benefit to making use of this discretization method is that a recovered (finite element) Hessian is a byproduct of the solution process. We build on the linear method and ultimately construct two different methodologies for the solution of second order fully nonlinear PDEs. Benchmark numerical results illustrate the convergence properties of the scheme for some test problems as well as the Monge–Amp`ere equation and the Pucci equation.

The composite-tendency RAW filter in semi-implicit integrations

The time discretization in weather and climate models introduces truncation errors that limit the accuracy of the simulations. Recent work has yielded a method for reducing the amplitude errors in leapfrog integrations from first-order to fifth-order. This improvement is achieved by replacing the Robert--Asselin filter with the RAW filter and using a linear combination of the unfiltered and filtered states to compute the tendency term. The purpose of the present paper is to apply the composite-tendency RAW-filtered leapfrog scheme to semi-implicit integrations. A theoretical analysis shows that the stability and accuracy are unaffected by the introduction of the implicitly treated mode. The scheme is tested in semi-implicit numerical integrations in both a simple nonlinear stiff system and a medium-complexity atmospheric general circulation model, and yields substantial improvements in both cases. We conclude that the composite-tendency RAW-filtered leapfrog scheme is suitable for use in semi-implicit integrations.