Numerical methods for boundary value problems with applications to bioengineering

Abstract

Boundary value problems (BVPs) modeling biomedical systems involve linearly or nonlinearly coupled systems of PDEs in quite non-trivial geometries. These BVPs can be broadly classified into exterior: modeling for instance applications of ultrasound techniques, and interior: modeling biological flows, tumor growth, source localization, and so on. The design and implementation of robust and efficient numerical methods for both types of BVPs in realistic three-dimensional geometries is the purpose of this thesis. Targeting applications in sonography, the direct scattering problem of acoustic waves by a prolate spheroidal scatterer was investigated by introducing a novel analytical approach accompanied by an extensive numerical study. Analytical methods in scattering theory lead to extremely ill-conditioned matrices in the linear systems involved. Thus we adopted arbitrary precision arithmetic in our computer implementation. This allowed us to handle previously unsolvable linear systems, and paved the way for the first extensive numerical study. In the first chapter we suggest and analyze a new method for obtaining the expansion coefficients of the scattered field. We demonstrate that our approach can successfully handle elongated scatterers and high frequency regimes. The modeling of biological flows and transport of biological molecules involves convection-diffusion-reaction PDE systems. The numerical treatment of such problems must overcome two main obstacles. The first obstacle is the inability of classical stabilization approaches to remove numerical ripples completely, especially in areas where the solution exhibits steep gradients. These shortcomings motivated for our second chapter the investigation of flux-corrected transport techniques for finite element methods (FEM-FCT). The latter are fully algebraic parameter-independent schemes that show superlinear convergence and (unlike classical stabilization methods) preserve positivity. A new semi-implicit FEM-FCT scheme is introduced, which is demonstrated to be significantly faster than the original algorithm, especially when combined with inexact Newton methods. The next major challenge in large-scale biological applications is the efficient solution of the linear systems arising from three-dimensional discretization of the underlying PDEs. This is the most demanding part of a computer simulation, as it consumes more than 95% of the total running time and is still considered an open problem. Motivated by the widespread availability of multi-processing computing systems, in our third chapter we consider domain decomposition approaches, which are the most attractive methods for parallel implementations. There, concepts from Functional Analysis, including Interpolation of Sobolev Spaces, are merged with tools and algorithms of Linear Algebra, such as the Lanczos algorithm for symmetric generalized eigenvalue problems. They target the sparse, efficient, and yet scalable evaluation of the action of a non-sparse representation of fractional Sobolev norms. The outcome is convergence of the solution approach independent of the problem size, and in many cases independent of the number of partitions. This prediction is tested on a range of discretizations of elliptic problems in both two and three dimensions. The performance of the resulting iterative solvers surpasses existing methodologies. Finally, the scalability of the parallel implementation of our preconditioning approach is extensively benchmarked on SMP multiprocessing platforms under several different operating systems.