Index of /rozprawy2/10850

Pełen tekst

(1)Akademia Górniczo-Hutnicza im. Stanisława Staszica w Krakowie Wydział Matematyki Stosowanej Katedra Analizy Matematycznej, Matematyki Obliczeniowej i Metod Probabilistycznych. Rozprawa doktorska. Efektywne ogólne metody liniowe dla równań różniczkowych zwyczajnych Michał Braś. Promotor: prof. dr hab. inż. Zdzisław Jackiewicz Kraków 2014.

(2) Mojemu Promotorowi, Profesorowi Zdzisławowi Jackiewiczowi.

(3) Spis treści. Spis rysunków . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. 4. Spis tablic . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. 6. Streszczenie . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. 7. Abstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. 8. Wstęp . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. 9. Rozdział 1. Wprowadzenie do ogólnych metod liniowych . . . . . . . . .. 13. 1.1. Ogólne metody liniowe . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.2. Stabilność, zgodność, zbieżność, warunki rzędu metody . . . . . . . . . . . . . 1.3. Stabilność ogólnych metod liniowych . . . . . . . . . . . . . . . . . . . . . . . . 1.4. Przykłady ogólnych metod liniowych . . . . . . . . . . . . . . . . . . . . . . . . 1.4.1. Metody Rungego–Kutty . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.4.2. Dwukrokowe metody Rungego-Kutty . . . . . . . . . . . . . . . . . . . . . 1.4.3. Metody liniowe wielokrokowe . . . . . . . . . . . . . . . . . . . . . . . . . . 1.4.4. Metody DIMSIM . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.5. Postać Nordsiecka ogólnych metod liniowych . . . . . . . . . . . . . . . . . . . 1.6. Ogólne metody liniowe z drugą pochodną - uogólnienie ogólnych metod liniowych. 13 15 18 20 20 21 22 23 24 26. Rozdział 2. Konstrukcja metod Nordsiecka z wbudowaną kwadratową funkcją stabilności . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. 28. 2.1. Metody Nordsiecka . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.2. Warunki rzędu i formuły reprezentacyjne . . . . . . . . . . . . . . . . . . . . 2.3. Wbudowana kwadratowa stabilność . . . . . . . . . . . . . . . . . . . . . . . 2.3.1. Charakteryzacja metod z wbudowaną kwadratową stabilnością . . . . . . 2.3.2. Postać macierzy X . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.4. Konstrukcja metod zamkniętych . . . . . . . . . . . . . . . . . . . . . . . . . 2.4.1. Konstrukcja metod zamkniętych o parametrach p = q = s = 1 oraz r = 2 2.4.2. Konstrukcja metod zamkniętych o parametrach p = q = s = 2 oraz r = 3 2.4.3. Konstrukcja metod zamkniętych o parametrach p = q = s = 3 oraz r = 4 2.4.4. Konstrukcja metod zamkniętych o parametrach p = q = s = 4 oraz r = 5 2.4.5. Konstrukcja metod zamkniętych o parametrach p = q = s = 5 oraz r = 6 2.4.6. Konstrukcja metod zamkniętych o parametrach p = q = s = 6 oraz r = 7. 29 30 33 33 35 37 38 38 41 43 45 47. . . . . . . . . . . . .. 1.

(4) Spis treści 2.5. Konstrukcja metod otwartych . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.6. Przykłady otwartych metod Nordsiecka o parametrach s = p = q oraz r = s + 1 2.6.1. Konstrukcja metod otwartych o parametrach p = q = s = 1, r = 2 . . . . . 2.6.2. Konstrukcja metod otwartych o parametrach p = q = s = 2, r = 3 . . . . . 2.6.3. Konstrukcja metod otwartych o parametrach p = q = s = 3, r = 4 . . . . . 2.6.4. Konstrukcja metod otwartych o parametrach p = q = s = 4, r = 5 . . . . . 2.6.5. Konstrukcja metod otwartych o parametrach p = q = s = 5, r = 6 . . . . . 2.6.6. Konstrukcja metod otwartych o parametrach p = q = s = 6, r = 7 . . . . . 2.6.7. Konstrukcja metod otwartych o parametrach p = q = s = 7, r = 8 . . . . . 2.6.8. Konstrukcja metod otwartych o parametrach p = q = s = 8, r = 9 . . . . . 2.6.9. Kwadratowa stabilność i wbudowana kwadratowa stabilność . . . . . . . . 2.7. Struktura macierzy A . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. 49 51 51 52 52 53 54 55 56 57 59 59. Rozdział 3. Implementacja metod Nordsiecka . . . . . . . . . . . . . . . . .. 61. 3.1. Eksperymenty ze stałym krokiem metody . . . . . . . . . . 3.2. Zmiennokrokowe sformułowanie metod Nordsiecka . . . . . 3.3. Elementy implementacji . . . . . . . . . . . . . . . . . . . . 3.3.1. Propagacja błędu . . . . . . . . . . . . . . . . . . . . . 3.3.2. Wybór długości kroku całkowania . . . . . . . . . . . . 3.3.3. Przeskalowanie wektora Nordsiecka . . . . . . . . . . . 3.3.4. Procedura startowa oraz kończąca . . . . . . . . . . . . 3.4. Eksperymenty numeryczne z jawnymi metodami Nordsiecka 3.4.1. Problemy testowe . . . . . . . . . . . . . . . . . . . . . 3.4.2. Wyniki eksperymentów numerycznych . . . . . . . . . .. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . zmiennym kroku . . . . . . . . . . . . . . . . . . . .. 61 63 64 64 69 71 71 72 72 73. Rozdział 4. Algebraicznie stabilne metody Nordsiecka . . . . . . . . . . .. 81. 4.1. 4.2. 4.3. 4.4. 4.5. 4.6.. . . . . . . . o . .. Kryteria algebraicznej stabilności . . . . . . . . . . . . . . . . . . Poszukiwanie algebraicznie stabilnych metod w oparciu o funkcję Konstrukcja metod o parametrach p = q = s = r − 1 = 1 . . . . Konstrukcja metod o parametrach p = q = s = r − 1 = 2 . . . . Konstrukcja metod o parametrach p = q = s = r − 1 = 3 . . . . Konstrukcja metod o parametrach p = q = s = r − 1 = 4 . . . .. . . . . . . Nyquista . . . . . . . . . . . . . . . . . . . . . . . .. . . . . . .. 82 84 85 87 92 95. Rozdział 5. Metody Nordsiecka z podwyższonym rzędem . . . . . . . . .. 98. 5.1. Warunki rzędu . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.2. Formuły reprezentacyjne . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.3. Propagacja błędu oraz estymacja błędu lokalnego . . . . . . . . . . . . . 5.4. Konstrukcja metod . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.4.1. Konstrukcja metod o parametrach p = r = 2, s = 1, q = 1 lub q = 2 5.4.2. Konstrukcja metod o parametrach p = r = 3, s = 2, q = 2 lub q = 3 5.4.3. Konstrukcja metod o parametrach p = r = 4, s = 3, q = 3 lub q = 4 5.5. Konstrukcja metod rzędu p = s + 1 z własnością IQS . . . . . . . . . . . 5.5.1. Konstrukcja metod o parametrach p = r = 2, q = s = 1 . . . . . . . 5.5.2. Konstrukcja metod o parametrach p = r = 3, q = s = 2 . . . . . . . 5.5.3. Konstrukcja metod o parametrach p = r = 4, q = s = 3 . . . . . . . 5.5.4. Konstrukcja metod o parametrach p = r = 5, q = s = 4 . . . . . . .. . . . . . . . . . . . .. . . . . . . . . . . . .. . . . . . .. . . . . . . . . . . . .. . . . . . . . . . . . .. 99 101 105 112 113 113 118 121 121 122 123 124. 2.

(5) Spis treści 5.5.5. Konstrukcja metod o parametrach p = r = 6, q = s = 5 . . . . . . . . . . . 126 5.5.6. Konstrukcja metod o parametrach p = r = 7, q = s = 6 . . . . . . . . . . . 128 5.5.7. Konstrukcja metod o parametrach p = r = 8, q = s = 7 . . . . . . . . . . . 130. Rozdział 6. Podsumowanie . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 134 Bibliografia . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 137. 3.

(6) Spis rysunków. [. 2.1. L- i A-stabilne metody Nordsiecka z c =. 2.2. L- i A-stabilne metody Nordsiecka z c =. 2.3. L- i A-stabilne metody Nordsiecka z c =. [. 1 2 3 3. 1. 1 1 3 4 2 4. [. ]T. 1. . . . . . . . . . . . . . . . . .. 42. . . . . . . . . . . . . . . . .. 44. ]T. 1 2 3 4 5 5 5 5. 1. ]T. , λ ∈ [0, 1], v12 ∈ [−3, 1],. v13 ∈ [−2.5, −4] . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . [. 2.4. L- i A-stabilne metody Nordsiecka z c =. 1 2 3 4 5 5 5 5. ]T. 1. przy ustalonym v14 = −2. oraz v13 = −3 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . [. 2.5. L- i A-stabilne metody Nordsiecka z c =. 1 1 1 2 5 6 3 2 3 6. 46. 1. ]T. 46. przy ustalonym. v13 = −3, v14 = 0 oraz v15 = 2 . . . . . . . . . . . . . . . . . . . . . . . . . . . .. 48. 2.6. Obliczanie pola powierzchni obszaru absolutnej stabilności . . . . . . . . . . . .. 50. 2.7. Obszar absolutnej stabilności metody RK rzędu p = 1 oraz metody Nordsiecka p = q = s = 1, r = 2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. 2.8. Obszar absolutnej stabilności metody RK rzędu p = 2 oraz metody Nordsiecka p = q = s = 2, r = 3 z IQS . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. 2.9. 52 53. Obszar absolutnej stabilności metody RK rzędu p = 3 oraz metody Nordsiecka p = q = s = 3, r = 4 z IQS . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. 53. 2.10 Obszar absolutnej stabilności metody RK rzędu p = 4 oraz metody Nordsiecka p = q = s = 4, r = 5 z IQS . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. 54. 2.11 Obszar absolutnej stabilności metody RK rzędu p = 5 oraz metody Nordsiecka p = q = s = 5, r = 6 z IQS . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. 55. 2.12 Obszar absolutnej stabilności metody RK rzędu p = 6 oraz metody Nordsiecka p = q = s = 6, r = 7 z IQS . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. 56. 2.13 Obszar absolutnej stabilności metody RK rzędu p = 7 oraz metody Nordsiecka p = q = s = 7, r = 8 z IQS . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. 57. 2.14 Obszar absolutnej stabilności metody RK rzędu p = 8 oraz metody Nordsiecka p = q = s = 8, r = 9 z IQS . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. 58. 4.

(7) 3.1. Metoda rzędu p = 3. Lewy rysunek - standardowy operator, prawy P I-operator; górny rysunek: odrzucone kroki (*) vs. długość kroku; dolny rysunek: błąd lokalny (◦) vs. jego aproksymacja (·). . . . . . . . . . . . . . . . .. 3.2. 75. Metoda rzędu p = 4. Lewy rysunek - standardowy operator, prawy P I-operator; górny rysunek: odrzucone kroki (*) vs. długość kroku; dolny rysunek: błąd lokalny (◦) vs. jego aproksymacja (·). . . . . . . . . . . . . . . . .. 75. 4.1. ˜ Wartości własne λ1 = λ1 (θ) i λ2 = λ2 (θ) macierzy He(DN(ξ)) dla ξ = exp(iθ), 88. 4.2. θ ∈ [0, 2π] . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ˜ Wartości własne λ1 = λ1 (θ) macierzy He(DN(ξ)) dla metody z zadanym. 89. 4.3. dla ξ = exp(iθ), θ ∈ [0, 2π] . . . . . . . . . . . . . . . . . . . . . . . . ˜ Wartości własne λ1 = λ1 (θ) macierzy He(DN(ξ)) dla L-stabilnej metody dla ξ = exp(iθ), θ ∈ [0, 2π] . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. 89. [. c=. 1 2. 1. ]T. [. 1 2. ]T. T. 4.5. (lewy rysunek) i c = [0 1] (prawy rysunek) 91 ˜ Wartości własne λ1 = λ1 (θ) macierzy He(DN(ξ)) dla metody rzędu 92. 4.6. p = q = s = 3 z wolnym wektorem c dla ξ = exp(iθ), θ ∈ [0, 2π] . . . . . . . . . ˜ Wartości własne λ1 = λ1 (θ) macierzy He(DN(ξ)) metody rzędu p = q = s = 3 o. 92. 4.7. wektorze c = [−1 0 1]T dla ξ = exp(iθ), θ ∈ [0, 2π] . . . . . . . . . . . . . . . . . ˜ Wartości własne λ1 = λ1 (θ) macierzy He(DN(ξ)) metody rzędu p = q = s = 3. 93. 4.8. typu 4 o wolnym wektorze c dla ξ = exp(iθ), θ ∈ [0, 2π] . . . . . . . . . . . . . . ˜ Wartości własne λ1 = λ1 (θ) macierzy He(DN(ξ)) metody rzędu p = q = s = 4 o wolnym wektorze c dla ξ = exp(iθ), θ ∈ [0, 2π] . . . . . . . . . . . . . . . . . . .. 95. 4.4. Residuum res versus λ dla c =. 1. 5.1. A-stabilne metody rzędu poziomów q = 2 oraz rzędu metody p = 3 . . . . . . . . 115. 5.2. A- i L-stabilne metody rzędu poziomów q = 4 oraz rzędu metody p = 5 . . . . . 125. 5.3. A- i L-stabilne metody rzędu poziomów q = 5 oraz rzędu metody p = 6 . . . . . 127. 5.4. A- i L-stabilne metody rzędu poziomów q = 6 oraz rzędu metody p = 7 . . . . . 129. 5.5. A- i L-stabilne metody rzędu poziomów q = 7 oraz rzędu metody p = 8 . . . . . 131.

(8) Spis tablic. 3.1. Wyniki eksperymentów numerycznych ze stałym krokiem dla metod rzędu p = 1, 2, . . . , 6 dla problemu (3.1) . . . . . . . . . . . . . . . . . . . . . . . . . . .. 3.2. Wyniki eksperymentów numerycznych ze stałym krokiem dla metod rzędu p = 6, 7, 8 dla problemu (3.1) . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. 3.3. . . . . . . . . . . . . . . . . . . . . . . .. . . . . . . . . . . . . . . . . . . . . . . .. 79. Porównanie metody CERK43 oraz jawnej metody Nordsiecka rzędu p = 4 dla problemu (3.17) ze standardowym operatorem wyboru długości kroku . . . . . .. 3.8. 78. Wyniki numeryczne dla problemu van der Pohla (3.18) dla t ∈ [0, 20] z P I-operatorem wyboru długości kroku. 3.7. 77. Wyniki numeryczne dla problemu van der Pohla (3.18) dla t ∈ [0, 20] ze standardowym operatorem wyboru długości kroku . . . . . . . . . . . . . . . . .. 3.6. 76. Wyniki numeryczne dla problemu Prothero–Robinsona (3.17) dla t ∈ [0, 100] z P I-operatorem wyboru długości kroku. 3.5. 63. Wyniki numeryczne dla problemu Prothero–Robinsona (3.17) dla t ∈ [0, 100] ze standardowym operatorem wyboru długości kroku . . . . . . . . . . . . . . . . .. 3.4. 62. 80. Porównanie metody CERK43 oraz jawnej metody Nordsiecka rzędu p = 4 dla problemu (3.18) ze standardowym operatorem wyboru długości kroku . . . . . .. 80. 6.

(9) Streszczenie W rozprawie zajmujemy się badaniem metod Nordsiecka, szerokiej podklasy ogólnych metod liniowych. W szczególności omówimy konstrukcję otwartych oraz zamkniętych metod Nordsiecka z tak zwaną wbudowaną kwadratową stabilnością. Zaprezentowane zostaną zagadnienia związane z implementacją skonstruowanych metod. Zaczniemy od wprowadzenia pojęcia wbudowanej kwadratowej stabilności. Są to pewne algebraiczne warunki, które gwarantują, że funkcja stabilności metody Nordsiecka ma dwa niezerowe pierwiastki. W związku z tym badanie stabilności można ograniczyć do kwadratowego czynnika wielomianu stabilności, co znacznie ułatwia analizę metod i umożliwia konstrukcję metod otwartych oraz zamkniętych o wysokim rzędzie metody. Wyniki te są istotnym rozszerzeniem znanych w literaturze metod z tak zwaną wbudowaną stabilnością Rungego-Kutty. Efektywność skonstruowanych otwartych metod Nordsiecka zostanie sprawdzona w testach numerycznych ze stałym oraz ze zmiennym krokiem metody. W szczególności porównamy standardowy operator wyboru długości kroku oraz wywodzący się z teorii optymalnego sterowania P I-operator. Kolejnym rozważanym przypadkiem są metody Nordsiecka z podwyższonym rzędem metody. Rząd metody jest tu o jeden wyższy niż liczba aproksymacji wewnętrznych, zatem skonstruowane metody są efektywniejsze w stosunku do metod zazwyczaj rozważanych w literaturze. Konstrukcja metod jest jednak znacznie trudniejsza. Ostatnim rozważanym zagadnieniem jest konstrukcja algebraicznie stabilnych metod Nordsiecka. Konstrukcję takich metod przeprowadzimy poprzez numeryczną optymalizację funkcji celu opartej na funkcji stabilności Nyquista. Własności algebraicznie stabilnych metod są przydatne przy rozwiązywaniu monotonicznych układów równań różniczkowych zwyczajnych.. Słowa kluczowe ogólne metody liniowe, postać Nordsiecka, warunki rzędu oraz rzędu poziomów metody, liniowa oraz nieliniowa stabilność, implementacja. 7.

(10) Abstract We study in this thesis the class of Nordsieck methods, a subclass of general linear methods for ordinary differential equations. We consider the construction of explicit and implicit Nordsieck methods with so-called inherent quadratic stability property. We also present some results of numerical experiments. We begin with discussing the inherent quadratic stability property. It means that some algebraic conditions that guarantee that stability function has only two non-zero roots. Hence, the analysis of linear stability properties can be restricted to quadratic part of stability polynomial. These results are substantial extension of well-known inherent Runge-Kutta stability property. The efficiency of constructed explicit methods is tested in several numerical experiments with constant and variable stepsize. In particular, we compare the standard stepsize control strategy and based on control theory arguments, so-called P I stepsize control. Next, we search for efficient Nordsieck methods, i.e., methods with order of the method higher then number of stages, so the constructed methods are more efficient then methods usually considered in the literature on the subject. However, the construction of such methods is much more complicated. Finally, we consider the construction of algebraically stable Nordsieck methods. The numerical search for such methods is based on Nyquist stability function. These methods are suitable for solving monotone systems of ordinary differential equations.. Key words general linear methods, Nordsieck form, order nad stage order conditions, linear and non-linear stability, implementation. 8.

(11) Wstęp. Otrzymywanie analitycznego rozwiązania problemu Cauchy’ego dla równań różniczkowych zwyczajnych jest zadaniem bardzo trudnym, a często niemożliwym. Pomocą w uzyskaniu informacji o rozwiązaniu są metody numeryczne. Najprostszą metodą numeryczną dla tego problemu jest schemat otwarty Eulera. Mając dane rozwiązanie w pewnym punkcie rozważanego przedziału przesuwamy się wzdłuż stycznej do rozwiązania w tym punkcie o pewną odległość h. Z kroku na krok rozwiązujemy jednak inny problem początkowy, co prowadzi do idei błędu lokalnego metody. Porównując otrzymane schematem Eulera rozwiązanie z rozwinięciem Taylora w danym punkcie okazuje się, że są one zgodne z dokładnością do wyrazów pierwszej potęgi kroku h, zatem schemat Eulera ma rząd równy 1. W celu poprawienia dokładności metod numerycznych stosuje się zazwyczaj dwa podejścia. Pierwsze z nich polega na wykorzystaniu informacji o rozwiązaniu uzyskanej w poprzednich krokach, co prowadzi do metod wielowartościowych, w tym metod wielokrokowych. Drugie podejście polega na obliczeniu dodatkowych aproksymacji rozwiązania w dodatkowych punktach, co prowadzi do metod wielopoziomowych, w tym metod Rungego-Kutty. Pierwszymi metodami o wielowartościowym charakterze były metody zaproponowana w 1883 roku przez Bashfortha oraz Adamsa [13], znane teraz pod nazwą metod Adamsa-Bashfortha. Są one szczególną podklasą rodziny liniowych metod wielokrokowych. Już w tej pracy, a później niezależnie w 1926 roku przez Moultona [91], zostały wprowadzone zamknięte metody wielokrokowe - metody Adamsa-Moultona. W porównaniu z metodami Adamsa-Bashfortha metody Adamsa-Moultona potrzebują mniejszej liczby informacji z poprzednich kroków, aby uzyskać taką samą dokładność, jednak. 9.

(12) Wstęp. obliczenie wartości rozwiązania wymaga rozwiązania nieliniowego równania algebraicznego, co często wykonywane jest przy użyciu metody Newtona. Inne podejście, zaproponowane przez Milne [90], polegaj na używaniu metod w parach predyktor-korektor: najpierw rozwiązywanie jest aproksymowane za pomocą otwartej metody (predyktor), a następnie jest poprawiane iteracyjnie za pomocą metody zamkniętej (korektor). W 1962 roku Nordsieck [92] zaobserwował, że wszystkie metody numerycznego całkowania są równoważne poszukiwaniu odpowiedniego wielomianu interpolacyjnego i zaproponował, by reprezentować ten wielomian za pomocą kolejnych skalowanych pochodnych, dzięki czemu zmiana długości kroku w metodach wielokrokowych staje się znacznie łatwiejsza, w przeciwieństwie do standardowej reprezentacji metod wielokrokowych. Ważny wkład w rozwój liniowych metod wielokrokowych wniósł Dahlquist, któremu zawdzięczamy tak zwaną pierwszą i drugą barierę Dahlquista. Pomimo że rząd metody k-krokowej może wynosić teoretycznie 2k, to Dahlquist [58] w 1956 roku pokazał, że rząd zbieżnej k-krokowej metody nie może być większy niż k+1, gdy k jest nieparzyste, większy niż k + 2, gdy k jest parzyste oraz większy niż k, gdy metoda jest jawna. Twierdzenie to nazywane jest pierwszą barierą Dahlquista. Druga bariera Dahlquista [59] z 1963 roku mówi, że A-stabilna liniowa metoda wielokrokowa nie może mieć rzędu wyższego niż 2. Innym problemem, jaki dotyczy liniowych metod wielokrokowych, jest to, że wraz ze wzrostem rzędu metody maleje jej obszar absolutnej stabilności. Spowodowane jest to między innymi tym, że w danym kroku obliczana jest tylko jedna nowa wartość funkji prawej strony rozważanego problemu. Pierwsza metoda wykorzystująca strategię wielopoziomową przypisywana jest Rungemu [95], 1895 rok. Dalszy rozwój to prace Heuna [72], który skonstruował metody niskiego rzędu oraz Kutty [87], który scharakteryzował wszystkie metody rzędu 4 oraz niektóre rzędu 5. Rozważane przez nich metody znane są jako schematy Rungego-Kutty. W tej klasie duże znaczenie mają zaproponowane przez Alexandra [7] w 1977 metody Rungego–Kutty, w których macierz współczynników A jest miacierzą trójkątną dolną. Dzięki takiemu założeniu nieliniowe równania wiążące ze sobą poziomy mogą być rozwiązywane sekwencyjnie, co znacznie redukuje koszt. Z drugiej strony rząd wewnętrznych aproksymacji jest niski, przez co metody te nie nadawały się do rozwiązywania sztywnych układów równań. Innym zaproponowanym rozwiązaniem było żądanie, by macierz A miała jednopunktowe sepktrum [93]. Metodom Rungego-Kutty poświęcone są, między innymi, monografie Butchera [27, 31, 33].. 10.

(13) Wstęp. Jedno z pierwszych połączeń strategii wielowartościowych i wielopoziomowych zostało zaproponowane w 1964 roku przez Gragga i Stettera [67] w celu ominięcia pierwszej bariery Dahlquista dla liniowych metod wielokrokowych. Autorzy rozważali metodę typu predyktor-korektor z dwoma predyktorami: pierwszy z nich był standardowym predyktorem wartości rozwiązania, drugi aproksymował wartość rozwiązania w pewnym punkcie pozasiatkowym, a następnie korektor używał obu wartości uzyskanych przez predyktory do obliczenia aproksymacji rozwiązania w punkcie siatkowym. Metody te autorzy nazwali uogólnionymi metodami predyktor-korektor. Podobne połączenia zaproponowali też Gear [65] - metody hybrydowe oraz Butcher [23] - zmodyfikowane liniowe metody wielokrokowe. Odwrotnie, wykorzystanie cech liniowych metod wielokrokowych w metodach Rungego-Kutty, zaproponowali między innymi Byrne i Lambert [50]. Ogólne metody liniowe zostały zaproponowane w 1966 roku przez Butchera [24] jako połączenie obu podejść oraz dostarczenie ogólnych narzędzi do badania tradycyjnych metod. Najczęściej używane obecnie ich sformułowanie pochodzi z 1980 roku z pracy Butchera i Burrage’a [22]. Monografia Jackiewicza [79] z 2009 roku jest w całości poświęcona ogólnym metodom liniowym oraz podsumowuje stan wiedzy w tej dziedzinie. Na podstawie tej monografii w rozdziale pierwszym tej rozprawy przedstawimy sformułowanie ogólnych metod liniowych wraz z podstawowymi pojęciami. Przedstawimy przykłady klasycznych metod zapisanych w postaci ogólnych metod liniowych, a także pewne nowe klasy metod. Ogólne metody liniowe stanowią szeroką klasę metod numerycznych dla problemu początkowego dla równań różniczkowych zwyczajnych, zatem znalezienie użytecznych metod jest zadaniem nietrywialnym. Dużo uwagi w literaturze poświęcono podklasie metod DIMSIM. Innym kierunkiem badań były ogólne metody liniowe z własnością stabilności Rungego–Kutty. W rozdziale drugim zajmiemy się badaniem metod z wbudowaną kwadratową funkcją stabilności, narzędziem, które umożliwia konstrukcję metod otwartych oraz zamkniętych o wysokim rzędzie metody i które jest mniej ograniczające niż badany wcześniej przypadek wbudowanej stabilności Rungego-Kutty. Zajmiemy się analizą przypadku p = q = s = r − 1, gdzie p jest rzędem metody, q jest rzędem poziomów metody, s jest liczbą poziomów (aproksymacji wewnętrznych) obliczanych w danym kroku, a r jest liczbą aproksymacji zewnętrznych, które przechodzą z jednego kroku do kolejnego. Opiszemy konstrukcję przy użyciu twierdzenia Schura L-stabilnych ogólnych metod liniowych zamkniętych rzędu p = 1, 2, . . . , 6 [14]. 11.

(14) Wstęp. oraz poszukiwania otwartych metod o maksymalnej powierzchni obszaru absolutnej stabilności p = 1, 2, . . . , 8 [15]. Wektor aproksymacji zewnętrznych metod skonstruowanych w rozdziale 2 jest wektorem Nordsiecka, przez co metody te w naturalny sposób można zaimplementować w algorytmach o zmiennej długości kroku. W rozdziale 3 zajmiemy się testowaniem metod jawnych [16] z rozdziału 2 w wybranych eksperymentach numerycznych o stałym oraz zmiennym kroku metody. Rozdział 4 poświęcony jest poszukiwaniu zamkniętych metod rzędu p = s+1, gdzie s jest liczbą poziomów. Metody te są efektywniejsze pod względem kosztu rozumianego jako liczba wywołań funkcji prawej strony równania w stosunku do metod rozważanych w rozdziale 2 oraz najczęściej rozważanych w literaturze. Warunki rzędu otrzymane dla tej rodziny metod są trudniejsze do rozwiązania, a otrzymywanie A- i L-stabilnych metod jest technicznie bardzo trudne. Bezpośrednio możliwe jest skonstruowanie metod rzędu p = 2, 3 oraz 4 [19]. Przy konstrukcji metod wyższego rzędu pomocne jest wykorzystanie warunków wbudowanej kwadratowej stabilności, dzięki czemu możemy skonstruować metody do rzędu p = 8 [18]. Niestety, odbywa się to kosztem ograniczenia rzędu poziomów do wartości q = p − 1. Ostatnim zagadnieniem tej rozprawy, którym zajmujemy się w rozdziale 5, jest konstrukcja algebraicznie stabilnych ogólnym metod liniowych [17]. Definicja algebraicznej stabilności jest bardzo techniczna, jednak metody posiadające tę własność są przydatne przy rozwiązywaniu monotonicznych układów równań różniczkowych zwyczajnych. Konstrukcja metod Nordsiecka opiera się na kryteriach podanych przez Hewitt i Hilla [73, 74] poprzez optymalizację odpowiedniej funkcji celu. Skonstruowane metody mają rząd p = 1, 2, 3 oraz 4. W rozdziale 6 zawarte jest podsumowanie wyników przedstawionych w rozprawie.. 12.

(15) Rozdział 1. Wprowadzenie do ogólnych metod liniowych. 1.1. Ogólne metody liniowe Zajmujemy się metodami numerycznymi dla problemu początkowego dla układu równań różniczkowych zwyczajnych pierwszego rzędu (ODEs) w autonomicznej postaci  ( )   y ′ (t) = f y(t) ,. t ∈ [t0 , T ],.   y(t0 ) = y0 ∈ Rm .. (1.1). Będziemy zakładali, że funkcja f : Rm → Rm jest dostatecznie gładka, a y0 ∈ Rm jest danym warunkiem początkowym. Założenie o autonomiczności problemu ułatwia analizę i upraszcza zapis. Nie jest ono przy tym ograniczające, gdyż dowolne nieautonomiczny układ m równań różniczkowych można sprowadzić do postaci autonomicznej poprzez dodanie trywialnego równania związanego z czasem. Rozwiązanie dokładne tego problemu zazwyczaj będziemy oznaczali przez y(t).. 13.

(16) 1.1. Ogólne metody liniowe. Ogólne metody liniowe (GLMs) zadawane są przez wektor współczynników c = [c1 , . . . , cs ]T , ci ̸= cj i cztery macierze współczynników A = [aij ], U = [uij ], B = [bij ], V = [vij ] o wymiarach A ∈ Rs×s ,. U ∈ Rs×r ,. B ∈ Rr×s ,. V ∈ Rr×r .. Na siatce jednostajnej tn = t0 + nh, n = 0, 1, . . . , N , gdzie krok h = (T − t0 )/N przy zadanej liczbie kroków N + 1, przyjmują one następującą postać   [n]   Yi      [n]     yi [n]. Wielkości Yi. = h. s ∑. [n]. aij f (Yj ) +. j=1. = h. s ∑. r ∑. [n−1]. uij yj. , i = 1, 2, . . . , s,. j=1 [n] bij f (Yj ). j=1. +. r ∑. (1.2) [n−1] vij yj ,. i = 1, 2, . . . , r.. j=1. nazywamy poziomami i pełnią one analogiczną funkcję do poziomów me-. tod Rungego-Kutty. Są one aproksymacjami tak zwanego rzędu poziomu q rozwiązania równania (1.1) w punktach określanych przez wektor c, czyli wielkości y(tn−1 + ci h), i = 1, 2, . . . , s. Poziomy obliczane i używane są jedynie w bieżącym kroku z tn−1 do tn . [n]. Wielkości yi. nazywamy aproksymacjami zewnętrznymi rzędu p kombinacji li[n−1]. niowych pochodnych rozwiązania y(t) problemu (1.1) w punkcie tn . Wielkości yi. obliczone w poprzednim kroku są wykorzystywane w kolejnym kroku do obliczenia [n]. wielkości yi w kroku bieżącym. W dalszej części rozprawy będziemy w szczególności [n]. zainteresowani metodami, gdzie wielkości yi. aproksymują kolejne współrzędne wek-. tora Nordsiecka. Pojęcie rzędu poziomów oraz rzędu metody będzie szczegółowo omawiane w dalszej części tego rozdziału. Stosując zapis . Y. [n]. [n]. Y  1  .. = . .    , . . (. [n].  f Y1 ( )  .. [n] f Y = .  . Ys[n]. f. (. Ys[n]. )    ,  ) . . y [n]. [n]. . y  1   ..  =  . , . . yr[n]. ogólną metodę liniową możemy zapisać w postaci macierzowej   Y [n]  y [n]. = h (A ⊗ Im ) f (Y [n] ) + (U ⊗ Im ) y [n−1] , = h (B ⊗ Im ) f (Y [n] ) + (V ⊗ Im ) y [n−1] ,. (1.3). n = 1, 2 . . . , N , gdzie Im jest macierzą jednostkową wymiaru m, to znaczy wymiaru układu (1.1), a ⊗ jest iloczynem Kroneckera macierzy [89]. Współczynniki metody zapisywane są często w postaci macierzy blokowej wymiary (s + r) × (s + r)  . A U B V. . .. 14.

(17) 1.2. Stabilność, zgodność, zbieżność, warunki rzędu metody. Macierz współczynników A pełni podobną rolę do macierzy A w metodach Rungego-Kutty i określa koszt implementacji metody. W przypadku, gdy macierz A [n]. jest silnie trójkątna dolna, wtedy każdy poziom Yi. zależy od znanych wartości, za-. tem mamy do czynienia z metodami jawnymi (otwartymi). W przeciwnym przypadku mamy do czynienia z metodami uwikłanymi (zamkniętymi), gdzie konieczne jest rozwiązanie układu nieliniowych równań algebraicznych. Ze względu na strukturę macierzy współczynników A możemy wyróżnić 4 typy ogólnych metod liniowych. Niech macierz A ma postać .   λ   a21  . A=  .  . . λ .. .. ... .. as1 as2 . . .. λ.     ,   . gdzie λ = 0 dla metod typu 1 lub λ > 0 dla metod typu 2. Metody 1 typu są przewidziane dla problemów niesztywnych, podczas gdy metody typu 2 przewidziane są dla problemów sztywnych. Oba typy wykorzystuje się przy obliczaniach sekwencyjnych. W przypadku obliczeń równoległych, macierz A jest macierzą diagonalną A = diag(λ, λ, . . . , λ) = λI, przy czym metody typu 3 to metody o współczynniku λ = 0, a metody typu 4 to metody o współczynniku λ > 0.. 1.2. Stabilność, zgodność, zbieżność, warunki rzędu metody Podobnie jak w [79] (porównaj także [31]), w celu znalezienia minimalnych warunków dokładności zakładamy, że istnieją takie wektory q0 oraz q1 , q0 =. [. ]. q1,0 q2,0 . . .. qr,0. q1 =. [. ]. q1,1 q2,1 . . .. qr,1. ,. że dane wejściowe y [n−1] w danym kroku spełniają [n−1]. yi. = qi,0 y(tn−1 ) + qi,1 hy ′ (tn−1 ) + O(h2 ),. i = 1, 2, . . . , r.. Od metody wymagamy, żeby wielkości Y [n] oraz y [n] otrzymane w danym kroku spełniały minimalne warunki dotyczące dokładności, to znaczy by [n]. Yi. = y(tn−1 + ci h) + O(h2 ),. i = 1, 2, . . . , s 15.

(18) 1.2. Stabilność, zgodność, zbieżność, warunki rzędu metody. oraz yi = qi,0 y(tn ) + qi,1 hy ′ (tn ) + O(h2 ), [n]. i = 1, 2, . . . , r.. Prowadzi to do następującej definicji. Definicja 1.1. Ogólna metoda liniowa (1.2) jest przedzgodna (lub prezgodna), jeżeli istnieje taki wektor q0 , że Uq0 = e,. Vq0 = q0 ,. (1.4). gdzie e = [1, . . . , 1]T . Wektor q0 nazywamy wektorem przedzgodności. Zauważmy, że z warunku w powyższej definicji wynika, że q0 jest wektorem własnym macierzy V odpowiadającym wartości własnej równej 1. Definicja 1.2. Ogólna metoda liniowa (1.2) jest zgodna, jeżeli jest przedzgodna z wektorem przedzgodności q0 oraz istnieje taki wektor q1 , że Be + Vq1 = q0 + q1 . Wektor q1 nazywamy wektorem zgodności. Definicja 1.3. Ogólna metoda liniowa (1.2) jest poziomowo zgodna, jeżeli Ae + Uq1 = c. Przykłady wektorów q0 i q1 dla wybranych podklas GLMs można znaleźć w [79]. Definicja 1.4. Mówimy, że ogólna metoda liniowa (1.2) jest zero-stabilna, gdy istnieje taka stała C, że ||Vn || ¬ C dla każdego n 0. Zero-stabilność jest minimalnym warunkiem stabilności, który gwarantuje, że błędy popełniane w danym kroku nie mają zbyt dużego wpływu na dalsze kroki. Związane jest to z równaniem trywialnym y ′ (t) = 0, którego rozwiązanie analityczne jest ograniczone. Rozwiązanie numeryczne otrzymane za pomocą (1.2) ma postać y [n] = Vy [n−1] = Vn y [0] i jest ono także ograniczone, gdy metoda jest zero-stabilna. Bardziej praktyczne jest sformułowanie warunku zerowej stabilności w poniższym twierdzeniu. 16.

(19) 1.2. Stabilność, zgodność, zbieżność, warunki rzędu metody. Twierdzenie 1.5. [27, 79] Ogólna metod liniowa (1.2) jest zero-stabilna, gdy wielomian minimalny macierzy współczynników V nie ma pierwiastków o module większym niż 1, a wszystkie jego pierwiastki o module równym 1 są pojedyncze. Dla metod, którymi zajmujemy się w rozprawie, zero-stabilność jest gwarantowana poprzez wybór postaci macierzy współczynników V. Definicja 1.6. Ogólna metod liniowa (1.2) jest zbieżna, jeżeli dla dowolnego problemu początkowego (1.1) spełniającego warunek Lipschitza ||f (y) − f (z)|| ¬ L||y − z|| istnieje taki niezerowy wektor q0 ∈ Rr oraz taka procedura startowa Sh : Rm → Rmr , lim Sh (y0 ) = lim y [0] = (q0 ⊗ Im )y(t0 ), że dla dowolnego t¯ > t0 ciąg wektorów y [n]. h→0. h→0. obliczonych używając n kroków ogólnej metody liniowej z krokiem h = (t¯ − t0 )/n jest zbieżny do q0 y(t¯), gdzie y(t) jest rozwiązaniem (1.1). Dowód zbieżności metody jest na ogół nietrywialny, jednak poniższe twierdzenie podaje związek między zbieżnością metody, a jej zgodnością i zero-stabilnością. Jest ono analogiczne do znanych wyników dla liniowych metod wielokrokowych oraz metod Rungego–Kutty [27]. Twierdzenie 1.7. [27, 79] Ogólna metod liniowa (1.2) jest zbieżna wtedy i tylko wtedy, gdy jest zgodna i zero-stabilna. Przejdziemy teraz do zdefiniowania pojęcia rzędu poziomów oraz rzędu ogólnej metody liniowej. W tym celu zakładamy, że składowe wektora aproksymacji zewnętrznych y [n−1] używane jako zbiór danych wejściowych w danym kroku z tn−1 do tn = tn−1 + h są kombinacjami liniowymi rzędu p skalowanych pochodnych rozwiązania równania (1.1) w punkcie tn−1 , czyli [n−1] yi. =. p ∑. qik hk y (k) (tn−1 ) + O(hp+1 ),. i = 1, 2, . . . , r,. (1.5). k=0. dla pewnych rzeczywistych współczynników qik , i = 1, 2, . . . , r, k = 0, 1, . . . , p. Definicja 1.8. Mówimy, że ogólna metoda liniowa ma rząd poziomów równy q, jeżeli z [n]. (1.5) wynika, że składowe wektora wewnętrznych aproksymacji Yi. są aproksymacjami. rzędu poziomów q rozwiązania (1.1) w punktach tn−1 + ci h, czyli [n]. Yi. = y(tn−1 + ci h) + O(hq+1 ),. i = 1, 2, . . . , s.. (1.6). 17.

(20) 1.3. Stabilność ogólnych metod liniowych. Definicja 1.9. Mówimy, że ogólna metoda liniowa ma rząd metody równy p, jeżeli z [n]. (1.5) wynika, że składowe wektora zewnętrznych aproksymacji yi po wykonaniu kroku z tn−1 do tn = tn−1 + h spełniają (1.5) z n − 1 zastąpionym przez n, to znaczy [n] yi. =. p ∑. qik hk y (k) (tn ) + O(hp+1 ),. i = 1, 2, . . . , r,. (1.7). k=0. z tymi samymi współczynnikami qik . Używając funkcji zmiennej zespolonej warunki rzędu metody dane przez (1.5), (1.6) i (1.7) można zapisać w postaci, która jest znacznie wygodniejsza do użycia w programach do obliczeń symbolicznych (Mathematica, Maple). Ograniczymy się tutaj do przedstawienia wyników uzyskanych przez Butchera [28] dla przypadku ogólnej metody liniowej z p = q. Wyniki te zostaną użyte w rozdziale 2, a technika dowodowa oraz analiza innych przypadków niż p = q zostanie zaprezentowana w rozdziale 5. Definiujemy zatem wektory qi = [ q1,i · · · oraz macierz W=. qr,i ]T ,. i = 0, 1, . . . , p. [. ]. q0 q1 . . . qp. .. (1.8). (1.9). Definiujemy też wspomnianą funkcję zmiennej zespolonej jako w = w(z) =. p ∑. qi z i ,. z ∈ C.. (1.10). i=0. Twierdzenie 1.10. [28, 79] Ogólna metoda liniowa ma rząd p = q wtedy i tylko wtedy, gdy spełnione są warunki ecz = zAecz + Uw(z) + O(z p+1 ) oraz ez w(z) = zBecz + Vw(z) + O(z p+1 ), gdzie ecz = [ ec1 z · · ·. e cs z ] T .. 1.3. Stabilność ogólnych metod liniowych Zajmiemy się teraz podstawowymi własnościami liniowej stabilności, jakie ogólna metoda liniowa powinna posiadać. Rozważamy skalarne liniowe równanie testowe zaproponowane przez Dahlquista [59] y ′ = ξy,. t 0,. (1.11) 18.

(21) 1.3. Stabilność ogólnych metod liniowych. gdzie ξ ∈ C oraz Re ξ ¬ 0. Rozwiązanie tego problemu jest ograniczone przy t dążącym do nieskończoności, zatem chcemy by rozwiązanie numeryczne otrzymane za pomocą ogólnej metody liniowej posiadało podobne własności. Wymaganie to prowadzi do pojęcia liniowej stabilności. Stosując (1.2) do (1.11), pierwsze równanie (1.3) wiążące ze sobą poziomy metody przyjmuje postać Y [n] = zAY [n] + Uy [n−1] , gdzie z = hξ, co jest równoważnie Y [n] = (I − zA)−1 Uy [n−1] , gdzie chcemy, by macierz (I − zA) była nieosobliwa. Drugie równanie ogólnej metody liniowej (1.3) daje nam, że wektor aproksymacji zewnętrznych w bieżącym kroku ma postać y [n] = zBY [n] + Vy [n−1] = (V + zB(I − zA)−1 U)y [n−1] . Definiujemy macierz stabilności M(z) = V + zB(I − zA)−1 U, w związku z czym wektory aproksymacji zewnętrznych w kolejnych krokach są związane równaniem rekurencyjnym y [n] = M(z)y [n−1] . Definiujemy też funkcję stabilności, jako wielomian charakterystyczny macierzy stabilności M(z) p(w, z) = det(wI − M(z)).. (1.12). Funkcja ta jest wielomianem stopnia r ze względu na w, a jego współczynniki są funkcjami wymiernymi zmiennej z. Oznaczmy przez wi (z), i = 1, 2, . . . , r pierwiastki funkcji stabilności. Definicja 1.11. Mówimy, że ogólna metoda liniowa (1.2) jest absolutnie stabilna w danym punkcie z ∈ C, jeżeli wszystkie pierwiastki wi (z), i = 1, 2, . . . , r funkcji (1.12) leżą wewnątrz okręgu jednostkowego. Definicja 1.12. Obszarem absolutnej stabilności A metody (1.2) nazywamy zbiór A = {z ∈ C : |wi (z)| < 1, i = 1, 2, . . . , r} .. 19.

(22) 1.4. Przykłady ogólnych metod liniowych. Definicja 1.13. Metoda (1.2) jest A-stabilna, jeżeli jej obszar absolutnej stabilności zawiera lewą półpłaszczyznę zespoloną, to znaczy {z ∈ C : Re(z) < 0} ⊂ A. Poszukiwanie A-stabilnych metod jest nietrywialnym zadaniem - otrzymywanie takich metod zostanie opisane w rozdziale 2 przy pewnych założeniach upraszczających strukturę funkcji stabilności oraz w rozdziale 5. Będziemy szukali także metod otwartych o możliwe dużej powierzchni obszaru absolutnej stabilności. W przypadku metody A-stabilnej nie ma ograniczeń długości kroku związanych ze stabilnością, co jest pożądane przy numerycznym całkowaniu sztywnym układów równań. Silniejszą własnością jest L-stabilność wymagana w celu stłumienia bardzo sztywnych komponentów w równaniu. Definicja 1.14. Ogólna metoda liniowa (1.2) jest L-stabilna, gdy jest A-stabilna oraz lim ρ(M(z)) = 0,. z→∞. gdzie ρ(M(z)) oznacza promień spektralny macierzy M(z). Obok liniowej stabilności rozważa się także nieliniową stabilności ogólnych metod liniowych, w szczególności algebraiczną stabilność. Konstrukcji metod algebraicznie stabilnych poświęcony jest rozdział 4.. 1.4. Przykłady ogólnych metod liniowych Ogólne metody liniowe łączą w sobie cechy zarówno metod wielokrokowych, jak i wielopoziomowych (metody Rungego–Kutty). Przedstawimy teraz wybrane klasyczne metody oraz ich zapis w postaci (1.2) oraz wprowadzimy ważną klasę ogólnych metod liniowych - rodzinę metod DIMSIM.. 1.4.1. Metody Rungego–Kutty s-poziomowa metoda Rungego–Kutty [27, 29, 31, 33, 69, 70] w klasycznej postaci dana jest przez.  s ∑  [n] [n]   Y = y + h aij f (Yj ) n−1   i j=1. s ∑  [n]   y = y + h bj f (Yj )  n n−1 . .. j=1. 20.

(23) 1.4. Przykłady ogólnych metod liniowych. W postaci macierzy Butchera zapisywana jest jako. c A. =. bT. c1 a11 . . . a1s .. .. . . . . .. . . cs as1 . . .. ass. b1. bs. .... .. W postaci ogólnej metody liniowej przyjmuje ona postać (z r = 1 i wektorem c) .   . A B. a11 . . . a1s .. . . . . .. ..    U =   a V  s1  . b1. .... ass. .... bs. 1 .. ..     . 1   . 1. Poziomy metody Rungego–Kutty odpowiadają tutaj poziomom ogólnej metody liniowej.. 1.4.2. Dwukrokowe metody Rungego-Kutty Dwukrokowe metody Rungego–Kutty (TSRK) zaproponowane zostały przez Jackiewicza i Tracognę w [80, 81, 99] jako uogólnienie schematów Rungego–Kutty i oparte są na wartościach aproksymacji rozwiązania w bieżącym oraz poprzednim kroku. Mają one postać   [n]   Y   i. = (1 − ui )yn−1 + ui yn−2 + h. s ( ∑. [n]. [n−1]. aij f (Yj ) + bij f (Yj. j=1.      yn. = (1 − ϑ)yn−1 + ϑyn−2 + h. s ( ∑. [n]. [n−1]. vj f (Yj ) + wj f (Yj. ). ) ,. ). ) .. j=1 [n]. W postaci macierzowej przedstawiane są jako wektor c = [c1 , . . . , cs ]T z Yi. będącym. aproksymacją y(tn−1 + ci h) oraz u A. B. ϑ vT wT. .. W postaci ogólnej metody liniowej mają one postać [79]        . Y [n] yn yn−1 hf (Y [n] ). . . A e−u u.     T   v =     0  . I. B.   .  1 − ϑ ϑ wT     1 0 0  . 0. 0. 0. (. hf Y [n]. ). yn−1 yn−2. (. hf Y [n−1]. .    ,   ) . 21.

(24) 1.4. Przykłady ogólnych metod liniowych. gdzie e oznacza wektor odpowiedniego wymiaru złożony z samych jedynek. Dwukrokowe metody Rungego–Kutty są szeroko znane w literaturze [8, 9, 11, 48, 53, 54, 55, 61, 79, 82, 83], w szczególności postać Nordsiecka była badana w [10].. 1.4.3. Metody liniowe wielokrokowe Metody liniowe wielokrokowe [27, 31, 33, 69, 70, 71, 88, 97, 98] mają postać yn =. k ∑. αj yn−j + h. j=1. k ∑. βj f (yn−j ) ,. (1.13). j=0. n = k, k + 1, . . . , N . Przyjmując w naturalny sposób Y [n] = yn oraz [. y [n] =. ]T. yn yn−1 . . .. y [n−1] =. yn−k+1 hf (yn ) hf (yn−1 ) . . .. hf (yn−k+1 ). [. ,. ]T. yn−1 yn−2 . . .. yn−k hf (yn−1 ) hf (yn−2 ) . . .. hf (yn−k ). ,. w postaci GLM dostajemy [22] .  . A B.             U =   V          . . β0 α1 . . .. αk−1 αk β1 . . .. β0 α1 . . .. αk−1 αk β1 . . .. βk−1 βk   βk−1 βk  . 0 .. .. 1 .. .. ... .. .. 0 .. .. 0 .. .. 0 .. .. ... .. .. 0 .. .. 0 .. .. 0. 0. .... 1. 0. 0. .... 0. 0. 1. 0. .... 0. 0. 0. .... 0. 0. 0 .. .. 0 .. .. ... .. .. 0 .. .. 0 .. .. 1 .. .. ... .. .. 0 .. .. 0 .. .. 0. 0. .... 0. 0. 0. .... 0. 0.         .           . (1.14). Zapisując (1.13) w postaci yn = hβ0 f (yn ) +. k ∑. (αj yn−j + hβj f (yn−j )). j=1. i przyjmując [n−1]. yi. =. k ∑. (αj yn+k−i−j + hβj f (yn+k−i−j )) ,. j=k−i+1. 22.

(25) 1.4. Przykłady ogólnych metod liniowych. j = 1, 2, . . . , k dostajemy następującą reprezentacje metody wielokrokowej w postaci ogólnej metody liniowej [36] . β0.    αk β0 + βk    αk−1 β0 + βk−1    αk−2 β0 + βk−2   ..  .     α2 β0 + β2 . 1. 0. .... 0. 0 0. 0. .... 0. 1 0. 0. .... 0 αk−1  . . 1 0 . . . 0 αk−2 .. . . . . .. .. . . . . . .. . 0 α2 0 0 0 0 0. 0. .... .  αk  . 0 .. .. α1 β0 + β1. . 0 0. 1.  .        . (1.15). α1. Przedstawienie (1.15) metody wielokrokowej w postaci ogólnej metody jest nieredukowalne, w przeciwieństwie do (1.14). W [36] udowodniono, że dla nieredukowalnej wielokrokowej metody A-stabilność metody jest równoważna jej algebraicznej stabilności. Przypomnijmy, że ogólna metoda liniowa jest nieredukowalna, gdy nie jest redukowalna. Redukowalność ogólnej metody liniowej oznacza, że równoważną metodę można zapisać w postaci macierzy rzadkiej  s1 s2 r1 r2. s1. s2. r1. r2. r3. A11. 0. U11. 0. U13.    A21    B11    B21 . 0. r3. A22 0 B22 0. U21 U22 V11. 0. .  0   .  V13   . V21 V22 V23   0. 0. V33. gdzie s = s1 + s2 , r = r1 + r2 + r3 oraz s2 + r2 + r3 > 0. Dwie ogólne metody liniowe . . ˜ U ˜ A   ˜ V ˜ B. . oraz. . A U.  . B V. są natomiast równoważne, gdy istnieją taka macierz permutacyjna P oraz taka nieosobliwa macierz Q, że . . . . ˜ U ˜ PT AP PT UP A =  . −1 −1 ˜ ˜ B V Q BQ Q VQ. 1.4.4. Metody DIMSIM Ważną podklasą ogólnych metod liniowych jest klasa DIMSIM (ang. diagonally implicit multistage integration methods), która została zaproponowana przez Butchera 23.

(26) 1.5. Postać Nordsiecka ogólnych metod liniowych. [28] w 1993 roku i dalej badana między innymi w [35, 38, 39, 40, 41, 42, 46, 84, 85, 86, 101], a także w [79]. Klasa ta jest określana przez następujące własności: • Macierz współczynników A jest macierzą trójkątną dolną z tym samym elementem λ 0 na głównej przekątnej. • Macierz współczynników V ma rząd równy jeden, a jej jedyną niezerową wartością własną jest liczba 1. • Rząd metody p, rząd etapów metody q, liczba poziomów s oraz liczba aproksymacji zewnętrznych r są w przybliżeniu równe. Najczęściej p = q lub p = q + 1 oraz r = s lub r = s + 1. W przypadku r = s przyjmuje się, że U = I oraz że macierz V ma postać V = evT , gdzie vT e = 1, e jest wektorem złożonym z samych jedynek, a sama metoda przyjmuje postać. .  . A B. gdzie λ 0, a. s ∑.            U =   V        . . λ a21 .. .. λ .. .. ... .. 1. 0. .... 0 .. .. 1 .. .. ... .. .. 0. 0. .... as1 as2 . . .. λ. b11 b12 . . .. b1s v1 v2 . . .. 0   0   ..  .   .  1  . ,. vs  . (1.16). . b21 b22 . . . b2s v1 v2 . . . vs   .. .. . . .. .. .. . . ..  . . . .  . . . .  bs1 bs2 . . .. bss v1 v2 . . .. vs. . vs = 1.. i=1. Metody DIMSIM można podzielić na cztery typy, podobnie jak to ma miejsce w przypadku GLM. W następnym podrozdziale przedstawimy postać Nordsiecka metod DIMSIM. Uogólnienie tej postaci prowadzi to klasy metod, którą tutaj nazywać będziemy metodami Nordsiecka i które są głównym tematem rozprawy.. 1.5. Postać Nordsiecka ogólnych metod liniowych Zaczniemy od przedstawienia metody DIMSIM o parametrach p = q = r = s danej przez (1.16) w postaci Nordsiecka. Następnie poprzez uogólnienie tej postaci wprowadzimy klasę metod Nordsiecka będącą głównym tematem badań w tej rozprawie. Postać Nordsiecka metody DIMSIM została uzyskana w [34]. Ma ona następującą formę na siatce jednostajnej z krokiem h 24.

(27) 1.5. Postać Nordsiecka ogólnych metod liniowych   Y [n] = h (A ⊗ I ) f (Y [n] ) + (P ⊗ I ) z [n−1] , m m  z [n] = h (G ⊗ I ) f (Y [n] ) + (Q ⊗ I ) z [n−1] , m. (1.17). m. gdzie A ∈ Rs×s ,. Q ∈ Rs×(s+1) ,. G ∈ R(s+1)×s ,. Q ∈ R(s+1)×(s+1) ,. a współrzędne wektora z [n] aproksymują kolejne współrzędne wektora Nordsiecka, to [n]. znaczy zi = hi−1 y (i−1) (tn ), i = 1, 2, . . . , r. Związek między macierzami postaci (1.16) oraz (1.17) jest następujący [34] P = UW, gdzie W zdefiniowana jest przez (1.9), Q = e1 gdzie e1 =. [. ]T. 1 0 .... [. 1 v T q1 . . .. ∈ Rs+1 , a v =. 0. ]. v T qp. ,. (1.18). [. ]T. v1 v2 . . .. vs. jest wierszem macie-. rzy współczynników V w (1.16), G = LC−1 , gdzie kolumny Lk macierzy L dane są przez . k ∑. . ej+1 Lk = (k − 1)!  − Qek+1  , (k − j)! j=0. k = 1, 2, . . . , p,. a macierz C ma postać C=. [. e c c2 · · ·. cs−1. ]. ∈ Rs×s ,. gdzie ci oznacza wektor złożony z i-tych potęg elementów wektora c, a e = [. ]T. 1 1 ... 1. ∈ Rs .. Z równania (1.18) wynika, że macierz Q, poza pierwszym wierszem, składa się z samych zer. Uogólnienie tej rodziny metod polega na dodaniu nad główną przekątną dodatkowych wolnych współczynników do macierzy Q. Daje to nam nowe możliwości, w szczególności umożliwia konstrukcję metod z wbudowaną kwadratową stabilnością, a przez to konstrukcję metod wysokiego rzędu z zadanymi własnościami stabilności [14, 15, 17, 18, 19]. Wcześniejsze rozważane uogólnienia postaci Nordsiecka metod DIMSIM to dodanie do macierzy Q jednego niezerowego elementu w drugim jej wierszu, porównaj [30]. Inne uogólnienie, to założenie w [43], że wszystkie potęgi macierzy Q powstałej przez. 25.

(28) 1.6. Ogólne metody liniowe z drugą pochodną - uogólnienie ogólnych metod liniowych. usunięcie pierwszego wiersza i pierwszej kolumny macierzy Q są wspólnie ograniczone oraz że element q11 = 1. Narzuca to jednak pewne dodatkowe warunki. Szczegółową analizę klasy metod Nordsiecka przeprowadzimy w dalszych rozdziałach rozprawy.. 1.6. Ogólne metody liniowe z drugą pochodną uogólnienie ogólnych metod liniowych Rozwinięciem klasy ogólnych metod liniowych są ogólne metody liniowe z drugą pochodną dla układów równań różniczkowych zwyczajnych pierwszego rzędu postaci (1.1). Metody te wykorzystują obok poziomów opartych na pierwszej pochodnej rozwiązania, czyli wartości f (Y [n] ), także poziomy oparte na drugiej pochodnej, czyli g(Y [n] ), gdzie g(·) = f ′ (·)f (·). Metody te, w skrócie oznaczane jako SGLM lub SD-GLM, zostały zaproponowane przez Butchera i Hojjati’ego w 2005 roku w [37], a następnie były rozwijane między innymi w [1, 2, 3, 4, 5]. Przedstawimy tutaj tylko ogólne sformułowanie metod oraz wybrane wyniki, po szczegóły odsyłamy do wymienionych prac. Na jednostajnej siatce tn = t0 + hn o ustalonym kroku h = (T − t0 )/N dla ustalonego, naturalnego N postać wektorowa SGLM jest następujaca   Y [n] = h(A ⊗ I )f (Y [n] ) + h2 (A ⊗ I )g(Y [n] ) + (U ⊗ I )y [n−1] , m m m  y [n] = h(B ⊗ I )f (Y [n] ) + h2 (B ⊗ I )g(Y [n] ) + (V ⊗ I )y [n−1] , m. m. m. gdzie n = 1, 2, · · · , N . Wymiary macierzy współczynników są następujące A, A ∈ Rs×s , U ∈ Rs×r , B, B ∈ Rr×s oraz V ∈ Rr×r . Metody SGLM możemy zapisać także w postaci macierzy blokowej wymiaru (s + r) × (2s + r) jako.  . A A U B B V.  .. Podobnie jak w przypadku ogólnych metod liniowych możemy wyróżnić 4 typy metod SGLM. Przyjmując . . .  λ   a21 A=  .  .  . .     ,   .  µ   a21 A=  .  .  . . λ .. .. ... .. as1 as2 · · ·. λ. . µ .. .. ... .. as1 as2 · · ·. µ.     ,   . 26.

(29) gdzie λ 0, µ ¬ 0 ze względu na stabilność, metody typu 1 oraz 2 to metody z, odpowiednio, λ = µ = 0 oraz λ > 0 i µ < 0. Dla aij = aij = 0 dostajemy metody typu 3, gdy λ = µ = 0 oraz typu 4, gdy λ > 0 i µ < 0. W [37] zostały sformułowane pojęcia przedzgodności, zgodności oraz stabilności SGLM analogicznie do GLM. Macierz stabilności metod SGLM ma postać (. )−1. ¯ Im − zA − z 2 A ¯ M(z) = V + z(B + z B). U,. a funkcja stabilności to p(w, z) = det(wI − M(z)). W [1] zaproponowano podklasę metod SDIMSIM jako rozwinięcie metod DIMSIM poprzez wykorzystanie drugiej pochodnej rozwiązania. Metody te określają następujące własności: ¯ są macierzami trójkątnymi dolnymi z tymi • Macierze współczynników A oraz A samymi elementami na głównej przekątnej: λ 0 dla macierzy A oraz µ ¬ 0 dla ¯ macierzy A. • Macierz współczynników V ma rząd równy jeden, a jej jedyną niezerową wartością własną jest liczba 1. • Rząd metody p, rząd etapów metody q, liczba poziomów s oraz liczba aproksymacji zewnętrznych r są w przybliżeniu równe. W [1] zostało także udowodnione twierdzenie o reprezentacji macierzy współczynników B. Jest to rozwinięcie twierdzenia Butchera [28] dla metod DIMSIM. Przy implementacji metod SGLM pojawia się dodatkowy kosz związany z obliczeniem wartości funkcji g(·), jednak wartości f ′ (·) wykorzystywane są też przy rozwiązywaniu metodą Newtona nieliniowego równania występującego przy metodach zamkniętych. Wstępne eksperymenty numeryczne wykazują, że metody SGLM cechuje wysoka dokładność [1]..

(30) Rozdział 2. Konstrukcja metod Nordsiecka z wbudowaną kwadratową funkcją stabilności. Metody Nordsiecka stanowią szeroką podklasę ogólnych metod liniowych. Ich duża liczba współczynników umożliwia konstrukcje metod wysokiego rzędu o zadanych własnościach stabilności. Wektor aproksymacji zewnętrznych jest aproksymacją wektora Nordsiecka, co umożliwia efektywną implementacje w algorytmach o zmiennym kroku metody. Z drugiej strony duża liczba współczynników, często powiązanych między sobą w sposób nieliniowy, utrudnia analizę oraz konstrukcję metod za pomocą obliczeń symbolicznych. W celu badania ogólnych metod liniowych wprowadzono pojęcie stabilności Rungego-Kutty (RKS). Funkcja stabilności takich metod ma jeden niezerowy pierwiastek, przez co mają one podobne własności stabilności jak metody Rungego-Kutty. W [45, 49, 79, 102, 103] wprowadzono pojęcie wbudowanej stabilności Rungego-Kutty (IRKS). Są to pewne algebraiczne warunki na macierze współczynników, których spełnienie gwarantuje, że funkcja stabilności ma jeden niezerowy pierwiastek, a metoda ma własności stabilności podobne do metod Rungego–Kutty. Warunki te są jednak. 28.

(31) 2.1. Metody Nordsiecka. mocno ograniczające, w szczególności dla metod jawnych pozostaje tylko jeden wolny współczynnik, który jednocześnie odpowiada za pole powierzchni obszaru absolutnej stabilności i stałą błędu. W tym rozdziale wprowadzimy pojęcie wbudowanej kwadratowej stabilności (IQS). Będziemy badać metody Nordsiecka, których funkcja stabilności ma dwa niezerowe pierwiastki - jest to zatem ogólniejszy przypadek, niż metody z własnością IRKS. Badanie własności stabilności może być w takim przypadku ograniczone do kwadratowej części wielomianu. Następnie wykorzystamy własność IQS do konstrukcji metod jawnych z parametrami p = q = s = r − 1 o maksymalnym obszarze absolutnej stabilności rzędu p = 1, 2, . . . , 8 oraz zamkniętych metod L-stabilnych rzędu p = 1, 2, . . . , 6. Wyniki zaprezentowane w tym rozdziale oparte są na pracach [14, 15]. Zagadnieniami związanymi z implementacją skonstruowanych metod zajmiemy się w rozdziale 3.. 2.1. Metody Nordsiecka W tym rozdziale będziemy zajmować się ogólnymi metodami liniowymi o parametrach p, q, r i s związanych ze sobą w następujący sposób: p = q = s oraz r = s + 1. Wybór takich samych wartości rzędu metody oraz rzędu poziomów metody ma na celu uniknięcie zjawiska redukcji rzędu. Zjawisko to polega na tym, że obserwowany w eksperymentach numerycznych ze sztywnymi układami równań różniczkowych rząd metody może być niższy, niż wynikałoby to z teoretycznych własności. O macierzy współczynników A będziemy zakładali, że ma ona postać trójkątną dolną z tym samym elementem λ 0 na głównej przekątnej, to znaczy . .  λ   a21  A=  ..  . . λ .. .. ... .. as1 as2 . . .. λ.     .   . Taka postać umożliwia efektywną implementacje metody. W przypadku, gdy λ = 0 metoda Nordsiecka jest jawna (otwarta), dla λ > 0 metoda jest uwikłana (zamknięta). Konstrukcją obu typów metod zajmiemy się później w tym rozdziale.. 29.

(32) 2.2. Warunki rzędu i formuły reprezentacyjne. Jak już wspomniano, macierz V ma następującą postać .  1 v12 v13 · · ·   0 0 v23 · · ·   . .. . . . ... V= . .  .    0 . 0. 0. ···. 0. 0. 0. ···. . v1,r   v2,r   ..  .  . . vr−1,r   . 0. Jedynymi wartościami własnymi macierzy V są zatem 1 o krotności równej 1 oraz 0 o krotności r − 1, z czego wynika, że każda metoda tej postaci jest automatycznie zero-stabilna. Przyjmując, że wektory qi , i = 0, 1, . . . , s zdefiniowane przez (1.8) są kolejnymi wektorami bazy kanonicznej Rr , czyli q0 = e 1 , . . . , qp = e r , dostajemy, że wektor y [n] aproksymacji zewnętrznych jest aproksymacją rzędu p = r−1 wektora Nordsiecka y [n] = z(tn , h) + O(hp+1 ), gdzie wektor Nordsiecka z(t, h) ma postać     z(t, h) =     . . y(t) ′. hy (t) .. . hr−1 y (r−1) (t).      ∈ Rrm .   . (2.1). Metody o powyższych własnościach będziemy nazywali metodami Nordsiecka. W rozdziale 5 rozszerzymy wprowadzoną tutaj klasę metod Nordsiecka na przypadek, gdy rząd p metody jest równy liczbie aproksymacji zewnętrznych r = s + 1.. 2.2. Warunki rzędu i formuły reprezentacyjne Konsekwencja wyboru wektorów qi jako kolejnych elementów bazy kanonicznej Rr jest szczególna postać funkcji wektorowej w = w(z) zdefiniowanej przez (1.10), którą tutaj będziemy dla uproszczenia oznaczać jako wektor Z      w(z) = Z =    . . 1 z .. . z r−1.     ,   . 30.

(33) 2.2. Warunki rzędu i formuły reprezentacyjne. z ∈ C. Dostajemy następującą wersje twierdzenia 1.10 (porównaj także [28]). Twierdzenie 2.1. [14] Metoda Nordsiecka ma rząd poziomów równy q oraz rząd p wtedy i tylko wtedy, gdy ecz = zAecz + UZ + O(z p+1 ),. (2.2). ez Z = zBecz + VZ + O(z q+1 ),. (2.3). gdzie ecz = [ ec1 z ec2 z · · ·. ecs z ]T .. Twierdzenie to jest konsekwencją twierdzenia 3.1 z pracy [28], dlatego też dowód pomijamy. Użyjemy teraz wyników powyższego twierdzenia w celu wyrażenia macierzy U oraz B jako funkcji A, V oraz c. W tym celu definiujemy macierze Kp ∈ R(p+1)×(p+1) oraz Ep ∈ R(p+1)×(p+1) przez   0 1   0 0   . . . . Kp =   . .. . 0 ··· 1 .. ..    0 0 0 . . 0   ··· 0   . . ..  . .  , ···. 0 0 0 ···.  1   . 1 1 2!1 · · ·    0 1 1 ···  . Ep = exp(Kp ) =   0  .  .  . . 0. 0 1 ··· .. .. . . . . .. 0 0 0 ···. 1 p! 1 (p−1)! 1 (p−2)!. .. ..       ,     . 1. oraz macierz Cp przez [. Cp =. 2 e c c2! · · ·. cp p!. ]. ∈ Rs×(p+1) ,. gdzie ci oznacza wektor złożony z i-tych potęg elementów wektora c. Macierz Kp nazywana jest macierzą przesunięcia w prawo, gdyż przemnożenie dowolnej macierzy z prawej strony przez macierz Kp powoduje przesunięcie kolumn macierzy w prawo o jedną kolumnę i wypełnienie pierwszej kolumny zerami. Podobnie przemnożenie lewostronne przez macierz Kp powoduje przesunięcie wierszy o jeden w górę i uzupełnienie ostatniego wiersza zerami. Zaczniemy od przedstawienia macierzy U oraz V w terminach pozostałych współczynników metody. Twierdzenie 2.2. [14] Załóżmy, że p = q = s oraz r = s + 1. Wtedy U = Cp − ACp Kp ,. (2.4). V = Ep − BCp Kp .. (2.5). 31.

(34) 2.2. Warunki rzędu i formuły reprezentacyjne. Dowód. Rozwijając ecz w szereg Taylora dostajemy  . ecz. . . p ∑ cj1 z j.  ec1 z  j=0 j!     ..  .. = . = .     p j j  ∑c z cp z p e  j=0 j!.       + O(z p+1 ) = Cp Z + O(z p+1 ).    . (2.6). Ponieważ z · Z = Kp Z + O(z p+1 ), zatem z równania (2.6) wynika, że zecz = Cp Kp Z + O(z p+1 ).. (2.7). Podstawiając (2.6) i (2.7) do (2.2) dostajemy (2.4). Podobnie biorąc pod uwagę, że ez Z = Ep Z + O(z p+1 ) z warunku (2.3) dostajemy (2.5), co kończy dowód. Ze względu na to, że macierz V jest macierzą trójkątną górną, wygodniejsze do bezpośredniego wykorzystania jest odwrócenie relacji (2.5) i wyrażenie macierzy B w terminach V i c. W tym celu podzielimy macierze B, V, Ep i Cp w następujący sposób . . bT , B= ˜ B. . . 1 v , V= ˜ 0 V. . Ep = . 1 eTp−1 0 Ep−1.  ,. [. Cp =. p Cp−1 cp!. ]. ,. gdzie bT jest pierwszym wierszem macierzy B oraz 0 jest wektorem odpowiedniego wymiaru złożonym z samych 0. Twierdzenie 2.3. [14, 79] Załóżmy, że p = q = s, r = s + 1 oraz że elementy wektora c są parami różne. Wtedy bT = (eTp−1 − v)C−1 p−1 ,. (2.8). ˜ = (Ep−1 − V)C ˜ −1 . B p−1. (2.9). Dowód. Z równania (2.5) oraz relacji Cp Kp = [0|Cp−1 ] wynika, że . . . . . . 1 eTp−1 1 v 0 bT Cp−1 = .  − ˜ ˜ 0 V 0 Ep−1 0 BCp−1 Porównując teraz odpowiednie elementy powyższego równania macierzowego oraz z założenia o odwracalności macierzy Cp−1 dostajemy (2.8) i (2.9).. 32.

(35) 2.3. Wbudowana kwadratowa stabilność. 2.3. Wbudowana kwadratowa stabilność Własności stabilności ogólnej metody liniowej określone są przez jej funkcję stabilności (1.12). W przypadku, gdy funkcja stabilności ma dwa niezerowe pierwiastki, badanie stabilności możemy ograniczyć do badania kwadratowego czynnika wielomianu stabilności. Funkcja stabilności p(w, z) dana przez (1.12) metody Nordsiecka jest wielomianem zmiennej w o współczynnikach będących funkcjami wymiernymi zmiennej z. W celu ułatwienia analizy przemnażamy funkcję p(w, z) przez czynnik (1 − λz)s , a otrzymany wielomian stabilności oznaczamy zazwyczaj tym samym symbolem p(w, z). Ma on postać p(w, z) = (1 − λz)s ws+1 − ps (z)ws + · · · + (−1)s p1 (z)w + (−1)s+1 p0 (z), gdzie ps (z) = 1 + ps1 z + · · · + pss z s , ps−1 (z) = ps−1,1 z + · · · + ps−1,s z s , .. . p0 (z) = p01 z + · · · + p0s z s . Bezpośrednie podejście poprzez wygenerowanie wielomianu stabilności za pomocą obliczeń symbolicznych staje się niemożliwe, gdy s i r przekraczają 4. Możliwe jest jednak numeryczne wygenerowanie wielomianu stabilności dla danej ogólnej metody liniowej przy pomocy algorytmu opartego na wariancie szeregów Fouriera [41] nawet dla metod o dużych wartościach parametrów s oraz r. W tym podrozdziale omówimy narzędzie, jakim jest wbudowana kwadratowa stabilność, oznaczana w skrócie jako IQS. Idea ta została zaproponowana najpierw dla dwukrokowych metod Rungego–Kutty w [54], a następnie przeniesiona na metody Nordsiecka [14, 15, 51].. 2.3.1. Charakteryzacja metod z wbudowaną kwadratową stabilnością W celu badania metod z kwadratową funkcją stabilności, czyli metod których wielomian stabilności posiada dwa niezerowe pierwiastki, wprowadzamy relację równoważności między macierzami tego samego wymiaru. Mówimy zatem, że dwie macierze D, E ∈ Rj×k są równoważne, co oznaczamy przez D ≡ E, wtedy i tylko wtedy, gdy 33.

(36) 2.3. Wbudowana kwadratowa stabilność. są sobie równe za wyjątkiem co najwyżej dwóch pierwszych wierszy. Relacja ta jest zachowywana przy prawostronnym mnożeniu przez dowolną macierz. Oznacza to, że jeżeli F ∈ Rk×l , to z faktu D ≡ E wynika równoważność DF ≡ EF. Mamy następującą definicję. Definicja 2.4 ([14, 54, 79] ). Metoda Nordsiecka o macierzach współczynników A, U, B, V oraz wektorze c posiada własność wbudowanej kwadratowej stabilności (oznaczane przez IQS, ang. inherent quadratic stability), jeżeli istnieje taka macierz X ∈ Rr×r , że BA ≡ XB. (2.10). BU ≡ XV − VX.. (2.11). Warunki (2.10) oraz (2.11) będziemy nazywali warunkami kwadratowej stabilności lub krócej - warunkami IQS. Znaczenie tej własności wynika z następującego twierdzenia. Twierdzenie 2.5. [14] Załóżmy, że metoda Nordsiecka posiada własność IQS. Wtedy jej wielomian stabilności p(w, z) przyjmuje kwadratową postać p(w, z) = ws−1 (w2 − p1 (z)w + p0 (z)),. (2.12). gdzie p1 (z) i p0 (z) są funkcjami wymiernymi zmiennej z (metody zamknięte) lub wielomianami zmiennej z (metody otwarte). Dowód. Zauważmy, że (2.10) jest równoważne B(I − zA) ≡ (I − zX)B. Przy założeniu, że I − zA jest nieosobliwa, otrzymujemy B ≡ (I − zX)B(I − zA)−1 .. (2.13). f Zamiast macierzy stabilności M(z) rozważamy macierz M(z) podobną do M(z) i. zdefiniowaną niżej. Korzystając z (2.11) i (2.13) oraz przy założeniu nieosobliwości I − zX dostajemy, że f M(z) = (I − zX)M(z)(I − zX)−1. = (I − zX)(V + zB(I − zA)−1 U)(I − zX)−1 = (V − zXV + z(I − zX)B(I − zA)−1 U)(I − zX)−1 ≡ (V − zXV + zBU)(I − zX)−1 ≡ (V − zXV + z(XV − VX))(I − zX)−1 = (V − zVX)(I − zX)−1 . 34.

(37) 2.3. Wbudowana kwadratowa stabilność. Zatem (I − zX)M(z)(I − zX)−1 ≡ V. Powyższą relację zapisujemy w postaci blokowej . (I − zX)M(z)(I − zX)−1. f (z) M 11    0 0   . . . . =  . .    . . f (z) M 12. 0 v34 . . . v3,r .. .. . . .. . . . .. 0 0. 0. 0. .... vs,r. 0 0. 0. 0. .... 0.      ,     . f (z) ∈ R2×2 i M f (z) ∈ R2×(s−1) . Funkcja stabilności p(w, z) o macierzach gdzie M 11 12 f stabilności M(z) lub M(z) przyjmuje zatem postać (2.12), co kończy dowód.. 2.3.2. Postać macierzy X Następujące twierdzenie dotyczy struktury macierzy X występującej w definicji 2.4. Twierdzenie 2.6. [15] Najogólniejsza postać macierzy X dla metody Nordsiecka o parametrach p = q = s = r − 1 i spełniająca warunki (2.10) oraz (2.11) to . x x1, x13 · · ·  11   x2,1 x2,2 x2,3 · · · .    X=      . . x1,r−1 x1,r   x2,r−1 x2,r   . 0. x3,r  . ··· .. .. 0 .. .. x4,r .. .. ···. 1. xr,r. 0. 1. 0. ···. 0 .. .. 0 .. .. 1 .. .. 0. 0. 0. .      . (2.14). Dowód. Mnożąc lewostronnie równanie (2.2) przez zB oraz równanie (2.3) przez I−zX, a następnie dodając do siebie otrzymane równości dostajemy ez (I − zX)Z = z 2 (BA − XB)ecz + VZ + z(BU − XV)Z + O(z p+1 ). Korzystając teraz z warunków IQS (2.10) i (2.11) dostajemy ez (I − zX)Z ≡ V(I − zX)Z + O(z p+1 ), a dalej (ez I − V)(I − zX)Z ≡ O(z p+1 ).. (2.15). 35.

(38) 2.3. Wbudowana kwadratowa stabilność. Podzielmy macierz V oraz wektor Z˜ := (I − zX)Z w następujący sposób . V1,1 V1,2. V= gdzie. 0. . V1,1 = . V2,2. . . ,. . 1 v12 0. 0. ,. . Z˜1 , Z˜ =  ˜ Z2 . V1,2 = . . v13 . . .. v1,r. v23 . . .. v2,r. ,. 0 jest macierzą zerową wymiaru (r − 2) × 2, a Z˜1 składa się z dwóch pierwszych ˜ elementów Z. Z równania (2.15) dostajemy . . (ez I2 − V1,1 )Z˜1 − V1,2 Z˜2  ≡ O(z p+1 ), (ez I − V)Z˜ =  z ˜ (e Ir−2 − V2,2 )Z2 gdzie Ik jest macierzą jednostkową wymiaru k. Macierz ez Ir−2 − V2,2 jest odwracalna dla |z| < ε dla pewnego ε > 0 (porównaj [51]). Zdefiniujmy macierz F jako . F= Mamy. I2. . 0 −1. 0 (e Ir−2 − V2,2 ) z. .. . . (ez I2 − V1,1 )Z˜1 − V1,2 Z˜2 z ˜   ≡ O(z p+1 ), F(e I − V)Z = Z˜2 co jest równoważne temu, że Z˜ = (I − zX)Z ≡ O(z p+1 ). Zauważmy, że zJZ ≡ Z, gdzie J = KTp jest macierzą przesunięcia w lewo przy mnożeniu lewostronnym lub w dół przy mnożeniu prawostronnym, zatem z(J − X)Z ≡ O(z p+1 ), co jest równoważne (J − X)Z ≡ O(z p ). W związku z powyższym macierz X − J musi być równa 0 za wyjątkiem 2 pierwszych wierszy oraz ostatniej kolumny, co należało pokazać.. 36.

(39) 2.4. Konstrukcja metod zamkniętych. 2.4. Konstrukcja metod zamkniętych W tym podrozdziale prezentujemy konstrukcje uwikłanych A- i L- stabilnych metod Nordsiecka posiadających własność IQS. Głównym narzędziem wykorzystywanym przy konstrukcji metod, obok IQS, jest kryterium (twierdzenie) Schura [96]. Dla metod o rzędzie s = p 3 będziemy zakładali, że wektor c jest z góry zadany. Metody rzędu p = 1, 2, 3 i 4 zostały skonstruowane w [14], metody rzędu p = 5 oraz p = 6 są nowymi wynikami. Konstrukcja metod uwikłanych jest technicznie trudniejsza ze względu na to, że macierz współczynników A ma niezerowy element λ > 0 na głównej przekątnej. Zaczynamy od wykorzystania warunków rzędu metody (2.4), (2.8) oraz (2.9). Układ równań IQS dany równaniami (2.10) oraz (2.11), przy uwzględnieniu postaci macierzy X danej przez (2.14), rozwiązujemy ze względu na wszystkie pozadiagonalne współczynniki macierzy A, współczynniki macierzy V z wierszy od numeru 2 do r − 1 oraz współczynniki xir , i = 3, 4, . . . , r macierzy X. Na mocy twierdzenia 2.5 wielomian charakterystyczny p(w, z) ma postać (2.12), którą tutaj zapiszemy w postaci przemnożonej przez czynnik (1 − λz)s p(w, z) = ws−1 ((1 − λz)s w2 − p1 (z)w + p0 (z)), gdzie p1 (z) i p0 (z) są teraz wielomianami zmiennej z: p1 (z) = 1 + p11 z + · · · + p1s z s oraz p0 (z) = p01 z + · · · + p0s z s . Na tym etapie konstrukcji pozostaje nam r wolnych współczynników: λ oraz v12 , . . . , v1r . Wykorzystujemy je do konstrukcji metod A-stabilnych oraz, o ile to możliwe, L-stabilnych. W celu zagwarantowania L-stabilności (przy założeniu, że możliwe będzie zapewnienie A-stabilności), przyrównujemy do zera współczynniki przy najwyższych potęgach wielomianów p1 oraz p0 , to znaczy rozwiązujemy układ równań p1s = 0,. p0s = 0.. (2.16). Układ ten możemy rozwiązać ze względu na współczynniki v1,r−1 oraz v1,r (poza przypadkiem s = 1, r = 2). Zredukowanie liczby parametrów ułatwia zastosowanie twierdzenia Schura [96] do kwadratowej części wielomianu stabilności p(w, z).. 37.

(40) 2.4. Konstrukcja metod zamkniętych. 2.4.1. Konstrukcja metod zamkniętych o parametrach p = q = s = 1 oraz r = 2 Używając formuł (2.4), (2.8) oraz (2.9) dla przypadku p = q = s = 1, r = 2 dostajemy trzy parametrową rodzinę metod zależnych od c1 , λ i v12 . Współczynniki tej metody mają postać  . λ.   = 1−v 12   V. A U B. . . 1. 1 c1 − λ 1. v12. 0. 0.    ,  . a wielomian stabilności p(w, z) = (1 − λz)w2 − p1 (z)w + p0 (z), gdzie p1 (z) = 1 + (1 + c1 − 2λ − v12 )z, p0 (z) = (c1 − λ − v12 )z. Wielomian ten ma stopień 2 ze względu na w, w związku z czym warunki IQS nie są potrzebne. Biorąc c1 = 1, λ = 1, v12 = 0 dostajemy metodę Nordsiecka równoważną zamkniętemu schematowi Eulera. Metoda dana jest przez  . . 1 1 0.   = 1   V. A U B. . .   1 0  . . 1 0 0 Wielomian stabilności tej metody ma postać p(w, z) = (1 − z)w2 − w. Metoda ta jest zatem A- i L-stabilna oraz, dodatkowo, jest algebraicznie stabilna (porównaj [17] oraz rozdział 4).. 2.4.2. Konstrukcja metod zamkniętych o parametrach p = q = s = 2 oraz r = 3 Przypadek p = q = s = 2 jest pierwszym nietrywialnym przypadkiem wśród metod uwikłanych, gdzie wykorzystamy warunki IQS do otrzymania wielomianu stabilności o 38.