MIMS EPrints
Not a member yet
2151 research outputs found
Sort by
Testing matrix function algorithms using identities
No full textAlgorithms for computing matrix functions are typically tested by comparing the forward error with the product of the condition number and the unit roundoff. The forward error is computed with the aid of a reference solution, typically computed at high precision. An alternative approach is to use functional identities such as the ``round trip tests'' and , as are currently employed in a SciPy test module.
We show how a linearized perturbation analysis for a functional identity allows the determination of a maximum residual consistent with backward stability of the constituent matrix function evaluations. Comparison of this maximum residual with a computed residual provides a necessary test for backward stability. We also show how the actual linearized backward error for these relations can be computed. Our approach makes use of Fr\'echet derivatives and estimates of their norms. Numerical experiments show that the proposed approaches are able both to detect instability and to confirm stability
Backward error analysis of the shift-and-invert Arnoldi algorithm
Get PDFWe perform a backward error analysis of the inexact shift-and-invert Arnoldi algorithm.
We consider inexactness in the solution of the arising linear systems, as well as in the orthonormalization steps, and
take the non-orthonormality of the computed Krylov basis into account.
We show that the computed basis and Hessenberg matrix satisfy an exact shift-and-invert Krylov relation for a perturbed matrix, and we give bounds for the perturbation.
We show that the shift-and-invert Arnoldi algorithm is backward stable if the condition number of the small Hessenberg matrix is not too large.
This condition is then relaxed using implicit restarts.
Moreover we give notes on the Hermitian case, considering Hermitian backward errors, and
finally, we use our analysis to derive a sensible breakdown condition
Generalized rational Krylov decompositions with an application to rational approximation
Get PDFGeneralized rational Krylov decompositions are matrix relations which, under certain conditions, are associated with rational Krylov spaces. We study the algebraic properties of such decompositions and present an implicit Q theorem for rational Krylov spaces. Transformations on rational Krylov decompositions allow for changing the poles of a rational Krylov space without recomputation, and two algorithms are presented for this task. Using such transformations we develop a rational Krylov method for rational least squares fitting. Numerical experiments indicate that the proposed method converges fast and robustly. A MATLAB toolbox with implementations of the presented algorithms and experiments is provided
Estimating the Condition Number of f(A)b
Get PDFNew algorithms are developed for estimating the condition number of , where is a matrix and is a vector. The condition number estimation algorithms for already available in the literature require the explicit computation of matrix functions and their Fr\'{e}chet derivatives and are therefore unsuitable for the large, sparse typically encountered in problems. The algorithms we propose here use only
matrix-vector multiplications. They are based on a modified version of the power iteration for estimating the norm of the Fr\'{e}chet derivative of a matrix function, and work in conjunction with any existing algorithm for computing . The number of matrix-vector multiplications required to estimate the condition number is proportional to the square of the number of matrix-vector multiplications required by the underlying
algorithm. We develop a specific version of our algorithm for estimating the
condition number of , based on the algorithm of Al-Mohy and Higham [SIAM J. Matrix Anal. Appl., 30(4):1639--1657, 2009].
Numerical experiments demonstrate that our condition estimates are reliable and of reasonable cost
Linear and weakly nonlinear instability of a premixed curved flame under the influence of its spontaneous acoustic field.
Get PDFThe stability of premixed flames in a duct is investigated using an asymptotic formulation, which is derived from first principles and based on high-activation-energy and low-Mach-number assumptions (Wu et al., J. Fluid Mech., vol. 497, 2003, pp. 23�53). The present approach takes into account the dynamic coupling between the flame and its spontaneous acoustic field, as well as the interactions between the hydrodynamic field and the flame. The focus is on the fundamental mechanisms of combustion instability. To this end, a linear stability analysis of some steady curved flames is undertaken. These steady flames are known to be stable when the spontaneous acoustic perturbations are ignored. However, we demonstrate that they are actually unstable when the latter effect is included. In order to corroborate this result, and also to provide a relatively simple model guiding active control, we derived an extended Michelson�Sivashinsky equation, which governs the linear and weakly nonlinear evolution of a perturbed flame under the influence of its spontaneous sound. Numerical solutions to the initial-value problem confirm the linear instability result, and show how the flame evolves nonlinearly with time. They also indicate that in certain parameter regimes the spontaneous sound can induce a strong secondary subharmonic parametric instability. This behaviour is explained and justified mathematically by resorting to Floquet theory. Finally we compare our theoretical results with experimental observations, showing that our model captures some of the observed behaviour of propagating flames
EIT Reconstruction Algorithms for Respiratory Intensive Care
Get PDFElectrical impedance tomography (EIT) is an emerging medical imaging technique that aims to reconstruct the internal conductivity distribution of a subject from electrical measurements obtained on the skin. In this thesis we explore the promising application of EIT to the respiratory monitoring of humans. We pay particular focus to the forward problem, highlighting the need to have an
accurately known external boundary shape and electrode positions on a reconstruction model. A theoretical study of uniqueness results of EIT with an unknown external
boundary shape is presented. A novel sensitivity study of the external boundary shape is presented as well as results from a reconstruction algorithm to account for errors in electrode position with simulated data in 3D. We also demonstrate results of a shape correction algorithm from a pilot study of lung EIT with data collected using the fEITER system, and MR images used to inform the external boundary shape of healthy subjects. After image co-registration of the resulting dynamic 3D EIT reconstruction images with the lung-segmented MR image, we outline a novel mutual information performance criterion to measure the quality of reconstructed images. We also outline the computation of the forward problem of the complete electrode model in 3D using high order polynomial finite elements and present convergence results in 2D for the continuum, point and complete electrode model. Our numerical study demonstrates that the convergence rate of the forward problem is independent of the polynomial approximation order for the complete electrode model and there is no global convergence for the point electrode model in the energy norm.
Reconstructed conductivity images can be difficult to interpret at the bedside. Moreover clinicians would like clinically meaningful indices, such as regional lung compliance, to determine the pathologies of patients in real time. By modelling the respiratory system as a coupled time dependent system of simple mechanical functional units, we propose a novel methodology to couple mechanical ventilation and EIT. The mechanical properties of the lungs are estimated through an inverse coefficient problem on coupled ODEs, with the measurable data being the time series of pressure at airway opening and interior air volume data. We present results with simulated data as well as a discussion on extensions and limitations to the mechanical models.
Finally we present a theoretical discussion of anisotropic EIT. It is well known that any diffeomorphism fixing points on the boundary gives rise to a conductivity with the same electrical measurements on the skin, generating a large class of conductivities that are electrically equivalent. We define novel classes of anisotropic media with constraints on their eigenspace: prescribed eigenvalues, prescribed orthogonal coordinates,
prescribed eigenvectors, fibrous and layered conductivities. By drawing analogies with elasticity theory, we discuss how these constraints on the eigenspace restrict the set of diffeomorphisms fixing points on the boundary, and present two uniqueness results for anisotropic conductivities with prescribed eigenvalues and prescribed eigenvectors
Restoring Definiteness via Shrinking, with an Application to Correlation Matrices with a Fixed Block
Get PDFIndefinite approximations of positive semidefinite matrices arise in many data analysis applications involving covariance matrices and correlation matrices. We propose a method for restoring positive semidefiniteness of an indefinite matrix that constructs a convex linear combination of and a positive semidefinite target matrix . In statistics, this construction for improving an estimate by combining it with new information in is known as shrinking. We make no statistical assumptions about and define the optimal shrinking parameter as \alpha_* = \min \{\alpha \in [0,1] : \mbox{S(\alpha) is positive semidefinite}\}. We describe three \alg s for computing . One algorithm is based on the bisection method, with the use of Cholesky factorization to test definiteness, a second employs Newton's method, and a third finds the smallest eigenvalue of a symmetric definite generalized eigenvalue problem. We show that weights that reflect confidence in the individual entries of can be used to construct a natural choice of the target matrix . We treat in detail a problem variant in which a positive semidefinite leading principal submatrix of remains fixed, showing how the fixed block can be exploited to reduce the cost of the bisection and generalized eigenvalue methods. Numerical experiments show that when applied to indefinite approximations of correlation matrices shrinking can be at least an order of magnitude faster than computing the nearest correlation matrix
A mathematical model of the colon crypt capturing compositional dynamic interactions between cell types
Get PDFModels of the development and early progression of colorectal cancer are based upon understanding the cycle of stem cell turnover, proliferation, differentiation and death. Existing crypt compartmental models feature a linear pathway of cell types, with little regulatory mechanism. Previous work has shown that there are perturbations in the enteroendocrine cell population of macroscopically normal crypts, a compartment not included in existing models. We show that existing models do not adequately recapitulate the dynamics of cell fate pathways in the crypt. We report the progressive development, iterative testing and fitting of a developed compartmental model with additional cell types, and which includes feedback mechanisms and cross-regulatory mechanisms between cell types. The fitting of the model to existing data sets suggests a need to invoke cross-talk between cell types as a feature of colon crypt cycle models
A new strain energy function for the hyperelastic modelling of ligaments and tendons based on fascicle microstructure
Get PDFA new strain energy function for the hyperelastic modelling of ligaments and tendons based on the geometrical arrangement of their fibrils is derived. The distribution of the crimp angles of the fibrils is used to determine the stress-strain response of a single fascicle, and this stress-strain response is used to determine the form of the strain energy function, the parameters of which can all potentially be directly measured via experiments - unlike those of commonly used strain energy functions such as the Holzapfel-Gasser-Ogden (HGO) model, whose parameters are phenomenological. We compare the new model with the HGO model and show that the new model gives a better match to existing stress-strain data for human patellar tendon than the HGO model, with the average relative error in matching this data when using the new model being 0.053 (compared with 0.57 when using the HGO model), and the average absolute error when using the new model being 0.12MPa (compared with 0.31MPa when using the HGO model)
Zolotarev quadrature rules and load balancing for the FEAST eigensolver
Get PDFThe FEAST method for solving large sparse eigenproblems is equivalent to subspace iteration with an approximate spectral projector and implicit orthogonalization. This relation allows to characterize the convergence of this method in terms of the error of a certain rational approximant to an indicator function. We propose improved rational approximants leading to FEAST variants with faster convergence, in particular, when using rational approximants based on the work of Zolotarev. Numerical experiments demonstrate the possible computational savings especially for pencils whose eigenvalues are not well separated and when the dimension of the search space is only slightly larger than the number of wanted eigenvalues. The new approach improves both convergence robustness and load balancing when FEAST runs on multiple search intervals in parallel