This thesis studies the so-called “criticality problem”, an important generalised eigenvalue problem arising in neutron transport theory. The smallest positive real eigenvalue of the problem contains valuable information about the status of the fission chain reaction in the nuclear reactor (i.e. the criticality of the reactor), and thus plays an important role in the design and safety of nuclear power stations. Because of the practical importance, efficient numerical methods to solve the criticality problem are needed, and these are the focus of this thesis. In the theory we consider the time-independent neutron transport equation in the monoenergetic homogeneous case with isotropic scattering and vacuum boundary conditions. This is an unsymmetric integro-differential equation in 5 independent variables, modelling transport, scattering, and fission, where the dependent variable is the neutron angular flux. We show that, before discretisation, the nonsymmetric eigenproblem for the angular flux is equivalent to a related eigenproblem for the scalar flux, involving a symmetric positive definite weakly singular integral operator(in space only). Furthermore, we prove the existence of a simple smallest positive real eigenvalue with a corresponding eigenfunction that is strictly positive in the interior of the reactor. We discuss approaches to discretise the problem and present discretisations that preserve the underlying symmetry in the finite dimensional form. The thesis then describes methods for computing the criticality in nuclear reactors, i.e. the smallest positive real eigenvalue, which are applicable for quite general geometries and physics. In engineering practice the criticality problem is often solved iteratively, using some variant of the inverse power method. Because of the high dimension, matrix representations for the operators are often not available and the inner solves needed for the eigenvalue iteration are implemented by matrix-free inneriterations. This leads to inexact iterative methods for criticality computations, for which there appears to be no rigorous convergence theory. The fact that, under appropriate assumptions, the integro-differential eigenvalue problem possesses an underlying symmetry (in a space of reduced dimension) allows us to perform a systematic convergence analysis for inexact inverse iteration and related methods. In particular, this theory provides rather precise criteria on how accurate the inner solves need to be in order for the whole iterative method to converge. The theory is illustrated with numerical examples on several test problems of physical relevance, using GMRES as the inner solver. We also illustrate the use of Monte Carlo methods for the solution of neutron transport source problems as well as for the criticality problem. Links between the steps in the Monte Carlo process and the underlying mathematics are emphasised and numerical examples are given. Finally, we introduce an iterative scheme (the so-called “method of perturbation”) that is based on computing the difference between the solution of the problem of interest and the known solution of a base problem. This situation is very common in the design stages for nuclear reactors when different materials are tested, or the material properties change due to the burn-up of fissile material. We explore the relation ofthe method of perturbation to some variants of inverse iteration, which allows us to give convergence results for the method of perturbation. The theory shows that the method is guaranteed to converge if the perturbations are not too large and the inner problems are solved with sufficiently small tolerances. This helps to explain the divergence of the method of perturbation in some situations which we give numerical examples of. We also identify situations, and present examples, in which the method of perturbation achieves the same convergence rate as standard shifted inverse iteration. Throughout the thesis further numerical results are provided to support the theory.

Very fine discretizations of differential operators often lead to large, sparse matrices A, where the condition number of A is large. Such ill-conditioning has well known effects on both solving linear systems and eigenvalue computations, and, in general, computing solutions with relative accuracy independent of the condition number is highly desirable. This dissertation is divided into two parts. In the first part, we discuss a method of preconditioning, developed by Ye, which allows solutions of Ax=b to be computed accurately. This, in turn, allows for accurate eigenvalue computations. We then use this method to develop discretizations that yield accurate computations of the smallest eigenvalue of the biharmonic operator across several domains. Numerical results from the various schemes are provided to demonstrate the performance of the methods. In the second part we address the role of the condition number of A in the context of multigrid algorithms. Under various assumptions, we use rigorous Fourier analysis on 2- and 3-grid iteration operators to analyze round off errors in floating point arithmetic. For better understanding of general results, we provide detailed bounds for a particular algorithm applied to the 1-dimensional Poisson equation. Numerical results are provided and compared with those obtained by the schemes discussed in part 1.

The solution of large and sparse eigenvalue problems arising from the discretization of the diffusion equation is considered. The multigroup diffusion equation is discretized by means of the Nodal expansion Method (NEM) [9, 10]. A new formulation of the higher order NEM variants revealing the true nature of the problem, that is, a generalized eigenvalue problem, is proposed. These generalized eigenvalue problems are solved using the Jacobi-Davidson (JD) method [26]. The most expensive part of the method consists of solving a linear system referred to as correction equation. It is solved using Krylov subspace methods in combination with aggregation-based Algebraic Multigrid (AMG) techniques. In that context, a particular aggregation technique used in combination with classical smoothers, referred to as oblique geometric coarsening, has been derived. Its particularity is that it aggregates unknowns that are not coupled, which has never been done to our knowledge. A modular code, combining JD with an AMG preconditioner, has been developed. The code comes with many options, that have been tested. In particular, the instability of the Rayleigh-Ritz [33] acceleration procedure in the non-symmetric case has been underlined. Our code has also been compared to an industrial code extracted from ARTEMIS.

Directional wave-field decomposition for heterogeneousanisotropic media with in-stantaneous response is establishedfor both the acoustic and the electromagnetic equations. We derive a sufficient condition for ellipticity of thesystem's matrix in the Laplace domain and show that theconstruction of the splitting matrix via a Dunford-Taylorintegral over the resolvent of the non-compact, non-normalsystem's matrix is well de ned. The splitting matrix also hasproperties that make it possible to construct the decompositionwith a generalized eigenvector procedure. The classical way ofobtaining the decomposition is equivalent to solving analgebraic Riccati operator equation. Hence the proceduredescribed above also provides a solution to the algebraicRiccati operator equation. The solution to the wave-field decomposition for theisotropic wave equation is expressed in terms of theDirichlet-to-Neumann map for a plane. The equivalence of thisDirichlet-to-Neumann map is the acoustic admittance, i.e. themapping between the pressure and the particle velocity. Theacoustic admittance, as well as the related impedance aresolutions to algebraic Riccati operator equations and are keyelements in the decomposition. In the electromagnetic case thecorresponding impedance and admittance mappings solve therespective algebraic Riccati operator equations and henceprovide solutions to the decomposition problem. The present research shows that it is advantageous toutilize the freedom implied by the generalized eigenvectorprocedure to obtain the solution to the decomposition problemin more general terms than the admittance/impedancemappings. The time-reversal approach to steer an acoustic wave eld inthe cavity and half space geometries are analyzed from aboundary control perspective. For the cavity it is shown thatwe can steer the field to a desired final configuration, withthe assumption of local energy decay. It is also shown that thetime-reversal algorithm minimizes a least square error forfinite times when the data are obtained by measurements. Forthe half space geometry, the boundary condition is expressedwith help of the wave-field decomposition. In the homogeneousmaterial case, the response of the time-reversal algorithm iscalculated analytically. This procedure uses the one-wayequations together with the decomposition operator.

<p>Directional wave-field decomposition for heterogeneousanisotropic media with in-stantaneous response is establishedfor both the acoustic and the electromagnetic equations.</p><p>We derive a sufficient condition for ellipticity of thesystem's matrix in the Laplace domain and show that theconstruction of the splitting matrix via a Dunford-Taylorintegral over the resolvent of the non-compact, non-normalsystem's matrix is well de ned. The splitting matrix also hasproperties that make it possible to construct the decompositionwith a generalized eigenvector procedure. The classical way ofobtaining the decomposition is equivalent to solving analgebraic Riccati operator equation. Hence the proceduredescribed above also provides a solution to the algebraicRiccati operator equation.</p><p>The solution to the wave-field decomposition for theisotropic wave equation is expressed in terms of theDirichlet-to-Neumann map for a plane. The equivalence of thisDirichlet-to-Neumann map is the acoustic admittance, i.e. themapping between the pressure and the particle velocity. Theacoustic admittance, as well as the related impedance aresolutions to algebraic Riccati operator equations and are keyelements in the decomposition. In the electromagnetic case thecorresponding impedance and admittance mappings solve therespective algebraic Riccati operator equations and henceprovide solutions to the decomposition problem.</p><p>The present research shows that it is advantageous toutilize the freedom implied by the generalized eigenvectorprocedure to obtain the solution to the decomposition problemin more general terms than the admittance/impedancemappings.</p><p>The time-reversal approach to steer an acoustic wave eld inthe cavity and half space geometries are analyzed from aboundary control perspective. For the cavity it is shown thatwe can steer the field to a desired final configuration, withthe assumption of local energy decay. It is also shown that thetime-reversal algorithm minimizes a least square error forfinite times when the data are obtained by measurements. Forthe half space geometry, the boundary condition is expressedwith help of the wave-field decomposition. In the homogeneousmaterial case, the response of the time-reversal algorithm iscalculated analytically. This procedure uses the one-wayequations together with the decomposition operator.</p>

The presence of structure in a computational problem can often be exploited and can lead to a more efficient numerical algorithm. In this dissertation, we look at structured numerical problems that arise from applications in wireless communications and machine learning that also impact other areas of scientific computing. In wireless communication system designs, certain structured matrices (frames) need to be generated. The design of such matrices is equivalent to a symmetric inverse eigenvalue problem where the values of the diagonal elements are prescribed. We present algorithms that are capable of generating a larger set of these constructions than previous algorithms. We also discuss the existence of equiangular tight frames---frames that satisfy additional structural properties. Kernel learning is an important class of problems in machine learning. It often relies on efficient numerical algorithms that solve underlying convex optimization problems. In our work, the objective functions to be minimized are the von Neumann and the LogDet Bregman matrix divergences. The algorithm that solves this optimization problem performs matrix updates based on repeated eigendecompositions of diagonal plus rank-one matrices in the case of von Neumann matrix divergence, and Cholesky updates in case of the LogDet Bregman matrix divergence. Our contribution exploits the low-rank representations and the structure of the constraint matrices, resulting in more efficient algorithms than previously known. We also present two specialized zero-finding algorithms where we exploit the structure through the shape and exact formulation of the objective function. The first zero-finding task arises during the matrix update step which is part of the above-mentioned kernel learning application. The second zero-finding problem is for the secular equation; it is equivalent to the computation of the eigenvalues of a diagonal plus rank-one matrix. The secular equation arises in various applications, the most well-known is the divide-and-conquer eigensolver. In our solutions, we build upon a somewhat forgotten zero-finding method by P. Jarratt, first described in 1966. The method employs first derivatives only and needs the same amount of evaluations as Newton's method, but converges faster. Our contributions are the more efficient specialized zero-finding algorithms. / text