Ćwiczenie 2b: Schematy iteracyjne dla równania Poissona: relaksacja lokalna, globalna i
wielosiatkowa.
26 kwietnia 2016
Zadanie polega na rozwiązaniu numerycznym równania Poissona
∇2V = ∂2V
∂x2 +∂2V
∂y2 =−ρ(x, y) (1)
na dwuwymiarowej siatce. Wprowadzamy siatkę indeksując węzły i, j = 0, 1, 2, . . . , n dla n = 128 oraz odległości międzywęzłowe ∆ = ∆x = ∆y = 0.01. Wartości x i y wyliczamy według wzorów: xi = i· ∆, yj = j· ∆. W naszym problemie wykorzystamy warunki brzegowe typu Dirichleta (patrz rysunek 1):
1. lewy brzeg V0,j = VD1, j = 1, 2, . . . , n− 1, VD1= +0.1 2. prawy brzeg Vn,j= VD2, j = 1, 2, . . . , n− 1, VD2=−0.1 3. dolny brzeg Vi,0= VD3, i = 0, 2, . . . , n, VD3= 0
4. górny brzeg Vi,n= VD4, i = 0, 2, . . . , n, VD4= 0 Gęstość definiujemy jako
ρ(x, y) =
+Q0 ⇐⇒ (x− x1)2+ (y− y1)2¬ R2
−Q0 ⇐⇒ (x− x2)2+ (y− y2)2¬ R2 0 w pozostałym obszarze
(2) Przyjmujemy: Q0= 80, x1= 0.4, y1= 0.64, x2= 0.88, y2= 0.64, R = 0.051.
1 Relaksacja globalna i lokalna
Stosując relakasację globalną, w jednej iteracji, wyznaczamy najpierw nowe war- tości potencjału w całej tablicy ( ˜V ):
V˜i,j=Vi+k,jold + Viold−k,j+ Vi,j+kold + Vi,jold−k+ (k· ∆)2ρi,j
4 , i, j = k, 2k, . . . , n− k (3)
V
D2V
D1V
D4V
D3(i,j) ( i,j+k)
i = 0 k 2k n
k 2k n
j = 0
Rysunek 1: Schemat rozmieszczenia węzłów na siatce przy przejściu z siatki rzadkiej (czerwone punkty) do siatki dwa razy gęstszej (czerwone+czarne). Ob- szary zaznaczone VD1, VD2, VD3, VD4 pokazują węzły brzegowe na które należy nałożyć warunki brzegowe Dirichleta. Żółty i biały kwadrat pokazują obszary wewnątrz których należy uśredniać gęstość na rzadkiej siatce (wartości średnie liczone są na najgęstszej siatce ale przypisywane są do węzłów czerwonych).
gdzie: k - krok na siatce, a następnie nowe przybliżenie mieszając potencjał ˜V ze starymi wartościami Vold:
Vi,jnew= (1− ωG)Vi,jold+ ωGV˜i,j, i, j = k, 2k, . . . , n− k (4) gdzie: ωG ∈ (0, 1) jest parametrem zbieżności dla relaksacji lokalnej.
W relaksacji lokalnej sprawa wygląda prościej gdyż zmiany wprowadzamy na bieżąco do tablicy na której operujemy:
Vi,j = (1−ωL)Vi,j+ωL
Vi+k,j+ Vi−k,j+ Vi,j+k+ Vi,j−k+ (k· ∆)2ρi,j
4 , i, j = k, 2k, . . . , n−k (5)
gdzie: ωL∈ (0, 2) jest parametrem zbieżności.
W celu kontroli jakości bieżącego rozwiązania korzystamy z funkcjonału energii:
S =
∫ ∫ dxdy
{ 1 2
[(∂V
∂x )2
+ (∂V
∂y )2]
− ρV }
(6)
który w wersji dyskretnej przybiera postać:
S =
n∑−k i=0
n∑−k j=0
(k· ∆)2 {
1 2
[(Vi+k,j− Vi,j
k· ∆ )2
+
(Vi,j+k− Vi,j
k· ∆
)2]
− ρi,jVi,j }
(7) przy czym należy pamiętać, że dla ustalonej wartości k, na siatce przechodzimy co k węzłów w kierunkach ’x’ i ’y’.
Relaksację prowadzimy aż do momentu gdy spełniony zostanie warunek Sit+1− Sit
Sit
< TOL (8)
gdzie: it jest numerem iteracji.
Zadania do wykonania:
1. Przyjąć parametry: T OL = 10−12, k = 1
2. Wykonać relaksację globalną (wzory 3 i 4) dla ωG = 0.1, 0.2, . . . , 0.9 aż do uzyskania zbieżności. Na jednym rysunku zamieścić zmiany wartości S (wzór 7) w funkcji numeru iteracji dla wszystkich wartości ωG. Po zre- laksowaniu potencjału należy wykonać wykresy: a) potencjału, b)gęstości oraz c) różnicy δ =∇2V + ρ (w celu sprawdzenia błędów) (25pkt) 3. Wykonać relaksację lokalną (wzór 5) dla ωL = 1.1, 1.2, . . . , 1.9 aż do uzy-
skania zbieżności. Na jednym rysunku zamieścić zmiany wartości S (wzór 7) w funkcji numeru iteracji dla wszystkich wartości ωL. Po zrelaksowa- niu potencjału należy wykonać wykresy: a) potencjału, b)gęstości oraz c) różnicy δ =∇2V + ρ (w celu sprawdzenia błędów) (25pkt)
2 Prosta relaksacja wielosiatkowa
Relaksację prowadzić będziemy startując od najrzadzszej siatki k = 16. Po uzyskaniu zbieżności przechodzimy na siatkę dwa razy gęstszą w kolejności k = 16 → 8 → 4 → 2 → 1. Przy przejściu z siatki k na siatkę k/2 w nowo dodanych węzłach potencjał przybliżamy licząc średnią arytmetyczną potencja- łów dla najbliższych sąsiadów (rys. 2).
V
i,jV
i+k,jV
i+k/2,j=(V
i,j+V
i+k,j)/2 V
i+k/2,j+k/2=(V
i,j+V
i+k,j+V
i,j+k+V
i+k,j+k)/4
Rysunek 2: Schemat uśredniania potencjału w nowo dodanych węzłach (kolor czerwony - stare węzły, kolor niebieski i czarny - nowe węzły).
W relaksacji wielosiatkowej należy jeszcze rozwiązać problem sposobu uwzględ- nienia gęstości. Mianowicie wykonując duże kroki na siatce (np. k = 14) może się okazać że gęstość, która jest skupiona w małych obszarach może zostać po- minięta. Aby temu zapobiec należy zmodyfikować tablicę gęstości. Przez ρ(1) oznaczmy tablicę gęstości na najgęstszej siatce (k = 1). Elementy tablicy gęsto- ści ρ(k) wyznaczamy uśredniając gęstość w zakresach i∈ [i − k/2, i + k/2] oraz j∈ [j − k/2, j + k/2]. Należy pamiętać że węzły skrajne znajdujące się w danej komórce wnoszą też wkłady do sąsiednich komórek więc musimy je skalować tak aby całka z gęstości po całym obszarze była zachowana. Skalowanie gęstości wykonujemy przed rozpoczęciem relaksacji licząc kolejno (w poniższych sumach i′ oraz j′ zmieniają się co 1 niezależnie od wartości k):
1. sumę wag( wi = wj = 1.0 dla węzłów w środku, oraz wi = wj = 0.5 dla węzłów brzegowych - żółty kwadrat na rysunku 1 ):
Wi,j(k)=
i+k/2∑
i′=i−k/2 j+k/2∑
j′=j−k/2
wi′· wj′ (9)
2. całkę z gęstości w komórce scentrowanej w punkcie (i, j) uwzględniając wagi:
Ci,j(k)=
i+k/2∑ j+k/2∑
ρ(1)i′,j′· wi′· wj′ (10)
3. uśredniając wartość gęstości w komórce (a dokładniej w węźle centralnym komórki):
ρ(k)i,j = Cij(k) Wij(k)
(11)
Zadania do wykonania:
Przeprowadzić relaksację wielosiatkową dla siatek k = 16, 8, 4, 2, 1. Dla każdego k wykonać relaksację według schematu lokalnego (wzór 5) przyjmując wartość parametru zbieżności ωL= 1.9, reszta parametrów jak w zadaniu 1. Na każdej z siatek w każdej iteracji wyliczać wartość funkcjonału Sit(k)(wzór 7). Na jednym rysunku proszę pokazać wykresy wartości funkcjonału energii w funkcji numeru iteracji dla k = 16, 8, 4, 2, 1 natomiast dla każdego k sporządzić oddzielny wykres zrelaksowanego potencjału. Za wyniki z zadania 2 można uzyskać (50 pkt).
