Przykład:
Dana jest linia długa o długości L
0
bez strat o stałych
kilometrycznych L,C.Na początku linii zostaje załączona
siła elektromotoryczna e(t), a koniec linii jest zwarty. Określić
rozkład prądów i napięć w linii.
Równania linii długiej bez strat są:
t
,
x
,
t
,
x
,
Cu
i
Li
u
Różniczkując pierwsze równanie po x, a drugie po t i eliminując
prąd otrzymujemy równanie hiperboliczne (falowe) dla napięcia:
tt
,
2
xx
,
u
v
1
u
tt
,
2
xx
,
u
v
1
u
gdzie
LC
1
v
- prędkość propagacji fali w linii.
Warunki brzegowe zapisane dla napięcia są:
0
L
x
0
x
0
t
,
x
u
t
e
t
,
x
u
Warunki początkowe wynikają z założenia, że linia była
rozładowana przed załączeniem, czyli:
0
t
t
,
0
t
0
t
,
x
u
0
t
,
x
u
Jeszcze raz stosujemy przekształcenie Laplace’a i mamy:
0
u
v
s
u
2
xx
,
którego rozwiązanie jest:
v
sx
sin
A
v
sx
cos
A
s
,
x
u
2
1
gdzie stałe A
1
i A
2
wyznaczamy z warunków brzegowych:
0
L
x
0
x
0
s
,
x
u
s
E
s
,
x
u
gdzie E(s)=L[e(t)].
Otrzymujemy:
v
sL
ctg
s
E
A
s
E
A
0
2
1
Rozwiązanie dla transformaty napięcia jest:
v
sL
sin
x
L
v
s
sin
s
E
s
,
x
u
0
0
Korzystając jak poprzednio z twierdzenia o splocie znajdujemy
transformatę odwrotną L
-1
[E(s)]=e(t) i
v
sL
sin
x
L
v
s
sin
L
0
0
1
którą obliczamy korzystając z twierdzenia Cauchy.
Całka po linii zamkniętej z funkcji analitycznej równa się
sumie residuów funkcji znajdujących się wewnątrz krzywej
całkowania pomnożonej przez 2i.
Funkcja
v
sL
sin
x
L
v
s
sin
0
0
ma bieguny w punktach s
k
zerowania się
mianownika, tj.
v
L
k
s
0
k
gdzie k=0,1, 2,.....
Obliczając residuum mamy:
k
k
0
0
0
0
0
0
1
L
vt
ik
exp
k
cos
v
L
L
x
L
k
sin
v
sL
sin
x
L
v
s
sin
L
Rozdzielając sumę dla dodatnich i ujemnych k i zmieniając
znak k otrzymujemy ostatecznie:
1
k
0
0
0
k
0
0
0
1
L
vt
k
cos
x
L
L
k
sin
1
v
L
2
v
sL
sin
x
L
v
s
sin
L
Rozkład napięcia w linii określa wzór:
1
k
t
0
0
0
0
k
0
d
L
v
k
cos
t
e
L
x
L
k
sin
1
L
v
2
t
,
x
u
Prąd w linii można otrzymać wykorzystując równania linii
np.:
x
t
,
u
L
1
,
i
i mamy:
1
k
t
0
0
0
0
k
0
d
L
v
k
sin
t
e
L
x
L
k
sin
1
LL
2
t
,
x
i
Jako ostatni przykład zastosowania metody rozdzielania zmiennych
rozważmy:
Dana jest cienka prostokątna axb membrana zamocowana na brzegu.
Na powierzchnię membrany działa ciśnienie p(x,y,t). Wyznaczyć
drgania membrany, jeżeli w chwili t=0 była ona nieodkształcona
i nieruchoma. Grubość membrany wynosi h, jej masa właściwa .
x
y
w
a
b
dy
dx
w
x
dx
dy
Nw
x
,
dy
dx
w
w
N
xx
,
x
,
dxdy
y
,
x
p
N
N
dxdy
w
h
dB
w
x
dx
dy
Nw
x
,
dy
dx
w
w
N
xx
,
x
,
dxdy
y
,
x
p
N
N
dxdy
w
h
dB
Bilansując siły na oś w otrzymujemy:
0
pdxdy
dB
dx
Nw
dx
dy
w
w
N
dy
Nw
dy
dx
w
w
N
y
,
yy
,
y
,
x
,
xx
,
x
,
Wykonując redukcję i dzieląc przez hdxdy otrzymujemy równanie
opisujące odkształcenie membrany:
t
,
y
,
x
p
h
1
w
w
2
2
gdzie
h
N
2
Warunki brzegowe wyrażające fakt umocowania brzegu mają
postać:
b
y
0
y
a
x
0
x
0
w
;
0
w
0
w
;
0
w
Warunki początkowe są:
0
t
0
t
0
w
0
w
Przyjmując rozwiązanie w postaci:
t
T
y
Y
x
X
t
,
y
,
x
w
otrzymujemy dla równania jednorodnego po podzieleniu przez XYT:
0
dy
Y
d
Y
1
dx
X
d
X
1
dt
T
d
T
1
2
2
2
2
2
2
2
Przyjmujemy:
2
2
2
2
2
2
dy
Y
d
Y
1
dx
X
d
X
1
Rozwiązując powyższe równania mamy:
y
cos
D
y
sin
C
y
Y
x
cos
B
x
sin
A
x
X
0
y
0
x
0
w
;
0
w
Na mocy warunków brzegowych:
mamy:
0
C
;
0
B
Natomiast stałe , określamy z warunków brzegowych:
b
y
a
x
0
w
;
0
w
i mamy:
b
k
;
a
n
k
n
gdzie n,k=0,1,2,...
Rozwiązanie spełniające warunki brzegowe:
b
y
0
y
a
x
0
x
0
w
;
0
w
0
w
;
0
w
ma postać:
1
n
1
k
nk
b
y
k
sin
a
x
n
sin
t
T
t
,
y
,
x
w
Podstawiając do równania
t
,
y
,
x
p
h
1
w
w
2
2
otrzymujemy:
t
,
y
,
x
p
h
1
b
y
k
sin
a
x
n
sin
b
k
a
n
T
T
1
k
1
n
2
2
nk
2
nk
Funkcje
b
y
k
sin
a
x
n
sin
tworzą ortogonalny ciąg na prostokącie [0,a]x[0,b]. Mnożąc
równanie:
t
,
y
,
x
p
h
1
b
y
k
sin
a
x
n
sin
b
k
a
n
T
T
1
k
1
n
2
2
nk
2
nk
obustronnie przez
b
y
k
sin
a
x
n
sin
i całkując po powierzchni prostokąta [0,a]x[0,b] otrzymujemy:
t
P
T
T
nk
nk
2
nk
nk
t
P
T
T
nk
nk
2
nk
nk
gdzie
2
2
2
2
nk
b
k
a
n
a
0
b
0
nk
dxdy
b
y
k
sin
a
x
n
sin
t
,
y
,
x
p
abh
1
t
P
Biorąc pod uwagę warunki początkowe:
0
t
0
t
0
w
0
w
rozwiązanie równania można zapisać w postaci:
t
0
nk
nk
nk
nk
d
sin
t
P
t
T
Ugięcie membrany w pod wpływem ciśnienia p jest opisane
wzorem:
1
n
1
k
t
0
nk
nk
nk
b
y
k
sin
a
x
n
sin
d
sin
t
P
t
,
y
,
x
w
gdzie
2
2
2
2
nk
b
k
a
n
a
0
b
0
nk
dxdy
b
y
k
sin
a
x
n
sin
t
,
y
,
x
p
abh
1
t
P
Metody numeryczne rozwiązywania równań
różniczkowych cząstkowych
Metoda różnic skończonych (siatek)
Uwagi ogólne
Dane równanie różniczkowe cząstkowe opisane operatorem L:
f
Lu
w obszarze i warunki brzegowe:
1,2,...p
i
dla
u
l
i
x
i
i
X
k
W metodach różnicowych poszukuje się tablicy wartości
przybliżonych u
h
rozwiązania dokładnego u na zbiorze
izolowanych punktów X
k
(k=1,2,...,N
h
)
zwanym siatką.
Punkty X
k
są nazywane węzłami siatki.
węzeł pomocniczy
węzeł podstawowy
Równania służące do wyznaczania wartości przybliżonych
nazywamy równaniami różnicowymi.
x
y
h
x
h
y
h=(h
x
,h
y
)
Parametr h
charakteryzuje
siatkę
h
Dla równania różniczkowego:
f
Lu
w obszarze z warunkami brzegowymi:
1,2,...p
i
dla
u
l
i
x
i
i
otrzymujemy jego odpowiednik różnicowy:
ih
h
ih
h
h
h
u
r
f
u
R
Zakładając, że problem opisany równaniem różniczkowym ma
jednoznaczne rozwiązanie, to równania różnicowe będą jego
odpowiednikiem jeżeli są spełnione następujące warunki:
ih
h
ih
h
h
h
u
r
f
u
R
1. Układ równań różnicowych posiada jednoznaczne rozwiązanie:
h
h
k
h
k
h
X
dla
X
u
dla każdego dopuszczalnego h.
2. Zbieżność do rozwiązania dokładnego u.
Oznacza to, że rozwiązanie u
h
powinno przy h 0 dążyć do
rozwiązania dokładnego u.
Dla określenia zbieżności jest koniecznym wprowadzenie
odpowiednich przestrzeni funkcyjnych i norm w nich.
Wprowadzamy przestrzeń funkcyjną U z normą || ||
U
, do której
należy rozwiązanie dokładne u.
oraz przestrzeń N
h
- wymiarową U
h
z normą || ||
Uh
,
której elementami są układy N
h
liczb i do której
należy rozwiązanie u
h
.
Normy || ||
U
i || ||
Uh
winny być zgodne, tzn. ponieważ funkcja u(x)
jest określona w węzłach podstawowych X
k
siatki, to mówimy,
że normy są zgodne jeżeli zachodzi:
U
Uh
0
h
x
u
x
u
lim
dla każdego uU.
Przykłady norm zgodnych:
d
f
f
2
L
2
0
h
h
2
j
,
i
2
ij
2
L
f
h
f
- zbiór węzłów wewnętrznych.
0
h
d
,
u
u
u
i
2
i
2
H
1
2
h
ij
1
h
2
j
,
i
2
ij
1
j
,
i
2
ij
j
,
1
i
2
2
H
h
h
u
u
h
u
u
u
h
u
Wielkość:
Uh
h
h
x
u
u
nazywamy błędem rozwiązania przybliżonego u
h
.
U
h
dąży do rozwiązania dokładnego u(x), jeżeli
0
x
u
u
lim
Uh
h
0
h
Jeżeli można znaleźć taką funkcję (h), że
h
x
u
u
Uh
h
to mówimy, że zostało znalezione oszacowanie błędu.
3. Stabilność
ih
h
ih
h
h
h
u
r
f
u
R
Różnicowe zagadnienie brzegowe jest stabilne, jeżeli istnieją
takie liczby >0 i h
0
>0, że dla dowolnych h<h
0
i f
h
F
h
, takich,
że ||f
h
||
Fh
< zadanie brzegowe:
ih
ih
h
ih
h
h
h
h
z
r
f
f
z
R
posiada jedno jednoznaczne rozwiązanie i zachodzi:
Fh
h
Uh
h
h
f
C
u
z
gdzie C stała niezależna od h.
Twierdzenie wiążące stabilność i zbieżność:
Jeżeli zadanie różnicowe
ih
h
ih
h
h
h
u
r
f
u
R
jest przybliżeniem różniczkowego problemu brzegowego:
f
Lu
1,2,...p
i
dla
u
l
i
x
i
i
i rozwiązanie u
h
jest stabilne wtedy zachodzi
0
x
u
u
lim
Uh
h
0
h
i rząd zbieżności w funkcji h jest taki sam jak rząd aproksymacji.
Zastępowanie pochodnych ilorazami różnicowymi
na siatce prostokątnej
i-1
i
i+1
k-1
k
k+1
24
h
x
u
6
h
x
u
2
h
x
u
h
x
u
u
u
4
i
ik
4
4
3
i
ik
3
3
2
i
ik
2
2
i
ik
ik
k
,
1
i
24
h
x
u
6
h
x
u
2
h
x
u
h
x
u
u
u
4
i
ik
4
4
3
i
ik
3
3
2
i
ik
2
2
i
ik
ik
k
,
1
i
h
i
h
k
2
k
2
k
1
k
,
i
1
k
,
i
2
i
2
i
k
,
1
i
k
,
1
i
h
O
h
2
u
u
y
u
h
O
h
2
u
u
x
u
i-1
i
i+1
k-1
k
k+1
h
i
h
k
Druga pochodna dodając
stronami:
24
h
x
u
6
h
x
u
2
h
x
u
h
x
u
u
u
4
i
ik
4
4
3
i
ik
3
3
2
i
ik
2
2
i
ik
ik
k
,
1
i
24
h
x
u
6
h
x
u
2
h
x
u
h
x
u
u
u
4
i
ik
4
4
3
i
ik
3
3
2
i
ik
2
2
i
ik
ik
k
,
1
i
2
k
2
k
1
k
,
i
ik
1
k
,
i
2
2
2
i
2
i
k
,
1
i
ik
k
,
1
i
2
2
h
O
h
u
u
2
u
y
u
h
O
h
u
u
2
u
x
u
i-1
i
i+1
k-1
k
k+1
h
i
h
k
Dla pochodnej mieszanej
i
3
k
k
3
i
k
i
1
k
,
1
i
1
k
,
1
i
1
k
,
1
i
1
k
,
1
i
2
h
h
,
h
h
max
O
h
h
4
u
u
u
u
y
x
u
Konstrukcja warunków brzegowych na siatce
i
k
k
,
i
h
i
k
i
ik
x
k
ik
x
1. Przeniesienie wartości:
Przyjmujemy:
i
k
k
ik
k
i
i
ik
h
dla
x
dla
x
k
,
i
r(i,k,)
lub
,
k
,
i
r
min
x
h
x
k
,
i
2. Interpolacja liniowa dla warunku brzegowego Dirichleta
i,k
i-1,k
i+1,k
h
Zakładając liniowy rozkład rozwiązania między sąsiednimi
węzłami mamy:
)
x
(
u
h
x
x
u
h
x
x
u
1
i
ik
i
k
,
1
i
i dla x=x
i
+ mamy:
i
ik
k
,
1
i
h
1
u
h
u
i
3. Interpolacja liniowa dla warunku brzegowego Neumanna.
n
i,k
i-1,k
i,k-1
h
k
h
i
k
,
i
sin
h
u
u
cos
h
u
u
k
1
k
,
i
ik
i
k
,
1
i
ik
y
,
x
y
,
x
u
y
,
x
d
n
u
y
,
x
c
Równania eliptyczne
Dla uproszczenia rozważań będą analizowane przypadki
dwuwymiarowe x=(x,y)
y
,
x
f
u
y
,
x
g
y
u
y
,
x
b
y
x
u
y
,
x
a
x
Warunki brzegowe:
1
2
3
4
5
6
7
8
9 x(i)
1
2
3
4
5
6
7
8
9
10 y(k)
h
y
h
x
(0,0)
j
j+1
j-1
l
p
y
,
x
f
u
y
,
x
g
y
u
y
,
x
b
y
x
u
y
,
x
a
x
Dla węzłów wewnętrznych będzie:
j
j+1
j-1
l
p
j
j
j
y
l
p
y
l
p
2
y
p
j
l
j
x
1
j
1
j
x
1
j
1
j
2
x
1
j
j
1
j
j
f
u
g
h
2
u
u
h
2
b
b
h
u
u
2
u
b
h
2
u
u
h
2
a
a
h
u
u
2
u
a
lub w formie macierzowej:
f
u
A
w
h
y
h
x
(0,0)
j
j+1
j-1
l
p
styczna
normalna
y
,
x
y
,
x
u
y
,
x
d
n
u
y
,
x
c
p
Przyjmując:
(xk
p
,yk
p
)
p
p
p
p
p
p
p
p
p
yk
,
xk
d
d
yk
,
xk
c
c
yk
,
xk
p-1
m
p
otrzymujemy:
p
p
p
p
y
p
m
p
x
1
p
p
p
u
d
sin
h
u
u
cos
h
u
u
c
Uwaga dotycząca błędu obliczeń.
Generalnie jeżeli węzły nie leżą na krzywej brzegowej
i liczymy metodą przeniesienia wartości, to dokładność
obliczeń jest rzędu h. Jeżeli węzły na krzywej brzegowej
bądź wyliczamy wartości funkcji brzegowej interpolując
liniowo dokładność wzrasta do h
2
.
W zagadnieniu Dirichleta oprócz trudności z wyznaczeniem
wartości brzegowych nie ma innych problemów i otrzymany
układ równań algebraicznych najczęściej można rozwiązać
bez kłopotów.
Sytuacja może się komplikować przy zagadnieniu Neumanna.
Siatki praktycznie nie stosowane w zagadnieniach
eliptycznych liniowych i nieliniowych.
Równania opisujące ewolucję układu w czasie
Równania paraboliczne
Dane jednowymiarowe równanie przewodnictwa:
T
0,
t
dla
t
f
t
1,
u
t
f
t
,
0
u
0,1
x
dla
x
0
,
x
u
T
0,
t
0,1
x
dla
t
,
x
f
,
u
,
u
2
1
0
xx
t
x
t
0
1
h
i-1 i i+1
k+1
k
k-1
x
t
0
N
h
i-1 i i+1
k+1
k
k-1
Oznaczamy operator różnicy II rzędu:
2
k
1
i
k
i
k
1
i
k
i
h
u
u
2
u
u
i wprowadzamy schematy różnicowe z wagą :
k
i
k
i
1
k
i
k
i
1
k
i
u
1
u
u
u
N
,
0
i
i
K
,
0
k
5
.
0
k
,
x
f
i
k
i
T
0,
t
dla
t
f
t
1,
u
t
f
t
,
0
u
0,1
x
dla
x
0
,
x
u
T
0,
t
0,1
x
dla
t
,
x
f
,
u
,
u
2
1
0
xx
t
k
i
k
i
1
k
i
k
i
1
k
i
u
1
u
u
u
Problem brzegowy jest aproksymowany przez
i
0
i
ih
u
dla i=1,2,...N
k
2
k
N
k
1
k
0
2
1
f
k
f
u
f
k
f
u
dla k=1,2,...,K.
Warunki zgodności
k+1
k
i-1 i i+1
Schemat sześciowęzłowy
k
i
k
i
1
k
i
k
i
1
k
i
u
1
u
u
u
Jeżeli =0 schemat jest
nazywany jawnym lub
explicite
k+1
k
i-1 i i+1
Jeżeli 0 schemat jest
nazywany niejawnym lub
implicite
k+1
k
i-1 i i+1
k+1
k
i-1 i i+1
Wartości w warstwie k+1
otrzymujemy rozwiązując
układ równań:
k
i
k
i
k
i
1
k
i
1
k
i
u
1
u
1
u
u
1
Czysto niejawny schemat:
k+1
k
i
1
k+1
k
i-1 i i+1
Schemat Cranka - Nicholsona:
5
.
0
Oszacowanie dokładności aproksymacji.
Rozwiązanie dokładne zagadnienia brzegowego
T
0,
t
dla
t
f
t
1,
u
t
f
t
,
0
u
0,1
x
dla
x
0
,
x
u
T
0,
t
0,1
x
dla
t
,
x
f
,
u
,
u
2
1
0
xx
t
jest u(x,t) i jego wartość w węzłach (x
i
,t
k
) siatki będzie oznaczana
u(i,k).
Rozwiązanie zagadnienia brzegowego sformułowanego dyskretnie
k
2
k
N
k
1
k
0
i
0
i
k
i
k
i
1
k
i
k
i
1
k
i
f
u
f
u
u
u
1
u
u
u
jest
k
i
u
Błąd aproksymacji
k
i
z
jest
)
k
,
i
(
u
u
z
k
i
k
i
Dla oceny błędu
k
i
z
w kroku k-tym wprowadza się normę, np.:
k
i
N
,
0
i
z
max
z
lub
1
N
1
i
2
k
i
z
h
z
Z
)
k
,
i
(
u
u
z
k
i
k
i
wynika, że
)
k
,
i
(
u
z
u
k
i
k
i
i podstawiając do
w miejsce
otrzymujemy
równoważne
zadanie różnicowe
dla
k
2
k
N
k
1
k
0
i
0
i
k
i
k
i
1
k
i
k
i
1
k
i
f
u
f
u
u
u
1
u
u
u
k
i
u
k
i
z
0
z
0
z
0
z
z
1
z
z
z
k
N
k
0
0
i
k
i
k
i
1
k
i
k
i
1
k
i
0
z
0
z
0
z
z
1
z
z
z
k
N
k
0
0
i
k
i
k
i
1
k
i
k
i
1
k
i
gdzie
k
i
k
i
k
,
i
u
1
k
,
i
u
k
,
i
u
1
1
k
,
i
u
błąd schematu różnicowego
w stosunku do rozwiązania
dokładnego u(x,t).
k
2
k
N
k
1
k
0
i
0
i
k
i
k
i
1
k
i
k
i
1
k
i
f
u
f
u
u
u
1
u
u
u
Mówimy, że
k
2
k
N
k
1
k
0
i
0
i
k
i
k
i
1
k
i
k
i
1
k
i
f
u
f
u
u
u
1
u
u
u
przybliża rozwiązanie problemu brzegowego
T
0,
t
dla
t
f
t
1,
u
t
f
t
,
0
u
0,1
x
dla
x
0
,
x
u
T
0,
t
0,1
x
dla
t
,
x
f
,
u
,
u
2
1
0
xx
t
z dokładnością rzędu (m,n) lub O(h
m
+
n
), jeżeli
k
i
k
i
k
,
i
u
1
k
,
i
u
k
,
i
u
1
1
k
,
i
u
spełnia nierówność:
n
m
h
M
gdzie M - stała.
Dla uproszczenia zapisu wprowadzamy:
u
k
,
i
u
1
k
,
i
u
i mamy:
u
5
.
0
k
,
i
u
1
k
,
i
u
5
.
0
k
,
i
u
1
k
,
i
u
5
.
0
k
,
i
u
1
k
,
i
u
5
.
0
1
k
,
i
u
u
5
.
0
k
,
i
u
1
k
,
i
u
5
.
0
k
,
i
u
u
5
.
0
k
,
i
u
1
k
,
i
u
5
.
0
k
,
i
u
u
k
,
i
u
k
,
i
u
1
k
,
i
u
k
,
i
u
1
1
k
,
i
u
Uwzględniając powyższe równości i podstawiając do
k
i
k
i
k
,
i
u
1
k
,
i
u
k
,
i
u
1
1
k
,
i
u
mamy
k
i
k
i
u
5
.
0
k
,
i
u
1
k
,
i
u
5
.
0
Rozwijając funkcje u(n,p) w szereg Taylora w otoczeniu
punktu: x
i
,t
k
+0.5 oraz wprowadzając oznaczenie:
u
5
.
0
t
,
x
u
k
i
Będziemy mieli:
2
t
3
tt
2
3
tt
2
t
3
tt
2
t
4
xxxx
2
xx
O
,
u
u
O
,
u
8
u
)
k
,
i
(
u
1
k
,
i
u
5
.
0
O
,
u
8
,
u
5
.
0
u
k
,
i
u
O
,
u
8
,
u
5
.
0
u
1
k
,
i
u
h
O
,
u
12
h
,
u
k
,
i
u
Uwzględniając powyższe zależności można
k
i
k
i
u
5
.
0
k
,
i
u
1
k
,
i
u
5
.
0
zapisać w postaci:
4
2
xxxx
2
txx
k
i
t
xx
k
i
h
O
,
u
12
h
,
u
5
.
0
,
u
,
u
Ale
0
f
,
u
,
u
t
xx
gdyż u jest rozwiązaniem dokładnym
T
0,
t
dla
t
f
t
1,
u
t
f
t
,
0
u
0,1
x
dla
x
0
,
x
u
T
0,
t
0,1
x
dla
t
,
x
f
,
u
,
u
2
1
0
xx
t
i w każdym punkcie obszaru spełnia równanie paraboliczne.
Uwzględniając, że
xx
txx
xxxx
,
f
,
u
,
u
mamy
2
4
xx
2
txx
2
k
i
k
i
h
O
,
f
12
h
,
u
12
h
5
.
0
f
2
4
xx
2
txx
2
k
i
k
i
h
O
,
f
12
h
,
u
12
h
5
.
0
f
Jeżeli wyrażenie w nawiasie kwadratowym jest równe zeru, tzn.:
12
h
5
.
0
0
12
h
5
.
0
2
2
oraz
xx
2
k
i
,
f
12
h
f
W obliczeniach numerycznych wygodniej przyjąć:
5
.
0
k
1
i
5
.
0
k
1
i
5
.
0
k
i
k
i
f
f
12
1
f
6
5
to schemat ma
dokładność O(h
4
+
2
)
k
2
k
N
k
1
k
0
i
0
i
k
i
k
i
1
k
i
k
i
1
k
i
f
u
f
u
u
u
1
u
u
u
Schemat:
z
5
.
0
k
1
i
5
.
0
k
1
i
5
.
0
k
i
k
i
f
f
12
1
f
6
5
jest schematem o podwyższonej dokładności wynoszącej:
2
4
h
O
Jeżeli =0.5 jak w schemacie Cranka-Nicholsona, to
2
4
xx
2
txx
2
k
i
k
i
h
O
,
f
12
h
,
u
12
h
5
.
0
f
2
4
xx
2
txx
2
k
i
k
i
h
O
,
f
12
h
,
u
12
h
5
.
0
f
i dla =0.5 mamy:
2
2
k
i
k
i
h
O
f
Dla zachowania oceny zbieżności O(h
2
+
2
) należy przyjąć:
5
.
0
k
i
k
i
f
lub
k
i
1
k
i
k
i
f
f
5
.
0
Jeżeli 0.5 i
*
, to dokładność obliczeń jest rzędu O(h
2
+).
Stabilność
T
0,
t
dla
0
t
1,
u
0
t
,
0
u
0,1
x
dla
0
0
,
x
u
T
0,
t
0,1
x
dla
,
u
,
u
xx
t
0
u
0
u
0
u
u
1
u
u
u
k
N
k
0
0
i
k
i
1
k
i
k
i
1
k
i
Zbadamy zachowanie się schematu:
jawnego tj. =0
0
u
0
u
0
u
u
u
2
u
h
u
u
k
N
k
0
0
i
k
1
i
k
i
k
1
i
2
k
i
1
k
i
Schemat jest
Przyjmijmy:
1
h
2
i załóżmy, że warunek początkowy w
punkcie i-tym jest dany z błędem . Badamy jak przenosi się błąd
na siatce.
x
t
0
N
0
0
u
0
u
i
n
i
n
0
u
u
u
u
u
k
N
k
0
0
n
k
1
i
k
i
k
1
i
1
k
i
2
3
2
3
6
7
6
3
4
10
16
19
16
10
4
5
15
30
45
51
45
30
15
5
0
0
6
21
50
90
126
141
126
90
50
20 0
Schemat jawny
z =0 i
1
h
2
x
t
0
N
0
0
u
0
u
i
n
i
n
0
u
u
1
.
0
u
8
.
0
u
1
.
0
u
k
N
k
0
0
n
k
1
i
k
i
k
1
i
1
k
i
1
.
0
8
.
0
1
.
0
01
.
0
16
.
0
66
.
0
16
.
0
01
.
0
01
.
0
02
.
0
2
.
0
56
.
0
2
.
0
02
.
0
01
.
0
0
01
.
0
04
.
0
22
.
0
49
.
0
22
.
0
04
.
0
01
.
0
0
0
0
01
.
0
06
.
0
23
.
0
44
.
0
23
.
0
06
.
0
01
.
0
0
0
0
0
0
01
.
0
07
.
0
23
.
0
4
.
0
07
.
0
01
.
0
0
0
Schemat jawny
z =0 i
1
.
0
h
2
Obliczenia z dokładnością do 2 miejsc
23
.
0
Analiza stabilności metodą spektralną
Niech rozwiązanie jednorodnego zagadnienia różnicowego:
0
u
0
u
u
u
1
u
u
u
k
N
k
0
n
0
n
k
n
1
k
n
k
n
1
k
n
ma postać
in
exp
u
k
k
n
gdzie
1
i
2
1
n
i
in
1
n
i
k
k
n
h
e
e
2
e
u
ale
2
sin
e
4
i
2
e
e
e
4
e
2
e
e
e
e
2
e
2
in
2
2
i
2
i
in
i
i
in
1
n
i
in
1
n
i
Po podstawieniu do
k
n
1
k
n
k
n
1
k
n
u
1
u
u
u
mamy
2
sin
e
1
h
4
2
sin
e
h
4
e
e
2
in
k
2
2
in
1
k
2
in
k
in
1
k
2
sin
h
4
1
2
sin
h
1
4
1
2
2
2
2
Po podzieleniu przez
in
k
e
otrzymujemy równanie:
Warunek konieczny stabilności Neumanna stwierdza, że schemat
różnicowy jest stabilny, jeżeli
1
Dla =0 mamy
2
sin
h
4
1
2
2
i na podstawie kryterium Neumanna otrzymujemy:
1
2
sin
h
4
1
2
2
czyli
0
2
sin
h
4
2
2
2
Dla = znajdujemy warunek na stosunek:
2
1
h
2
2
1
h
2
Warunkiem zbieżności schematu jawnego jest spełnienie
powyższego warunku. Warunek jest również prawdziwy
dla schematu jawnego w przypadku wielowymiarowego
równania parabolicznego.
Dowolne >0.
Ocenę prowadzimy przy =.
1
h
4
1
h
1
4
1
2
2
czyli
2
2
2
h
4
1
h
1
4
1
h
4
1
Prawa nierówność jest spełniona dla dowolnych , a z lewej
mamy:
4
h
2
1
2
4
h
2
1
2
Dla spełniających nierówność:
1
h
4
1
h
1
4
1
2
2
warunek:
jest spełniony dla dowolnego stosunku
2
h
W szczególności schemat Cranka-Nicholsona =0.5 jest stabilny
dla dowolnego stosunku kroków
2
h
Dla schematu o podwyższonej dokładności
12
h
5
.
0
2
mamy
4
h
2
1
2
bo
4
h
12
h
2
2
Przedstawione rozważania można rozszerzyć na przypadki
wielowymiarowe jak również na równania o zmiennych
współczynnikach.
W przypadku równań wielowymiarowych ocena zbieżności
zależy również od sposobu aproksymacji warunków
brzegowych podobnie jak w przypadku równań eliptycznych.
t
e
L
,
t
u
t
e
0
,
t
u
x
w
x
,
0
,
u
x
w
x
,
0
u
,
u
v
,
u
L
P
1
t
0
xx
2
tt
Równania hiperboliczne
Jako przykład zostanie rozpatrzone równanie linii długiej
bez strat o długości L. Dla napięcia u mamy:
Wprowadzamy siatkę prostokątną:
x
t
0
N
h
i-1 i i+1
k+1
k
k-1
i funkcję węzłową oznaczamy:
ih
x
,
k
t
u
u
i
k
k
i
Przyjmujemy aproksymację pochodnych:
2
1
k
i
k
i
1
k
i
k
i
tt
v
u
u
2
u
Du
,
u
i
2
k
1
i
k
i
k
1
i
k
i
xx
h
u
u
2
u
u
,
u
i rozpatrujemy następujący schemat trójwarstwowy
z parametrem >0:
i
1
0
i
1
i
i
0
0
i
k
L
k
N
k
P
k
0
1
k
i
k
i
1
k
i
k
i
w
~
u
u
w
u
e
u
e
u
u
u
2
1
u
Du
gdzie warunek początkowy
i
1
w
~ konstruujemy tak, aby
zachować rząd aproksymacji O(
2
).
Mamy
3
0
t
tt
2
0
t
t
O
t
,
x
,
u
5
.
0
t
,
x
,
u
0
,
x
u
,
x
u
czyli
2
0
t
tt
i
0
t
t
i
0
i
1
i
O
,
u
5
.
0
,
u
u
u
Z równania falowego mamy:
2
i
2
tt
i
h
O
u
v
,
u
Warunek początkowy dla pierwszej pochodnej będzie określony
z dokładnością O(h
2
+
2
), jeżeli przyjąć, że
2
2
i
0
2
i
1
0
i
1
i
h
O
w
v
5
.
0
w
u
u
czyli
i
0
2
i
1
i
1
w
v
5
.
0
w
w
~
Ostatecznie schemat różnicowy dla rozwiązania równania
falowego jest
i
0
2
i
1
i
0
1
i
i
0
0
i
k
L
k
N
k
P
k
0
1
k
i
k
i
1
k
i
k
i
w
v
5
.
0
w
w
u
w
u
e
u
e
u
u
u
2
1
u
Du
Ocena dokładności aproksymacji
Postępujemy podobnie jak poprzednio, a więc niech
k
i
k
i
k
i
t
,
x
u
u
z
gdzie
k
i
u
jest rozwiązaniem różnicowego zagadnienia
i
0
2
i
1
i
0
1
i
i
0
0
i
k
L
k
N
k
P
k
0
1
k
i
k
i
1
k
i
k
i
w
v
5
.
0
w
w
u
w
u
e
u
e
u
u
u
2
1
u
Du
a u(x
i
,t
k
) jest rozwiązaniem problemu brzegowego:
t
e
L
,
t
u
t
e
0
,
t
u
x
w
x
,
0
,
u
x
w
x
,
0
u
,
u
v
,
u
L
P
1
t
0
xx
2
tt
w punkcie x
i
,t
k
.
Pisząc
k
i
k
i
k
i
t
,
x
u
z
u
otrzymujemy:
i
1
i
0
i
k
N
k
0
k
i
1
k
i
k
i
1
k
i
k
i
z
0
z
0
z
0
z
z
z
2
1
z
Dz
gdzie
i
0
2
i
1
i
i
0
i
k
i
1
k
i
k
i
1
k
i
k
i
w
v
5
.
0
w
1
,
x
u
w
t
,
x
Du
t
,
x
u
t
,
x
u
2
1
t
,
x
u
Z konstrukcji warunku początkowego dla pochodnej wynika, że
2
2
i
h
O
Rozwijając w szereg Taylora mamy:
2
t
,
x
tt
k
i
1
k
i
1
k
i
k
i
,
u
t
,
x
u
2
t
,
x
u
t
,
x
u
Korzystając z otrzymanego wyniku mamy:
k
i
t
,
x
tt
2
k
i
k
i
1
k
i
1
k
i
,
u
t
,
x
u
t
,
x
u
2
1
t
,
x
u
t
,
x
u
Stąd otrzymujemy z dokładnością do małych 4-go rzędu:
k
i
2
2
4
t
,
x
t
x
2
x
2
xx
1
k
i
k
i
1
k
i
,
u
,
u
12
h
,
u
t
,
x
u
t
,
x
u
2
1
t
,
x
u
Podobnie
k
i
4
t
,
x
t
2
2
tt
2
k
i
,
u
v
12
,
u
v
1
t
,
x
Du
z dokładnością do małych 4-go rzędu. Ostatecznie otrzymujemy
ocenę błędu:
k
i
4
2
2
4
t
,
x
t
2
2
t
x
2
x
2
tt
2
xx
k
i
,
u
v
12
,
u
,
u
12
h
,
u
v
1
,
u
Funkcja u spełnia równanie falowe, a więc
tt
2
xx
,
u
v
1
,
u
stąd
2
2
k
i
h
O
niezależnie od
!!!
Oznacza to również zależność odwrotną, a mianowicie
dobór zależy tylko i wyłącznie od stabilności, a nie ma
wpływu na dokładność obliczeń.
Stabilność
n
0
2
n
1
n
0
1
n
n
0
0
n
k
N
k
0
1
k
n
k
n
1
k
n
k
n
w
v
5
.
0
w
w
u
w
u
0
u
0
u
u
u
2
1
u
Du
Analizujemy stabilność schematu różnicowego przy jednorodnych
warunkach brzegowych:
Przyjmujemy rozwiązanie w postaci:
in
exp
u
k
k
n
1
k
n
k
n
1
k
n
k
n
u
u
2
1
u
Du
i podstawiając do równania różnicowego:
po wykonaniu kilku przekształceń otrzymujemy:
2
sin
2
1
h
v
2
2
2
1
k
k
1
k
2
1
k
k
1
k
Dzieląc równanie przez
k-1
i grupując wyrazy otrzymujemy
równanie określające :
0
1
2
sin
h
v
2
1
2
sin
h
v
2
2
1
2
1
1
2
2
2
2
0
1
r
4
1
r
2
1
2
2
2
2
Badamy pierwiastki równania:
gdzie
h
v
r
przy =.
Wyróżnik:
2
2
2
2
r
1
4
1
r
4
1
r
4
Jeżeli <0, to równanie ma dwa sprzężone pierwiastki
zespolone
1
,
2
o module równym 1 co wynika ze wzoru
Viety:
1
2
1
2
1
Dla =0 otrzymujemy warunek Couranta:
1
h
v
r
który oznacza, że prędkość wędrówki fali na siatce h/
jest większa od prędkości fazowej.
1
h
v
r
Przypadek 0 prowadzi do pierwiastków większych co
do modułu od jedności i dlatego należy te przypadki
odrzucić.
Analizę można rozszerzyć na przypadki wielowymiarowe.
Wady i zalety metody różnicowej
Zalety:
1. Proste konstruowanie siatki podziałowej.
2. Prosta konstrukcja układu równań różnicowych szczególnie
w środowiskach izotropowych.
3. Opracowane oceny błędów metody i warunki stabilności.
Wady:
1. Duże trudności z dobrą aproksymacją brzegu lub wymuszony
mały krok siatki.
2. Trudności z utrzymaniem rzędu aproksymacji przy interpolacji
warunków brzegowych.
3. Praktycznie konieczność obliczania całego obszaru z tym
samym krokiem podziałowym.