Robusne numeričke metode za nelinearne probleme svojstvenih vrijednosti

Abstract

In this thesis we study numerical methods for solving nonlinear eigenvalue problems of polynomial type, i.e. $P(\lambda)x \equiv (\Sigma_{\ell=0}^k \lambda^{\ell} A_{\ell})x = 0$, where $A_{\ell} \in \mathbb{C}^{n \times n}, \lambda \in \mathbb{C}, 0 \neq x \in \mathbb{C}^n$. In particular, we are interested in the quadratic $(k = 2)$ and the quartic $(k = 4)$ eigenvalue problems. The methods are based on the corresponding linearization – the nonlinear problem is replaced with an equivalent linear problem of the type $(A - \lambda B)y = 0$, of dimension $kn$. We propose several modifications and improvements of the existing methods for both the complete and partial solution; this results in new numerical algorithms that are a substantial improvement over the existing ones. In particular, as an improvement of the state of the art quadeig method of Hammarling, Munro and Tisseur, we develop a scheme to deflate all zero and infinite eigenvalues before calling the QZ algorithm for the linear problem. This provides numerically more robust procedure, which we illustrate by numerical examples. Further, we supplement the parameter scaling (designed to equilibrate the norms of the coefficient matrices) with a two–sided diagonal scaling to nearly equilibrate (in modulus) the nonzero matrix entries. In addition, we analyze the fine details of the rank revealing factorization used in the deflation process. We advocate to use complete pivoting in the QR factorization, and we also propose a LU based approach, which is shown to be competitive, or even better than the one based on the QR factorization. The new method is extended to the quartic problem. For the partial quadratic eigenvalue problem (computing only a part of the spectrum), the iterative Arnoldi–like methods are studied, especially the implicitly restarted two level orthogonal Arnoldi algorithm (TOAR). We propose several improvements of the method. In particular, new shift selection strategy is proposed for the implicit restart for the class of overdamped quadratic eigenvalue problems. Also, we show the benefit of choosing the starting vector for TOAR, based on spectral information of a nearby proportionally damped pencil. Finally, we provide some new ideas for the development of a Krylov–Schur like methods that is capable of using arbitrary polynomial filters in the implicit restarting.Nelinearni problemi svojstvenih vrijednosti se javljaju u mnogim primjenama kako u prirodnim znanostima, tako i u inženjerstvu. Jedna od najpoznatijih klasa nelinearnih svojstvenih problema su polinomni svojstveni problemi. Tako se, na primjer, kvadratični svojstveni problem $(\lambda^2 M + \lambda C + K)x = 0$ pojavljuje u dinamičkoj analizi mehaničkih i električnih struktura, u vibro–akustici, mehanici fluida, obradi signala. S druge strane, polinomni se problem četvrtog reda $(\lambda^4 A + \lambda^3 B + \lambda^2 C + \lambda D + K)x = 0$ pojavljuje u analizi stabilnosti Poiseuilleovog toka u cijevi. Za razliku od linearnih problema svojstvenih vrijednosti, numeričke metode za nelinearne probleme još uvijek nisu dovoljno razrađene, niti numerički pouzdane, iako je algebarska teorija za polinomne probleme svojstvenih vrijednosti dobro razvijena. Naglasak ove disertacije je na numeričkom rješavanju kvadratičnog svojstvenog problema. Cilj je razviti nove, robusnije numeričke metode koje se mogu koristiti u praksi kao pouzdan numerički softver. U disertaciji se proučavaju dvije vrste metoda: direktne i iterativne. Direktne metode se razvijaju za računanje svih svojstvenih vrijednosti i odgovarajućih svojstvenih vektora zadanog problema. Kada nas zanima samo dio spektra, recimo one svojstvene vrijednosti koje su najveće po modulu ili one koje se nalaze u lijevoj kompleksnoj poluravnini, tada koristimo iterativne metode. Ovdje je najšešće slučaj da je dimenzija originalnog problema mnogo veća od broja svojstvenih vrijednosti koje želimo izračunati. Ideja iterativnih metoda je konstruirati potprostor mnogo manje dimenzije od originalnog problema koji sadrži informaciju o traženom dijelu spektra, a aproksimacija traženog dijela spektra se onda izračuna koristeći projekciju problema na nađeni potprostor. Osnova većine metoda za rješavanje polinomnih svojstvenih problema je linearizacija, to jest polinomni problem se zamijeni ekvivalentnim linearnim problemom koji se onda rješava koristeći već razvijene metode za linearne probleme. Međutim, naivno direktno korištenje linearnih metoda ne garantira zadovoljavajuće rezultate za originalni problem. Čak i ako izračunati svojstveni par ima malu grešku unazad za odgovarajuću linearizaciju, greška unazad za rekonstruirani svojstveni par originalnog problema može biti puno veća. Prije razvijanja metoda, u Poglavlju 2 je predstavljena analiza grešaka unazad za polinomni svojstveni problem, bazirana na radu F. Tisseur [66]. Ideja analize grešaka unazad je da se izračunate aproksimacije interpretiraju kao egzaktna rješenja problema koji je blizu originalnom problemu, i čiji matrični koeficijenti su definirani kao $A_\ell + \Delta A_\ell$ pri čemu je $\Delta A_\ell$ malo. Međutim, u mnogim primjenama matrice $A_\ell$ imaju određenu strukturu, npr. hermitske su, ili anti hermitske. Prema tome, bilo bi prirodno zahtijevati da greška unazad $\Delta A_\ell$ čuva ovu strukturu. U slučaju kad je ta struktura hermitska i anti hermitska, postojeći rezultati za realne svojstvene vrijednosti su prošireni na općenite svojstvene vrijednosti. U poglavlju 3 se proučavaju direktne metode za rješavanje kvadratičnog svojstvenog problema. Standardni pristup je korištenje QZ algoritma na odgovarajućoj linearizaciji. Međutim, ako originalni problem ima svojstvene vrijednosti koje su nula ili beskonačno, ovakav pristup je sklon numeričkim poteškoćama. 2011. Hammarling, Munro i Tisseur [37] su razvili quadeig algoritam koji prije korištenja QZ metode za linearni problem skalira originalni problem kako bi norme matričnih koeficijenata bile ujednačene te pokuša detektirati postojanje svojstvenih vrijednosti nula i beskonačno koje ona procesom deflacije ukloni iz linearizacije. Deflacija se temelji na određivanju ranga matrica M i K. Kod quadeiga se koristi QR faktorizacija pivotiranjem stupaca. Koristeći ortogonalne transformacije $n-rank(M)$ beskonačnih i $n-rank(K)$ svojstvenih vrijednosti nula je uklonjeno iz odgovarajuće linearizacije. Glavni doprinos ovog poglavlja je novi algoritam za nalaženje svih svojstvenih vrijednosti kvadratično problema kojeg zovemo KVADeig. Kao motivacija za potrebu poboljšanja quadeiga je predstavljen primjer kod kojeg quadeig nije uspio detektirati sve beskonačne svojstvene vrijednosti. Štoviše, nakon što je uklonjen određen broj ovih svojstvenih vrijednosti, preostale izračunate svojstvene vrijednosti koje su konačne čak nemaju ni veliku apsolutnu vrijednost koja bi nas možda mogla nagnati na zaključak da bi one trebale biti proglašene beskonačnim. Problem nastane kada postoji više od jednog Jordanovog bloka za svojstvene vrijednosti nula i beskonačno. Naime, deflacija u quadeigu ukloni samo jedan Jordanov blok. Kako bismo riješili ovaj problem razvili smo test koji služi za provjeru postoji li više od jednog Jordanovog bloka za svojstvene vrijednosti nula i beskonačno. On je baziran na Van Doorenovom algoritmu za određivanje Kroneckerove strukture generaliziranog svojstvenog problema. Dodatno se analizira utjecaj metoda koje se koriste kao faktorizacije za određivanje ranga te utjecaj kriterija po kojem se rang određuje. Pored skaliranja koje je predloženo u quadeigu uvodimo i dvostrano dijagonalno balansiranje čiji je cilj ujednačavanje elemenata u matricama koje definiraju problem. Na kraju razvijamo metodu baziranu na LU faktorizaciji potpunim pivotiranjem za određivanje ranga. Numerički eksperimenti u Sekciji 3.7 ilustriraju prednosti predložene metode. U poglavlju 4 je razvijen novi algoritam KVARTeig za rješavanje polinomnog svojstvenog problema stupnja četiri. Umjesto direktne linearizacije koristimo kvadratifikaciju koja je uvedena u [17], tj. definiramo ekvivalentan kvadratični problem. Novi algoritam je baziran na KVADeigu, s tim da je skaliranje definirano na matricama originalnog problema i proces deflacije je prilagođen tako da što više iskoristi strukturu originalnog problema. Kao i za kvadratični problem, i ovdje je razvijen test za provjeru postojanja više od jednog Jordanovog bloka za svojstvene vrijednosti nula i beskonačno. Numerički primjeri u Sekciji 4.5 prikazuju prednost nove metode nad quadeigom i polyeigom koji je implementiran u MATLABu. U Poglavlju 5 se proučavaju iterativne metode Arnoldijevog tipa za kvadratični svojstveni problem. Bai i Su [3] su prvi primijetili da je u slučaju iterativnih metoda Arnoldijevog tipa bolje primijeniti Rayleigh–Ritzovu projekciju direktno na originalni kvadratični problem. U tu svrhu su definirani Krilovljev potprostor drugog reda i odgovarajući algoritam SOAR (Second Order Arnoldi) za računanje odgovarajuće baze. Ovaj algoritam je dodatno modificiran te je razvijen takozvani TOAR (Two level orthogonal Arnoldi) algoritam [49]. U ovom poglavlju predlažemo nekoliko modifikacija implicitno restartanog TOAR algoritma koje su temeljene na činjenici da algoritam koristimo za rješavanje kvadratičnog problema svojstvenih vrijednosti. Pod implicitnim restartanjem se misli na korištenje polinomih filtera kako bi se definirao novi početni vektor koji uvelike utječe na konvergenciju metode. Za posebnu klasu pregušenih problema svojstvenih vrijednosti predlažemo novi način definiranja polinomnih filtera. Također, za općenite probleme, predlažemo novi izbor početnog vektora koji se temelji na aproksimaciji kvadratičnog svojstvenog problema problemom čije je gušenje linearno. Numerički primjeri pokazuju da predložene modifikacije rezultiraju manjim brojem restartanja potrebnih za nalaženje svojstvenih parova sa zadovoljavajućom greškom unatrag. U drugom dijelu Poglavlja 5 dajemo pregled implicitno restartanog Krylov–Schurovog algoritma kojeg je uveo Stewart [64]. Ideja ovog algoritma je da se definira faktorizacija koja ne zahtijeva posebnu strukturu kao Arnoldijeva, i na koju će se lakše primijeniti implicitno restartanje. Međutim, prilikom ovakvog restartanja moguće je koristiti samo egzaktne pomake za definiranje polinomnog filtera. Drmač i Bujanović su razvili metodu koja omogućava korištenje proizvoljnih pomaka kod implicitno restartanog Krylov–Schurovog algoritma. U ovom poglavlju generaliziramo predloženi proces u svrhu korištenja Krylov–Schurovog algoritma za rješavanje kvadratičnog svojstvenog problema