Linear programming and the calculation of maximum norm approximations

This thesis examines methods for the solution of a number of approximation problems which include multivariate approximation and the approximate solution of differential equations. After a brief discussion of the literature in Chapter 1, the best linear Chebyshev approximation prob l em is discussed in Chapter 2. A brief introduc ti on to the genera l theory is given, which includes a number of characterization theorems . The di screte linear Chebys hev approximation problem is then discussed and it is shown how this may be solved by linear programm ing techniques. In Chapter 3, the problem of determining, numerically, best linear Chebyshev approximations is considered, and two algorithms are given. Chapter 4 generalizes the methods of Chapters 2 and 3 to enable the numerical calculation of approximate solutions of differential equations to be made. A method is developed which can be used when the approximating family does not satisfy the boundary conditions. In §4 .2 linear elliptic partial differential equations are considered in detail, and an error bound is derived. Use can be made of this to produce an approximation with the best error bound. Non-linear problems are considered in Chapters 5 and 6. In §5 .3, two methods are presented for non-linear continuous approximation which I are analogues of the methods of Chapter 3. The methods depend on an algorithm for the solution of a discrete problem. Two such algorithms are presented in §5.l, §5.2. In Chapter 6, an investigation is made into the question of improving the efficiency of these algorithms. An algorithm is presented which is similar to those of Chapter 5, and which depends on a parameter, k. Conditions are given under which the algorithm has an order of convergence of at least k+l . It is then shown how k may be chosen in an attempt to maximize efficiency. Because many of the algorithms in this thesis require the solution of discrete Chebyshev problems, a set of linear programming routines is included in the appendix.· The code is written to take advantage of the special structure of these problems

Refresher course in maths and a project on numerical modeling done in twos

These lecture notes accompany a refresher course in applied mathematics with a focus on numerical concepts (Part I), numerical linear algebra (Part II), numerical analysis, Fourier series and Fourier transforms (Part III), and differential equations (Part IV). Several numerical projects for group work are provided in Part V. In these projects, the tasks are threefold: mathematical modeling, algorithmic design, and implementation. Therein, it is important to draw interpretations of the obtained results and provide measures (Parts I-IV) how to build confidence into numerical findings such intuition, error analysis, convergence analysis, and comparison to manufactured solutions. Both authors have been jointly teaching over several years this class and bring in a unique mixture of their respective teaching and research fields

Refresher course in maths and a project on numerical modeling done in twos

These lecture notes accompany a refresher course in applied mathematics with a focus on numerical concepts (Part I), numerical linear algebra (Part II), numerical analysis, Fourier series and Fourier transforms (Part III), and differential equations (Part IV). Several numerical projects for group work are provided in Part V. In these projects, the tasks are threefold: mathematical modeling, algorithmic design, and implementation. Therein, it is important to draw interpretations of the obtained results and provide measures (Parts I-IV) how to build confidence into numerical findings such intuition, error analysis, convergence analysis, and comparison to manufactured solutions. Both authors have been jointly teaching over several years this class and bring in a unique mixture of their respective teaching and research fields

Numerical Methods for Partial Differential Equations

These lecture notes are devoted to the numerical solution of partial differential equations (PDEs). PDEs arise in many fields and are extremely important in modeling of technical processes with applications in physics, biology, chemisty, economics, mechanical engineering, and so forth. In these notes, not only classical topics for linear PDEs such as finite differences, finite elements, error estimation, and numerical solution schemes are addressed, but also schemes for nonlinear PDEs and coupled problems up to current state-of-the-art techniques are covered. In the Winter 2020/2021 an International Class with additional funding from DAAD (German Academic Exchange Service) and local funding from the Leibniz University Hannover, has led to additional online materials such as links to youtube videos, which complement these lecture notes. This is the updated and extended Version 2. The first version was published under the DOI: https://doi.org/10.15488/9248

Advanced Treatment of Fission Yield Effects and Method Development for Improved Reactor Depletion Calculations

Fission product yield data play an important role in simulations of nuclear fission reactors, aimed at fuel cycle and safety analyses. The respective evaluated data libraries still have shortcomings regarding the treatment of energy dependencies and uncertainty information. This work has been aimed at the development of a fission model for future fission product yield evaluations as well as its validation on the levels of cross-sections, fission product yields and time dependent decay radiation

Advanced Treatment of Fission Yield Effects and Method Development for Improved Reactor Depletion Calculations

Daten über Spaltproduktausbeuten spielen eine wichtige Rolle in Simulationen von Kernreaktoren, welche auf Brennstoffzyklus- und Sicherheitsanalysen abzielen. Zusammen mit den Daten über Wirkungsquerschnitte und radioaktiven Zerfall bestimmen sie das Nuklidinventar, die neutronischen Eigenschaften und die Entstehung von Nachzerfallswärmeleistung von bestrahltem Kernbrennstoff. In bestehenden evaluierten Kerndatenbibliotheken sind die Spaltproduktausbeuten in einer Struktur von nicht mehr als drei oder vier Energiegruppen bzw. Datenpunkten gegeben. Die zugehörige Unsicherheitsinformation ist bislang auf Varianzen beschränkt, was die Quantifizierung der Unsicherheit beispielsweise der Kritikalität bei hohem Abbrand einschränkt. Die vorstehende Arbeit zielte darauf ab, einen Modellcode für Kernspaltung zur Anwendung in künftigen Evaluationen der Spaltproduktausbeuten zu entwickeln, mit dem Ziel einer besseren Berücksichtigung der den Kernspaltungsprozess beherrschenden physikalischen Effekte. Der GEF-Code, welcher von K.-H. Schmidt und B. Jurado im Auftrag der OECD NEA entwickelt wurde, basiert auf einem detaillierten und physikalisch fundierten Kernspaltungsmodell. Um eine mögliche Abregung des Compoundkerns vor der Spaltung adäquat zu modellieren, wurde dessen Version GEF-2013/2.2 mit dem Modellcode TALYS-1.4 gekoppelt. Die Modellrechnung der Spaltproduktausbeuten wurde in 77 Energiegruppen bis hinauf zu 20 MeV ausgeführt, was eine detaillierte Analyse der Energieabhängigkeiten ermöglicht. Validierungen der Ergebnisse wurden auf verschiedenen Ebenen vorgenommen, darunter primäre und kumulative Spaltproduktausbeuten, die Ausbeute verzögerter Neutronen und die zeitabhängige Emission von verzögerter Neutronen-, Beta- und Gammastrahlung. Das zeitabhängige Nuklidinventar wurde mithilfe der rationalen Tschebyscheff-Näherung berechnet. In vielen Fällen ergab sich eine zufriedenstellende Übereinstimmung mit den experimentellen Daten. Die Vorhersage der verzögerten Neutronenemission erwies sich jedoch als größte Herausforderung. Defizite des Kernspaltungsmodells, welche in diesem Bereich noch größere Abweichungen verursachten, wurden identifiziert. Weiterhin wurde die Anwendung des Modells zur Beschreibung der mit der Kernspaltung zusammenhängenden Observablen im aufgelösten Resonanzbereich diskutiert und für das Target U-235 demonstriert. Schließlich wurden die von den gekoppelten Codes erzeugten Spaltproduktausbeuten in der Simulation eines während der 1970er Jahre im DWR bei Obrigheim (Deutschland) durchgeführten Abbrandexperiments angewandt. Diese Simulation wurde in einer gemeinsamen Arbeit mit dem hauseigenen KAPROS-KANEXT Codesystem ausgeführt. Auf der Ebene des Abbrandcodes wurde die Berücksichtigung der Energieabhängigkeiten der Spaltproduktausbeuten durch den Einbau neu entwickelter Module, welche der Verarbeitung der Daten über Spaltproduktausbeuten dienen, verbessert. Die unter Anwendung der modellbasierten Daten erhaltenen Ergebnisse stimmen gut mit den meisten im Abbrandexperiment gemessenen Werten überein und erwiesen sich als konkurrenzfähig zu den Ergebnissen aus bestehenden Standardverfahren

Gluing Spaces and Analysis

In the first part of this work we study the gluing of metric measure spaces. The gluing is defined by a bilipschitz map which identifies the gluing sets. The new metric is given by the minimal length of all possible paths on the glued space. On each of these spaces a strongly local, regular Dirichlet form is defined. Additionally, each space satisfies a doubling property and a strong scaling invariant Poincaré inequality for all balls holds true. We derive the doubling property and the scaling invariant Poincaré inequality on the glued space provided a lower bound on the "heat transmission coefficient" for certain sets holds true. For the proof only assumptions on the Dirichlet forms on the separate pieces are used. These results imply upper and lower Gaussian estimates on the heat kernel, short-time asymptotics and the Feller property of the associated process on the glued space. In the second part we give some generalizations of results by Charles Amick. These are characterizations of the validity of the Poincaré inequality and of Rellichs compact embedding theorem on a domain in terms of a quantity extracted from the boundary. We prove characterizations of this kind for strongly local regular Dirichlet forms on metric measure spaces which satisfy a scaling invariant Poincaré inequality

Improved Numerical Methods for Elliptic Problems in Astrophysics

Les ePDEs (elliptic partial differential equations, en anglès) apareixen en una àmplia varietat d'àrees de les matemàtiques, la física i l'enginyeria. Són de particular interès en Astrofísica on apareixen, per exemple, quan es calcula el potencial gravitacional, en la solució de l'equació de Grad-Shafranov per magnetosferes lliures de forces, o d'imposar lligadures de divergència zero en la integració numèrica de les equacions MHD (magnetohydrodynamics, en anglès). En general, les ePDEs s'han de resoldre numèricament, establint una demanda cada vegada més gran d'algoritmes eficients i altament paral·lels per abordar la seua resolució computacional. El SRJ (scheduled relaxation Jacobi, en anglès) pertany a una prometedora classe de mètodes, atípic per la combinació de senzillesa i eficàcia, que s'ha introduït recentment per resoldre ePDEs lineals de tipus Poisson. És una extensió del mètode iteratiu clàssic de Jacobi utilitzat per resoldre sistemes d'equacions lineals del tipus Au = b. Hereta, d'entre altres, la seua robustesa. La seua metodologia es basa en el càlcul d'uns paràmetres apropiats per a una aproximació multinivell amb l'objectiu de minimitzar el nombre d'iteracions necessàries per a reduir el residual per davall d'una tolerància especificada. L'eficiència en la reducció del residual augmenta amb el nombre de nivells emprats en l'algoritme. Tanmateix, l'aplicació de la metodologia original per calcular els paràmetres d'estos esquemes SRJ òptims més enllà de 5 nivells és enormement dificultosa. Això és degut fonamentalment a la presència d'un sistema mixt algebraic-diferencial (no lineal) d'equacions el qual es torna cada vegada més rígid a mesura que augmenta el nombre de nivells. D'una banda, hem trobat una nova metodologia per a l'obtenció dels paràmetres dels esquemes òptims de l'algoritme SRJ que supera les limitacions de la metodologia original i proporciona els paràmetres per a estos esquemes amb un nombre elevat de nivells, fóra bo fins a 15, i per a resolucions de fins a 215 punts per dimensió. Això dóna lloc a factors d'acceleració de diversos centenars respecte del mètode de Jacobi en el cas de resolucions típiques i de milers en alguns casos amb altes resolucions. La major part de l'èxit en la recerca d'estos esquemes òptims amb més de 10 nivells es basa en una reducció analítica de la complexitat del sistema d'equacions abans esmentat. A més, s'estén l'algoritme original per aplicar-lo a certs sistemes d'equacions el·líptiques no lineals. D'altra banda, en un esquema típic SRJ, s'empra l'anterior conjunt de paràmetres en cicles de M iteracions consecutives fins que s'arriba a la tolerància prescrita. Presentem la forma analítica del conjunt òptim de factors de relaxació per al cas en què tots ells són estrictament diferents, i veiem que l'algoritme resultant és equivalent al mètode no estacionari de Richardson generalitzat, en el que es precondiciona la matriu del sistema d'equacions multiplicant per D = diag(A). El nostre mètode per estimar els pesos té l'avantatge que el càlcul explícit dels valors propis mínim i màxim de la matriu A (o la matriu d'iteració corresponent de l'esquema de Jacobi amb pes subjacent) es substitueix pel càlcul (molt més fàcil) de les freqüències mínima i màxima derivades de l'anàlisi d'estabilitat de von Neumann de l'operador el·líptic continu. Este conjunt de pesos també és l'òptim per al problema general, la qual cosa ens dóna la convergència més ràpida de tots els possibles esquemes SRJ per una estructura de malla donada. Ens referirem a ell com el mètode de Chebyshev-Jacobi. El factor d'amplificació del mètode es pot trobar analíticament i permet l'estimació exacta del nombre d'iteracions necessàries per a assolir la tolerància desitjada. També mostrem que a partir del conjunt de pesos calculats per l'esquema SRJ òptim per a una mida de cicle fix és possible calcular numèricament el valor òptim del factor de relaxació del mètode SOR (successive overrelaxation, en anglès) en alguns casos. Demostrem amb exemples pràctics, d'aplicació en Astrofísica, que el nostre mètode també funciona molt bé per als problemes de tipus Poisson en els que es fa servir una discretització d'alt ordre de l'operador Laplacià (per exemple, discretitzacions de 9- o 17- punts). Això té molt d'interès, ja que estes discretitzacions no produeixen matrius CO (consistently ordered, en anglès) i, per tant, la teoria de Young no es pot utilitzar per calcular el valor òptim del paràmetre de relaxació òptim de SOR. D'altra banda, els esquemes SRJ òptims deduïts ací són avantatjoses respecte a les implementacions existents per SOR pel que fa a discretitzacions d'alt ordre de l'operador Laplacià en la mesura que no cal recórrer als esquemes multicolors per a la seua execució en paral·lel. Presentem el mètode de Chebyshev-Jacobi fent servir una implementació purament MPI i una implementació híbrida OpenMP/MPI, ambdues sobre màquines de memòria compartida i de memòria distribuïda. Mostrem el seu rendiment i com escalen. També mostrem com arribar a velocitats de convergència notables amb execucions en paral·lel sobre GPUs quan la resolució d'equacions en derivades parcials el·líptiques amb diferències finites es fa utilitzant de manera conjunta el mètode de Chebyshev-Jacobi i les discretitzacions d'alt ordre. Finalment, tractar d'aplicar els nostres mètodes més enllà de l'àmbit de l'Astrofísica. En particular, abordem el problema de trobar els modes normals de vibració de l'ull humà. Este problema es pot resoldre amb una variant millorada de la metodologia que ací es presenta. La millora consisteix a estendre el càlcul del conjunt òptim de paràmetres al cas de matrius no definides positives. Les nostres idees sobre com procedir en este camp s'esbossen en el treball futur d'esta tesi.Elliptic partial differential equations appear in a wide variety of areas of mathematics, physics and engineering. They are of particular interest in Astrophysics and appear, e.g., when computing the gravitational potential, in the solution of the Grad-Shafranov equation for force-free magnetospheres, to impose divergence free constraints in the numerical integration of MHD equations or when solving the constraint equations in General Relativity. Typically, elliptic equations must be solved numerically, which sets an ever-growing demand for efficient and highly parallel algorithms to tackle their computational solution. The Scheduled Relaxation Jacobi is a promising class of methods, atypical for combining simplicity and efficiency, that has been recently introduced for solving linear Poisson-like elliptic equations. It is an extension of the classical Jacobi iterative method to solve linear systems of equations (Au=b). It also inherits its robustness. Its methodology relies on computing the appropriate parameters of a multilevel approach with the goal of minimizing the number of iterations needed to cut down the residuals below specified tolerances. The efficiency in the reduction of the residual increases with the number of levels employed in the algorithm. Applying the original methodology to compute the algorithm parameters with more than 5 levels notably hinders obtaining optimal schemes, as the mixed (non-linear) algebraic-differential system of equations from which they result become notably stiff. On one hand, we present a new methodology for obtaining the parameters of Scheduled Relaxation Jacobi schemes that overcomes the limitations of the original algorithm and provides parameters for these schemes with up to 15 levels and resolutions of up to 215 points per dimension, allowing for acceleration factors larger than several hundreds with respect to the Jacobi method for typical resolutions and, in some high resolution cases, close to 1000. Most of the success in finding these optimal schemes with more than 10 levels is based on an analytic reduction of the complexity of the previously mentioned system of equations. Furthermore, we extend the original algorithm to apply it to certain systems of non-linear elliptic equations. On the other hand, in a typical Scheduled Relaxation Jacobi scheme, the former set of factors is employed in cycles of M consecutive iterations until a prescribed tolerance is reached. We present the analytic form for the optimal set of relaxation factors for the case in which all of them are strictly different, and find that the resulting algorithm is equivalent to a non-stationary generalized Richardson's method where the matrix of the system of equations is preconditioned multiplying it by D=diag(A). Our method to estimate the weights has the advantage that the explicit computation of the maximum and minimum eigenvalues of the matrix A (or the corresponding iteration matrix of the underlying weighted Jacobi scheme) is replaced by the (much easier) calculation of the maximum and minimum frequencies derived from a von Neumann analysis of the continuous elliptic operator. This set of weights is also the optimal one for the general problem, resulting in the fastest convergence of all possible Scheduled Relaxation Jacobi schemes for a given grid structure. We refer to it as the Chebyshev-Jacobi method. The amplification factor of the method can be found analytically and allows for the exact estimation of the number of iterations needed to achieve a desired tolerance. We also show that with the set of weights computed for the optimal SRJ scheme for a fixed cycle size it is possible to estimate numerically the optimal value of the relaxation factor in the successive overrelaxation method in some cases. We demonstrate with practical examples, with application in Astrophysics, that our method also works very well for Poisson-like problems in which a high-order discretization of the Laplacian operator is employed (e.g., a 9- or 17-points discretization). This is of interest since the former discretizations do not yield consistently ordered A matrices and, hence, the theory of Young cannot be used to predict the optimal value of the SOR parameter. Furthermore, the optimal SRJ schemes deduced here are advantageous over existing SOR implementations for high-order discretizations of the Laplacian operator in as much as they do not need to resort to multi-coloring schemes for their parallel implementation. We present the implementation of the Chebyshev-Jacobi method using a purely MPI implementation, an openMP / MPI hybrid implementation over both shared memory machines and distributed memory machines. They show ideal speed up. We also show how to reach a remarkable speed up when solving elliptic partial differential equations with finite differences thanks to the joint use of the Chebyshev-Jacobi method with high order discretizations and its parallel implementation over GPUs. Finally, we have tried to apply our methods beyond the realm of Astrophysics with limited success though. Specially, we have addressed the problem of finding normal modes of human eyeballs. This problem is ready for being solved with an improved variant of the methodology here presented. The improvement consists on extending the calculation of the optimal set of parameters to non positive-definite matrices. Our ideas on how to proceed in this field are sketched in the outlook of this thesis