A direct method for the numerical solution of optimization problems with time-periodic PDE constraints

In der vorliegenden Dissertation entwickeln wir auf der Basis der Direkten Mehrzielmethode eine neue numerische Methode für Optimalsteuerungsprobleme (OCPs) mit zeitperiodischen partiellen Differentialgleichungen (PDEs). Die vorgeschlagene Methode zeichnet sich durch asymptotisch optimale Skalierung des numerischen Aufwandes in der Zahl der örtlichen Diskretisierungspunkte aus. Sie besteht aus einem Linearen Iterativen Splitting Ansatz (LISA) innerhalb einer Newton-Typ Iteration zusammen mit einer Globalisierungsstrategie, die auf natürlichen Niveaufunktionen basiert. Wir untersuchen die LISA-Newton Methode im Rahmen von Bocks kappa-Theorie und entwickeln zuverlässige a-posteriori kappa-Schätzer. Im Folgenden erweitern wir die LISA-Newton Methode auf den Fall von inexakter Sequentieller Quadratischer Programmierung (SQP) für ungleichungsbeschränke Probleme und untersuchen das lokale Konvergenzverhalten. Zusätzlich entwickeln wir klassische und Zweigitter Newton-Picard Vorkonditionierer für LISA und beweisen gitterunabhängige Konvergenz der klassischen Variante auf einem Modellproblem. Anhand numerischer Ergebnisse können wir belegen, dass im Vergleich zur klassichen Variante die Zweigittervariante sogar noch effizienter ist für typische Anwendungsprobleme. Des Weiteren entwickeln wir eine Zweigitterapproximation der Lagrange-Hessematrix, welche gut in den Rahmen des Zweigitter Newton-Picard Ansatzes passt und die im Vergleich zur exakten Hessematrix zu einer Laufzeitreduktion von 68% auf einem nichtlinearen Benchmarkproblem führt. Wir zeigen weiterhin, dass die Qualität des Feingitters die Genauigkeit der Lösung bestimmt, während die Qualität des Grobgitters die asymptotische lineare Konvergenzrate, d.h., das Bocksche kappa, festlegt. Zuverlässige kappa-Schätzer ermöglichen die automatische Steuerung der Grobgitterverfeinerung für schnelle Konvergenz. Für die Lösung der auftretenden, großen Probleme der Quadratischen Programmierung (QPs) wählen wir einen strukturausnutzenden zweistufigen Ansatz. In der ersten Stufe nutzen wir die durch den Mehrzielansatz und die Newton-Picard Vorkonditionierer bedingten Strukturen aus, um die großen QPs auf äquivalente QPs zu reduzieren, deren Größe von der Zahl der örtlichen Diskretisierungspunkte unabhängig ist. Für die zweite Stufe entwickeln wir Erweiterungen für eine Parametrische Aktive Mengen Methode (PASM), die zu einem zuverlässigen und effizienten Löser für die resultierenden, möglicherweise nichtkonvexen QPs führen. Weiterhin konstruieren wir drei anschauliche, contra-intuitive Probleme, die aufzeigen, dass die Konvergenz einer one-shot one-step Optimierungsmethode weder notwendig noch hinreichend für die Konvergenz der entsprechenden Methode für das Vorwärtsproblem ist. Unsere Analyse von drei Regularisierungsansätzen zeigt, dass de-facto Verlust von Konvergenz selbst mit diesen Ansätzen nicht verhindert werden kann. Des Weiteren haben wir die vorgestellten Methoden in einem Computercode mit Namen MUSCOP implementiert, der automatische Ableitungserzeugung erster und zweiter Ordnung von Modellfunktionen und Lösungen der dynamischen Systeme, Parallelisierung auf der Mehrzielstruktur und ein Hybrid Language Programming Paradigma zur Verfügung stellt, um die benötigte Zeit für das Aufstellen und Lösen neuer Anwendungsprobleme zu minimieren. Wir demonstrieren die Anwendbarkeit, Zuverlässigkeit und Effektivität von MUSCOP und damit der vorgeschlagenen numerischen Methoden anhand einer Reihe von PDE OCPs von steigender Schwierigkeit, angefangen bei linearen akademischen Problemen über hochgradig nichtlineare akademische Probleme der mathematischen Biologie bis hin zu einem hochgradig nichtlinearen Anwendungsproblem der chemischen Verfahrenstechnik im Bereich der präparativen Chromatographie auf Basis realer Daten: Dem Simulated Moving Bed (SMB) Prozess

Semilinear geometric optics with boundary amplification

We study weakly stable semilinear hyperbolic boundary value problems with highly oscillatory data. Here weak stability means that exponentially growing modes are absent, but the so-called uniform Lopatinskii condition fails at some boundary frequency $\beta$ in the hyperbolic region. As a consequence of this degeneracy there is an amplification phenomenon: outgoing waves of amplitude O(\eps^2) and wavelength \eps give rise to reflected waves of amplitude O(\eps), so the overall solution has amplitude O(\eps). Moreover, the reflecting waves emanate from a radiating wave that propagates in the boundary along a characteristic of the Lopatinskii determinant. An approximate solution that displays the qualitative behavior just described is constructed by solving suitable profile equations that exhibit a loss of derivatives, so we solve the profile equations by a Nash-Moser iteration. The exact solution is constructed by solving an associated singular problem involving singular derivatives of the form \partial_{x'}+\beta\frac{\partial_{\theta_0}}{\eps}, $x'$ being the tangential variables with respect to the boundary. Tame estimates for the linearization of that problem are proved using a first-order calculus of singular pseudodifferential operators constructed in the companion article \cite{CGW2}. These estimates exhibit a loss of one singular derivative and force us to construct the exact solution by a separate Nash-Moser iteration. The same estimates are used in the error analysis, which shows that the exact and approximate solutions are close in $L^\infty$ on a fixed time interval independent of the (small) wavelength \eps. The approach using singular systems allows us to avoid constructing high order expansions and making small divisor assumptions

Numerical methods for systems of highly oscillatory ordinary differential equations

Current research made contribution to the numerical analysis of highly oscillatory ordinary differential equations. Highly oscillatory functions appear to be at the forefront of the research in numerical analysis. In this work we developed efficient numerical algorithms for solving highly oscillatory differential equations. The main important achievements are: to the contrary of classical methods, our numerical methods share the feature that asymptotically the approximation to the exact solution improves as the frequency of oscillation grows; also our methods are computationally feasible and as such do not require fine partition of the integration interval. In this work we show that our methods introduce better accuracy of approximation as compared with the state of the art solvers in Matlab and Maple.This thesis presents methods for efficient numerical approximation of linear and non-linear systems of highly oscillatory ordinary differential equations. Phenomena of high oscillation is considered a major computational problem occurring in Fourier analysis, computational harmonic analysis, quantum mechanics, electrodynamics and fluid dynamics. Classical methods based on Gaussian quadrature fail to approximate oscillatory integrals. In this work we introduce numerical methods which share the remarkable feature that the accuracy of approximation improves as the frequency of oscillation increases. Asymptotically, our methods depend on inverse powers of the frequency of oscillation, turning the major computational problem into an advantage. Evolving ideas from the stationary phase method, we first apply the asymptotic method to solve highly oscillatory linear systems of differential equations. The asymptotic method provides a background for our next, the Filon-type method, which is highly accurate and requires computation of moments. We also introduce two novel methods. The first method, we call it the FM method, is a combination of Magnus approach and the Filon-type method, to solve matrix exponential. The second method, we call it the WRF method, a combination of the Filon-type method and the waveform relaxation methods, for solving highly oscillatory non-linear systems. Finally, completing the theory, we show that the Filon-type method can be replaced by a less accurate but moment free Levin-type method.The work was supported by Trinity College, University of Cambridge

Efficient method for detection of periodic orbits in chaotic maps and flows

An algorithm for detecting unstable periodic orbits in chaotic systems [Phys. Rev. E, 60 (1999), pp. 6172-6175] which combines the set of stabilising transformations proposed by Schmelcher and Diakonos [Phys. Rev. Lett., 78 (1997), pp. 4733-4736] with a modified semi-implicit Euler iterative scheme and seeding with periodic orbits of neighbouring periods, has been shown to be highly efficient when applied to low-dimensional system. The difficulty in applying the algorithm to higher dimensional systems is mainly due to the fact that the number of stabilising transformations grows extremely fast with increasing system dimension. In this thesis, we construct stabilising transformations based on the knowledge of the stability matrices of already detected periodic orbits (used as seeds). The advantage of our approach is in a substantial reduction of the number of transformations, which increases the efficiency of the detection algorithm, especially in the case of high-dimensional systems. The performance of the new approach is illustrated by its application to the four-dimensional kicked double rotor map, a six-dimensional system of three coupled H\'enon maps and to the Kuramoto-Sivashinsky system in the weakly turbulent regime.Comment: PhD thesis, 119 pages. Due to restrictions on the size of files uploaded, some of the figures are of rather poor quality. If necessary a quality copy may be obtained (approximately 1MB in pdf) by emailing me at [email protected]

Hybridizable compatible finite element discretizations for numerical weather prediction: implementation and analysis

There is a current explosion of interest in new numerical methods for atmospheric modeling. A driving force behind this is the need to be able to simulate, with high efficiency, large-scale geophysical flows on increasingly more parallel computer systems. Many current operational models, including that of the UK Met Office, depend on orthogonal meshes, such as the latitude-longitude grid. This facilitates the development of finite difference discretizations with favorable numerical properties. However, such methods suffer from the ``pole problem," which prohibits the model to make efficient use of a large number of computing processors due to excessive concentration of grid-points at the poles. Recently developed finite element discretizations, known as ``compatible" finite elements, avoid this issue while maintaining the key numerical properties essential for accurate geophysical simulations. Moreover, these properties can be obtained on arbitrary, non-orthogonal meshes. However, the efficient solution of the resulting discrete systems depend on transforming the mixed velocity-pressure (or velocity-pressure-buoyancy) system into an elliptic problem for the pressure. This is not so straightforward within the compatible finite element framework due to inter-element coupling. This thesis supports the proposition that systems arising from compatible finite element discretizations can be solved efficiently using a technique known as ``hybridization." Hybridization removes inter-element coupling while maintaining the desired numerical properties. This permits the construction of sparse, elliptic problems, for which fast solver algorithms are known, using localized algebra. We first introduce the technique for compatible finite element discretizations of simplified atmospheric models. We then develop a general software abstraction for the rapid implementation and composition of hybridization methods, with an emphasis on preconditioning. Finally, we extend the technique for a new compatible method for the full, compressible atmospheric equations used in operational models.Open Acces

Semilinear geometric optics with boundary amplification

International audienceWe study weakly stable semilinear hyperbolic boundary value problems with highly oscillatory data. Here weak stability means that exponentially growing modes are absent, but the so-called uniform Lopatinskii condition fails at some boundary frequency $\beta$ in the hyperbolic region. As a consequence of this degeneracy there is an amplification phenomenon: outgoing waves of amplitude $O(\varepsilon^2)$ and wavelength $\varepsilon$ give rise to reflected waves of amplitude $O(\varepsilon)$, so the overall solution has amplitude $O(\varepsilon)$. Moreover, the reflecting waves emanate from a radiating wave that propagates in the boundary along a characteristic of the Lopatinskii determinant. An approximate solution that displays the qualitative behavior just described is constructed by solving suitable profile equations that exhibit a loss of derivatives, so we solve the profile equations by a Nash-Moser iteration. The exact solution is constructed by solving an associated singular problem involving singular derivatives of the form $\partial_{x'}+\beta\frac{\partial_{\theta_0}}{\varepsilon}$, $x'$ being the tangential variables with respect to the boundary. Tame estimates for the linearization of that problem are proved using a first-order calculus of singular pseudodifferential operators constructed in the companion article \cite{CGW2}. These estimates exhibit a loss of one singular derivative and force us to construct the exact solution by a separate Nash-Moser iteration. The same estimates are used in the error analysis, which shows that the exact and approximate solutions are close in $L^\infty$ on a fixed time interval independent of the (small) wavelength $\varepsilon$. The approach using singular systems allows us to avoid constructing high order expansions and making small divisor assumptions

SciCADE 95: International conference on scientific computation and differential equations

This report consists of abstracts from the conference. Topics include algorithms, computer codes, and numerical solutions for differential equations. Linear and nonlinear as well as boundary-value and initial-value problems are covered. Various applications of these problems are also included

Modified Chebyshev-Picard Iteration Methods for Solution of Initial Value and Boundary Value Problems

The solution of initial value problems (IVPs) provides the evolution of dynamic system state history for given initial conditions. Solving boundary value problems (BVPs) requires finding the system behavior where elements of the states are defined at different times. This dissertation presents a unified framework that applies modified Chebyshev-Picard iteration (MCPI) methods for solving both IVPs and BVPs. Existing methods for solving IVPs and BVPs have not been very successful in exploiting parallel computation architectures. One important reason is that most of the integration methods implemented on parallel machines are only modified versions of forward integration approaches, which are typically poorly suited for parallel computation. The proposed MCPI methods are inherently parallel algorithms. Using Chebyshev polynomials, it is straightforward to distribute the computation of force functions and polynomial coefficients to different processors. Combining Chebyshev polynomials with Picard iteration, MCPI methods iteratively refine estimates of the solutions until the iteration converges. The developed vector-matrix form makes MCPI methods computationally efficient. The power of MCPI methods for solving IVPs is illustrated through a small perturbation from the sinusoid motion problem and satellite motion propagation problems. Compared with a Runge-Kutta 4-5 forward integration method implemented in MATLAB, MCPI methods generate solutions with better accuracy as well as orders of magnitude speedups, prior to parallel implementation. Modifying the algorithm to do double integration for second order systems, and using orthogonal polynomials to approximate position states lead to additional speedups. Finally, introducing perturbation motions relative to a reference motion results in further speedups. The advantages of using MCPI methods to solve BVPs are demonstrated by addressing the classical Lambert’s problem and an optimal trajectory design problem. MCPI methods generate solutions that satisfy both dynamic equation constraints and boundary conditions with high accuracy. Although the convergence of MCPI methods in solving BVPs is not guaranteed, using the proposed nonlinear transformations, linearization approach, or correction control methods enlarge the convergence domain. Parallel realization of MCPI methods is implemented using a graphics card that provides a parallel computation architecture. The benefit from the parallel implementation is demonstrated using several example problems. Larger speedups are achieved when either force functions become more complicated or higher order polynomials are used to approximate the solutions