1 Wykład no 9

2 Wyznaczanie obszaru stabilności bezwzględnejPrzekształćmy równanie charakterystyczne do postaci:

3 i stawiamy pytanie odwrotnie:Jeżeli z leży wewnątrz lub na kole jednostkowym to gdzie będzie leżało  płaszczyzna  płaszczyzna z Im() iy x Re()

4 Przykłady metoda jawna Adamsa - Bashforthaschemat Eulera Przyjmując, że z leży na okręgu jednostkowym i podstawiając do równania:

5 dla metody Eulera otrzymujemy:Wnętrze koła i dla rzeczywistych  mamy:

6 Algorytm Adamsa – Bashfortha II-go rzędui po odpowiednich przekształceniach mamy:

7 Odwzorowanie koła jednostkowego jestJeżeli  rzeczywiste

8 Dla rzędu III-go Dla rzeczywistych 

9 Dla czwartego rzędu h jeszcze mniejsze niż dla metody rzędu trzeciego

10 Metody niejawne Adamsa – Moultonametoda niejawna Eulera Obszar stabilności określa równanie: Jest to równanie okręgu o środku w punkcie –1 i promieniu 1. Obszar stabilności na zewnątrz tego okręgu.

11

12 niejawna metoda Eulera:Oznacza to, że niejawna metoda Eulera: jest stabilna dla dowolnej długości kroku h. Oczywiście nie możemy zapomnieć o fizyce procesu opisywanego równaniem czy układem równań różniczkowych

13 Algorytm II-go rzędu (algorytm trapezów)Dla odwzorowania koła jednostkowego mamy: Badamy jak odwzorowuje się koło czyli

14

15 Wnętrze koła odwzorowuje się na prawą półpłaszczyznę co oznacza, że jeżeli rozwiązanie równania jest stabilne, czyli Re()  0, to krok całkowania h może być dowolny. Zawsze nie możemy zapomnieć o fizyce procesu opisywanego równaniem czy układem równań różniczkowych i to jest główne ograniczenie wielkości kroku.

16 Algorytm III-go rzędu Adamsa - Moultonai dla koła jednostkowego czyli

17 dla rzeczywistych  mamy warunek:

18 Algorytm VI-go rzędu Adamsa - Moultonai mamy odwzorowanie:

19

20 i Zbierając wyniki: Metody jawne Adamsa – Bashfortha: i algorytm niejawny Adamsa –Moultona: h – dowolny h  – dowolny h  < 6 h  < 3

21 Otrzymane wyniki pokazują wyższość metod niejawnychnad metodami jawnymi. W metodzie predyktor – korektor, gdzie metoda jawna służy tylko i wyłącznie do otrzymania zerowego przybliżenia rozwiązania równania metody niejawnej, o stabilności decyduje tylko metoda niejawna. Jak widać również z podanych ocen z punktu widzenia stabilności zbyt wysoki rząd metody nie jest korzystny

22 na kondensatorze jest:Równania sztywne Dany jest obwód elektryczny: Równanie różniczkowe dla napięcia uC na kondensatorze jest: uC Warunki początkowe są: i Równanie charakterystyczne: które ma pierwiastki: i

23 r2=-10000 r1=-10

24 Zapiszmy równanie: w postaci normalnej: lub podstawiając dane: warunki początkowe: Wybierzmy metodę jawną Eulera, krok h=10-5

25

26 stała czasowa 0.1 więc możnaPonieważ dla czasów zanikła składowa u2(t)~exp(-10000t) weźmy dokładne wartości startowe w punkcie t=0.005 i(t=0.005)= uC(t=0.005)= stała czasowa 0.1 więc można przyjąć krok h=0.01

27 krok h=

28 krytyczny krok hkr<0.0002start w punkcie t=0.005 z dokładnych wartości początkowych

29 Rozważamy układ n równań różniczkowych:i=1,2,...,n i może być liczbą zespoloną postaci: Jedno z rozwiązań układu równań różniczkowych będzie postaci: gdzie Ci jest stałą całkowania. Przypadek 1. ponieważ =h, więc czyli jest to prawa półpłaszczyzna

30 rozwiązanie powinno być i dokładne i stabilne, bo składowa przejściowa Im liczba  dobrana tak, że po jednym kroku praktycznie: III stabilny W obszarze:  Re rozwiązanie powinno być i dokładne i stabilne, bo składowa przejściowa nie znikła i mamy oscylacje Dla uzyskania dokładności w fazie początkowej powinno być N<1/8. Oznaczając =ih Liczba oscylacji:

31 rozwiązanie numeryczneIm mamy: czyli III  stabilny II W obszarze II rozwiązanie numeryczne musi być: dokładne i stabilne  Re - Przypadek 2. jest narastające i można liczyć tylko odpowiednio małym krokiem Rozwiązanie

32 dokładny i względnie stabilny III Im czyli algorytm powinien być dokładny i względnie stabilny III  stabilny II I Algorytmy spełniające warunki I, II, III nazywamy sztywno stabilnymi -  Re - Twierdzenie Dahlquista: Algorytm wielokrokowy, stabilny bezwzględnie w obszarze nie może być rzędu wyższego niż 2. Najlepszy jest algorytm trapezów.

33 Trapezy prawa półpłaszczyzna =0

34 =0 trzeci rząd dla rzeczywistych  mamy warunek:

35 =0 czwarty rząd

36 Algorytmy sztywno stabilne GearaPierwszego

37 drugiego:

38 trzeciego:

39 czwartego:

40 Metoda Runge - Kutty Równanie Rozwiązujemy stosując szereg Taylora ale

41 czyli

42 ale a Podstawiając do i porządkując mamy

43 Metoda Runge -Kutty

44 Sposób wyznaczania współczynników na przykładziemetody drugiego rzędu (p=2): Drugi składnik rozwijamy w szereg Taylora w otoczeniu punktu xn, tn Podstawiając i porządkując mamy:

45 a porównując z szeregiem Tayloraprzy tych samych potęgach h otrzymujemy:

46 lub w1=w2=w i rozwiązując otrzymujemy: Przyjmując w2=1 mamy: w1=0, b21=1/2 i stąd algorytm: lub w1=w2=w i rozwiązując otrzymujemy: w=0.5, b21=1 i stąd inny algorytm:

47 Przykład: Dany jest dławik o charakterystce: Rezystancja dławika wynosi 0.5. Obliczyć prąd płynący w obwodzie zasilanym sem e(t)=100sin314t. Schemat obwodu możemy przyjąć w postaci:

48 Suma spadków napięć pozwala zapisać równanie:Biorąc pod uwagę krzywą magnesowania: Podstawiając do równania obwodu i porządkując:

49 Warunek początkowy jest i0=i(t=0)=0.Wybór kroku całkowania: Stała czasowa liniowej części obwodu wynosi 0.1/0.5=0.2s. Krok czasowy można przyjąć 0.2/10=20ms. Okres wymuszenia T=20ms krok należy przyjąć rzędu T/20=1ms. Prawdopodobnie będzie trzecia harmoniczna więc przyjmujemy krok h=0.2ms.

50 Obliczenia metodą Runge – Kutty według schematu:x=i; Start: i(t=0)=i0=0 h=0.0002

51 i mamy: t=h=0.0002 Metoda Runge – Kutty pozwala zmienić krok na każdym etapie. Zwiększamy krok dwukrotnie. h=0.0004

52 i2= K2 i2=

53 Jak ocenić czy wolno zmienić długość kroku?Czy zmniejszyć czy zwiększyć? Ocena błędu metodą Rungego: Dla metody rzędu p-go mamy:

54 stąd ocena błędu: Znając ocenę błędu można poprawić rozwiązanie podstawiając do

55 lub dokładniej z równania:

56 W obliczanym przypadku musimy powtórzyćobliczenia z krokiem i mamy dla t=0.0004: K1= K2=0.0188 i1+1/2= Dla t= mamy: K1= K2= i1+2*1/2=

57 i1+2*1/2= i2= Obliczone z krokiem h= było: W tym przypadku p=2 i ze wzoru: mamy oceną błędu:

58 Rozwiązanie poprawione ze wzoru:Na wykonanie jednego kroku należało policzyć funkcję f(in,tn) 2 – h=0.0004 1+2 – h=0.0002 razem 5 - razy

59 Metoda IV –go rzędu Przy ocenie dokładności obliczeń metodą Rungego wymaga 11-krotnego obliczenia f(x,t).

60 Metoda Mersona

61 tylko 5-cio krotne obliczanie f(x,t).Przykład Równanie wahadła: Niech =1s-2 Warunki początkowe: około 86°

62 Sprowadzamy do układu równań I-go rzęduWarunki początkowe: Obliczenia chcemy prowadzić z dokładnością 0.001 Startujemy z krokiem h=0.1. Krok wybrano jako 0.1 okresu wahadła liniowego.

63

64 Błąd:

65 Dokładność założona została osiągnięta.W następnym kroku można zwiększyć krok. Rozwiązanie w chwili t=0.1 i do następnego kroku możemy wystartować z nową wartością kroku h

66 Metody włożone lub Metody Fehelberga – Runge -Kutty Stosujemy metodę Runge – Kutty rzędu p i rzędu p+1 i aby zmniejszyć liczbę obliczanych współczynników wybieramy je tak, że w obu metodach jest pierwszych p współczynników K jednakowe, czyli i=2,3,..,p+1

67 i mamy dla metody rzędu p-goa dla metody rzędu (p+1)-go Ocenę błędu można zrobić stosunkowo prosto

68 Po odjęciu stronami otrzymujemy:gdzie

69 Znając błąd możemy postępować jak w metodzieMersona i rozwiązanie przyjmować z dokładniejszej metody rzędu p+1. Najczęściej stosowana metoda RKF45 ma współczynniki

70 Błąd

71 Rozwiązanie wykorzystując metodę dokładniejszą jestMetoda gwarantuje obliczenia z błędem rzędu h4.