Rational Krylov methods for f ( A ) b Michael Eiermann, Oliver G. - PowerPoint PPT Presentation

Rational Krylov methods for f ( A ) b Michael Eiermann, Oliver G. Ernst, and Stefan Güttel DWCAA09 September 8th, 2009

Pr♦❜❧❡♠ ❲❡ ❝♦♥s✐❞❡r t❤❡ ✈❡❝t♦r f ( A ) b ✱ ✇❤❡r❡ ◮ A ✐s ❛ ❧❛r❣❡ N ✲❜②✲ N ♠❛tr✐①✱ ◮ b ✐s ❛ ✈❡❝t♦r ♦❢ ❧❡♥❣t❤ N ✱ ◮ f ✐s ❛ s✉✐t❛❜❧❡ ❢✉♥❝t✐♦♥✳ ❈♦♠♣✉t❡ ❛♣♣r♦①✐♠❛t✐♦♥ f m ≈ f ( A ) b ❢r♦♠ ❛ r❛t✐♦♥❛❧ ❑r②❧♦✈ s♣❛❝❡✳

❲❤❛t ✐s ❛ r❛t✐♦♥❛❧ ❑r②❧♦✈ s♣❛❝❡❄ ▲❡t { ξ 1 , ξ 2 , . . . } ⊆  ⊂ C ❜❡ ❛ ❣✐✈❡♥ s❡q✉❡♥❝❡ ♦❢ ♣♦❧❡s ✳ ❉❡✜♥❡ t❤❡ ♣♦❧②♥♦♠✐❛❧s � m q m ( z ) = ( z − ξ j ) ∈ P m . j = 1 ξ j � = ∞ ❆ss✉♠❡ t❤❛t q m ( A ) − 1 ❡①✐sts✳ ❚❤❡♥ Q m + 1 ( A, b ) = K m + 1 ( A, q m ( A ) − 1 b ) ✐s t❤❡ r❛t✐♦♥❛❧ ❑r②❧♦✈ s♣❛❝❡ ❛ss♦❝✐❛t❡❞ ✇✐t❤ ( A, b , q m ) ✳

❙♣❡❝✐❛❧ ❝❛s❡s ◮  = {∞} ⇒ ♣♦❧②♥♦♠✐❛❧ ❑r②❧♦✈ Q m + 1 = K m + 1 ❬◆❛✉ts ✫ ❲②❛tt ✽✸❪ ❬✈❛♥ ❞❡r ❱♦rst ✽✼❪ ❬❉r✉s❦✐♥ ✫ ❑♥✐③❤♥❡r♠❛♥ ✽✽❪ ❬●❛❧❧♦♣♦✉❧♦s ✫ ❙❛❛❞ ✾✷❪ ❬❍♦❝❤❜r✉❝❦ ✫ ▲✉❜✐❝❤ ✾✼❪ ❬❊✐❡r♠❛♥♥ ✫ ❊r♥st ✵✻❪ ◮  = { ξ } ⇒ s❤✐❢t✲✐♥✈❡rt ❑r②❧♦✈ ❬▼♦r❡t ✫ ◆♦✈❛t✐ ✵✹❪ ❬✈❛♥ ❞❡♥ ❊s❤♦❢ ✫ ❍♦❝❤❜r✉❝❦ ✵✻❪ ◮  = {0 , ∞} ⇒ ❡①t❡♥❞❡❞ ❑r②❧♦✈ ❬❉r✉s❦✐♥ ✫ ❑♥✐③❤♥❡r♠❛♥ ✾✽❪ ❬❑♥✐③❤♥❡r♠❛♥ ✫ ❙✐♠♦♥❝✐♥✐ ✵✽❪ ◮  ❛r❜✐tr❛r② ⇒ r❛t✐♦♥❛❧ ❑r②❧♦✈ ❬❘✉❤❡ ✽✹❪ ❬❇❡❛tt✐❡ ✵✹❪ ❬❇❡❝❦❡r♠❛♥♥ ✫ ❘❡✐❝❤❡❧ ✵✽❪ ❬❑♥✐③❤♥❡r♠❛♥ ❡t ❛❧ ✵✽❪

❘❛t✐♦♥❛❧ ❆r♥♦❧❞✐ ❛❧❣♦r✐t❤♠ ❬❘✉❤❡ ✽✹✴✾✹❪ Input ✿ A, b , { ξ 1 , ξ 2 , . . . , ξ m } v 1 : = b / � b � for j = 1 , 2 , . . . , m do x : = ( I − A/ξ j ) − 1 A v j H ( 1 : j, j ) : = [ v 1 , . . . , v j ] ∗ x x : = x − [ v 1 , . . . , v j ] H ( 1 : j, j ) H ( j + 1 , j ) : = � x � v j + 1 : = x /H ( j + 1 , j ) end ❨✐❡❧❞s ❞❡❝♦♠♣♦s✐t✐♦♥ AV m + 1 ( I m + H m X − 1 m ) = V m + 1 H m ✳

❚❤❡♥ ✐s ♦❢ r❛♥❦ ✳ ❢♦r ❛ ✈❡❝t♦r ✳ ❋♦r ❡✈❡r② ✈❡❝t♦r t❤❡r❡ ❡①✐sts ❛ ✉♥✐q✉❡ ♣♦❧②♥♦♠✐❛❧ ✱ ✱ s✉❝❤ t❤❛t ✳ ❍❡♥❝❡✱ ✐❢ ✐s ✐♥✈❡rt✐❜❧❡✱ ✳ ❘❛t✐♦♥❛❧ ❑r②❧♦✈ ❞❡❝♦♠♣♦s✐t✐♦♥s ❚❤❡♦r❡♠ ✭●✳✱ ✷✵✵✾✮✿ ▲❡t ❛ ❣❡♥❡r❛❧ ❞❡❝♦♠♣♦s✐t✐♦♥ AV m + 1 K m = V m + 1 H m ❜❡ ❣✐✈❡♥✱ ✇❤❡r❡ V m + 1 ❤❛s m + 1 ❧✐♥❡❛r❧② ✐♥❞❡♣❡♥❞❡♥t ❝♦❧✉♠♥s✱ K m ∈ C ( m + 1 ) × m ✱ H m ∈ C ( m + 1 ) × m ✱ ❛♥❞ H m ✐s ♦❢ r❛♥❦ m ✳

❘❛t✐♦♥❛❧ ❑r②❧♦✈ ❞❡❝♦♠♣♦s✐t✐♦♥s ❚❤❡♦r❡♠ ✭●✳✱ ✷✵✵✾✮✿ ▲❡t ❛ ❣❡♥❡r❛❧ ❞❡❝♦♠♣♦s✐t✐♦♥ AV m + 1 K m = V m + 1 H m ❜❡ ❣✐✈❡♥✱ ✇❤❡r❡ V m + 1 ❤❛s m + 1 ❧✐♥❡❛r❧② ✐♥❞❡♣❡♥❞❡♥t ❝♦❧✉♠♥s✱ K m ∈ C ( m + 1 ) × m ✱ H m ∈ C ( m + 1 ) × m ✱ ❛♥❞ H m ✐s ♦❢ r❛♥❦ m ✳ ❚❤❡♥ 1. K m ✐s ♦❢ r❛♥❦ m ✳ 2. colspan ( V m + 1 ) = K m + 1 ( A, q ) ❢♦r ❛ ✈❡❝t♦r q ✳ 3. ❋♦r ❡✈❡r② ✈❡❝t♦r b ∈ colspan ( V m + 1 ) t❤❡r❡ ❡①✐sts ❛ ✉♥✐q✉❡ ♣♦❧②♥♦♠✐❛❧ q m ✱ deg ( q m ) ≤ m ✱ s✉❝❤ t❤❛t b = q m ( A ) q ✳ ❍❡♥❝❡✱ ✐❢ q m ( A ) ✐s ✐♥✈❡rt✐❜❧❡✱ colspan ( V m + 1 ) = Q m + 1 ( A, b ) ✳

❘❛t✐♦♥❛❧ ❑r②❧♦✈ ❛♣♣r♦①✐♠❛t✐♦♥s ❆ s♣❡❝✐❛❧ ❝❛s❡ ✐s t❤❡ ✭r❡❞✉❝❡❞✮ ❞❡❝♦♠♣♦s✐t✐♦♥ AV m K m = V m + 1 H m . ❆s ❛♥ ❛♣♣r♦①✐♠❛t✐♦♥ t♦ f ( A ) b ✇❡ ❝♦♥s✐❞❡r f m : = V m f ( H m K − 1 m ) V † m b . ❚❤❡♦r❡♠ ✭■♥t❡r♣♦❧❛t✐♦♥✮✿ ❚❤❡r❡ ❤♦❧❞s p m − 1 f m = r m ( A ) b = ( A ) b , q m − 1 ✇❤❡r❡ r m ❍❡r♠✐t❡✲✐♥t❡r♣♦❧❛t❡s f ❛t Λ ( H m K − 1 m ) ✳

❊①❛♠♣❧❡ ❚❤❡ ✐t❡r❛t✐♦♥ v 1 = b , ( I − A/ξ j ) − 1 ( A − α j I ) v j , β j v j + 1 = j = 1 , . . . , m, ②✐❡❧❞s ❛ ❞❡❝♦♠♣♦s✐t✐♦♥ AV m + 1 K m = V m + 1 H m ✇✐t❤ V m + 1 = [ v 1 , . . . , v m + 1 ] ✱     1 α 1          β 1 / ξ 1 1   β 1 α 2          ✳✳✳ ✳✳✳     K m = β 2 / ξ 2 ❛♥❞ H m = β 2 .         ✳✳✳ ✳✳✳      1   α m      β m / ξ m β m

P A I N Preassigned Poles and Interpolation Nodes m e t h o d ❊①❛♠♣❧❡ ❚❤❡ ✐t❡r❛t✐♦♥ v 1 = b , ( I − A/ξ j ) − 1 ( A − α j I ) v j , β j v j + 1 = j = 1 , . . . , m, ❝❛♥ ❜❡ ✉s❡❞ ❢♦r ❡①♣❧✐❝✐t r❛t✐♦♥❛❧ ✐♥t❡r♣♦❧❛t✐♦♥✿ ❇② ❚❤❡♦r❡♠ ✭■♥t❡r♣♦❧❛t✐♦♥✮ ✇❡ ❦♥♦✇ t❤❛t p m − 1 f m = V m f ( H m K − 1 m ) e 1 = r m ( A ) b = ( A ) b , q m − 1 ✇❤❡r❡ r m ❍❡r♠✐t❡✲✐♥t❡r♣♦❧❛t❡s f ❛t Λ ( H m K − 1 m ) = { α 1 , . . . , α m } ✳

P ❊①❛♠♣❧❡ A I N Preassigned Poles and Interpolation Nodes ❚❤❡ ✐t❡r❛t✐♦♥ m e t h o d v 1 = b , ( I − A/ξ j ) − 1 ( A − α j I ) v j , β j v j + 1 = j = 1 , . . . , m, ❝❛♥ ❜❡ ✉s❡❞ ❢♦r ❡①♣❧✐❝✐t r❛t✐♦♥❛❧ ✐♥t❡r♣♦❧❛t✐♦♥✿ ❇② ❚❤❡♦r❡♠ ✭■♥t❡r♣♦❧❛t✐♦♥✮ ✇❡ ❦♥♦✇ t❤❛t p m − 1 f m = V m f ( H m K − 1 m ) e 1 = r m ( A ) b = ( A ) b , q m − 1 ✇❤❡r❡ r m ❍❡r♠✐t❡✲✐♥t❡r♣♦❧❛t❡s f ❛t Λ ( H m K − 1 m ) = { α 1 , . . . , α m } ✳

❘❡♠❛r❦s ◮ 2 ✈❡❝t♦rs st♦r❛❣❡ ♥❡❡❞✱ 0 ✐♥♥❡r✲♣r♦❞✉❝ts ◮ ■❢ ❛❧❧ ξ j = ∞ ⇒ ♣♦❧②♥♦♠✐❛❧ ✐♥t❡r♣♦❧❛t✐♦♥ ❛t { α 1 , . . . , α m } ◮ P♦❧②♥♦♠✐❛❧ ✐♥t❡r♣♦❧❛t✐♦♥ ♠❡t❤♦❞s ❤❛✈❡ ❜❡❡♥ ❝♦♥s✐❞❡r❡❞ ❜❡❢♦r❡ ❬❍✉✐s✐♥❣❛ ❡t ❛❧ ✾✾❪ ❬❇❡r❣❛♠❛s❝❤✐✱ ❈❛❧✐❛r✐ ✫ ❱✐❛♥❡❧❧♦ ✵✹❪ ◮ ❋♦r { α 1 , . . . , α m } ✉s❡ ▲❡❥❛ ♣♦✐♥ts✱ s❝❛❧❡❞ t♦ ❛ s❡t ♦❢ ✉♥✐t ❝❛♣❛❝✐t② ❢♦r st❛❜✐❧✐t② ❬❘❡✐❝❤❡❧ ✾✵❪ ✳ ◮ ◆♦ s✉❝❤ s❝❛❧✐♥❣ ✐s ♥❡❝❡ss❛r② ✇✐t❤ t❤❡ P❆■◆ ♠❡t❤♦❞✿ s✐♠♣❧② ❝❤♦♦s❡ β j s✉❝❤ t❤❛t � v j + 1 � = 1 ✱ j = 1 , . . . , m ✳ ◮ ❋♦r r❛t✐♦♥❛❧ ✐♥t❡r♣♦❧❛t✐♦♥ ✉s❡ ▲❡❥❛✲❇❛❣❜② ♣♦✐♥ts✳

� A b ✱ A = diag ( 1 , . . . , 1000 ) ✱ b = [ 1 , . . . , 1 ] T ✳ ❈♦♠♣✉t❡✿ f ( A ) b =

� A b ✱ A = diag ( 1 , . . . , 1000 ) ✱ b = [ 1 , . . . , 1 ] T ✳ ❈♦♠♣✉t❡✿ f ( A ) b =

0 10 � κ − 1 R = � κ + 1 ≈ 0 . 94 −5 10 error −10 10 R ≈ 0 . 36 −15 10 0 20 40 60 80 100 120 140 160 order m

0 10 −5 10 error −10 10 ❩♦❧♦t❛r❡✈✬s ♣♦❧❡s −15 10 0 20 40 60 80 100 120 140 160 order m