• Nie Znaleziono Wyników

Ćwiczenie 1b: Schematy jawne i niejawne dla równań różniczkowych zwyczajnych ∗

N/A
N/A
Protected

Academic year: 2021

Share "Ćwiczenie 1b: Schematy jawne i niejawne dla równań różniczkowych zwyczajnych ∗"

Copied!
8
0
0

Pełen tekst

(1)

Ćwiczenie 1b: Schematy jawne i niejawne dla równań różniczkowych zwyczajnych

12 kwietnia 2016

Równanie różniczkowe o ogólnej postaci du

dt = f (t, u) (1)

rozwiążemy prostymi schematami różnicowymi i porównamy dokładność uzy- skanych rozwiązań.

Zadanie 1

Należy znaleźć numerycznie rozwiązanie rówanania:

dy

dt = f (t) = −2sin(t)

y − t + 1 (2)

przy użyciu metod

a) jawnej Eulera (10pkt) b) niejawnej Eulera (10pkt)

c) trapezów (10pkt)

d) Rungego-Kutty 2 rzędu (10pkt) e) Rungego-Kutty 4 rzędu (10pkt)

Rozwiązanie należy znaleźć dla następujących parametrów: t ∈ [0, 5 · π], ∆t = 0.1, 0.01 i warunku poczatkowego y(0) = 25 4 . Rozwiązanie dokładne:

y(t) = (

cos(t) + 3 2

) 2

+ t (3)

Laboratorium z inżynierskich metod numerycznych, Wydział Fizyki i Informatyki Sto-

sowanej AGH 2014/2015. Bartłomiej Szafran (bszafran@agh.edu.pl), Krzysztof Kolasiński

(kolasinski@fis.agh.edu.pl), Elżbieta Wach (Elzbieta.Wach@fis.agh.edu.pl), Dariusz Żebrowski

(Dariusz.Zebrowski@fis.agh.edu.pl)

(2)

Jako wyniki należy przedstawić wykresy: a) rozwiązania numerycznego (y num (t)) i b) błędu metody liczonego jako δ(t) = y dok (t) − y num (t) dla każdej z metod.

Dla danej wartości ∆t należy: i) zrobić jeden wykres zbiorczy przedstwiający uzyskane rozwiązania numeryczne (wszystkie metody), ii) jeden wykres błędu dla jawnej i niejawnej metody Eulera (rząd zbieżności p = 1), iii) jeden wykres błędu dla metody trapezów i RK2 (ten sam rząd zbieżności p = 2) oraz iv) oddzielny wykres błędu dla RK4 (rząd zbieżności p = 4). W sumie będzie 8 wykresów.

Uwaga 1: wszystkie metody najlepiej zaimplementować w tej samej pętli i zapisywać wyniki (rozwiązania + błąd) do jednego pliku (dla określonego ∆t) - pozwoli to na łatwą kontrolę/podgląd uzyskanych danych.

Uwaga 2: wszystkie wykresy można przedstawić na jednym rysunku (1 plik zamiast 8) korzystając z opcji multiplot w Gnuplocie. Szkielet skryptu reali- zującego takie zadanie (wykresy są umieszczane kolejno od lewej do prawej ):

set term postscript color enhanced solid size 20cm,40cm font 12 set out ’z1.eps’

set multiplot layout 4,2 rowsfirst

plot ... #rozwiazania dla dt=0.1 plot ... #rozwiazania dla dt=0.01 plot ... #błąd metod Eulera dt=0.1 plot ... #błąd metod Eulera dt=0.01

plot ... #błąd metod trapezów i RK2 dla dt=0.1 plot ... #błąd metod trapezów i RK2 dla dt=0.01 plot ... #błąd metody RK4 dla dt=0.1

plot ... #błąd metody RK4 dla dt=0.01 unset multiplot

jawna metoda Eulera

Po dyskretyzacji zmiennej czasowej t n = n · ∆t, n = 0, 1, 2, . . . , t max /∆t ak- tualne rozwiązanie oznaczymy jako u n a rozwiązanie w chwili następnej (t n+1 ) przez u n+1 . Do znalezienia u n+1 wykorzystujemy proste podstawienie (metoda jawna):

u n+1 = u n + ∆tf (t n , u n ) (4) Powyższy wzór stosujemy aż osiagniemy warunek t n = t max .

niejawna metoda Eulera

Do znalezienia u n+1 wykorzystujemy formułę:

u n+1 = u n + ∆tf (t n+1 , u n+1 ) (5)

w której wartość f (t n+1 , u n+1 ) jest nieznana (dlatego metoda niejawna). Aby

znaleźć u n+1 wykorzystujemy metodę Newtona poszukiwania pierwiastków rów-

(3)

Wykorzystujemy wzór iteracyjny:

u k+1 n+1 = u k n+1 F (u k n+1 )

∂F (u

kn+1

)

∂u

kn+1

= u k n+1 u k n+1 − u n − ∆t [

−2sin(t n+1 )

u k n+1 − t n+1 + 1 ]

1 + ∆tsin(t

n+1

)

u

kn+1

−t

n+1

(7) Wzór (7) iterujemy K = 20 razy startując od przybliżenia u 1 n+1 = u n , a ostatnie przybliżenie zachowujemy do dalszych obliczeń u n+1 = u K n+1 .

metoda trapezów

W metodzie trapezów rozwiązanie w kolejnym kroku wyznaczamy według wzo- ru:

u n+1 = u n + ∆t

2 [f (t n , u n ) + f (t n+1 , u n+1 )] (8) Ponieważ wartości f (t n , u n+1 ) nie znamy (metoda niejawna) więc u n+1 znaj- dujemy w identyczny sposób jak w metodzie niejawnej Eulera. Rozwiązanie w chwili t n+1 znajdujemy korzystając z metody Newtona. Czyli najpierw definiu- jemy równanie nieliniowe:

F (u n+1 ) = u n+1 − u n ∆t

2 [f (t n , u n ) + f (t n+1 , u n+1 )] (9) a następnie znajdujemy jego miejsce zerowe korzystając z formuły iteracyjnej:

u k+1 n+1 = u k n+1 u k n+1 − u n ∆t 2 [

( −2sin(t n )

u n − t n + 1) + ( −2sin(t n+1 )

u k n+1 − t n+1 + 1) ]

1 + ∆t 2sin(t

n+1

)

u

kn+1

−t

n+1

(10)

metoda Rungego-Kutty 2 rzędu (jawna)

W jawnej metodzie RK2 rozwiązanie w chwili t n+1 wyznaczamy korzystając z dwóch etapów pośrednich:

u n+1 = u n + ∆t

2 [k 1 + k 2 ] (11)

gdzie wielkości k 1 i k 2 są zdefiniowane w następujący sposób:

k 1 = f (t n , u n ) (12)

k 2 = f (t n+1 , u n + ∆tk 1 ) (13)

(14)

(4)

metoda Rungego-Kutty 4 rzędu (jawna)

W jawnej metodzie RK4 rozwiązanie w chwili t n+1 wyznaczamy korzystając z czterech etapów pośrednich:

u n+1 = u n + ∆t

6 [k 1 + 2k 2 + 2k 3 + k 4 ] (15) gdzie wielkości k i są zdefiniowane w następujący sposób:

k 1 = f (t n , u n ) (16)

k 2 = f (t n+1/2 , u n + ∆t

2 k 1 ) (17)

k 3 = f (t n+1/2 , u n + ∆t

2 k 2 ) (18)

k 4 = f (t n+1 , u n + ∆tk 3 ) (19) (20) oraz t n+1/2 = t n + ∆t 2

Zadanie 2 - problem sztywny, ekstrapolacja Ri- chardsona

Należy rozwiązać numerycznie problem oscylatora van der Pola tj. znaleźć roz- wiązanie jego równania ruchu:

d 2 u

dt 2 − λ(1 − u 2 ) du

dt + u = 0 (21)

Powyższe RRZ 2 rzędu zamieniamy na układ 2 RRZ 1 rzędu:

du

dt = f (t, u, v) = v (22)

dv

dt = g(t, u, v) = λ(1 − u 2 )v − u (23) Ponieważ dla λ = 100 problem jest sztywny (czyli mamy do czynienia z róż- nymi skalami czasowymi) zastosujemy automatyczną kontrolę kroku czasowego.

Sposób postępowania jest następujący. Jeśli znamy rozwiązanie w chwili t n to rozwiązanie w chwili t n+2 możemy uzyskać na dwa sposoby. Pierwszy sposób to oczywiście wykonanie jednego długiego kroku o długości 2∆t które daje nam rozwiązanie u (1) n+2 różnące od dokładnego o pewien błąd lokalny:

u dok (t n+2 ) = u (1) n+2 + C(2∆t) N + O(∆t N +1 ) (24)

Możemy też wykonać dwa krótsze kroki co ∆t i otrzymać lepsze rozwiązanie

u (2) co da nam inny błąd lokalny:

(5)

Stałą błędu C wyznaczamy odejmując od (25) rozwiązanie (24). Następnie mo- żemy określić błąd rozwiązania uzyskanego w dwóch krokach:

E = 2C∆t n = u (2) n+2 − u (1) n+2

2 N −1 − 1 (26)

i wykorzystać go do modyfikacji kroku czasowego:

∆t nowy =

( S · tol

|E|

) 1/N

∆t (27)

N oznacza błąd lokalny metody.

Algorytm numerycznego rozwiązywania równania różniczkowego z doborem kro- ku czasowego:

1. dysponując rozwiązaniem u n znajdujemy u (1) n+2 i u (2) n+2 2. obliczamy błąd |E| zgodnie z wzorem (26)

3. zmieniamy krok czasowy wg (27)

4. jeśli jest spełniony warunek |E| < tol to akceptujemy ”lepsze” rozwią- zanie u n+2 = u (2) n+2 i wykonujemy kolejny krok, w przeciwnym wypadku wykonujemy algorytm dla u n+2 zaczynając od punktu 1 ze zmienionym krokiem ∆t nowy (który automatycznie się zmniejszy)

Uwaga 1. Do wyboru są 2 warianty tego zadania - wybór wariantu należy do prowadzącego ćwiczenia laboratoryjne.

Uwaga 2. Dla każdego wariantu należy przyjąć parametry: λ = 100, startowy krok czasowy ∆t = 0.001, warunki początkowe u(0) = 1 i v(0) = 0, parametr antystagnacyjny S = 0.75, tolerancja błędu tol = 10 −4 , 10 −9 , t ∈ [0, t max ], t max = 400. Dla obu wartości parametru tol należy sporządzić jeden wykres przedstawiający oba rozwiązania oraz dwa wykresy pokazujące zmiany czasowe kroku ∆t(t). Bonus za komplet wyników (50pkt).

Wariant 1 - jawna metoda RK4

Rozwiązać problem przy użyciu jawnej metody RK4. Ponieważ mamy do czynie- nia z układem dwóch równań różniczkowych więc musimy zmodyfikować wzory podane w zadaniu 1.

Aby znaleźć rozwiązania u n+1 i v n+1 obliczamy ściśle w podanej kolejności:

(6)

k u1 = f (u n , v n ) (28)

k v1 = g(u n , v n ) (29)

k u2 = f (u n + ∆t

2 k u1 , v n + ∆t

2 k v1 ) (30)

k v2 = g(u n + ∆t

2 k u1 , v n + ∆t

2 k v1 ) (31)

k u3 = f (u n + ∆t

2 k u2 , v n + ∆t

2 k v2 ) (32)

k v3 = g(u n + ∆t

2 k u2 , v n + ∆t

2 k v2 ) (33)

k u4 = f (u n + ∆tk u3 , v n + ∆tk v3 ) (34) k v4 = g(u n + ∆tk u3 , v n + ∆tk v3 ) (35)

u n+1 = u n + ∆t

6 (k u1 + 2k u2 + 2k u3 + k u4 ) (36) v n+1 = v n + ∆t

6 (k v1 + 2k v2 + 2k v3 + k v4 ) (37) Przyjąć N = 5 przy obliczaniu (26) i (27) - rząd błędu lokalnego jest większy o 1 w stosunku do błędu globalnego metody (N = p + 1).

Wariant 2 - metoda trapezów

Rozwiązać problem stosując metodę trapezów. Ponieważ metoda jest niejawna (patrz zadanie 1) więc znalezienie rozwiązania w kolejnej chwili czasowej wyma- ga zastosowania metody Newtona. Funkcje f i g nie zależą w sposób jawny od czasu więc pomijamy go w dalszych równaniach.

Chcemy wyznaczyć 2 rozwiązania tj. u n+1 i v n+1 dla t n+1 : u n+1 = u n + ∆t

2 [f (u n , v n ) + f (u n+1 , v n+1 )] (38) v n+1 = v n + ∆t

2 [g(u n , v n ) + g(u n+1 , v n+1 )] (39) Problem jest dwuwymiarowy więc definiujemy dwie funkcje nieliniowe:

F = u n+1 − u n ∆t

2 [f (u n , v n ) + f (u n+1 , v n+1 )] (40) G = v n+1 − v n ∆t

[g(u n , v n ) + g(u n+1 , v n+1 )] (41)

(7)

szereg Taylora z dokładnością do wyrazów liniowych uzyskamy:

F (u n+1 + ∆u, v n+1 + ∆v) = F (u n+1 , v n+1 ) + ∂F

∂u ∆u + ∂F

∂v ∆v (42) G(u n+1 + ∆u, v n+1 + ∆v) = G(u n+1 , v n+1 ) + ∂G

∂u ∆u + ∂G

∂v ∆v (43) Z założenia metody Newtona wynika że F (u n+1 + ∆u, v n+1 + ∆v) = 0 oraz G(u n+1 + ∆u, v n+1 + ∆v) = 0, rozwinięcie możemy zapisać w postaci macierzo- wej:

[ −F

−G ]

=

[ ∂F/∂u n+1 ∂F/∂v n+1

∂G/∂u n+1 ∂G/∂v n+1

] [ ∆u

∆v ]

(44) To jest układ równań którego rozwiązanie ([∆u, ∆V ]) pozwoli zmienić u n+1 i v n+1 , robimy to oczywiście w sposób iteracyjny: u k+1 n+1 = u k n+1 + ∆u oraz v n+1 k+1 = v k n+1 + ∆v, starując od przybliżenia u 0 n+1 = u n i v n+1 0 = v n . W k + 1 iteracji wartości wektora [ −F, −G] (wzory 40 i 41) oraz macierzy układu wyznaczamy korzystając poprzedniego przybliżenia [u k n+1 , v k n+1 ].

Elementy macierzy układu A:

a 11 = ∂F/∂u n+1 = 1 (45)

a 12 = ∂F/∂v n+1 = ∆t

2 (46)

a 21 = ∂G/∂u n+1 = ∆t

2 ( −2λu k n+1 v n+1 k − 1) (47) a 22 = ∂G/∂v n+1 = 1 ∆t

2 (λ(1 − (u k n+1 ) 2 )) (48) (49) Rozwiązanie układu 2 × 2 znajdujemy metodą wyznacznikową:

∆u = ( −F )a 22 − (−G)a 12

a 11 a 22 − a 12 a 21 (50)

∆v = a 11 ( −G) − a 21 ( −F ) a 11 a 22 − a 12 a 21

(51) a następnie poprawiamy rozwiązanie.

Algorytm iteracyjnego wyznaczania u n+1 i v n+1 (czyli próba przejścia od t n do t n+1 - pamiętamy że uzyskane rozwiązanie może nie zostać zaakceptowane ze względu na błąd):

1. ustalamy przybliżenia startowe u 0 n+1 = u n i v n+1 0 = v n

2. dla danych u k n+1 i v k n+1 wyznaczamy wartości elementów wektora [ −F, −G]

oraz macierzy układu [a ij ]

3. rozwiązujemy układ równań i wyznaczamy nowe przybliżenia u k+1 n+1 i v k+1 n+1

(8)

4. punkt 2 i 3 powtarzamy K = 20 razy lub kończymy jeśli |∆u|, |∆v| <

10 −12

5. po zakończeniu iteracji Newtona wyznaczamy E u i E v wg wzoru (26) przyjmując N = 3

6. wyznaczamy |E| = max(|E u |, |E v |) 7. wyznaczamy ∆t nowy wg wzoru (27)

8. sprawdzamy warunek czy |E| < tol jeśli tak to zachowujemy u n+1 i v n+1 , w przeciwnym wypadku powtarzamy iterację Newtona ale z mniejszym

∆t nowy

0.1 Wariant 3 - niejawna metoda RK2

W przygotowaniu

Cytaty

Powiązane dokumenty

Ekstrapolacji Richardsona można użyć również do kontroli (zmiany w trakcie obliczeń) kroku czasowego ∆t, tak żeby błędy obcięcia nie przekraczały pewnej zadanej wartości

Na jednym rysunku pokazać trzy rozwiązania numeryczne i rozwiązanie analityczne... Na jednym rysunku umieścić wykresy

Często rozwiązanie zagadnienia brzegowego jest równocześnie roz- wiązaniem pewnego zagadnienia wariacyjnego, tzn... Aby sprawdzić czy rozwiązania są stabilne, porównać

We węzłach brzegowych u jest równa zeru jak w warunkach, więc nie trzeba

Wówczas, aby rozwiązać równanie wystarczy podać wszystkie jego rozwiązania integralne, gdyż każde inne rozwiązanie jest obcięciem pewnego rozwiązania integralnego do

Jak pokazaliśmy w przykładzie 1.3.1., każde rozwią- zanie tego równania określone jest na pewnym przedziale zawartym w dziedzinie jednego z powyższych rozwiązań, więc

Przykład: Funkcja obliczająca

[r]