Projekt 3: RRZ - kontrola kroku czasowego w problemach sztywnych.

(1)

Projekt 3: RRZ - kontrola kroku czasowego w problemach sztywnych.

Tomasz Chwiej 4 listopada 2020

1 Wstęp

Na zajęciach rozwiążemy równanie różniczkowe 2 rzędu oscylatora Van der Pola d

²

x

dt

²

− α(1 − x

²

) dx

dt + x = 0 (1)

RRZ 2 rzędu zamieniamy na układ 2 RRZ 1 rzędu (wprowadzając nową zmienną dx/dt = v):

dx

dt = f (t, x, v) = v (2)

dv

dt = g(t, x, v) = α(1 − x

²

)v − x (3)

Taki układ RRZ 1 rzędu rozwiążemy stosując metodę jawną RK2 i metodę niejawną trapezów.

Dla α = 0 problem sprowadza się do problemu oscylatora harmonicznego, natomiast dla α >> 1 problem staje się sztywny. Z tego powodu problem rozwiążemy kontrolując krok czasowy.

-2.5 -2 -1.5 -1 -0.5 0 0.5 1 1.5 2 2.5

-2 -1.5 -1 -0.5 0 0.5 1 1.5 2

v(x)

x alfa = 0.2

-3 -2 -1 0 1 2 3

-2.5 -2 -1.5 -1 -0.5 0 0.5 1 1.5 2 2.5

v(x)

x alfa = 1.0

-8 -6 -4 -2 0 2 4 6 8

-2.5 -2 -1.5 -1 -0.5 0 0.5 1 1.5 2 2.5

v(x)

x alfa = 5.0

Rysunek 1: Rozwiązanie w przestrzeni fazowej V (x) dla parametru α = 0.2; 1; 5. Warunki początkowe:

x = 0.01 i V = 0. Niezależnie od WP, rozwiązanie zawsze będzie dążyć do ustalonej trajektorii (istnieje atraktor).

2 Kontrola kroku czasowego

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 x

⁽¹⁾_n+2

różniące od dokładnego o pewien błąd lokalny:

(1) p+1 p+2

(2)

gdzie: p oznacza rząd metody (dla RK2 p = 2, dla trapezów p = 2). Możemy też wykonać dwa krótsze kroki co ∆t i otrzymać lepsze rozwiązanie x

⁽²⁾_n+2

co da nam inny błąd lokalny:

x

dok

(t

n+2

) = x

⁽²⁾_n+2

+ 2C

x

(∆t)

^p+1

+ O(∆t

^p+2

) (5) Stałą błędu C

_x

wyznaczamy odejmując od (5) równanie (4). Następnie możemy określić błąd rozwią- zania uzyskanego w dwóch krokach (dla x i V):

E

_x

= 2C

_x

∆t

^p+1

= x

⁽²⁾_n+2

− x

⁽¹⁾_n+2

2

^p

− 1 (6)

E

_v

= 2C

_v

∆t

^p+1

= v

⁽²⁾_n+2

− v

_n+2⁽¹⁾

2

^p

− 1 (7)

i wykorzystać go do modyﬁkacji kroku czasowego:

∆t

_nowy

=

( S · tol max ( |E

x

|, |E

v

|)

)

¹

p+1

· ∆t

stary

(8)

gdzie funkcja max(a, b) zwraca większą z dwóch porównywanych liczb.

Ogólny algorytm numerycznego rozwiązywania równania różniczkowego z doborem kroku czasowego:

i n i c j a l i z a c j a : t =0 ,∆t = ∆t

0

, x

n

= x

0

, v

n

= v

0

, t

max

DO

// s t a w i a m y dwa k r o k i ∆t

(x

⁽²⁾_n+1

, v

⁽²⁾_n+1

) ← schemat numeryczny(x

n

, v

_n

, ∆t, α) (x

⁽²⁾_n+2

, v

⁽²⁾_n+2

) ← schemat numeryczny(x

⁽²⁾_n+1

, v

_n+1⁽²⁾

, ∆t, α) // s t a w i a m y j e d e n k r o k 2 · ∆t

(x

⁽¹⁾_n+2

, v

⁽¹⁾_n+2

) ← schemat numeryczny(x

n

, v

_n

, 2 · ∆t, α) l i c z y m y : E

_x

, E

_v

// s p r a w d z a m y czy w y n i k j e s t a k c e p t o w a n y IF (max( |E

x

|, |E

v

|)<TOL ) THEN

t = t + 2 · ∆t x

n

= x

⁽²⁾_n+2

v

n

= v

_n+2⁽²⁾

z a p i s d a n y c h do p l i k u : t ,∆t, x

_n

, v

_n

END IF

// z m i a n a k r o k u n a s t ę p u j e z a w s z e

∆t =

(

S·T OL max(|Ex|,|Ev|)

)

¹

p+1

· ∆t

E N D D O W H I L E ( t <t

max

)

(3)

2.1 schemat numeryczny: metoda trapezów

Chcemy wyznaczyć 2 rozwiązania tj. x

_n+1

i v

_n+1

dla t

_n+1

: x

n+1

= x

n

+ ∆t

2 [f (x

n

, v

n

) + f (x

n+1

, v

n+1

)] (9) v

_n+1

= v

_n

+ ∆t

2 [g(x

_n

, v

_n

) + g(x

_n+1

, v

_n+1

)] (10) Problem jest dwuwymiarowy więc deﬁniujemy dwie funkcje nieliniowe:

F = x

_n+1

− x

n

− ∆t

2 [f (x

_n

, v

_n

) + f (x

_n+1

, v

_n+1

)] (11) G = v

n+1

− v

n

− ∆t

2 [g(x

n

, v

n

) + g(x

n+1

, v

n+1

)] (12) Rozwijając funkcje F (x

n+1

+∆x, v

n+1

+∆v) i G(x

n+1

+∆x, v

n+1

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

F (x

_n+1

+ ∆x, v

_n+1

+ ∆v) = F (x

_n+1

, v

_n+1

) + ∂F

∂x ∆x + ∂F

∂v ∆v (13)

G(x

n+1

+ ∆x, v

n+1

+ ∆v) = G(x

n+1

, v

n+1

) + ∂G

∂x ∆x + ∂G

∂v ∆v (14)

Z założenia metody Newtona wynika że F (x

n+1

+∆x, v

n+1

+∆v) = 0 oraz G(x

n+1

+∆x, v

n+1

+∆v) = 0, więc rozwinięcie możemy zapisać w postaci macierzowej

[ −F

−G ]

= [

_∂F

∂xn+1

∂F

∂vn+1

∂G

∂xn+1

∂G

∂vn+1

] [

∆x

∆v ]

(15) Elementy macierzy układu A (górny indeks k to numer iteracji):

a

₁₁

= ∂F

∂x

_n+1

= 1 (16)

a

12

= ∂F

∂v

n+1

= − ∆t

2 (17)

a

21

= ∂G

∂x

_n+1

= − ∆t

2 ( −2αx

^kn+1

v

_n+1^k

− 1) (18) a

₂₂

= ∂G

∂v

_n+1

= 1 − ∆t 2 α

[

1 − (x

^kn+1

)

²

]

(19) (20) Rozwiązanie układu 2 × 2 znajdujemy metodą wyznacznikową:

∆x = (−F )a

22

− (−G)a

12

a

11

a

22

− a

12

a

21

(21)

∆v = a

11

(−G) − a

21

(−F ) a

11

a

22

− a

12

a

21

(22) a następnie poprawiamy rozwiązanie

x

^k+1_n+1

= x

^k_n+1

+ ∆x (23)

v

_n+1^k+1

= v

^k_n+1

+ ∆v (24)

Na starcie przyjmujemy: x

⁰_n+1

= x

_n

i v

_n+1⁰

= v

_n

, a następnie iterację prowadzimy aż do spełnienia

|∆x| < δ i |∆v| < δ dla δ = 10

⁻¹⁰

(4)

2.2 schemat numeryczny: metoda RK2 Obliczenia prowadzimy sekwencyjnie

k

1x

= f (t

n

, x

n

, v

n

) = v

n

(25)

k

1v

= g(t

n

, x

n

, v

n

) = α(1 − x

²_n

)v

n

− x

n

(26)

k

_2x

= f (t

_n

+ ∆t, x

_n

+ ∆tk

_1x

, v

_n

+ ∆tk

_1v

) = v

_n

+ ∆tk

_1v

(27) k

_2v

= g(t

_n

+ ∆t, x

_n

+ ∆tk

_1x

, v

_n

+ ∆tk

_1v

) (28)

= α [

1 − (x

n

+ ∆tk

1x

)

²

]

(v

n

+ ∆tk

1v

) − (x

n

+ ∆tk

1x

) (29)

x

n+1

= x

n

+ ∆t

2 (k

1x

+ k

2x

) (30)

v

_n+1

= v

_n

+ ∆t

2 (k

_1v

+ k

_2v

) (31)

3 Zadania do wykonania

1. Zaprogramować metody: trapezów i RK2 jako dwie osobne procedury. Do procedur przekazuje- my: x

n

, v

n

, ∆t, α; a mają zwracać: x

n+1

, v

n+1

.

2. Zaimplementować algorytm kontroli kroku czasowego.

3. Przyjąć parametry startowe: x

₀

= 0.01, v

₀

= 0, ∆t

₀

= 1, S = 0.75, p = 2 (rząd dokładności obu metod), t

max

= 40, α = 5.

4. Rozwiązać równanie (1) metodą RK2 z kontrolą kroku czasowego dla parametru T OL = 10

⁻²

; 10

⁻⁵

. Wykonać rysunki: x(t), v(t), ∆t(t) i v(x). Wykresy tej samej wielkości dla obu wartości T OL umieścić na jednym rysunku (będą 4 rysunki). (50 pkt)

5. Rozwiązać równanie (1) metodą trapezów z kontrolą kroku czasowego dla parametru T OL = 10

⁻²

; 10

⁻⁵

. Wykonać rysunki: x(t), v(t), ∆t(t) i v(x). Wykresy tej samej wielkości dla obu wartości T OL umieścić na jednym rysunku (będą 4 rysunki). (50 pkt)

4 Przykładowe wyniki

-8 -6 -4 -2 0 2 4 6 8

0 5 10 15 20 25 30 35 40 45 v(t)

tol=10^-5 tol=10^-2

example

-2.5 -2 -1.5 -1 -0.5 0 0.5 1 1.5 2 2.5

0 5 10 15 20 25 30 35 40 45 x(t)

tol=10^-5 tol=10^-2

example

0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1

0 5 10 15 20 25 30 35 40 45 dt(t)

tol=10^-5 tol=10^-2

example

-8 -6 -4 -2 0 2 4 6 8

-2.5 -2 -1.5 -1 -0.5 0 0.5 1 1.5 2 2.5 v(x)

tol=10^-5 tol=10^-2

example

Rysunek 2: Wyniki dla metody trapezów

(5)

-8 -6 -4 -2 0 2 4 6 8

0 5 10 15 20 25 30 35 40 45 v(t)

tol=10^-5 tol=10^-2

example

-2.5 -2 -1.5 -1 -0.5 0 0.5 1 1.5 2 2.5

0 5 10 15 20 25 30 35 40 45 x(t)

tol=10^-5 tol=10^-2

example

0 0.1 0.2 0.3 0.4 0.5 0.6 0.7

0 5 10 15 20 25 30 35 40 45 dt(t)

tol=10^-5 tol=10^-2

example

-8 -6 -4 -2 0 2 4 6 8

-2.5 -2 -1.5 -1 -0.5 0 0.5 1 1.5 2 2.5 v(x)

tol=10^-5 tol=10^-2

example

Rysunek 3: Wyniki dla metody RK2