Dr Piotr Fronczak
Metody numeryczne
Wykład nr 9
Równania róŜniczkowe zwyczajne - problemy brzegowe (BVP)
Dotychczas omawialiśmy problemy początkowe – rówania róŜniczkowe, w których dane były wartości zmiennych zaleŜnych (lub ich pochodne) dla pewnej szczególnej wartości zmiennej niezaleŜnej.
Teraz naszym zadaniem będzie wyznaczenie spośród funkcji spełniających
dane równanie róŜniczkowe zwyczajne, zdefiniowanych w rozwaŜanym obszarze, tych, które spełniają dodatkowe warunki na brzegu tego obszaru. Warunki takie nazywane są warunkami brzegowymi i są nałoŜone na wartości funkcji i jej
pochodnych w więcej niŜ jednym punkcie tego obszaru.
BVP są zwykle szczególnym przypadkiem równań róŜniczkowych cząstkowych, których rozwiązaniem są funkcje czasu i połoŜenia, np. pole elektryczne,
rozkład temperatury, prędkość przepływu itp.
2 2
x y t
y
∂
= ∂
∂
∂
równanie dyfuzji2 2 2
2
x y t
y
∂
= ∂
∂
∂
równanie faloweJeśli y = y(x) (nie zaleŜy od czasu – stany ustalone, równowagi), to otrzymujemy ogólną postać BVP (drugiego rzędu):
) ,
2
(
2
dx
y
dydx f y
d =
Dane jest równanie
w dziedzinie
a ≤ x ≤ b
oraz określone są w pewien sposób warunki brzegowe.
Typowe formy warunków brzegowych
Warunki brzegowe Dirichleta
Dwie wartości y(x) są dane – jedna dla x = a, druga dla x = b.
y(a) = Ya oraz y(b) = Yb
Warunki brzegowe Neumanna
Dwie wartości dy/dx są dane – jedna dla x = a, druga dla x = b.
a a
x
dx D
dy =
=
b b
x
dx D
dy =
=
oraz
Mieszane warunki brzegowe Robina
a a
x
C a
y dx c
c dy + =
=
)
2
(
1 b
b x
C b
y dx c
c dy + =
=
)
4
(
oraz 3 4
3 2
1
, c , c , c
c
- stałeMetoda strzałów
b x a x F y
x dx k
y
d
2+ ( ) = ( ); ≤ ≤
2
Rozpatrzmy ogólne równanie
W celu jego zdyskretyzowania przyjmijmy
N a
b
h = ( − ) /
gdzie N jest liczbą punktów, na które dzielimy przedział [a,b].
Dyskretyzując drugą pochodną
)
2 (
22
1 '' 1
h h O
y y
y
iy
i−
i+
i+
=
+ −i i
i i
i
i
k y F
h
y y
y
+− +
−+ =
2
1
1
2
otrzymujemy równanie
i i
i i
i
i
k y F
h
y y
y
+− +
−+ =
2
1
1
2
Chcemy scałkować to równanie od x0 = a do xN = b, więc przedstawmy je w postaci
i i
i i
i
i
y y h k y h F
y
+1= −
−1+ 2 −
2+
21 2 1
1 2 1
0
2
y 2 y h k y h F
y = − + − +
Czyli dla i = 1
Mamy dane y0 = y(a) = 0, ale y1 jest nieznane.
Znajomość y1 jest równoznaczna ze znajomością y’ dla x = 0:
h y y
iy
i−
i≈
+1'
0 0
1
y h y
y ≈ + ′
RozwaŜmy równanie struny zaczepionej na obu końcach.
Jedyną siłą działającą na element struny jest siła napręŜeń T. Jej wartość w kierunku pionowym
( sin( ) sin( ) )
) sin(
)
sin( T T
i 1 iT
F = θ + ∆ θ − θ = θ
+− θ
Zakładając, Ŝe kąty θ są małe
x y tg
iy
i ii
∆
= −
≈
+ ++
1 1
1
) ( )
sin( θ θ
oraz, Ŝe przyśpieszenie elementu struny jest proporcjonalne do wychylenia
y a = − ω
2dostajemy
∆ +
≈ −
−
⋅
∆
=
+ −x
y y
T y y
x
ma
2 i 12
i i 1) ( ω ρ
) 0 (
2
22
1
1
+ =
∆
+
−
−+
y
x
y y
T y
i i iρω
Dla
∆x → 0 Ty '' + ρω
2y = 0 0 '' + y =
y λ
T ρω
2λ =
22 2
2
x y t
y
∂
= ∂
∂
∂ λ
2 2
dt y a = d
Zakładając Ŝe
otrzymalibyśmy
0 )
( , 0 )
0 ( 0
'' + y = y = y L =
y λ
Gdy λ > 0, to istnieje rozwiązanie postaci
) sin(
)
cos( x B x
A
y = α + α
Uwzględniając warunki brzegowe otrzymujemy nieskończenie wiele rozwiązań
...
, 3 , 2 , 1 ,
sin )
( =
= n
L x x n
y
nπ
Gdy λ = 0
B Ax y = +
Uwzględniając warunki brzegowe otrzymujemy rozwiązanie trywialne y = 0.
z wartościami własnymi
...
, 3 , 2 , 1
2
,
2 2
=
= n
L n
n
λ π
PoniewaŜ problem fizyczny nie ma ujemnych wartości własnych, nie musimy analizować przypadku λ < 0.
Zatem mamy rozwiązanie
...
, 3 , 2 , 1 ,
sin )
( =
= n
L x x n
y
nπ
równania
0 ) ( , 0 ) 0 ( 0
''
22 2
=
=
=
+ y y y L
L y n π
Wybierzmy n = 4 i L = 1.
( 4 ) , ( 0 ) 0 , ( 1 ) 0
sin )
( x = x y = y =
y π
( x )
x
y ' ( ) = 4 π cos 4 π
Pochodna wynosi
=0 =0
) 16
2 (
4
2 22
π h h π
y = −
i i
i i
i
i
y y h k y h F
y
+1= −
−1+ 2 −
2+
2Oraz dalsze kroki zgodnie ze wzorem ze slajdu nr 6
) 4
cos(
4
00
1
y h x
y = + π π
1 2 1
1 2 1
0
2
y 2 y h k y h F
y = − + − +
A kolejne kroki rozwiązania
x yteor ynum 0 0.0000 0.0000 0.1 0.9511 1.2566 0.2 0.5878 0.5289 0.3 -0.5878 -1.0341 0.4 -0.9511 -0.9641 0.5 0.0000 0.6283 0.6 0.9511 1.2285 0.7 0.5878 -0.1113 0.8 -0.5878 -1.2753 0.9 -0.9511 -0.4255 1 0.0000 1.0963
Podzielmy obszar rozwiązań [0,1] na 10 równych części (h = 0.1)
x yteor ynum
0 0.0000 0.0000 0.05 0.5878 0.6283 0.1 0.9511 1.0086 0.15 0.9511 0.9907 0.2 0.5878 0.5817 0.25 0.0000 -0.0570 0.3 -0.5878 -0.6731 0.35 -0.9511 -1.0235 0.4 -0.9511 -0.9699 0.45 -0.5878 -0.5333 0.5 0.0000 0.1138 0.55 0.5878 0.7160 0.6 0.9511 1.0355 0.65 0.9511 0.9462 0.7 0.5878 0.4834 0.75 0.0000 -0.1703 0.8 -0.5878 -0.7567 0.85 -0.9511 -1.0444 0.9 -0.9511 -0.9198 0.95 -0.5878 -0.4321 1 0.0000 0.2262
-1.5000 -1.0000 -0.5000 0.0000 0.5000 1.0000 1.5000
0 0.2 0.4 0.6 0.8 1 1.2
yteor ynum -1.5000
-1.0000 -0.5000 0.0000 0.5000 1.0000 1.5000
0 0.2 0.4 0.6 0.8 1 1.2
yteor ynum
Podzielmy obszar rozwiązań [0,1] na 20 równych części (h = 0.05)
Zwykle nie znamy wartości i wektorów własnych układu. Trzeba je zgadnąć.
Przekształćmy ogólny problem brzegowy drugiego rzędu
b dx a
dy
a x b y a Y y b Y
y x dx f
y
d
2= ( , , ) ≤ ≤ ( ) = , ( ) =
2
z dla
w układ dwóch równań pierwszego rzędu:
Y
aa y dx w
dy = z war. pocz. ( ) =
) , ,
( x y w dx f
dw =
BEZ WARUNKU POCZĄTKOWEGOMusimy znaleźć warunek początkowy
α
=
=
=a
dx
xa dy w ( )
Innymi słowy – musimy znaleźć nachylenie αOK krzywej y w punkcie a.
Ya Yb
x=a x=b
α2 < αOK < α1
0 ) ( , 0 ) 0 ( 0
16
'' +
2y = y = y L =
y π
Wróćmy do naszego przykładu (dla n=4, L=1)
Przepiszmy to równanie w postaci układu dwóch równań pierwszego rzędu:
=
−
=
=
=
α ) 0 ( 16
0 ) 0 (
2
y w
π w'
y w
y'
α - parametr układu.
Musimy znaleźć miejsce zerowe funkcji błędu
0 )
1 ( )
1 ( )
( = y − y =
E α
αZwykle powyŜszy układ równań będziemy rozwiązywać jedną z metod podanych na poprzednim wykładzie (np. Rungego-Kutty), ale na razie, korzystając z metod analitycznych, zauwaŜmy Ŝe
Problem początkowy!
π π 4
) 4 ) sin(
( x
C x
y =
A zatem nie znajdziemy stałej C z warunku y(1) = 0.
0.2 0.4 0.6 0.8 1.0
-1.0 -0.5 0.5 1.0
Wybierzmy zatem na potrzeby dydaktyki inny warunek brzegowy:
587 .
0 )
7 . 0
( =
Metoda bisekcji.
y
100 100
1 0
=
−
= α α
-100 -50 50 100
-6 -4 -2 2 4
0 )
2 ( 1
1
0
+ =
= α α
α
śrśr
E
śrE ( α
0) ⋅ ( α ) < 0 ⇒ α
1= α 50
) 2 (
1
1
0
+ = −
= α α
α
śrśr
E
śrE ( α
0) ⋅ ( α ) > 0 ⇒ α
0= α
…
-10.5 -10.0 -9.5 -9.0 -0.06
-0.04 -0.02 0.02 0.04
2 .
− 10
śr
= α
Sprawdźmy:
17 . 10 )
7 . 0 4 cos(
4 )' 7 . 0 4
sin( π = π π = −
Metoda siecznych.
1 1
1
0
+
−
−
−
= −
−
−
i i
i i
i
i
i
E E
E
α α
α α
i i i
i i
i
i
E
E E −
− −
=
−
− +
1 1 1
)
( α α
α α
E=0
α E
E
i-1E
iα
i-1α
i369 . 6 100
194 . 5
100 1
1
−
=
=
=
−
= −
−
i i
i i
E E
α
α
16 . 10 )
369 . 6 369 (
. 6 194 . 5
) 100 100
100 (
1
− = −
+
−
− −
+
= α
i17 . 10 )
7 . 0 4 cos(
4 )' 7 . 0 4
sin( π = π π = −
Przypomnienie:
Wniosek: szybka zbieŜność juŜ po pierwszej iteracji.
Metoda róŜnic skończonych
W metodzie tej, pochodne w równaniu róŜniczkowym zastępujemy róŜnicami skończonymi.
Dziedzinę rozwiązań [a,b] dzielimy na N przedziałów o długości h = (b-a) / N.
Mamy N+1 punktów. Dla kaŜdego z nich zapisujemy równanie róŜnicowe – czyli równanie algebraiczne.
Mamy zatem układ równań algebraicznych, który rozwiązujemy jedną z metod omówionych na wykładzie nr 3.
x =a1 xi-1 xi xi+1 x =bN+1 h
x y
Przykład:
0 ) 1 ( , 0 ) 0 ( ,
2
1
2
=
=
= y y
dx y d
Drugą pochodną moŜna przybliŜyć za pomocą trójpunktowych róŜnic centralnych (zwykle, choć niekoniecznie). Mamy zatem
Rozwiązanie analityczne:
2 ) 2
(
2
x
x x
y = −
By otrzymać rozwiązanie numeryczne najpierw dyskretyzujemy równanie dla punktów [x0, x1, …, xN], gdzie x0 = 0, xN = 1, xi = ih.
Wybierzmy h = 0.2 (N = 5).
4 , 3 , 2 , 1 ,
2 1
2
1
1
− +
+= =
−
i
h
y y
y
i i idla
Zwróćmy uwagę, Ŝe są to równania tylko dla punktów wewnętrznych.
Układ nasz ma 4 niewiadome i 4 równania:
. h y
y - y
, h
y
y - y
, h
y
y - y
, h
y y
- y
5 2 4
3
4 2 3
2
3 2 2
1
2 2 1
0
2 2
2 2
= +
= +
= +
= +
Macierz powyŜsza jest przykładem macierzy trójdiagonalnej.
Choć powyŜszy układ moŜna rozwiązać jedną z metod omówionych na
wykładzie 3 (np. metodą Gaussa), to szczególna postać tej macierzy pozwala zredukować liczbę obliczeń z n3 do n.
−
−
=
−
−
−
−
5 2
2 2
0 2
4 3 2 1
2 1
1 2
1
1 2
1
1 2
y h
h h
y h
y y y y
=
04 . 0
04 . 0
04 . 0
04
.
0
. y
y
, y
y y
,
y y
y
,
y
y
100 4 4
3
100 4 4
3 2
100 3 4
2 1
100 2 4
1
2 2
2 2
=
−
= +
−
= +
−
= +
−
. y
y
, y
y y
,
y y
,
y
y
100 4 4
3
100 4 4
3 2
100 3 6
2 2 3
100 2 2
2 1 1
2 2
=
−
= +
−
= +
−
−
=
−
. y
y
, y
y
,
y y
,
y
y
100 4 4
3
100 4 8
3 3 4
100 3 4
3 2 2
100 2 2
2 1 1
2 =
−
= +
−
−
=
−
−
=
−
. y
, y
y
,
y y
,
y
y
10 4 1
4 5
100 4 6
4 3 3
100 3 4
3 2 2
100 2 2
2 1 1
=
−
−
=
−
−
=
−
−
=
−
. )
(
) (
) (
100 8 100
12 2
1 100 1 2
100 12 100
12 3
2 100 2 4
100 12 100
8 4
3 100 3 6
100 4 8
−
=
− +
−
=
−
=
− +
−
=
−
=
− +
−
=
−
=
y
, y
, y
, y
0.0 0.2 0.4 0.6 0.8 1.0 -0.14
-0.12 -0.10 -0.08 -0.06 -0.04 -0.02 0.00
y
x
Algorytm Thomasa dla układów z macierzą trójdiagonalną
Dany jest układ równań:
n i
Y y
c b
y
a
i i−1+
i+
i i+1=
idla = 1 , 2 , ...,
=
−
n n
n n
n
Y Y Y
y y y
b a
c c
b a
c b
a
c b
. . .
. 0
...
...
0
...
...
...
...
...
0 ...
0
0 ...
0
0 ...
...
0
2 1 2
1
1 3
3 3
2 2
2
1 1
lub w postaci macierzowej
a
1= 0 , c
n= 0
Algorytm składa się z dwóch faz:
Faza eliminacji wprzód – wykonując dla równań od i = 1 do n eliminację
niewiadomych uzyskujemy ostatnie równanie (i = n) z jedną tylko niewiadomą, którą moŜemy wyznaczyć.
Faza eliminacji wstecz – korzystając z wyznaczonej w równaniu i+1 niewiadomej yi+1
wyznaczamy z równania i niewiadomą yi, aŜ do otrzymania wartości y1.
(∗) (∗) (∗)
(∗)
A zatem szukamy schematu postaci:
i i
i
i
y
y
−1= γ + β
Podstawiając ten schemat do równania (∗) otrzymujemy:
(
i i i)
i i i i ii
y b y c y Y
a γ + β + +
+1=
Wyznaczając yi
i i
i
i i i
i i i
i
i
i
a b
a y Y
b a
y c
+ + −
+
= −
+γ
β
γ
1Porównując z otrzymujemy
i i
i
i i i
i i
i i
i
i
a b
a Y
b a
c
+
= − +
= −
++
γ
β β
γ
1γ ,
1Otrzymaliśmy równanie rekurencyjne na poszukiwane współczynniki γ i β.
Współczynniki początkowe γ1 i β1 nie mają znaczenia, bo mnoŜone są przez a1=0.
Musimy jeszcze znać yn, by móc rozpocząć iteracyjne obliczanie niewiadomych.
Podstawiając pierwsze równanie schematu
do równania (∗) otrzymujemy:
n n
n
n
y
y
−1= γ + β
(
n n n)
n n nn
y b y Y
a γ + β + =
Czyli
n n
n
n n n
n
a b
a y Y
+
= −
γ
β
i i
i
i i i
i
a b
a Y
+
= −
+
γ
β
1β
Warto zauwaŜyć, Ŝe skoro to wystarczy przyjąć
(∗∗)
1
= 0
+
y
n(∗∗∗)
, by móc bezpośrednio skorzystać ze wzoru (∗∗).
Podsumowując:
i i
i 1
- i
1 n
i i
i
i i i
1 i i
i i
i 1
i 1 1
x y
1..2 n
i for
0;
y
b a
a , Y
b a
c 1..n
i for
przykład //na
0;
β γ
γ β β
γ γ β γ
+
=
+
=
=
+
= − +
= −
=
=
=
+
+ +
Algorytm Thomasa jest niezawodny, gdy macierz jest diagonalnie dominująca
n i
c a
b
i≥
i+
i= 1 , ...,
Przykład: radiator prętowy
Równanie opisujące rozkład temperatury wzdłuŜ długości pręta:
TA
TB TS
L
x
L x T
kA T P h dx
T d
S C
C
− = ≤ ≤
− ( ) 0 , 0
2 2
hc – współczynnik wnikania ciepła P – obwód pręta
k – współczynnik przewodzenia ciepła
Ac – pole poprzecznego przekroju pręta T 293K K 293 T(L)
K 473 T(0)
m 0.1 L
m 10 1.6 A
W/m/K 240
k
m 0.016 P
/K W/m 40 h
S
2 5 c
2 c
=
=
=
=
×
=
=
=
=
−
0 )
2
(
2
=
−
− T T
Sdx
T
d β
−1− 2
2+
+1− (
i−
S) = 0
ii
i
T T
h
T T
T β
S i
i
i
h T T h T
T
−1− ( 2 +
2β ) +
+1= −
2β
Podzielmy domenę rozwiązań na 5 części (h = L / 5 = 2 cm)
T
Sh T
T h
T
i = 2
1− ( 2 +
2β )
2+
3= −
2β ) (
) 2
( + h
2T
2+ T
3= − h
2T
S+ T
1− β β
czyli
T
Sh T
T h
T
i = 3
2− ( 2 +
2β )
3+
4= −
2β
T
Sh T
T h
T
i = 4
3− ( 2 +
2β )
4+
5= −
2β
T
Sh T
T h
T
i = 5
4− ( 2 +
2β )
5+
6= −
2β ) (
) 2
( + h
2T
5+ T
4= − h
2T
S+ T
6− β β
czyli
+
−
−
−
+
−
=
+
− +
− +
− +
−
) (
) (
) 2
( 1
0 0
1 )
2 ( 1
0
0 1
) 2
( 1
0 0
1 )
2 (
6 2
2 2
1 2
5 4 3 2
2 2
2 2
T T
h
T h
T h
T T
h
T T T T
h h
h h
S S S S
β β β β
β β
β β
Mamy zatem układ 4 równań z czterema niewiadomymi.
0.00 0.02 0.04 0.06 0.08 0.10
300 350 400 450 500
Rozwiązując ten układ
za pomocą algorytmu Thomasa otrzymujemy:
Metoda róŜnic skończonych dla nieliniowych równań róŜniczkowych
Najciekawsze problemy współczesnej fizyki są nieliniowe.
Nieliniowe problemy brzegowe dyskretyzujemy w podobny sposób.
Wynikiem jest jednak układ nieliniowych równań algebraicznych.
Metody rozwiązywania takich równań nieliniowych omówiliśmy na wykładzie nr 2.
Najbardziej wydajne obliczeniowo są w tym przypadku metody iteracyjne.
Istnieje jednak potencjalny problem związany ze zbieŜnością schematu iteracyjnego.
Metoda punktu stałego
Układ równań nieliniowych moŜna zapisać w postaci
] [ ] [ ] ][
[ a y + Φ = b
[a] – macierz współczynników
[Φ] – wektor nieliniowych wyrazów będących funkcją niewiadomych yi [b] - wektor znanych wielkości stałych
Spośród wielu sposobów konstruowania procedury iteracyjnej wybierzmy najprostszy
] [ ] [ ] ][
[ a y + Φ = b
k
k
b
y
a ][ ] [ ] [ ]
[
+1= − Φ
obliczone na podstawie wcześniejszego kroku k
Jeśli liczba punktów jest mała moŜemy macierz [a] odwrócić. Jeśli nie, moŜemy
skorzystać z metody eliminacji Gaussa albo Thomasa (dla macierzy trójdiagonalnej)
(
k)
k
a b
y ] [ ] [ ] [ ]
[
+1=
−1− Φ
Przykład: radiator prętowy
Gdy uwzględnimy wyraz odpowiedzialny na wypromieniowywanie ciepła,
równanie opisujące rozkład temperatury wzdłuŜ długości pręta uzyska postać:
L x T
kA T T P
kA T P h dx
T d
S C
S C
C
− − − = ≤ ≤
− ( ) (
4 4) 0 , 0
2
2
εσ
ε - względna zdolność emisyjna σ - stała Stefana-Boltzmanna
0 ) (
)
2 (
4 42
1
1
− +
+− − − − =
−
S i
B S
i A i
i
i
T T T T
h
T T
T β β
Równanie zdyskretyzowane:
) (
) 2
(
2 2 4 1 2 41 A i B i i A S B S
i
h T h T T h T T
T
−− + β − β +
+= − β + β
=
−
−
−
− +
+
− +
− +
− +
−
) (
) (
) (
) (
) 2
( 1
0 0
1 )
2 ( 1
0
0 1
) 2
( 1
0 0
1 )
2 (
4 5 2
4 4 2
4 3 2
4 2 2
5 4 3 2
2 2
2 2
T h
T h
T h
T h
T T T T
h h
h h
B B B B
A A
A A
β β β β
β β
β β
− +
−
+
−
+
−
− +
−
=
6 4
2
4 2
4 2
1 4
2
) (
) (
) (
) (
T T
T h
T T
h
T T
h
T T
T h
S B S
A
S B S
A
S B S
A
S B S
A
β β
β β
β β
β β
(
k)
k
a b
T ] [ ] [ ] [ ] [
+1=
−1− Φ
Teraz stosujemy procedurę iteracyjną:
0.00 0.02 0.04 0.06 0.08 0.10 280
320 360 400 440 480
T [K]
x [m]
krok 0 krok 1 krok 2 krok 3
0 0.02 0.04 0.06 0.08 0.1
krok 0 473 400 400 400 400 293
krok 1 473 423.2293 382.8297 349.1078 319.8155 293 krok 2 473 423.3492 383.3225 349.8507 320.4519 293 krok 3 473 423.344 383.3132 349.8409 320.4456 293
• Zakładamy hx = hy= h [siatka kwadratowa]
• u(xi,yj) = ui,j
• u(xi+h, yj+h) = ui+1,j+1,
x
ix
i+1x
i-1y
jy
j-1y
j+1(
i j)
x[ u
i ju
i ju
i j]
y h x
u 1
2 1,2
, 1,,
'' ≈
−− +
+2
0
2 2
2
∂ = + ∂
∂
∂
y u x
u
u
i,jhx
hy
Dwuwymiarowe zagadnienie brzegowe
Równanie Laplace’a
( ) 1
2[
, 12
, , 1]
,
'' ≈
i j−−
i j+
i j+j y
i
u u u
y h x u
( ) ( )
{ } { }
[ ]
[ 4 ] 0
1
2 1 2
, '' ,
''
1 , ,
1 ,
1 , ,
2 1
1 , ,
1 , ,
1 ,
, 2 1
= +
+
− +
= +
− +
+
−
≈ +
+ +
−
−
+
− +
−
j i j
i j
i j
i j
i
j i j
i j
i j
i j
i j
i
j y x i
j i
u u
u u
h u
u u
u u
u h u
y x u y
x
u
Sumując wyrazy wyznaczamy ui,j: ,
{
1, , 1 1, , 1}
4 1
+ +
−
−
+ + +
=
i j i j i j i jj
i
u u u u
u
Przyjmujemy początkowe przybliŜenie
1 2
3 4
u = 0
u = 0 u = 0
u = 1
Przykład: siatka 3x3
) 1
0 (
) 0 1 (
) 0
0 (
) 0 0
(
3 4 2
1 1 4
4 4 1
1 1 3
1 4 4
1 1 2
3 4 2
1 1 1
m m
m
m m
m
m m
m
m m
m
u u
u
u u
u
u u
u
u u
u
+ + +
=
+ + +
=
+ +
+
=
+ +
+
=
+ + + +
1 0
4 1
3 1
2 1
1 =u =u =u = u
Korzystamy z metody iteracyjnej Jacobiego (lub Gaussa-Seidla – szybsza zbieŜność)
1 2
3
4 S1
S2 S3
S4
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
Korzystając z metod dokładnych (np. dekompozycja LU) musimy ułoŜyć macierz o rozmiarze liniowym n x n.
Wniosek: Musimy wprowadzić indeksowanie równań odpowiadających punktom dwuwymiarowej siatki n x n:
...,
2, 2 ,
1 n
P =
(na przykład wierszami)
) ,
2
(
2 2
2
y x y f
u x
u = −
∂ + ∂
∂
Przykład: równanie Poissona
∂
[ u
i ju
i ju
i ju
i ju
i j] f
i jh 1
2 1, , 14
, 1, , 1 ,=
−
− +
−
−
− − + +P k
j , ) → (
j n k
P = ( − 1 ) ⋅ +
[ u
Pu
P nu
Pu
Pu
P n] f
Ph 1
2−
−1−
−+ 4 −
+1−
+= ε
0ϕ = ρ
∆
potencjał
gęstość ładunku
Przykład:
Powierzchnia potencjału przy losowo rozmieszczonej gęstości ładunku