Numerical methods for solving linearly damped Hamiltonian ordinary differential equations are analyzed and compared. The methods are constructed from the well-known Störmer-Verlet and implicit midpoint methods. The structure preservation properties of each method are shown analytically and numerically. Each method is shown to preserve a symplectic form up to a constant and are therefore conformal symplectic integrators, with each method shown to accurately preserve the rate of momentum dissipation. An analytical linear stability analysis is completed for each method, establishing thresholds between the value of the damping coefficient and the step-size that ensure stability. The methods are all second order and the preservation of the rate of energy dissipation is compared to that of a third order Runge-Kutta method that does not preserve conformal properties. Numerical experiments will include the damped harmonic oscillator and the damped nonlinear pendulum.

This work presents a novel method to model immiscible incompressible fluids in a stable manner. Here, the immiscible behavior of the flow is described by the incompressible Navier-Stokes-Cahn-Hilliard model, which is based on a diffuse interface method. We introduce buoyancy effects in the model through the Boussinesq approximation in a consistent manner. A structure-preserving discretization is used to guarantee the linear stability of the discrete problem and to satisfy the incompressibility of the discrete solution at every point in space by construction. For the solution of the model, we developed the Portable Extensible Toolkit for Isogeometric Analysis with Multi-Field discretizations (PetIGA-MF), a high-performance framework that supports structure-preserving spaces. PetIGA-MF is built on top of PetIGA and the Portable Extensible Toolkit for Scientific Computation (PETSc), sharing all their user-friendly, performance, and flexibility features. Herein, we describe the implementation of our model in PetIGA-MF and the details of the numerical solution. With several numerical tests, we verify the convergence, scalability, and validity of our approach. We use highly-resolved numerical simulations to analyze the merging and rising of droplets. From these simulations, we detailed the energy exchanges in the system to evaluate quantitatively the quality of our simulations. The good agreement of our results when compared against theoretical descriptions of the merging, and the small errors found in the energy analysis, allow us to validate our approach. Additionally, we present the development of an unconditionally energy-stable generalized-alpha method for the Swift-Hohenberg model that offers control over the numerical dissipation. A pattern formation example demonstrates the energy-stability and convergence of our method.

Generally, large-scale dynamical systems pose tremendous computational difficulties when applied in numerical simulations. In order to overcome these challenges we use several model reduction techniques. For stable linear models these techniques work very well and provide good approximations for the full model. However, large-scale unstable systems arise in many applications. Many of the known model reduction methods are not very robust, or in some cases, may not even work if we are dealing with unstable systems. When approximating an unstable system by a reduced order model, accuracy is not the only concern. We also need to consider the structure of the reduced order model. Often, it is important that the number of unstable poles in the reduced system is the same as the number of unstable poles in the original system. The Iterative Rational Krylov Algorithm (IRKA) is a robust model reduction technique which is used to locally reduce stable linear dynamical systems optimally in the ℋ₂-norm. While we cannot guarantee that IRKA reduces an unstable model optimally, there are no numerical obstacles to the reduction of an unstable model via IRKA. In this thesis, we investigate IRKA's behavior when it is used to reduce unstable models. We also consider systems for which we cannot obtain a first order realization of the transfer function. We can use Realization-independent IRKA to obtain a reduced order model which does not preserve the structure of the original model. In this paper, we implement a structure preserving algorithm for systems with nonlinear frequency dependency. / Master of Science

Published Article / The role of the Projection Principle within Chomsky's Government-Binding (GB) Theory is to preserve the subcategorisation properties of lexical items at all levels of syntactic representation, viz. D-structure, S-structure, and Lexical Form. Arguments have been made that the Projection Principle is a new concept that is simply an extension of theTransformational Component (XFM) and Emonds' Structure-Preserving Constraint (SPC), and that it does not deserve the high status it has been accorded in GB theory. This paper provides evidence, based on sentences involving movement operations, that the Projection Principle is innovative and that it convincingly addresses what theXFMandSPChave failed to address.

This dissertation studies the model reduction and distributed control problems for interconnected systems, i.e., systems that consist of multiple interacting agents/subsystems. The study of the analysis and synthesis problems for interconnected systems is motivated by the multiple applications that can benefit from the design and implementation of distributed controllers. These applications include automated highway systems and formation flight of unmanned aircraft systems. The systems of interest are modeled using arbitrary directed graphs, where the subsystems correspond to the nodes, and the interconnections between the subsystems are described using the directed edges. In addition to the states of the subsystems, the adopted frameworks also model the interconnections between the subsystems as spatial states. Each agent/subsystem is assumed to have its own actuating and sensing capabilities. These capabilities are leveraged in order to design a controller subsystem for each plant subsystem. In the distributed control paradigm, the controller subsystems interact over the same interconnection structure as the plant subsystems. The models assumed for the subsystems are linear time-varying or linear parameter-varying. Linear time-varying models are useful for describing nonlinear equations that are linearized about prespecified trajectories, and linear parameter-varying models allow for capturing the nonlinearities of the agents, while still being amenable to control using linear techniques. It is clear from the above description that the size of the model for an interconnected system increases with the number of subsystems and the complexity of the interconnection structure. This motivates the development of model reduction techniques to rigorously reduce the size of the given model. In particular, this dissertation presents structure-preserving techniques for model reduction, i.e., techniques that guarantee that the interpretation of each state is retained in the reduced order system. Namely, the sought reduced order system is an interconnected system formed by reduced order subsystems that are interconnected over the same interconnection structure as that of the full order system. Model reduction is important for reducing the computational complexity of the system analysis and control synthesis problems. In this dissertation, interior point methods are extensively used for solving the semidefinite programming problems that arise in analysis and synthesis. / Ph. D. / The work in this dissertation is motivated by the numerous applications in which multiple agents interact and cooperate to perform a coordinated task. Examples of such applications include automated highway systems and formation flight of unmanned aircraft systems. For instance, one can think of the hazardous conditions created by a fire in a building and the benefits of using multiple interacting multirotors to deal with this emergency situation and reduce the risks on humans. This dissertation develops mathematical tools for studying and dealing with these complex systems. Namely, it is shown how controllers can be designed to ensure that such systems perform in the desired way, and how the models that describe the systems of interest can be systematically simplified to facilitate performing the tasks of mathematical analysis and control design.

This dissertation develops a variety of structure-preserving algorithms for mechanical systems with external forcing and also extends those methods to systems that evolve on non-Euclidean manifolds. The dissertation is focused on numerical schemes derived from variational principles – schemes that are general enough to apply to a large class of engineering problems. A theoretical framework that encapsulates variational integration for mechanical systems with external forcing and time-dependence and which supports the extension of these methods to systems that evolve on non-Euclidean manifolds is developed. An adaptive time step, energy-preserving variational integrator is developed for mechanical systems with external forcing. It is shown that these methods track the change in energy more accurately than their fixed time step counterparts. This approach is also extended to rigid body systems evolving on Lie groups where the resulting algorithms preserve the geometry of the configuration space in addition to being symplectic as well as energy and momentum-preserving. The advantages of structure-preservation in the numerical simulation are illustrated by various representative examples from engineering applications, which include limit cycle oscillations of an aeroelastic system, dynamics of a neutrally buoyant underwater vehicle, and optimization for spherical shape correlation and matching. / Doctor of Philosophy / Accurate numerical simulation of dynamical systems over long time horizons is essential in applications ranging from particle physics to geophysical fluid flow to space hazard analysis. In many of these applications, the governing physical equations derive from a variational principle and their solutions exhibit physically meaningful invariants such as momentum, energy, or vorticity. Unfortunately, most traditional numerical methods do not account for the underlying geometric structure of the physical system, leading to simulation results that may suggest nonphysical behavior. In this dissertation, tools from geometric mechanics and computational methods are used to develop numerical integrators that respect the qualitative features of the physical system. The research presented here focuses on numerical schemes derived from variational principles– schemes that are general enough to apply to a large class of engineering problems. Energy-preserving algorithms are developed for mechanical systems by exploiting the underlying geometric properties. Numerical performance comparisons demonstrate that these algorithms provide almost exact energy preservation and lead to more accurate prediction. The advantages of these methods in the numerical simulation are illustrated by various representative examples from engineering applications, which include limit cycle oscillations of an aeroelastic system, dynamics of a neutrally buoyant underwater vehicle, and optimization for spherical shape correlation and matching.

This dissertation embodies a comprehensive approach towards resilient monitoring of frid events using Structure Preserving Energy Functions (SPEFs) and introduces a novel control architecture in Multi Terminal Direct Current (MTDC) grids, for inter-area oscillation damping and achieving robustness to AC as well as DC side, post-contingency events in the modern power grid. This work is presented as a collection of several publications which investigate and address the proposed research topics. At first, SPEFs are derived for multi-machine IEEE benchmark models with the help of the Wide-Area Measurements (WAMs). A physics-based hybrid approach to develop one-to-one mapping between properties of energy function components with respect to the type of fault in the system is introduced. The proposed method is tested offline on a IEEE-39 bus, New England Test System (NETS), with particular interest in monitoring the most sensitive energy functions during relay misoperations. Such events can be precipitated by zone 3 trips in distance relays due to load encroachment during stressed conditions. These might include a genuine misoperation, a false trip due to cyber-attacks, or a load encroachment, all of which are undesirable under normal operating circumstances. An online monitoring scheme is introduced in an actual blackout simulation in the Western Electricity Coordinating Council (WECC) to examine what further indications these energy function components can provide, especially during a cascading sequence, and how they could supervise critical tripping decisions by distance relays. Next, a futuristic grid comprised of Voltage Source Converter (VSC) based AC-MTDC is considered due to its recent proliferation in integrating offshore wind farms to onshore grids, and additionally improving strength of weak AC grids. A robust control is designed using the converter station poles as actuators to provide damping support to the surrounding AC grid. Further, a design problem is envisioned and implemented that introduces disturbance rejection into control architecture by designing a novel explicitly modeled disturbance plant in the Linear Matrix Inequality (LMI) framework. Finally, a novel robust inter-area oscillation damping controller is designed that proves its effectiveness in inter-area mode settling times, and provides robustness to (n-1) contingencies in the AC as well as the DC side of the meshed AC-MTDC grid.

A new spline space, the hierarchical structure-preserving B-spline space, is introduced and implemented in the analysis of Stokes flow. The space, when properly constrained, is shown to be stable and to have at least optimal convergence rates in the velocity field and suboptimal convergence rates in the pressure field. However, results show that superoptimal convergence can often be expected in the pressure field, likely due to error reduction in the velocity field. Like other hierarchical spline spaces, these splines are shown to greatly increase accuracy and to drastically lower computation times for analyses on domains whose solution fields have singularities or could otherwise benefit from local refinement. With the advent of this adaptive, locally-refineable, high-fidelity technology, isogeometric methods can become more feasible for use in fluid analyses.

<p>This thesis consists of the structure-preserving numerical methods for PNP-NS equation and dynamic liquid crystal systems in Oseen-Frank energy. </p> <p>In Chapter 1, we give a brief introduction of the Poisson-Nernst-Planck-Navier-Stokes (PNP-NS) system, and the dynamical liquid system in Oseen-Frank energy in one-constant approximation case and a special non-one-constant case. Each of those systems has a special structure and properties we want to keep at the discrete level when designing numerical methods.</p> <p>In Chapter 2, we introduce a first-order numerical scheme for the PNP-NS system that is decoupled, positivity preserving, mass conserving, and unconditionally energy stable. The numerical scheme is designed in the context of Wasserstein gradient flow based on the form ∇ · (c∇ ln c). The mobility terms are treated explicitly, and the chemical potential terms are treated implicitly so that the solution of the numerical scheme is the minimizer of a convex functional, which is the key to the unique solvability and positivity preserving of the numer-<br> ical scheme. Proper boundary conditions for chemical potentials are chosen to guarantee the mass-conservation property. The convection term in Poisson-Nernst-Planck(PNP) equation part is treated explicitly with an O(∆t) term introduced so that the numerical scheme is decoupled and unconditionally energy stable. Pressure correction methods are used for the Navier-Stokes(NS) equation part. And we proved the optimal convergence rate with an irregular high-order asymptotic expansion technique.</p> <p>In Chapter 3, we propose a first-order implicit numerical method for a dynamic liquid crystal system in a one-constant-approximation case(which is also known as heat flow of harmonic maps to S2). The solution is the minimizer of a convex functional under the unit length constraint, and from this point, the weak convergence of the numerical scheme could be proved. The numerical scheme is solved in an iterative procedure. This procedure could be proved to be energy decreasing and this implies the convergence of the algorithm.</p> <p>In Chapter 4, we study the dynamic liquid crystal system in a more generalized Oseen- Frank energy compared to Chapter 3. We are assuming K2 = K3 = −K4, the domain Ω is a rectangular region in R3, and d satisfies the periodic boundary condition on ∂Ω. And we propose a class of numerical schemes for this system that preserve the unit length constraint. The convergence of the numerical scheme has been proved under necessary assumptions. And numerical experiments are presented to validate the accuracy and demonstrate the performance of the proposed numerical scheme.</p>

Les méthodes d'algèbres linéaire classiques, pour le calcul de valeurs et vecteurs propres d'une matrice, ou des approximations de rangs inférieurs (low-rank approximations) d'une solution, etc..., ne tiennent pas compte des structures de matrices. Ces dernières sont généralement détruites durant le procédé du calcul. Des méthodes alternatives préservant ces structures font l'objet d'un intérêt important par la communauté. Cette thèse constitue une contribution dans ce domaine. La décomposition SR peut être calculé via l'algorithme de Gram-Schmidt symplectique. Comme dans le cas classique, une perte d'orthogonalité peut se produire. Pour y remédier, nous avons proposé deux algorithmes RSGSi et RMSGSi qui consistent à ré-orthogonaliser deux fois les vecteurs à calculer. La perte de la J-orthogonalité s'est améliorée de manière très significative. L'étude directe de la propagation des erreurs d'arrondis dans les algorithmes de Gram-Schmidt symplectique est très difficile à effectuer. Nous avons réussi à contourner cette difficulté et donner des majorations pour la perte de la J-orthogonalité et de l'erreur de factorisation. Une autre façon de calculer la décomposition SR est basée sur les transformations de Householder symplectique. Un choix optimal a abouti à l'algorithme SROSH. Cependant, ce dernier peut être sujet à une instabilité numérique. Nous avons proposé une version modifiée nouvelle SRMSH, qui a l'avantage d'être aussi stable que possible. Une étude approfondie a été faite, présentant les différentes versions : SRMSH et SRMSH2. Dans le but de construire un algorithme SR, d'une complexité d'ordre O(n³) où 2n est la taille de la matrice, une réduction (appropriée) de la matrice à une forme condensée (J(Hessenberg forme) via des similarités adéquates, est cruciale. Cette réduction peut être effectuée via l'algorithme JHESS. Nous avons montré qu'il est possible de réduire une matrice sous la forme J-Hessenberg, en se basant exclusivement sur les transformations de Householder symplectiques. Le nouvel algorithme, appelé JHSJ, est basé sur une adaptation de l'algorithme SRSH. Nous avons réussi à proposer deux nouvelles variantes, aussi stables que possible : JHMSH et JHMSH2. Nous avons constaté que ces algorithmes se comportent d'une manière similaire à l'algorithme JHESS. Une caractéristique importante de tous ces algorithmes est qu'ils peuvent rencontrer un breakdown fatal ou un "near breakdown" rendant impossible la suite des calculs, ou débouchant sur une instabilité numérique, privant le résultat final de toute signification. Ce phénomène n'a pas d'équivalent dans le cas Euclidien. Nous avons réussi à élaborer une stratégie très efficace pour "guérir" le breakdown fatal et traîter le near breakdown. Les nouveaux algorithmes intégrant cette stratégie sont désignés par MJHESS, MJHSH, JHM²SH et JHM²SH2. Ces stratégies ont été ensuite intégrées dans la version implicite de l'algorithme SR lui permettant de surmonter les difficultés rencontrées lors du fatal breakdown ou du near breakdown. Rappelons que, sans ces stratégies, l'algorithme SR s'arrête. Finalement, et dans un autre cadre de matrices structurées, nous avons présenté un algorithme robuste via FFT et la matrice de Hankel, basé sur le calcul approché de plus grand diviseur commun (PGCD) de deux polynômes, pour résoudre le problème de la déconvolution d'images. Plus précisément, nous avons conçu un algorithme pour le calcul du PGCD de deux polynômes bivariés. La nouvelle approche est basée sur un algorithme rapide, de complexité quadratique O(n²), pour le calcul du PGCD des polynômes unidimensionnels. La complexité de notre algorithme est O(n²log(n)) où la taille des images floues est n x n. Les résultats expérimentaux avec des images synthétiquement floues illustrent l'efficacité de notre approche. / The classical linear algebra methods, for calculating eigenvalues and eigenvectors of a matrix, or lower-rank approximations of a solution, etc....do not consider the structures of matrices. Such structures are usually destroyed in the numerical process. Alternative structure-preserving methods are the subject of an important interest mattering to the community. This thesis establishes a contribution in this field. The SR decomposition is usually implemented via the symplectic Gram-Schmidt algorithm. As in the classical case, a loss of orthogonality can occur. To remedy this, we have proposed two algorithms RSGSi and RMSGSi, where the reorthogonalization of a current set of vectors against the previously computed set is performed twice. The loss of J-orthogonality has significantly improved. A direct rounding error analysis of symplectic Gram-Schmidt algorithm is very hard to accomplish. We managed to get around this difficulty and give the error bounds on the loss of the J-orthogonality and on the factorization. Another way to implement the SR decomposition is based on symplectic Householder transformations. An optimal choice of free parameters provided an optimal version of the algorithm SROSH. However, the latter may be subject to numerical instability. We have proposed a new modified version SRMSH, which has the advantage of being numerically more stable. By a detailes study, we are led to two new variants numerically more stables : SRMSH and SRMSH2. In order to build a SR algorithm of complexity O(n³), where 2n is the size of the matrix, a reduction to the condensed matrix form (upper J-Hessenberg form) via adequate similarities is crucial. This reduction may be handled via the algorithm JHESS. We have shown that it is possible to perform a reduction of a general matrix, to an upper J-Hessenberg form, based only on the use of symplectic Householder transformations. The new algorithm, which will be called JHSH algorithm, is based on an adaptation of SRSH algorithm. We are led to two news variants algorithms JHMSH and JHMSH2 which are significantly more stable numerically. We found that these algortihms behave quite similarly to JHESS algorithm. The main drawback of all these algorithms (JHESS, JHMSH, JHMSH2) is that they may encounter fatal breakdowns or may suffer from a severe form of near-breakdowns, causing a brutal stop of the computations, the algorithm breaks down, or leading to a serious numerical instability. This phenomenon has no equivalent in the Euclidean case. We sketch out a very efficient strategy for curing fatal breakdowns and treating near breakdowns. Thus, the new algorithms incorporating this modification will be referred to as MJHESS, MJHSH, JHM²SH and JHM²SH2. These strategies were then incorporated into the implicit version of the SR algorithm to overcome the difficulties encountered by the fatal breakdown or near-breakdown. We recall that without these strategies, the SR algorithms breaks. Finally ans in another framework of structured matrices, we presented a robust algorithm via FFT and a Hankel matrix, based on computing approximate greatest common divisors (GCD) of polynomials, for solving the problem pf blind image deconvolution. Specifically, we designe a specialized algorithm for computing the GCD of bivariate polynomials. The new algorithm is based on the fast GCD algorithm for univariate polynomials , of quadratic complexity O(n²) flops. The complexitiy of our algorithm is O(n²log(n)) where the size of blurred images is n x n. The experimental results with synthetically burred images are included to illustrate the effectiveness of our approach