Title: Numerical modeling of free oscillations applied to superconducting-gravimeter data in a low-frequency seismic range Author: Eliška Zábranová Department: Department of Geophysics Supervisor: Doc. RNDr. Ctirad Matyska, DrSc. Abstract: Deformations and changes of the gravitational potential of prestressed selfgravitating elastic bodies caused by free oscillations are described by means of the momentum and Poisson equations and the constitutive relation. For spheri- cally symmetric bodies we transform the equations and boundary conditions into ordinary differential equations of the second order by the spherical harmonic de- composition and further discretize the equations by highly accurate pseudospectral difference schemes on Chebyshev grids. We thus receive a series of matrix eigenvalue problems for eigenfrequencies and eigenfunctions of the free oscillations. Since elas- tic parameters are frequency dependent, we solve the problem for several fiducial frequencies and interpolate the results. Both the mode frequencies and the eigen- functions are benchmarked against the output from the Mineos software package based on Runge-Kutta integration techniques. Subsequently, we use our method to calculate low-frequency synthetic accelerograms of the recent megathrust events and compare them with the observed...

Optical testing in adverse environments, ophthalmology and applications where characterization by curvature is leveraged all have a common goal: accurately estimate wavefront shape. This dissertation investigates wavefront sensing techniques as applied to optical testing based on gradient and curvature measurements. Wavefront sensing involves the ability to accurately estimate shape over any aperture geometry, which requires establishing a sampling grid and estimation scheme, quantifying estimation errors caused by measurement noise propagation, and designing an instrument with sufficient accuracy and sensitivity for the application. Starting with gradient-based wavefront sensing, a zonal least-squares wavefront estimation algorithm for any irregular pupil shape and size is presented, for which the normal matrix equation sets share a pre-defined matrix. A Gerchberg–Saxton iterative method is employed to reduce the deviation errors in the estimated wavefront caused by the pre-defined matrix across discontinuous boundary. The results show that the RMS deviation error of the estimated wavefront from the original wavefront can be less than λ/130~ λ/150 (for λ equals 632.8nm) after about twelve iterations and less than λ/100 after as few as four iterations. The presented approach to handling irregular pupil shapes applies equally well to wavefront estimation from curvature data. A defining characteristic for a wavefront estimation algorithm is its error propagation behavior. The error propagation coefficient can be formulated as a function of the eigenvalues of the wavefront estimation-related matrices, and such functions are established for each of the basic estimation geometries (i.e. Fried, Hudgin and Southwell) with a serial numbering scheme, where a square sampling grid array is sequentially indexed row by row. The results show that with the wavefront piston-value fixed, the odd-number grid sizes yield lower error propagation than the even-number grid sizes for all geometries. The Fried geometry either allows sub-sized wavefront estimations within the testing domain or yields a two-rank deficient estimation matrix over the full aperture; but the latter usually suffers from high error propagation and the waffle mode problem. Hudgin geometry offers an error propagator between those of the Southwell and the Fried geometries. For both wavefront gradient-based and wavefront difference-based estimations, the Southwell geometry is shown to offer the lowest error propagation with the minimum-norm least-squares solution. Noll’s theoretical result, which was extensively used as a reference in the previous literature for error propagation estimate, corresponds to the Southwell geometry with an odd-number grid size. For curvature-based wavefront sensing, a concept for a differential Shack-Hartmann (DSH) curvature sensor is proposed. This curvature sensor is derived from the basic Shack-Hartmann sensor with the collimated beam split into three output channels, along each of which a lenslet array is located. Three Hartmann grid arrays are generated by three lenslet arrays. Two of the lenslets shear in two perpendicular directions relative to the third one. By quantitatively comparing the Shack-Hartmann grid coordinates of the three channels, the differentials of the wavefront slope at each Shack-Hartmann grid point can be obtained, so the Laplacian curvatures and twist terms will be available. The acquisition of the twist terms using a Hartmann-based sensor allows us to uniquely determine the principal curvatures and directions more accurately than prior methods. Measurement of local curvatures as opposed to slopes is unique because curvature is intrinsic to the wavefront under test, and it is an absolute as opposed to a relative measurement. A zonal least-squares-based wavefront estimation algorithm was developed to estimate the wavefront shape from the Laplacian curvature data, and validated. An implementation of the DSH curvature sensor is proposed and an experimental system for this implementation was initiated. The DSH curvature sensor shares the important features of both the Shack-Hartmann slope sensor and Roddier’s curvature sensor. It is a two-dimensional parallel curvature sensor. Because it is a curvature sensor, it provides absolute measurements which are thus insensitive to vibrations, tip/tilts, and whole body movements. Because it is a two-dimensional sensor, it does not suffer from other sources of errors, such as scanning noise. Combined with sufficient sampling and a zonal wavefront estimation algorithm, both low and mid frequencies of the wavefront may be recovered. Notice that the DSH curvature sensor operates at the pupil of the system under test, therefore the difficulty associated with operation close to the caustic zone is avoided. Finally, the DSH-curvature-sensor-based wavefront estimation does not suffer from the 2π-ambiguity problem, so potentially both small and large aberrations may be measured.

This thesis is on the numerical computation of eigenvalues of symmetric hierarchical matrices. The numerical algorithms used for this computation are derivations of the LR Cholesky algorithm, the preconditioned inverse iteration, and a bisection method based on LDLT factorizations. The investigation of QR decompositions for H-matrices leads to a new QR decomposition. It has some properties that are superior to the existing ones, which is shown by experiments using the HQR decompositions to build a QR (eigenvalue) algorithm for H-matrices does not progress to a more efficient algorithm than the LR Cholesky algorithm. The implementation of the LR Cholesky algorithm for hierarchical matrices together with deflation and shift strategies yields an algorithm that require O(n) iterations to find all eigenvalues. Unfortunately, the local ranks of the iterates show a strong growth in the first steps. These H-fill-ins makes the computation expensive, so that O(n³) flops and O(n²) storage are required. Theorem 4.3.1 explains this behavior and shows that the LR Cholesky algorithm is efficient for the simple structured Hl-matrices. There is an exact LDLT factorization for Hl-matrices and an approximate LDLT factorization for H-matrices in linear-polylogarithmic complexity. This factorizations can be used to compute the inertia of an H-matrix. With the knowledge of the inertia for arbitrary shifts, one can compute an eigenvalue by bisectioning. The slicing the spectrum algorithm can compute all eigenvalues of an Hl-matrix in linear-polylogarithmic complexity. A single eigenvalue can be computed in O(k²n log^4 n). Since the LDLT factorization for general H-matrices is only approximative, the accuracy of the LDLT slicing algorithm is limited. The local ranks of the LDLT factorization for indefinite matrices are generally unknown, so that there is no statement on the complexity of the algorithm besides the numerical results in Table 5.7. The preconditioned inverse iteration computes the smallest eigenvalue and the corresponding eigenvector. This method is efficient, since the number of iterations is independent of the matrix dimension. If other eigenvalues than the smallest are searched, then preconditioned inverse iteration can not be simply applied to the shifted matrix, since positive definiteness is necessary. The squared and shifted matrix (M-mu I)² is positive definite. Inner eigenvalues can be computed by the combination of folded spectrum method and PINVIT. Numerical experiments show that the approximate inversion of (M-mu I)² is more expensive than the approximate inversion of M, so that the computation of the inner eigenvalues is more expensive. We compare the different eigenvalue algorithms. The preconditioned inverse iteration for hierarchical matrices is better than the LDLT slicing algorithm for the computation of the smallest eigenvalues, especially if the inverse is already available. The computation of inner eigenvalues with the folded spectrum method and preconditioned inverse iteration is more expensive. The LDLT slicing algorithm is competitive to H-PINVIT for the computation of inner eigenvalues. In the case of large, sparse matrices, specially tailored algorithms for sparse matrices, like the MATLAB function eigs, are more efficient. If one wants to compute all eigenvalues, then the LDLT slicing algorithm seems to be better than the LR Cholesky algorithm. If the matrix is small enough to be handled in dense arithmetic (and is not an Hl(1)-matrix), then dense eigensolvers, like the LAPACK function dsyev, are superior. The H-PINVIT and the LDLT slicing algorithm require only an almost linear amount of storage. They can handle larger matrices than eigenvalue algorithms for dense matrices. For Hl-matrices of local rank 1, the LDLT slicing algorithm and the LR Cholesky algorithm need almost the same time for the computation of all eigenvalues. For large matrices, both algorithms are faster than the dense LAPACK function dsyev.:List of Figures xi List of Tables xiii List of Algorithms xv List of Acronyms xvii List of Symbols xix Publications xxi 1 Introduction 1 1.1 Notation 2 1.2 Structure of this Thesis 3 2 Basics 5 2.1 Linear Algebra and Eigenvalues 6 2.1.1 The Eigenvalue Problem 7 2.1.2 Dense Matrix Algorithms 9 2.2 Integral Operators and Integral Equations 14 2.2.1 Definitions 14 2.2.2 Example - BEM 16 2.3 Introduction to Hierarchical Arithmetic 17 2.3.1 Main Idea 17 2.3.2 Definitions 19 2.3.3 Hierarchical Arithmetic 24 2.3.4 Simple Hierarchical Matrices (Hl-Matrices) 30 2.4 Examples 33 2.4.1 FEM Example 33 2.4.2 BEM Example 36 2.4.3 Randomly Generated Examples 37 2.4.4 Application Based Examples 38 2.4.5 One-Dimensional Integral Equation 38 2.5 Related Matrix Formats 39 2.5.1 H2-Matrices 40 2.5.2 Diagonal plus Semiseparable Matrices 40 2.5.3 Hierarchically Semiseparable Matrices 42 2.6 Review of Existing Eigenvalue Algorithms 44 2.6.1 Projection Method 44 2.6.2 Divide-and-Conquer for Hl(1)-Matrices 45 2.6.3 Transforming Hierarchical into Semiseparable Matrices 46 2.7 Compute Cluster Otto 47 3 QR Decomposition of Hierarchical Matrices 49 3.1 Introduction 49 3.2 Review of Known QR Decompositions for H-Matrices 50 3.2.1 Lintner’s H-QR Decomposition 50 3.2.2 Bebendorf’s H-QR Decomposition 52 3.3 A new Method for Computing the H-QR Decomposition 54 3.3.1 Leaf Block-Column 54 3.3.2 Non-Leaf Block Column 56 3.3.3 Complexity 57 3.3.4 Orthogonality 60 3.3.5 Comparison to QR Decompositions for Sparse Matrices 61 3.4 Numerical Results 62 3.4.1 Lintner’s H-QR decomposition 62 3.4.2 Bebendorf’s H-QR decomposition 66 3.4.3 The new H-QR decomposition 66 3.5 Conclusions 67 4 QR-like Algorithms for Hierarchical Matrices 69 4.1 Introduction 70 4.1.1 LR Cholesky Algorithm 70 4.1.2 QR Algorithm 70 4.1.3 Complexity 71 4.2 LR Cholesky Algorithm for Hierarchical Matrices 72 4.2.1 Algorithm 72 4.2.2 Shift Strategy 72 4.2.3 Deflation 73 4.2.4 Numerical Results 73 4.3 LR Cholesky Algorithm for Diagonal plus Semiseparable Matrices 75 4.3.1 Theorem 75 4.3.2 Application to Tridiagonal and Band Matrices 79 4.3.3 Application to Matrices with Rank Structure 79 4.3.4 Application to H-Matrices 80 4.3.5 Application to Hl-Matrices 82 4.3.6 Application to H2-Matrices 83 4.4 Numerical Examples 84 4.5 The Unsymmetric Case 84 4.6 Conclusions 88 5 Slicing the Spectrum of Hierarchical Matrices 89 5.1 Introduction 89 5.2 Slicing the Spectrum by LDLT Factorization 91 5.2.1 The Function nu(M − µI) 91 5.2.2 LDLT Factorization of Hl-Matrices 92 5.2.3 Start-Interval [a, b] 96 5.2.4 Complexity 96 5.3 Numerical Results 97 5.4 Possible Extensions 100 5.4.1 LDLT Slicing Algorithm for HSS Matrices 103 5.4.2 LDLT Slicing Algorithm for H-Matrices 103 5.4.3 Parallelization 105 5.4.4 Eigenvectors 107 5.5 Conclusions 107 6 Computing Eigenvalues by Vector Iterations 109 6.1 Power Iteration 109 6.1.1 Power Iteration for Hierarchical Matrices 110 6.1.2 Inverse Iteration 111 6.2 Preconditioned Inverse Iteration for Hierarchical Matrices 111 6.2.1 Preconditioned Inverse Iteration 113 6.2.2 The Approximate Inverse of an H-Matrix 115 6.2.3 The Approximate Cholesky Decomposition of an H-Matrix 116 6.2.4 PINVIT for H-Matrices 117 6.2.5 The Interior of the Spectrum 120 6.2.6 Numerical Results 123 6.2.7 Conclusions 130 7 Comparison of the Algorithms and Numerical Results 133 7.1 Theoretical Comparison 133 7.2 Numerical Comparison 135 8 Conclusions 141 Theses 143 Bibliography 145 Index 153