Linear Solvers for Singularly Perturbed Problems Numerical Analysis for Singularly Perturbed Problems (dedicated to the 60 th birthday of Martin Stynes) 16 th – 18 th Nov, 2011 Niall Madden, NUI Galway Joint work with Scott MacLachlan, Tufts University, MA. 6 3 1.5 0.8 0.4 4 2 0.6 0.3 1 2 1 0.4 0.2 0.5 0 0 0.2 0.1 −2 −1 0 0 0 15 15 15 15 15 10 15 10 15 10 15 10 15 10 15 10 10 10 10 10 5 5 5 5 5 5 5 5 5 5 0 0 0 0 0 0 0 0 0 0 0.2 0.1 0.06 0.03 0.15 0.04 0.02 0.1 0.05 0.02 0.01 0.05 0 0 0 0 15 15 15 15 15 15 15 15 10 10 10 10 10 10 10 10 5 5 5 5 5 5 5 5 0 0 0 0 0 0 0 0 −3 −3 x 10 x 10 0.015 8 4 6 3 0.01 4 2 0.005 2 1 0 0 0 15 15 15 10 15 10 15 10 15 10 10 10 5 5 5 5 5 5 0 0 0 0 0 0 Numerical Analysis for Singularly Perturbed Problems, 16–18 Nov 2011 1/25
The abstract The abstract for this talk was: There has been great interest in the design and analysis of numerical schemes that yield robust solutions to singularly perturbed problems. However, there have been only a small number of studies concerning the solution to the associated linear systems. In this talk we discuss some difficulties that arise when applying standard direct solvers to finite difference discretizations of reaction-diffusion equations. Furthermore, we will present some recent developments of multigrid approaches for these problems. Particular emphasis will be given to coupled systems of reaction-diffusion problems... That there are been relatively few studies on linear solvers for SPPs, does not suggest there have been none, e.g.., [Roos, 1996, Ansari and Hegarty, 2003] and by the Zaragoza group. The focus here is on reaction-diffusion problems in one and two dimensions. In particular addressing some issues that have an effect on direct solvers. Systems are less of an issue for this talk than was planned when the abstract was written. Numerical Analysis for Singularly Perturbed Problems, 16–18 Nov 2011 2/25
Overview of the talk This work began (in earnest) Oct-Dec of last year at the Institute for Mathematics and Its Applications , the University of Minnesota, Minneapolis. It is joint work with Scott MacLachlan , Tufts University. Our collaboration has been supported by the NSF under grant DMS-0811022. Numerical Analysis for Singularly Perturbed Problems, 16–18 Nov 2011 3/25
Overview of the talk In this talk I will describe some work on solvers for discretizations of singularly perturbed reaction-diffusion equations. Particularly: Some observations on the use of direct solvers; Towards a multigrid preconditioner. Moreover... 1. I’ll explain why it is Martin Stynes’ fault that I’m interested in this problem. 2. We reconsider some difficulties associated with applying a direct solver to a standard discretization, and offer an explanation for its poor performance. 3. Demonstrate that an Algebraic Multigrid (AMG) solver can be applied with some success. 4. Describe some progress in the development of a Geometric Multigrid (GMG) solver for the problem. Numerical Analysis for Singularly Perturbed Problems, 16–18 Nov 2011 4/25
Background My interest in these differential equations stems from my time as a graduate student in University College Cork, working under Martin’s supervision, as part of a project led by Gareth Thomas. The main thesis topic was supposed to be the development of methods for problems in wave-current interactions. The equations we had to solve included a k - ǫ turbulence model : a coupled system of two nonlinear reaction-diffusion problems with leading coefficients with different orders of magnitude. Numerical Analysis for Singularly Perturbed Problems, 16–18 Nov 2011 5/25
Background My interest in these differential equations stems from my time as a graduate student in University College Cork, working under Martin’s supervision, as part of a project led by Gareth Thomas. The main thesis topic was supposed to be the development of methods for problems in wave-current interactions. The equations we had to solve included a k - ǫ turbulence model : a coupled system of two nonlinear reaction-diffusion problems with leading coefficients with different orders of magnitude. Earlier numerical analyses of fitted meshes for coupled systems, particularly [Shishkin, 1995] and [Matthews et al., 2002] led us apply piecewise uniform meshes to this applied problem. An analysis of the approach for the linear case was given in [Madden and Stynes, 2003]. Numerical Analysis for Singularly Perturbed Problems, 16–18 Nov 2011 5/25
Background Some years later we returned to the topic of analysing fitted mesh methods coupled systems. In [Clavero et al., 2005] an analysis is given for a Shishkin mesh for an (uncoupled) reaction-diffusion problem in two dimensions. This was extended to a system of M coupled equations in [Kellogg et al., 2008]. ∆ 0 . . . 0 0 . . . 0 ∆ − ε 2 . . u + A u = f . ... . . . . 0 0 0 . . . ∆ When generating the numerical results for that paper, we noticed that the time required by the linear solver (“backslash” in Matlab) did not appear to be uniform in the perturbation parameter. Numerical Analysis for Singularly Perturbed Problems, 16–18 Nov 2011 6/25
Background More recently, I was fortunate to be involved in a project with Martin on the development of a two-scale sparse grid finite element method for (uncoupled and coupled) reaction-diffusion problems in two dimensions [Liu et al., 2009]. The algorithm reduces the number of unknowns from O ( N 2 ) to O ( N 3 / 2 ) without compromising the accuracy of the underlying scheme. The price paid for this includes increased bandwidth in the linear system, and an increase in the number of non-zero entries. Justifying the assertion that the two-scale approach is more efficient than the standard Galerkin approach required a systematic study of the CPU time taken by the linear solver. Galerkin Mesh Sparse Grid Mesh(es) Numerical Analysis for Singularly Perturbed Problems, 16–18 Nov 2011 7/25
Background More recently, I was fortunate to be involved in a project with Martin on the development of a two-scale sparse grid finite element method for (uncoupled and coupled) reaction-diffusion problems in two dimensions [Liu et al., 2009]. The algorithm reduces the number of unknowns from O ( N 2 ) to O ( N 3 / 2 ) without compromising the accuracy of the underlying scheme. The price paid for this includes increased bandwidth in the linear system, and an increase in the number of non-zero entries. Justifying the assertion that the two-scale approach is more efficient than the standard Galerkin approach required a systematic study of the CPU time taken by the linear solver. Galerkin sparsity pattern Sparse Grid pattern Numerical Analysis for Singularly Perturbed Problems, 16–18 Nov 2011 7/25
Background More recently, I was fortunate to be involved in a project with Martin on the development of a two-scale sparse grid finite element method for (uncoupled and coupled) reaction-diffusion problems in two dimensions [Liu et al., 2009]. The algorithm reduces the number of unknowns from O ( N 2 ) to O ( N 3 / 2 ) without compromising the accuracy of the underlying scheme. The price paid for this includes increased bandwidth in the linear system, and an increase in the number of non-zero entries. Justifying the assertion that the two-scale approach is more efficient than the standard Galerkin approach required a systematic study of the CPU time taken by the linear solver. Again it was noticed that the standard Matlab solver did not behave in a uniform manner for our test problem... Numerical Analysis for Singularly Perturbed Problems, 16–18 Nov 2011 7/25
Background Galerkin sparsity pattern Sparse Grid pattern Numerical Analysis for Singularly Perturbed Problems, 16–18 Nov 2011 8/25
Example A two-dimensional model problem − ε 2 ∆u + bu = f on ( 0, 1 ) × ( 0, 1 ) + BCs. 3 2 1 0 1 1 0.5 0.5 0 0 The solution exhibits boundary and interior layers, in this example resolved using a piecewise uniform mesh. Numerical Analysis for Singularly Perturbed Problems, 16–18 Nov 2011 9/25
Example We’ll solve this problem using a finite difference scheme on a tensor-product grid. Let. ¯ x and ¯ Ω N Ω N y be arbitrary meshes with N intervals on [ 0, 1 ] . Set Ω N = { ( x i , y j ) } N i , j = 0 to be the tensor product of ¯ x and ¯ Ω N Ω N y . Set h i = x i − x i − 1 , k j = y j − y j − 1 , ¯ ¯ and h i = ( x i + 1 − x i − 1 ) / 2, k j = ( y j + 1 − y j − 1 ) / 2. Define the scaled 5-point second-order central difference operator: ¯ h i k j + 1 ¯ ¯ � 1 � 1 k j 1 1 k j ∆ N := − ¯ − ¯ � � k j + h i + . h i h i h i + 1 k j k j + 1 h i + 1 ¯ h i k j Then (neglecting boundary conditions), the numerical scheme is: − ε 2 ∆ N + ¯ h ¯ U = ¯ h ¯ � � kF . kB The system matrix is symmetric positive definite, and the system of equations should be easy to solve. Numerical Analysis for Singularly Perturbed Problems, 16–18 Nov 2011 10/25
Recommend
More recommend
