Metoda elementów skończonych (FEM)
Marcin Orchel
Spis treści
1 Wstęp 1
1.1 Rachunek wariacyjny . . . 1
1.2 Sprowadzanie problemu brzegowego do zagadnienia wariacyjnego . . . 5
1.2.1 Warunki brzegowe dla równania Poissona . . . 7
1.2.2 Jednowymiarowy problem brzegowy . . . 7
1.3 Metoda Galerkina . . . 10
1.4 Metoda elementów skończonych . . . 11
2 Zadania 13 2.1 Przydatne polecenia . . . 13
2.2 Zadania na 3.0 . . . 13
2.3 Zadania na 4.0 . . . 14
2.4 Zadania na 5.0 . . . 14
1 Wstęp
1.1 Rachunek wariacyjny
W rachunku wariacyjnym szukamy takiej funkcji y (x), dla której wyrażenie całkowe I [y] =
Z b a
Fx, y (x) , y0(x) , . . . , y(n)(x)dx (1) osiąga ekstremum, gdzie y (x) należy do ściśle określonej klasy funkcji, I to funkcjonał mapujący funkcję do wartości skalarnej. Na funkcję y (x) i jej pochodne nakłada się warunki brzegowe i warunki dodatkowe (poboczne).
Przykłady zagadnień wariacyjnych:
Przykład 1 (Zagadnienie izoperymetryczne). Znaleźć w kartezjańskim układzie współ- rzędnych taką krzywą y = y (x) łączącą punkty A (a, 0) i B (b, 0), że jej długość wynosi l i razem z odcinkiem AB odcina powierzchnię o możliwie największym polu. W rachunku wariacyjnym: znaleźć dwukrotnie różniczkowalną funkcję y (x), taką, że:
maxy I [y] = Z b
a
y (x) dx (2)
z warunkiem pobocznym:
G [y] = Z b
a
q
1 + y02(x)dx = l (3)
i z warunkami brzegowymi:
y (a) = y (b) = 0 (4)
Wyprowadzenie wzoru na długość łuku krzywej y(x):
ds = q
dx2+ dy2 (5)
l = Z b
a
ds = Z b
a
s
dx dt
2
+
dy dt
2
dt (6)
Jeśli t = x:
l = Z b
a
s 1 +
dy dx
2
dx (7)
Rozwiązaniem tego zagadnienia jest część okręgu.
Przykład 2 (Zagadnienie brachistochrony). Punkt P0(x0, y0), leżący w umieszczonej pionowo płaszczyźnie Oxy, należy tak połączyć krzywą y = y (x) z początkiem ukła- du współrzędnych, aby punkt materialny poruszający się bez tarcia wzdłuż tej krzywej, wyłącznie pod wpływem siły ciężkości, znalazł się w początku układu w możliwie najkrót- szym czasie. W rachunku wariacyjnym: szukamy funkcji y = y (x), ciągłej wraz ze swoją pierwszą pochodną, dla której:
T [y] = Z te
0
dt = Z x0
0
p1 + y02
p2g (y0− y)dx = min (8) z warunkami brzegowymi:
y (0) = 0 (9)
y (x0) = y0 (10)
Wyprowadzenie z mv2/2 = mg(y0− y), v =p2g(y0− y), a także v = ds/dt, dt = ds/v, ds =p1 + y02dx. Rozwiązaniem jest brachistochrona, czyli część cykloidy.
Proste zagadnienie wariacyjne: wyznaczyć wartości ekstremalne funkcjonału postaci:
I [y] = Z b
a
F x, y (x) , y0(x)dx (11) gdzie y (x) jest zbiorem funkcji ciągłych wraz ze swoimi pochodnymi do drugiego rzędu włącznie, spełniających warunki brzegowe:
y (a) = A (12)
y (b) = B (13)
Funkcjonał I [y] osiąga np. dla funkcji y0(x) maksimum lokalne, wtedy dla pewnego otoczenia tej funkcji zawartego w dziedzinie funkcyjnej zachodzi nierówność:
I [y0] ≥ I [y] (14)
Krzywą y = y0(x) nazywamy ekstremalą funkcjonału lub funkcją ekstremalną.
Warunek konieczny rozwiązania prostego zagadnienia wariacyjnego. Dla rozwiązania y0(x) spełniającego (14) konstruujemy tzw. funkcję porównawczą:
y (x) = y0(x) + εη (x) (15)
gdzie η (x) ma ciągłe pochodne aż do rzędu drugiego włącznie oraz spełnia specjalne warunki brzegowe
η (a) = η (b) = 0 (16)
Podstawiając (15) do (11) otrzymujemy:
I (ε) = Z b
a
F x, y0+ εη, y00 + εη0dx (17) Zadanie aby funkcjonał I [y] osiągał w punkcie y0(x) ekstremum, przechodzi w zadanie, aby funkcja I (ε) miała ekstremum w punkcie ε = 0. Warunek konieczny na ekstremum brzmi:
dI
dε = 0 (18)
dla
ε = 0 (19)
Wzór Taylora dla funkcji wielu zmiennych. Istnieje takie 0 < φ < 1, że:
(∆f )(x
0,y0) = (df )(x
0,y0)+1 2
d2f
(x0,y0)+ . . . + 1 (n − 1)!
dn−1f
(x0,y0) (20) + 1
(n)!(dnf )(x0+φdx,y0+φdy) (21)
gdzie
(∆)(x
0,y0)= f (x0+ dx, y0+ dy) − f (x0, y0) (22) Różniczka zupełna rzędu n dla 3 zmiennych ma postać:
dnu =
∂
∂xdx + ∂
∂ydy + ∂
∂zdz
n
u (23)
Dla n = 1
du = ∂u
∂xdx + ∂u
∂ydy +∂u
∂zdz (24)
Dla n = 2 otrzymujemy d (du):
d2u =∂2u
∂x2dx2+ ∂u2
∂xydxdy + ∂u2
∂xzdxdz (25)
+∂2u
∂y2dy2+ ∂u2
∂xydxdy + ∂u2
∂yzdydz (26)
+∂2u
∂z2dz2+ ∂u2
∂xzdxdz +∂u2
∂yzdydz = (27)
∂2u
∂x2dx2+ ∂2u
∂y2dy2+ ∂2u
∂z2dz2+ 2∂u2
∂xydxdy + 2∂u2
∂xzdxdz + 2∂u2
∂yzdydz (28) Dla funkcji trzech zmiennych wzór można zapisać w postaci:
f (x0+ h, y0+ k, z0+ l) = f (x0, y0, z0)+fx0 (x0, y0, z0) h+fy0 (x0, y0, z0) k+fz0(x0, y0, z0) l (29) +1
2
fxx00 h2+ fyy00 k2+ fzz00l2+ 2fxy00hk + 2fxz00hl + 2fyz00kl
(x0+φh,y0+φk,z0+φl) (30) Zakładając różniczkowalność funkcji F po rozwinięciu jej w szereg Taylora dla n = 1 w punkcie (x, y0, y00) otrzymujemy:
I (ε) = Z b
a
F x, y0, y00+∂F
∂y x, y0, y00εη +∂F
∂y0 x, y0, y00εη0+ Oε2dx (31) Różniczkując powyższe po i przyjmując po różniczkowaniu, że = 0 warunek konieczny przyjmuje postać:
Z b a
η∂F
∂ydx + Z b
a
η0∂F
∂y0dx = 0 (32)
Po całkowaniu przez części Z b
a
η∂F
∂ydx +
η∂F
∂y0
b a
− Z b
a
η d dx
∂F
∂y0
dx = 0 (33)
i zauważeniu, że η jest zerem na brzegu otrzymujemy Z b
a
η
∂F
∂y − d dx
∂F
∂y0
dx = 0 (34)
Równanie powyższe powinno być spełnione dla dowolnej funkcji η, dlatego musi zacho- dzić:
∂F
∂y − d dx
∂F
∂y0
= 0 (35)
Równanie to jest warunkiem koniecznym na rozwiązanie prostego zagadnienia wariacyj- nego i nazywane jest równaniem Eulera.
Przykład 3. Wyznaczyć najkrótszą linię, która łączy dwa punkty P1(a, A) i P2(b, B) leżące w płaszczyźnie Oxy:
miny I [y] = Z b
a
q
1 + y02dx (36)
F x, y, y0= F y0= q
1 + y02 (37)
∂F
∂y0 = y0
p1 + y02 (38)
d dx
y0 p1 + y02
!
=
y00p1 + y02− √y02y00
1+y02
1 + y02 = y00
p1 + y02 (39)
Warunek Eulera jest następujący:
y00
p1 + y02 = 0 (40)
Stad
y00= 0 (41)
A więc najkrótszą linią łączącą dwa punkty jest prosta.
Zadanie. Dla funkcjonału
I [y] = Z b
a
(x − y)2dx (42)
znaleźć równanie Eulera, rozwiązać je, oraz obliczyć wartość funkcjonału dla otrzymanej ekstremali.
1.2 Sprowadzanie problemu brzegowego do zagadnienia wariacyjnego Przykład dwuwymiarowy: Dane jest zagadnienie brzegowe z operatorem Laplace’a (la- plasjanem) ∆u:
∆u = ∂2u
∂x2 + ∂2u
∂y2 = f (x, y) (43)
we wnętrzu obszaru G, gdzie
u = 0 (44)
na brzegu obszaru G. Powyższe równanie różniczkowe znane jest jako równanie Poissona.
Będziemy wykorzystywać wystarczająco gładką funkcję v (x, y), zwaną funkcją testową, która znika na brzegu obszaru G (v = 0).
Stosujemy twierdzenie Greena opisujące zależność między całką po powierzchni ogra- niczonej i całką po krzywej zamkniętej postaci:
Z Z
(B)
∂Q (x, y)
∂x −∂P (x, y)
∂y
dxdy = I
(K)
[P (x, y) dx + Q (x, y) dy] (45)
gdzie B jest powierzchnią o brzegu K, P i Q są funkcjami ciągłymi i mają ciągłe pierwsze pochodne cząstkowe. Przyjmując
P (x, y) = −v (x, y) uy (46)
Q (x, y) = v (x, y) ux (47)
otrzymujemy:
Z Z
(G)
∂vux
∂x + ∂vuy
∂y
dxdy = I
(K)
[−vuydx + vuxdy] (48)
Po przekształceniu:
Z Z
(G)
"
∂u
∂x
∂v
∂x+ v∂2u
∂x2 +∂u
∂y
∂v
∂y + v∂2u
∂y2
#
dxdy = I
(K)
[−vuydx + vuxdy] (49)
Z Z
(G)
∂u
∂x
∂v
∂x +∂u
∂y
∂v
∂y
dxdy + Z Z
(G)
"
v∂2u
∂x2 + v∂2u
∂y2
#
dxdy = I
(K)
[−vuydx + vuxdy] (50)
Następnie po pomnożeniu (43) przez funkcję v otrzymujemy
∂2u
∂x2 +∂2u
∂y2
!
v = f v (51)
i po podstawieniu otrzymujemy:
Z Z
(G)
∂u
∂x
∂v
∂x+∂u
∂y
∂v
∂y
dxdy + Z Z
(G)
f vdxdy = I
(K)
[−vuydx + vuxdy] (52)
Prawa strona jest równa zeru, ponieważ v = 0. A więc mamy znaleźć takie u, że dla każdego v zachodzi:
− Z Z
(G)
∂u
∂x
∂v
∂x+∂u
∂y
∂v
∂y
dxdy = Z Z
(G)
f vdxdy (53)
Wprowadzając oznaczenia a na funkcjonał biliniowy, a (u, v) = −
Z Z
(G)
∂u
∂x
∂v
∂x+ ∂u
∂y
∂v
∂y
dxdy (54)
oraz b na funkcjonał liniowy
b (v) = Z Z
(G)
f vdxdy (55)
otrzymujemy
a (u, v) = b (v) (56)
Postać ta jest nazywana zagadnieniem wariacyjnym dla równania Poissona.
1.2.1 Warunki brzegowe dla równania Poissona
• warunki Dirichleta: specyfikacja wartości funkcji na brzegu danego obszaru
• warunki Neumanna: specyfikacja wartości pochodnej normalnej ∂u∂n na brzegu da- nego obszaru, pochodna normalna jest zdefiniowana jako
∂y
∂n(x) = ∇y (x) · n (x) (57)
gdzie ∇f to operator gradientu, przykładowo zdefiniowany jako
∇f = ∂f
∂x~i +∂f
∂y~j +∂f
∂z~k (58)
gdzie ~i, ~j, ~k to wektory jednostkowe.
• warunki Cauchy’ego – specyfikacja zarówno wartości funkcji jak i pochodnej nor- malnej na brzegu obszaru
• warunki Robina: specyfikacja wartości wyrażeń postaci:
αU + β∂u
∂n (59)
gdzie α, β = const, oraz
α2+ β2 6= 0 (60)
1.2.2 Jednowymiarowy problem brzegowy Dane jest równanie różniczkowe postaci:
d2u
dx2 = f (x) (61)
dla
u (a) = u (b) = 0 (62)
G = (a, b) (63)
Sprowadzamy powyższe równanie do zagadnienia wariacyjnego podobnie jak dla przy- padku dwuwymiarowego, z tą różnicą, że stosujemy wzór na całkowanie przez części zamiast wzoru całkowego Gaussa. Najpierw mnożymy obie strony przez funkcję v(x) i całkujemy otrzymując
Z
G
d2u dx2vdx =
Z
G
f vdx (64)
Z
G
d dx
du dx
vdx =
Z
G
f vdx (65)
Wzór na całkowanie przez części wygląda następująco:
b
Z
a
u0vdx = uv|ba−
b
Z
a
uv0dx (66)
Stosując wzór na całkowanie przez części do (65) otrzymujemy:
du dxv
b a
− Z
G
du dx
dv dxdx =
Z
G
f vdx (67)
Ponieważ v (a) = 0 oraz v (b) = 0, to:
− Z
G
du dx
dv dxdx =
Z
G
f vdx (68)
inaczej zapisane
− Z
G
u0v0dx = Z
G
f vdx (69)
Podobnie jak dla dwóch wymiarów, możemy powyższe zapisać w postaci
a (u, v) = b (v) (70)
gdzie
a (u, v) = − Z
G
u0v0dx (71)
b (v) = Z
G
f vdx (72)
Zagadnienie wariacyjne: znajdź u dla którego powyższe równanie jest spełnione dla każ- dego v. To równanie jest nazywane również sformułowaniem słabym. Możemy to sformu- łowanie słabe wyprowadzić również dla innego równania z dodatkowymi składnikami z pochodnymi pierwszego rzędu i z samą funkcją u. Wtedy dokonujemy całkowania przez części tylko do pierwszego składnika z drugimi pochodnymi. Równanie to zachodzi dla odpowiedniej klasy funkcji u i testowych.
Jaki jest związek tego równania z zagadnieniem wariacyjnym? Rozważmy zagadnienie wariacyjne
I (y) = min
y
1 2
Z b a
hy02(x) + 2f (x) y (x)idx (73) Zapiszmy warunek konieczny w formie (32) oznaczając η jako v i otrzymujemy
Z b
a
v∂F
∂ydx + Z b
a
v0∂F
∂y0dx = 0 (74)
1 2
Z b a
2f vdx +1 2
Z b a
2v0y0dx = 0 (75)
Z b a
f vdx + Z b
a
v0y0dx = 0 (76)
Przykład 4. Wyprowadzić postać słabą (wariacyjną) dla problemu brzegowego
− u00+ u0+ 2u = 3x2− x (77)
dla x ∈ (−1, 1) oraz
u (−1) = 1 (78)
u (1) = 0 (79)
Najpierw należy sprowadzić zagadnienie do warunków brzegowych zerowych. W tym celu u zapisujemy jako sumę u = w + ˜u, gdzie w musi spełniać powyższe warunki brzegowe, zaś ˜u warunki zerowe.
− ˜u00+ ˜u0+ 2˜u = 3x2− x + w00− w0− 2w (80) dla x ∈ (−1, 1) oraz
u (−1) = 0˜ (81)
u (1) = 0˜ (82)
A następnie podobnie jak wcześniej, stosując całkowanie przez części również dla w00. Możemy łatwo wybrać funkcję w tak aby były spełnione warunki brzegowe. Przykładowo
w (x) = α +β − α
b − a (x − a) (83)
w (x) = 1 + −1/2 (x + 1) (84)
Przykład 5. Wyprowadź postać słabą następującego problemu brzegowego
− u00+ x2u0+ 4u = x2− 2x3 (85) dla x ∈ (0, 2) z warunkami
− u0(0) + 2u (0) = −1 (86)
u (2) = 1 (87)
Rozwiązaniem pierwszego warunku jest funkcja u postaci u (x) = Ce2x−1
2 (88)
A więc
u (0) = C −1
2 (89)
Na początku włączamy pierwszy warunek Robina. Postać tego warunku dla przedzia- łu jednowymiarowego jest tutajhttp: // en. wikipedia. org/ wiki/ Robin_ boundary_
condition, a sposób włączenia omówiony jest tutaj http: // math. stackexchange.
com/ questions/ 361423/ variational-formulation-of-robin-boundary-value-problem-
for-poisson-equation-in. Po pomnożeniu przez funkcję v(x) i scałkowaniu otrzymu- jemy
− Z 2
0
u00vdx + Z 2
0
x2u0vdx + Z 2
0
4uvdx = Z 2
0
f vdx (90)
Następnie stosujemy wzór na całkowanie przez części dla pierwszego składnika i dla niego otrzymujemy
− Z 2
0
u00vdx = −u0v20+ Z 2
0
u0v0dx (91)
= −u0(2) v (2) + u0(0) v (0) + Z 2
0
u0v0dx (92)
Możemy podstawić z pierwszego warunku u0(0) = 2u(0) + 1 i otrzymujemy
− u0(2) v (2) + 2u (0) v (0) + v (0) + Z 2
0
u0v0dx (93)
Następnie włączamy drugi warunek Dirichleta sprowadzając go do warunku brzegowe- go zerowego za pomocą podstawienia u = w + ˆu, gdzie ˆu(2) = 0, w(2) = 1 i otrzymujemy
2w (0) v (0) + 2ˆu (0) v (0) + v (0) + Z 2
0
w0v0dx + Z 2
0
uˆ0v0dx (94) A więc całe równanie
2w (0) v (0) + 2ˆu (0) v (0) + v (0) + Z 2
0
w0v0dx + Z 2
0
uˆ0v0dx (95)
+ Z 2
0
x2w0vdx + Z 2
0
x2uˆ0vdx + Z 2
0
4wvdx + Z 2
0
4ˆuvdx = Z 2
0
f vdx (96) które można przedstawić w postaci (56). Należy mieć na uwadze przestrzenie do których należą funkcje ˆu i v, ˆu, v, w ∈ H1(Ω) z dodatkowymi warunkami ˆu(2) = 0, w(2) = 1.
1.3 Metoda Galerkina
Możemy rozwiązać zagadnienie wariacyjne za pomocą metody Galerkina, a więc w jed- nym z wariantów przybliżając funkcję u wyrażeniem:
u = ˜ˆ u +
n
X
i=1
αiφi (97)
Metoda jest przedstawiona dla przypadku jednowymiarowego. Funkcja ˜u odpowiada za spełnienie niezerowego warunku brzegowego Dirichleta. Funkcje φi powinny przyjmować wartość zero na brzegu obszaru na którym występują warunki Dirichleta. Funkcje φi nazywane są funkcjami bazowymi. Ponadto żądamy spełnienia równania wariacyjnego dla n funkcji v równych:
vi= φi (98)
dla i = 1, 2, . . . , n W ten sposób otrzymujemy układ n równań liniowych o n niewiado- mych αi postaci:
− Z
G
d u +˜ Pn
j=1
αjφj
!
dx
dφi
dxdx = Z
G
f φidx (99)
Po przekształceniu:
−
n
X
j=1
αj
Z
G
dφj dx
dφi dxdx −
Z
G
d˜u dx
dφi dxdx =
Z
G
f φidx (100)
Zauważmy, że zapis (64) jest tak naprawdę zapisem wymuszenia ortogonalności błę- du (defektu) dla równania (61), więc widzimy, że sformułowanie metody Galerkina z konspektu z metodami numerycznymi jest równoważne powyższemu.
Przykład 6.
−d2u
dx2 = 1 (101)
w przedziale (0, 1)
u (0) = u (1) = 0 (102)
Możemy wybrać funkcję ˜u = 0, ponieważ dla niej będą spełnione warunki brzegowe.
Ustalmy n = 3, wybierzmy funkcje bazowe
φ1 = x (103)
φ2= x2 (104)
φ3= x3 (105)
Konstruujemy układ równań (100).
1.4 Metoda elementów skończonych
W metodzie elementów skończonych dzielimy obszar na podobszary i znajdujemy funkcję u dla każdego podobszaru.
Przykład 7.
−d2u
dx2 = 1 (106)
w przedziale (0, 1)
u (0) = u (1) = 0 (107)
Do przybliżenia wybieramy takie funkcje φi, aby każda z nich przyjmowała wartość 1 w punkcie xi, a w pozostałych punktach podziału przyjmowała wartość 0, przykładowo dla x0 = 0, x1= 0, 5, x2 = 1:
φ1(x) =
( 1 − 2x; x ∈ (0, 0.5)
0; x ∈ (0.5, 1) (108)
φ2(x) =
( 2x; x ∈ (0, 0.5)
2 − 2x; x ∈ (0.5, 1) (109)
φ3(x) =
( 0; x ∈ (0, 0.5)
2x − 1; x ∈ (0.5, 1) (110)
u =ˆ
3
X
i=1
αiφi (111)
W każdym przedziale xk, xk+1 u będzie przybliżane linią prostą. W ogólności:
φk(x) =
x−xk−1
xk−xk−1; x ∈ (xk−1, xk)
xk+1−x
xk+1−xk; x ∈ (xk, xk+1) 0
(112)
np.: x0= 0, x1 = 0, 5, x2= 1,
φ0(x) =
( 1 − 2x x ∈ (0, 0.5)
0 (113)
φ1(x) =
( 2x; x ∈ (0, 0.5)
2 − 2x; x ∈ (0.5, 1) (114)
φ2(x) =
( 2x − 1; x ∈ (0.5, 1)
0 (115)
Przykład 8. Napisać układ Galerkina dla problemu brzegowego
− 2u00+ u = x + 2 (116)
dla x ∈ (0, 2) z warunkami
u (0) = 1 (117)
u (2) = 0 (118)
dla funkcji bazowych jednowymiarowej metody elementów skończonych v1, v2, v3 zwią- zanych z siatką punktów xk = k/2, k = 0, 1, 2, 3, 4, tzn.
vk(x) = 1 − |x − xk|
2 (119)
jeśli x ∈ [xk−1, xk+1], oraz
vk(x) = 0 (120)
w przeciwnym przypadku dla k = 1, 2, 3.
Należy najpierw naszkicować funkcje bazowe. Rozwiązanie będzie składało się z roz- wiązań dla podprzedziałów. Dla każdego podprzedziału musimy określić rozwiązanie postu- lowane które znajdujemy za pomocą odpowiedniego układu równań (100). Układ równań
musi być tak skonstruowany aby liczba równań była równa liczbie zmiennych i aby nie występowały równania nieokreślone postaci 0 = 0 lub sprzeczne. Dla pierwszego podprze- działu musi być spełniony warunek (117), więc dla tego podprzedziału rozwiązanie będzie miało postać
u = ˜u +
1
X
k=0
αkvk(x) (121)
Zauważmy, że dla tego podprzedziału nie możemy postulować rozwiązania z większą licz- bą funkcji bazowych, ponieważ otrzymujemy wtenczas jedno z równań z układu równań (100) postaci 0 = 0. Dlatego też rozpatrujemy tylko te funkcje bazowe, które są nieze- rowe w danym podprzedziale. W kolejnych podprzedziałach nie musimy dodawać funkcji u z wyjątkiem ostatniego podprzedziału dla którego musi być spełniony warunek (118).˜ Zauważmy, że ˜u zarówno dla pierwszego jak i ostatniego warunku wystarczy, że speł- nia tylko jeden odpowiedni warunek, więc postać tej funkcji może być odpowiednią stałą.
W jaki sposób wyznaczyć liczbę podprzedziałów? Na pewno będzie ona określona w ten sposób aby w każdym podprzedziale zestaw funkcji bazowych był ten sam, ponadto taki podprzedział może być dalej podzielony na kolejne podprzedziały ze względu na postać od- cinkową funkcji. Należy również sprawdzić, czy rzeczywiście dla pierwszego i ostatniego podprzedziału są spełnione warunki brzegowe, a więc przy obecności funkcji ˜u czy zerują się wartości funkcji bazowych vk(x).
2 Zadania
2.1 Zadania na 3.0
1. Wybrać dowolny problem brzegowy i rozwiązać.
2. Rozwiązać zagadnienie brzegowe metodą elementów skończonych:
− d dx
kdu
dx
= 0 (122)
z warunkiem Dirichleta w x = 0:
u (0) = uw (123)
i z warunkiem Cauchy’ego w x = 1 fdu
dx+ hu (1) = huz (124)
Wskazówki:
Wskazówki do R:
• pakiet fdaPDE Wskazówki do Matlaba:
• partial differential equation toolbox
• pdetool
• http://www.mathworks.com/help/pde/index.html Przykłady w R:
• przykład 1 triangulacja data(MeuseData)
mesh <- create.MESH.2D(nodes = MeuseData[,c(2,3)], order = 1) plot(mesh)
• przykład 2
data("mesh.2D.rectangular") plot(mesh.2D.rectangular)
FEMbasis = create.FEM.basis(mesh.2D.rectangular) coeff <-
sin(mesh.2D.rectangular$nodes[,1])*cos(mesh.2D.rectangular$nodes[,2]) FEM_object<- FEM(coeff, FEMbasis)
plot(FEM_object)
• przykład 3
data("mesh.2D.rectangular") plot(mesh.2D.rectangular)
FEMbasis = create.FEM.basis(mesh.2D.rectangular) coeff <-
sin(mesh.2D.rectangular$nodes[,1])*cos(mesh.2D.rectangular$nodes[,2]) FEM_object<- FEM(coeff, FEMbasis)
image(FEM_object)
• przykład 4
data(MeuseData) data(MeuseBorder)
mesh <- create.MESH.2D(nodes = MeuseData[,c(2,3)], segments = MeuseBorder, order = 1)
plot(mesh)
• przykład 5 Refine the triangulation data(MeuseData)
data(MeuseBorder)
mesh <- create.MESH.2D(nodes = MeuseData[,c(2,3)], segments = MeuseBorder, order = 1)
plot(mesh)
mesh_refine <- refine.MESH.2D(mesh, minimum_angle = 30, maximum_area = 10000)
plot(mesh_refine)
• przykład 6
library(fdaPDE) data(MeuseData) data(MeuseBorder) order=1
mesh <- create.MESH.2D(nodes = MeuseData[,c(2,3)], segments = MeuseBorder, order = order)
plot(mesh)
FEMbasis = create.FEM.basis(mesh) data = log(MeuseData[,"zinc"]) lambda = 10^3.5
ZincMeuse = smooth.FEM.basis(observations = data, FEMbasis = FEMbasis, lambda = lambda)
plot(ZincMeuse$fit.FEM)
desmat = matrix(1,nrow=nrow(MeuseData),ncol=2) desmat[,1] = sqrt(MeuseData[,"dist.log(m)"]) desmat[,2] = MeuseData[,"elev"]
smooth.FEM.PDE.basis
ZincMeuseCovar = smooth.FEM.basis(observations = data, covariates
= desmat,
FEMbasis = FEMbasis, lambda = lambda)
# Plot of the non parametric part (f) of the regression model y_i
= beta_1 x_i1 + beta_2 x_i2 + f plot(ZincMeuseCovar$fit.FEM)
# Print covariates’ regression coefficients print(ZincMeuseCovar$beta)
• przykład 7
# Load the mesh and plot it data(mesh.2D.simple)
plot(mesh.2D.simple)
# Create a vector of noisy samples of an underlying spatial field,
# located over the nodes of the mesh
observations = sin(pi*mesh.2D.simple$nodes[,1]) + rnorm(n = nrow(mesh.2D.simple$nodes), sd = 0.1)
# Create the FEM basis object
FEMbasis = create.FEM.basis(mesh.2D.simple)
# Set a vector of smoothing coefficients
lambda = c(10^-4, 1, 10^4)
# Set the anysotropic smoothing matrix K
PDE_parameters_anys = list(K = matrix(c(0.01,0,0,1), nrow = 2), b
= c(0,0), c = 0)
# Estimate one field for each smoothing parameter and plot these FEM_CPP_PDE = smooth.FEM.PDE.basis(observations = observations, FEMbasis = FEMbasis, lambda = lambda,
PDE_parameters = PDE_parameters_anys) plot(FEM_CPP_PDE$fit.FEM)
# Evaluate solution in three points
eval.FEM(FEM_CPP_PDE$fit.FEM, locations = rbind(c(0,0),c(0.5,0),c(-2,-2)))
• przykład 8
# Loading the mesh
data(mesh.2D.rectangular)
# Create the FEM basis object
FEMbasis = create.FEM.basis(mesh.2D.rectangular)
# Create a vector of noisy samples of an underlying spatial field,
# located over the nodes of the mesh
observations = sin(0.2*pi*mesh.2D.rectangular$nodes[,1]) + rnorm(n = nrow(mesh.2D.rectangular$nodes), sd = 0.1)
# Set the smoothing coefficient lambda = c(10^-2)
#Set the space vriant coefficients of the penalizying PDE K_func<-function(points)
{
mat<-c(0.01,0,0,1)
output = array(0, c(2, 2, nrow(points))) for (i in 1:nrow(points))
output[,,i] = 0.5*mat %*% t(points[i,1]^2) output
}
b_func<-function(points) {
output = array(0, c(2, nrow(points))) for (i in 1:nrow(points))
output[,i] = 0 output
}
c_func<-function(points) {
rep(c(0), nrow(points))
}
u_func<-function(points) {
rep(c(0), nrow(points)) }
# Assemble the parameters in one object
PDE_parameters = list(K = K_func, b = b_func, c = c_func, u = u_func)
# Estimate the underlying spatial field and plot these
FEM_CPP_PDE = smooth.FEM.PDE.sv.basis(observations = observations, FEMbasis = FEMbasis, lambda = lambda, PDE_parameters =
PDE_parameters)
plot(FEM_CPP_PDE$fit.FEM)
2.2 Zadania na 4.0
Rozwiązać metodą elementów skończonych równanie postaci:
− a (x) u0(x)0 = 0 (125)
dla x ∈ [0, 1] z warunkami:
u (0) = 0 (126)
u (1) = 5 (127)
oraz
a (x) = 1; x ∈ [0, 0.5] (128)
a (x) = 2; x ∈ [0.5, 1.0] (129)
2.3 Zadania na 5.0
Rozwiązać metodą elementów skończonych równanie postaci:
− d2u (x)
dx2 = 0 (130)
dla x ∈ [0, 2] z warunkami:
u (2) = 1 (131)
du (0)
dx + u (0) = 20 (132)
