Simulation, design, and control of dynamical systems play an important role in numerous scientific and industrial tasks. The need for detailed models leads to large-scale dynamical systems, posing tremendous computational difficulties when employed in numerical simulations. In order to overcome these challenges, we perform model reduction, replacing the large-scale dynamics with high-fidelity reduced representations. There exist a plethora of methods for reduced order modeling of linear systems, including the Iterative Rational Krylov Algorithm (IRKA), Balanced Truncation (BT), and Hankel Norm Approximation. However, these methods generally target stable systems and the approximation is performed over an infinite time horizon. If we are interested in a finite horizon reduced model, we utilize techniques such as Time-limited Balanced Truncation (TLBT) and Proper Orthogonal Decomposition (POD). In this dissertation we establish interpolation-based optimality conditions over a finite horizon and develop an algorithm, Finite Horizon IRKA (FHIRKA), that produces a locally optimal reduced model on a specified time-interval. Nonetheless, the quantities being interpolated and the interpolant are not the same as in the infinite horizon case. Numerical experiments comparing FHIRKA to other algorithms further support our theoretical results. Next, we discuss model reduction for nonlinear dynamical systems. For models with unstructured nonlinearities, POD is the method of choice. However, POD is input dependent and not optimal with respect to the output. Thus, we use operator splitting to integrate the best features of system theoretic approaches with trajectory based methods such as POD in order to mitigate the effect of the control inputs for the approximation of nonlinear dynamical systems. We reduce the linear terms with system theoretic methods and the nonlinear terms terms via POD. Evolving the linear and nonlinear terms separately yields the reduced operator splitting solution. We present an error analysis for this method, as well as numerical results that illustrate the effectiveness of our approach. While in this dissertation we only pursue the splitting of linear and nonlinear terms, this approach can be implemented with Quadratic Bilinear IRKA or Balanced Truncation for Quadratic Bilinear systems to further diminish the input dependence of the reduced order modeling. / Doctor of Philosophy / Simulation, design, and control of dynamical systems play an important role in numerous scientific and industrial tasks such as signal propagation in the nervous system, heat dissipation, electrical circuits and semiconductor devices, synthesis of interconnects, prediction of major weather events, spread of fires, fluid dynamics, machine learning, and many other applications. The need for detailed models leads to large-scale dynamical systems, posing tremendous computational difficulties when applied in numerical simulations. In order to overcome these challenges, we perform model reduction, replacing the large-scale dynamics with high-fidelity reduced representations. Reduced order modeling helps us to avoid the outstanding burden on computational resources. Numerous model reduction techniques exist for linear models over an infinite horizon. However, in practice we usually are interested in reducing a model over a specific time interval. In this dissertation, given a reduced order, we present a method that finds the best local approximation of a dynamical system over a finite horizon. We present both theoretical and numerical evidence that supports the proposed method. We also develop an algorithm that integrates operator splitting with model reduction to solve nonlinear models more efficiently while preserving a high level of accuracy.

In this paper we propose an algorithm for solving systems of coupled monotone inclusions in Hilbert spaces. The operators arising in each of the inclusions of the system are processed in each iteration separately, namely, the single-valued are evaluated explicitly (forward steps), while the set-valued ones via their resolvents (backward steps). In addition, most of the steps in the iterative scheme can be executed simultaneously, this making the method applicable to a variety of convex minimization problems. The numerical performances of the proposed splitting algorithm are emphasized through applications in average consensus on colored networks and image classification via support vector machines.

No description available.

Nous abordons le développement d'une nouvelle génération de méthodes numériques pour la résolution des EDP évolutives qui modélisent des phénomènes multi-échelles en temps et en espace issus de divers domaines applicatifs. La raideur associée à ce type de problème, que ce soit via le terme source chimique qui présente un large spectre d'échelles de temps caractéristiques ou encore via la présence de fort gradients très localisés associés aux fronts de réaction, implique en général de sévères difficultés numériques. En conséquence, il s'agit de développer des méthodes qui garantissent la précision des résultats en présence de forte raideur en s'appuyant sur des outils théoriques solides, tout en permettant une implémentation aussi efficace. Même si nous étendons ces idées à des systèmes plus généraux par la suite, ce travail se focalise sur les systèmes de réaction-diffusion raides. La base de la stratégie numérique s'appuie sur une décomposition d'opérateur spécifique, dont le pas de temps est choisi de manière à respecter un niveau de précision donné par la physique du problème, et pour laquelle chaque sous-pas utilise un intégrateur temporel d'ordre élevé dédié. Ce schéma numérique est ensuite couplé à une approche de multirésolution spatiale adaptative permettant une représentation de la solution sur un maillage dynamique adapté. L'ensemble de cette stratégie a conduit au développement du code de simulation générique 1D/2D/3D académique MBARETE de manière à évaluer les développements théoriques et numériques dans le contexte de configurations pratiques raides issue de plusieurs domaines d'application. L'efficacité algorithmique de la méthode est démontrée par la simulation d'ondes de réaction raides dans le domaine de la dynamique chimique non-linéaire et dans celui de l'ingénierie biomédicale pour la simulation des accidents vasculaires cérébraux caractérisée par un terme source "chimique complexe''. Pour étendre l'approche à des applications plus complexes et plus fortement instationnaires, nous introduisons pour la première fois une technique de séparation d'opérateur avec pas de temps adaptatif qui permet d'atteindre une précision donnée garantie malgré la raideur des EDP. La méthode de résolution adaptative en temps et en espace qui en résulte, étendue au cas convectif, permet une description consistante de problèmes impliquant une très large palette d'échelles de temps et d'espace et des scénarios physiques très différents, que ce soit la propagation des décharges répétitives pulsées nanoseconde dans le domaine des plasmas ou bien l'allumage et la propagation de flammes dans celui de la combustion. L'objectif de la thèse est l'obtention d'un solveur numérique qui permet la résolution des EDP raides avec contrôle de la précision du calcul en se basant sur des outils d'analyse numérique rigoureux, et en utilisant des moyens de calculs standard. Quelques études complémentaires sont aussi présentées comme la parallélisation temporelle, des techniques de parallélisation à mémoire partagée et des outils de caractérisation mathématique des schémas de type séparation d'opérateur. / We tackle the development of a new generation of numerical methods for the solution of time dependent PDEs modeling general time/space multi-scale phenomena issued from various application fields. This type of problem induces well-known numerical restrictions and potentially large stiffness, which stem from the broad spectrum of time scales in the nonlinear chemical terms as well as from steep, spatially very localized, spatial gradients in the reaction fronts. Therefore, dedicated numerical strategies are needed to ensure the accuracy of the numerical approximations from a theoretical point of view, taking also into account adequate practical implementations to reduce computational costs. In order to cope with these problems, this study introduces a few mathematical and numerical elements for the solution of stiff reaction-diffusion systems, extensible in practice to more general configurations. The core of the numerical strategy is thus based on a specially conceived operator splitting method with dedicated high order time integration schemes for each subproblem. An appropriate choice of splitting time steps allows us the simulation of the solution within a prescribed accuracy, according to the overall physics of the problem. The resulting numerical scheme is properly coupled with an adaptive multiresolution technique for dynamic spatial mesh representations of the solution. Such an approach has led to the conception of the academic, generic 1D/2D/3D MBARETE code in order to evaluate the proposed theoretical and numerical developments in practical stiff configurations arising in several research fields. The algorithmic efficiency of the method is assessed by the simulation of propagating stiff reaction waves issued from nonlinear chemical dynamics and from biomedical engineering applications for a brain stroke model with "detailed chemical mechanisms''. Moreover, in order to extend the applicability of the method to more complex and unsteady problems, we consider for the first time a time adaptive splitting scheme for stiff PDEs, that yields dynamic time stepping within the prescribed accuracy. The fully time/space adaptive method allows us then a consistent description of reaction-diffusion-convection problems disclosing a broad spectrum of time/space scales as well as different physical scenarios, such as highly nanosecond repetitively pulsed discharges or self-ignition and propagation of flames for, respectively, plasma and combustion applications. The main goal of this work is hence to numerically solve stiff PDEs with reasonable, standard computational resources and based on a mathematical background that ensures robust, general and accurate numerical schemes. Further studies are also presented that include time parallelization strategies, parallel computing techniques for shared memory architectures and complementary mathematical characterization of splitting schemes.

A non-iterative pressure based algorithm which consists of splitting the solution of momentum energy and species equations into a sequence of predictor-corrector stages was developed for the simulation of transient reacting radiating flows. A semi-discrete approach called the Method of Lines (MOL) which enables implicit time-integration at all splitting stages was used for the solution of conservation equations. The solution of elliptic pressure equation for the determination of pressure field was performed by a multi-grid solver (MUDPACK package). Radiation calculations were carried out by coupling previously developed gray and non-gray radiation models with the algorithm. A first order (global) reaction mechanism was employed to account for the chemistry. The predictions of the algorithm for the following test cases: i) non-isothermal turbulent pipe flow and ii) laminar methane-air diffusion flame / were benchmarked against experimental data and numerical solutions available in the literature and the capability of the code to predict transient solutions was demonstrated on these test cases. Favorable agreements were obtained for both test cases. The effect of radiation and non-gray treatment of the radiative properties were investigated on the second test case. It was found that incorporation of radiation has significant effect on Temeprature and velocity fields but its effect is limited in species predictions. Executions with both radiation models revealed that the non-gray radiation model considered in the present study produces similar results with the gray model at a considerably higher computational cost. The algorithm developed was found to be an efficient and versatile tool for the timedependent simulation of different flow scenarios constitutes the initial steps towards the computation of transient turbulent combustion.

This thesis studies numerical methods for integrating the master equations describing Markov chain models of cardiac ion channels. Such models describe the time evolution of the probability that ion channels are in a particular state. Numerical simulations of such models are often computationally demanding because many solvers require relatively small time steps to ensure numerical stability. The aim of this project is to analyse selected Markov chains and develop more efficient and accurate solvers. We separate a Markov chain model into fast and slow time-scales based on the speed of transitions between states. Eliminating the fast transitions, we find an asymptotic reduction of zeroth-order and first-order in a small parameter describing the time-scales separation. We apply the theory to a Markov chain model of the fast sodium channel INa. We consider several variants for classifying some transitions as fast in order to find reduced systems that yield a good accuracy. However, the time step size is still restricted by numerical instabilities. We adapt the Rush-Larsen technique originally developed for gate models. Assuming that a transition matrix can be considered constant during each time step, we solve the Markov chain model analytically. The solution provides a recipe for a stable exponential solver, which we call "Matrix Rush-Larsen" (MRL). Using operator splitting we design an even more flexible "hybrid" method that combines the MRL with other solvers. The resulting improvement in stability allows a large increase in the time step size. In some models, we obtain reasonably accurate results 27 times faster using a hybrid method than with the forward Euler method, even with the maximal time step allowed by the stability constraint. Finally, we extend the cardiac simulation package BeatBox by the developed exponential solvers. We upgrade a format of "ionic" modules which describe a cardiac cell, in order to allow for a specific definition of Markov chain models. We also modify a particular integrator for ionic modules to include the MRL and the hybrid method. To test the functionality of the code, we have converted a number of cellular models into the ionic format. The documented code is available in the official BeatBox package distribution.

Ce travail décrit le développement d’un outil de simulation du transport réactif, nommé SpeCTr (Spéciation Cinétique Transport), intégrant le fractionnement isotopique. Ce modèle est obtenu à travers le couplage d’un module décrivant le transport et d’un module décrivant les principales réactions chimiques (approche de séparation d’opérateur). Une grande partie du travail est dédiée à l’amélioration des algorithmes du module décrivant les réactions chimiques pour la résolution des équations de l’équilibre thermodynamique (méthode de Newton Raphson modifiée à travers les techniques du scaling et des Fractions Continues Positives) et du mélange de réactions cinétiques et à l’équilibre (étude sur la formulation et résolution des systèmes d’équations différentielles et différentielles-algébriques). L’outil est validé à travers la résolution de plusieurs tests (batch et transport réactif) et appliqué pour la simulation d’expériences de laboratoire en 1D, 2D et 3D portant sur la dissolution des cristaux de calcite dans une colonne de milieu poreux / The work describes the development of a reactive transport code named SpeCTr (Spéciation Cinétique Transport in French). The code, able to describe isotopic fractionation, is obtained through the coupling of a transport module and a reaction module that describes the main chemical reactions (operator splitting approach). A consistent portion of the work is dedicated to the improvement of the numerical methods employed in the reaction module for solving thermodynamic equilibrium (Newton Raphson method modified with scaling and Positive Continuous Fractions) and mixed equilibrium and kinetic reactions (formulation and solution of systems of differential and differential-algebraic equations). The code was verified through the solution of different benchmarks (batch and reactive transport simulations) and applied to perform 1D, 2D and 3D simulations of laboratory experiments dedicated to calcite crystals dissolution in a column of porous medium.

Die Anwenung von Operator-Splitting-Methoden auf gewöhnliche Differentialgleichungen ist gut etabliert. Für Differential-algebraische Gleichungen und partielle Differential-algebraische Gleichungen unterliegt sie jedoch vielen Einschränkungen aufgrund des Vorhandenseins von Nebenbedingungen. Die räumliche Diskretisierung reduziert PDAEs und lenkt unseren Fokus auf das Konzept der DAEs. Um eine reibungslose Übertragung des Operator-Splittings von ODEs auf DAEs durchzuführen, ist es wichtig, eine geeignete entkoppelte Struktur für das gewünschte Differential-algebraische System zu haben. In dieser Arbeit betrachten wir ein Modell, das partielle Differentialgleichungen für elektromagnetische Bauelemente - modelliert durch die Maxwell-Gleichungen - mit Differential-algebraischen Gleichungen koppelt, die die elementaren Schaltungselemente beschreiben. Nach der räumlichen Diskretisierung der klassischen Formulierung der Maxwell-Gleichungen mit Hilfe der finiten Integrationstechnik formulieren wir das resultierende gekoppelte System als Differential-algebraische Gleichung. Um eine geeignete Entkopplung zu bekommen, verwenden wir den zweigorientierten Loop-Cutset-Ansatz für die Schaltungsmodellierung. Daraus folgt, dass wir in der Lage sind, eine geeignete Operatorzerlegung so zu konstruieren, dass wir eine natürliche topologisch entkoppelte Port-Hamiltonsche DAE-Struktur erhalten. Wir schlagen einen Operator-Splitting-Ansatz für die Schaltungs-DAEs und gekoppelten Feld-Schaltungs-DAEs in entkoppelter Form vor und analysieren seine numerischen Eigenschaften. Darüber hinaus nutzen wir das Hamiltonsche Verhalten der inhärenten gewöhnlichen Differentialgleichung durch die Verwendung expliziter und energieerhaltender Zeitintegrations-methoden. Schließlich führen wir numerische Tests, um das mathematische Modell zu illustrieren und die Konvergenzergebnisse für das vorgeschlagene DAE-Operator-Splitting zu demonstrieren. / Le equazioni algebriche differenziali e algebriche alle derivate parziali hanno avuto un enorme successo come modelli di sistemi dinamici vincolati. Nella modellazione matem- atica, spesso si desidera catturare diversi aspetti di una situazione come le leggi di conservazione della fisica, il trasporto convettivo o la diffusione. Queste aspetti si riflettono nel sistema di equazioni del modello come operatori diversi. La tecnica dell’Operator Splitting si è rivelata una strategia di successo per affrontare problemi così complicati. L’applicazione dei metodi di Operator Splitting alle equazioni differenziali ordinarie (ODE) è ormai una tecnologia ben consolidata. Tuttavia, per equazioni algebriche differenziali (DAE) e algebriche differenziali parziali (PDAE), l’approccio è soggetto a molte restrizioni dovute alla presenza di vincoli e alla proprietà di indice. La discretizzazione spaziale riduce le PDAE e indirizza la nostra attenzione al concetto di DAE. Le DAE emergono in problemi dinamici vincolati come circuiti elettrici o reti di trasporto di energia. Al fine di generalizzare agevolmente la tecnica dell’Operator Splitting dalle ODE alle DAE, è importante avere una struttura disaccoppiata adeguata per il sistema algebrico differenziale desiderato. In questa tesi, consideriamo un modello che accoppia equazioni differenziali alle derivate parziali per dispositivi elettromagnetici -modellati dalle equazioni di Maxwell- con equazioni algebriche differenziali che descrivono gli elementi base del circuito. Dopo aver discretizzato spazialmente la formulazione classica delle equazioni di Maxwell usando la tecnica di integrazione finita, formuliamo il sistema accoppiato risultante come una equazione algebrica differenziale. Interpretando il dispositivo elettromagnetico come un elemento capacitivo, l’indice dell’intero sistema di circuito e campo accoppiato può essere specificato utilizzando le proprietà topologiche del circuito e non supera il valore di due. Per eseguire un disaccoppiamento appropriato, utilizziamo l’approccio loop-cutset per la modellazione dei circuiti. In tal modo siamo in grado di costruire una opportuna decomposizione dell’operatore tale da ottenere una naturale struttura disaccoppiata port-Hamiltonian DAE. Proponiamo un approccio di suddivisione dell’operatore per i DAE a circuito disaccoppiato e a circuito di campo accoppiato utilizzando gli algoritmi di divisione Lie-Trotter e Strang e per analizzare le proprietà numeriche di questi sistemi. Inoltre, sfruttiamo il comportamento hamiltoniano del sistema di equazioni differenziali ordinarie mediante l’utilizzo di metodi di integrazione temporale con esatta conservazione dell’energia. Poggiando sull’analisi di convergenza del metodo di suddivisione dell’operatore ODE, deriviamo i risultati di convergenza per l’approccio proposto che dipendono dall’indice delsistema e quindi dalla sua struttura topologica. Infine, eseguiamo prove numeriche di sistemi circuitali, nonchè sistemi accoppiati a circuito di campo, per testare il modello matematico e dimostrare i risultati di convergenza per la proposta Operator Splitting DAE. / The application of operator splitting methods to ordinary differential equations (ODEs) is well established. However, for differential-algebraic equations (DAEs) and partial differential-algebraic equations (PDAEs), it is subjected to many restrictions due to the presence of constraints. In constrained dynamical problems as electrical circuits or energy transport networks, DAEs arise. In order to perform a smooth transfer of the operator splitting from ODEs to DAEs, it is important to have a suitable decoupled structure for the desired differential-algebraic system. In this thesis, we consider a model which couples partial differential equations for electro- magnetic devices -modeled by Maxwell’s equations- with differential-algebraic equations describing the basic circuit elements. After spatially discretizing the classical formulation of Maxwell’s equations using the finite integration technique, we formulate the resulting coupled system as a differential-algebraic equation. To perform an appropriate decoupling, we use the branch oriented loop-cutset approach for circuit modeling. It follows that we are able to construct a suitable operator decomposition such that we obtain a natural topologically decoupled port-Hamiltonian DAE structure. We propose an operator splitting approach for the decoupled circuit and coupled field- circuit DAEs using the Lie-Trotter and Strang splitting algorithms and analyze its numerical properties. Furthermore, we exploit the Hamiltonian behavior of the system’s inherent ordinary differential equation by the utilization of explicit and energy-preserving time integration methods. Based on the convergence analysis of the ODE operator splitting method, we derive convergence results for the proposed approach that depends on the index of the system and thus on its topological structure. Finally, we perform numerical tests, to underline the mathematical model and to demonstrate the convergence results for the proposed DAE operator splitting.

The heart is one the most important organs in the human body and many other live creatures. The electrical activity in the heart controls the heart function, and many heart diseases are linked to the abnormalities in the electrical activity in the heart. Mathematical equations and computer simulation can be used to model the electrical activity in the heart. The heart models are challenging to solve because of the complexity of the models and the huge size of the problems. Several cell models have been proposed to model the electrical activity in a single heart cell. These models must be coupled with a heart model to model the electrical activity in the entire heart. The bidomain model is a popular model to simulate the propagation of electricity in myocardial tissue. It is a continuum-based model consisting of non-linear ordinary differential equations (ODEs) describing the electrical activity at the cellular scale and a system of partial differential equations (PDEs) describing propagation of electricity at the tissue scale. Because of this multi-scale, ODE/PDE structure of the model, splitting methods that treat the ODEs and PDEs in separate steps are natural candidates as numerical methods. First, we need to solve the problem at the cellular scale using ODE solvers. One of the most popular methods to solve the ODEs is known as the Rush-Larsen (RL) method. Its popularity stems from its improved stability over integrators such as the forward Euler (FE) method along with its easy implementation. The RL method partitions the ODEs into two sets: one for the gating variables, which are treated by an exponential integrator, and another for the remaining equations, which are treated by the FE method. The success of the RL method can be understood in terms of its relatively good stability when treating the gating variables. However, this feature would not be expected to be of benefit on cell models for which the stiffness is not captured by the gating equations. We demonstrate that this is indeed the case on a number of stiff cell models. We further propose a new partitioned method based on the combination of a first-order generalization of the RL method with the FE method. This new method leads to simulations of stiff cell models that are often one or two orders of magnitude faster than the original RL method. After solving the ODEs, we need to use bidomain solvers to solve the bidomain model. Two well-known, first-order time-integration methods for solving the bidomain model are the semi-implicit method and the Godunov operator-splitting method. Both methods decouple the numerical procedure at the cellular scale from that at the tissue scale but in slightly different ways. The methods are analyzed in terms of their accuracy, and their relative performance is compared on one-, two-, and three-dimensional test cases. As suggested by the analysis, the test cases show that the Godunov method is significantly faster than the semi-implicit method for the same level of accuracy, specifically, between 5 and 15 times in the cases presented. Second-order bidomain solvers can generally be expected to be more effective than first-order bidomain solvers under normal accuracy requirements. However, the simplest and the most commonly applied second-order method for the PDE step, the Crank-Nicolson (CN) method, may generate unphysical oscillations. We investigate the performance of a two-stage, L-stable singly diagonally implicit Runge-Kutta method for solving the PDEs of the bidomain model and present a stability analysis. Numerical experiments show that the enhanced stability property of this method leads to more physically realistic numerical simulations compared to both the CN and Backward Euler (BE) methods.

Distributed optimization over multi-agent networks has become an increasingly popular research topic as it incorporates many applications from various areas such as consensus optimization, distributed control, network resource allocation, large scale machine learning, etc. Parallel distributed solution algorithms are highly desirable as they are more scalable, more robust against agent failure, align more naturally with either underlying agent network topology or big-data parallel computing framework. In this dissertation, we consider a multi-agent optimization formulation where the global objective function is the summation of individual local objective functions with respect to local agents' decision variables of different dimensions, and the constraints include both local private constraints and shared coupling constraints. Employing and extending tools from the monotone operator theory (including resolvent operator, operator splitting, etc.) and fixed point iteration of nonexpansive, averaged operators, a series of distributed solution approaches are proposed, which are all iterative algorithms that rely on parallel agent level local updates and inter-agent coordination. Some of the algorithms require synchronizations across all agents for information exchange during each iteration while others allow asynchrony and delays. The algorithms' convergence to an optimal solution if one exists are established by first characterizing them as fixed point iterations of certain averaged operators under certain carefully designed norms, then showing that the fixed point sets of these averaged operators are exactly the optimal solution set of the original multi-agent optimization problem. The effectiveness and performances of the proposed algorithms are demonstrated and compared through several numerical examples.<br>