Stability of Krylov Subspace Spectral Methods James V. Lambers Department of Energy Resources Engineering Stanford University includes joint work with Patrick Guidotti and Knut Sølna, UC Irvine Margot Gerritsen, Stanford University Harrachov ‘07 August 23, 2007
Model Variable-coefficient Problem where We assume � � is positive semi-definite, � � � � � � � � � � � � � � � � � � � � � � and that the coefficients are smooth August 23, 2007 Harrachov '07 2/30
Quick-and-dirty Solution Let { φ φ ν ν } be a set of orthonormal � π -periodic φ φ ν ν functions. Then, an approximate solution is where This works if the φ φ ν ν are nearly eigenfunctions, φ φ ν ν but if not, how can we compute � � φ � � � � � � � � � � � � � � � � � � � � � � � � � � � � � � � � � � � � � � � � � as φ ν φ φ ν �� �� � �� �� � � � � ν ν accurately and efficiently as possible? August 23, 2007 Harrachov '07 3/30
Elements of Functions of Matrices � × × � × × � � � � � � � �� �� � If � � is � � and symmetric, then � � � � � �� �� � is � � � � � � � � � � � � given by a Riemann-Stieltjes integral provided the measure α α α � α � λ λ λ λ ��� ��� which is based on � � ��� ��� the spectral decomposition of � � , is positive and � � increasing This is the case if � � , or if � � is a small � � � � � � � � � perturbation of � � � � August 23, 2007 Harrachov '07 4/30
The � � ≠ ≠ � � Case ≠ ≠ � � � � � � � � � � � � � �� � � �� � �� �� For general � � and � � , the bilinear form � � � � � � � � � � � � can be obtained by writing it as the difference quotient where δ δ δ δ is a small constant. Both forms lead to Riemann-Stieltjes integrals with positive, increasing measures How can we approximate these integrals? August 23, 2007 Harrachov '07 5/30
Gaussian Quadrature • These two integrals can be approximated using Gaussian quadrature rules (G. Golub and G. Meurant, '94) • The nodes and weights are obtained by applying the Lanczos algorithm to � � with � � starting vectors � � and ������������� � , the ������������� � � � ������������� ������������� � � tridiagonal matrix of recursion coefficients • The nodes are the eigenvalues of � � , and the � � weights are the products of the first components of the left and right eigenvectors • We want rules for each Fourier component August 23, 2007 Harrachov '07 6/30
What if � � is a differential operator? � � � � � ω � � ω � ω ω � , then recursion coefficients � � • If � � � � � � � � � � � � ω �� � � �� � �� �� � � α � α α α � , β β β � β � are functions of ω ω ω ω � � � � • Let � � from before, with � � �� � � �� �� �� � � � � � � � � � � � � � � � � � � � � � � � � � � � � � � � � ������ � � ������ ������ ������ The recursion coefficients for a 2-node Gaussian rule are where � � ��� � � ��� � ��� ��� � � � • Similar formulas apply in higher dimensions August 23, 2007 Harrachov '07 7/30
Updating Coefficients for � ≠ ≠ � ≠ ≠ To produce modified recursion coefficients generated by � ��� and � � + δ δ � δ δ � : � � � � � ��� ��� ��� � � � � � � � � � � This is particularly useful if the � � , α α � α α � , β β β � β � are � � � � � � � � � � parameterized families and � � is a fixed vector � � August 23, 2007 Harrachov '07 8/30
Krylov Subspace Spectral Methods � � � � � � � �� � � � � � To compute Fourier components of � � : � � � � � � � � � � � � � � � � � � � � • Apply symmetric Lanczos algorithm to � � with � � � � � ω ω ω � ω � starting vector � � � � � � � • Use fast updating to obtain modified recursion � � � ω ω � � ω ω � � � � � ω ω � � � , � � � � � ω ω � � � � � � � coefficients for starting vectors � � δ δ δ δ � � � � � � � � � • Approximate, by Gaussian quadrature, • Finally, August 23, 2007 Harrachov '07 9/30
Why Do it This Way? • Other, more general Krylov subspace methods (e.g. M. Hochbruck and C. Lubich, '96) use a single Krylov subspace for each � as the starting vector, to � � � � time step, with � � � � � � approximate ���� � � ����� ���� ���� � � � �� �� � �� �� � � � � � � • KSS methods obtain Fourier components from derivatives of frequency-dependent Krylov � � � subspace bases in the direction of � � � � � • Thus, each component receives individual attention, enhancing accuracy and stability August 23, 2007 Harrachov '07 10/30
Properties • They're High-Order Accurate! � � � �� �� � Each Fourier component of � � � �� �� � is � � � � � � � � � � � � � � � � � � � � � � � � � � computed with local accuracy of � � ∆ ∆ � ∆ ∆ � � � � � , � � � � � � � � � � where � � is the number of nodes in the � � Gaussian rule • They're Explicit but Very Stable! If � � is constant and � � is bandlimited, then � � � � � � � � � � � � � � � � � � � � � � � � � � � � for �� ��� , method is unconditionally stable! �� �� ��� �� ��� ��� For �� ��� , solution operator is bounded �� ��� �� �� ��� ��� independently of ∆ ∆ � � and � ∆ ∆ � � � � � August 23, 2007 Harrachov '07 11/30
Demonstrating Stability They do not experience the same difficulties with stiffness as other Krylov subspace methods, or the same weak instability as the unsmoothed Fourier method Contrasting Krylov subspace methods Fourier, KSS methods applied to � � � � �� �� � � �� �� � � applied to parabolic problem on different ����� � ����� � � � � � � � � � , smooth initial data until � � ��� ��� ����� ����� � � � � � � � � � � ��� ��� � � grids August 23, 2007 Harrachov '07 12/30
Properties, cont’d • They’re efficient and scalable! – Performance of MATLAB implementation comparable or superior to that of built-in ODE solvers – Accuracy and efficiency scale to finer grids or higher spatial dimension August 23, 2007 Harrachov '07 13/30
In the Limit: Derivatives of Moments! • Each Fourier component approximates the Gâteaux derivative • � � -node approximate solution has the form � � August 23, 2007 Harrachov '07 14/30
The Splitting Perspective • Derivatives of nodes and weights w.r.t. δ δ δ δ are Fourier components of applications of � � � � pseudodifferential operators applied to � � � � � � � � � � � � � � � � � � � � � � � � • � � -node approximate solution has form � � where and each � � is a constant-coefficient � � � � � � pseudo-differential operator, of the same order as � � , and positive semi-definite � � � � � � � � � � � � � � � � � � � � � � August 23, 2007 Harrachov '07 15/30
Simplest Splittings • The case � � reduces to solving the � � � � � � � � � � averaged-coefficient problem exactly, after applying forward Euler with the residual operator • In the case � � , with � � constant, the � � � � � � � � � � � � � � � � � � � � � � � � operators � � are second-order, and yet the � � � � � � approximate solution operator � � ∆ ∆ � ∆ ∆ � satisfies � � � � � � � � � � � � � � � � where � � is a bounded operator! � � August 23, 2007 Harrachov '07 16/30
Recommend
More recommend