METODY NUMERYCZNE
W ELEKTROTECHNICE
METODY NUMERYCZNE
W ELEKTROTECHNICE
2
Metody numeryczne – dział matematyki stosowanej, zajmujący się opracowywaniem i
badaniem metod przybliżonego rozwiązywania problemów obliczeniowych w modelach
matematycznych innych dziedzin nauki
Przykładowe zastosowania:
elektrotechnika – obliczanie parametrów obwodów elektrycznych
medycyna – tomografia komputerowa, opracowywanie nowych leków
chemia – konstruowanie nowych cząsteczek
inżynieria – przemysł samochodowy i lotniczy
informatyka – konstruowanie nowych procesorów
W obrębie klasycznych metod numerycznych możemy wyróżnić m.in. takie zagadnienia jak:
analiza błędów zaokrągleń
interpolacja
aproksymacja
rozwiązywanie równań i układów równań nieliniowych
całkowanie i różniczkowanie numeryczne
rozwiązywanie układów równań liniowych
obliczanie wartości własnych i wektorów własnych macierzy
rozwiązywanie zagadnień dla równań różniczkowych zwyczajnych i cząstkowych
rozwiązywanie równań całkowych i układów równań całkowych
(każdy wynik musi podlegać weryfikacji!)
Reprezentacja liczb w maszynie cyfrowej
Liczby całkowite (stałopozycyjne = stałoprzecinkowe):
i
n
0
i
i
2
e
znak
l
, gdzie e
i
= 0 lub 1
Jeżeli rejestr ma d-bitów, wówczas liczba całkowita n może zawierać się w przedziale
-2
d-1
< n < 2
d-1
-1
przykład: Zapisać liczbę 18 w systemie dwójkowym:
0
1
2
3
4
2
0
2
1
2
0
2
0
2
1
l
, może być zatem przechowywana w słowie o długości D
+ 1 > 5 bitów jako
a
rozwinieci
cyfr
d
znak
0010010
...
00
0
Liczby rzeczywiste (zmiennopozycyjne):
c
2
m
zn
x
,
gdzie ½ < m < 1, m-mantysa, i-cecha
i
1
i
i
2
e
m
3
W maszynie cyfrowej mantysa zapisywana jest w t-bitach m
t
i
t
1
i
i
t
)
2
t
(
t
2
e
2
e
m
w wyniku zaokrąglenia do t-bitów mantysy
t
t
2
2
1
m
m
Zmiennopozycyjna reprezentacja liczby rzeczywistej x oznaczana jest symbolem rd(x) i jest
równa
t
c
m
2
zn
)
x
(
rd
t
2
x
x
)
x
(
rd
!
Gdy elementy macierzy są rzędu
lub n musimy dokonać przekształcenia całego
układu. Uciekamy się od operacji na małych liczbach.
Nie prowadzi się operacji na małych liczbach, co oznacza, że:
1
x
)
x
(
rd
,
gdzie
t
2
liczbę 2
-t
nazywamy dokładnością maszynową
przykład: Liczba 18,5 daje się przedstawić w postaci:
6
5
4
3
2
1
)
2
1
2
0
2
1
(
2
1
2
0
2
1
2
0
2
0
2
1
2
x
0
1
2
może być więc przechowywana w słowie o długości d + 1 > 11 bitów o t > 6 bitach mantysy:
mantysy
bitów
t
cechy
bitów
t
d
c
znak
znak
0
...
1001010
00101
...
000
0
0
0,5781525
W sytuacji, gdy elementy macierzy są < 1, np. mV, musimy wtedy dokonać przeskalowania
całego układu (np. poprzez pomnożenie macierzy przez jakąś liczbę).
Błędy obliczeń
Przy obliczeniach wykonywanych na maszynach cyfrowych spotyka się trzy podstawowe
rodzaje błędów:
błędy wejściowe (danych wejściowych) – występujące, gdy dane wprowadzane do
pamięci maszyny cyfrowej odbiegają od dokładnych wartości tych danych. Źródłem
4
tych błędów jest najczęściej skończona długość słów binarnych reprezentujących
liczby i w związku z tym nieuniknione jest wstępne ich zaokrąglenie
błędy odcięcia – powstają podczas obliczeń na skutek zmniejszenia liczby działań, na
przykład przy obliczaniu sum nieskończonych. Chcąc obliczyć wartość wyrażenia e
x
równego szeregowi:
...
x
!
N
1
...
x
2
1
x
1
N
2
zastępuje je sumą częściową o
odpowiednio dobranej wartości N:
N
x
N
x
x
!
1
...
2
1
1
2
Jeżeli liczba N będzie niedostatecznie duża, to uzyskana w ten sposób wartość liczby e
x
będzie obliczona niedokładnie, a popełniony w ten sposób błąd jest właśnie błędem odcięcia.
Błędy tego typu powstają też często przy obliczaniu wielkości będących granicami. Podobnie
zastąpienie wartości pochodnej funkcji jej ilorazem różnicowym powoduje powstanie błędu
odcięcia. W wielu przypadkach daje się uniknąć błędu wejściowego i odcięcia przez
ograniczenie danych.
Lemat Wilkinsona:
Błędy zaokrągleń powstające podczas wykonywania działań zmiennopozycyjnych są
równoważne zastępczemu zaburzeniu liczb, na których wykonujemy działania. W przypadku
działań arytmetycznych otrzymujemy:
)
1
(
x
)
1
(
x
)
x
x
(
fl
2
2
1
1
2
1
2
1
3
3
2
1
2
1
)
1
(
x
x
)
x
x
(
fl
1
2
5
2
1
4
5
2
1
2
4
1
2
1
)
1
(
x
x
x
)
1
(
x
x
x
fl
gdzie
ε
i
są co do modułu niewiększe niż dokładność maszyny
ε
. Przedrostek fl oznacza, że są
to wyniki działań wyznaczone przy użyciu maszyny cyfrowej.
Oszacowanie błędów w obliczeniach iteracyjnych:
Wiele algorytmów obliczeniowych polega na wyznaczeniu ciągu liczb zbieżnego do
poszukiwanej wartości.
Przykład:
wyznaczyć wartość funkcji
x
y
, przekształcając do równania iteracyjnego. Można to
uczynić obliczając ciąg y
0
, y
1
, y
2
..., gdzie y
0
jest dowolnie wybraną liczbą, natomiast każdy
element ciągu y
i
jest dany wzorem:
1
i
1
i
i
y
x
y
2
1
y
, dla i = 1, 2, 3, ...
(niewiadoma musi być po obu stronach równania)
2
x
y
x
y
2
x
y
y
y
x
1
y
, współczynnik przed funkcją (tutaj 1) musi być < 1 dlatego:
y
x
y
y
2
5
2
:
y
y
x
y
2
1
i
1
i
i
y
y
x
2
1
y
- iteracja i-ta
np.
1
i
1
i
i
2
y
y
x
7
2
1
y
y
x
7
)
y
y
2
(
x
7
y
takiego równania nie da się rozwiązać
bo 7/2
ale można zrobić tak:
1
i
1
i
i
y
7
y
x
7
8
1
y
y
x
7
)
y
7
y
8
(
jeśli założyć, że działania wykonywane są dokładnie, to ciąg y
0
, y
1
, y
2
, ... będzie zbieżny do
liczby x, ponieważ
x
p
y
i
i
, to
x
p
1
p
2
1
y
i
i
1
i
, a ciąg p
o
> 0
i
i
1
i
p
1
p
2
1
p
, przy i = 1, 2, 3, ... jest zbieżny do 1
Przykład:
rozwiązanie iteracyjne
5
4
)
4
5
(
2
3
7
7
8
5
2
4
3
7
x
y
y
y
x
x
y
x
y
x
5
1
2
5
4
2
8
1
3
7
7
1
1
i
i
i
i
i
i
y
x
y
y
x
x
wszystkie składniki przy niewiadomych są < 1 (7/8, -7/8, 3/5)
przyjmujemy dowolnie x
i
= 1, y
i
= 2
Działania potrzebne do wyznaczania kolejnych wartości ciągu y
0
, y
1
, y
2
, ... nazywamy
iteracją. Informację o tym czy w kolejnej iteracji nastąpiło przybliżenie do rozwiązania
uzyskujemy obliczając wyrażenie:
x
y
x
y
q
i
i
i
1
W przypadku dokładnych działań, gdy
x
p
y
i
i
i
i
i
i
p
2
1
p
q
, więc q
i
< 0 dla i = 1, 2, 3, ...
jeśli p
i
>1/3. Oznacza to, że maleje wartość bezwzględna różnicy
x
y
1
odpowiednio
dużego numeru iteracji i.
6
Dla obliczeń wykonywanych na maszynie cyfrowej, mamy:
yi
1
1
x
1
y
2
1
)
y
(
fl
i
3
i
2
i
1
i
1
i
, gdzie
3
i
k
, k = 1, 2, 3, ...
Dla liczby
)
y
(
fl
y
1
i
'
1
i
uzyskamy
x
y
x
y
q
i
'
'
i
1
i
Jeżeli y
i
jest już dostatecznie dobrym przybliżeniem liczby
x , a więc p
i
jest bliskie
jedności, to pierwszy składnik powyższego wyrażenia ma pomijalnie małą wartość, natomiast
drugi ze składników wynika z błędów zaokrągleń i może mieć tym większą wartość im
bliższe jedności jest p
i
.
Przykład: dla liczby 4 x=4
000609
,
2
05
,
2
4
05
,
2
2
1
y
05
,
2
5
,
2
4
5
,
2
2
1
y
5
,
2
4
4
4
2
1
y
1
y
3
2
1
0
Algorytm:
zmienne rzeczywiste x, y, eps
y:= x/4;
{pierwsze przybliżenie}
repeat
y:= 0,5*(y+(x/y));
until abs (x-sqr(y))<=eps
{eps jest zakładaną dokładnością obliczeń)
Uwarunkowanie zadania i stabilność algorytmów.
Załóżmy, że mamy skończoną liczbę danych rzeczywistych, x = ( x
1
, x
2
, ... x
n
), na ich
podstawie chcemy obliczyć skończenie wiele wyników rzeczywistych y = ( y
1
, y
2
, ... y
m
).
Będziemy więc chcieli określić wartość y według odwzorowania y = φ (x), gdzie
: D → R
m
jest ciągłym odwzorowaniem i D
R
m
.
Algorytm to jednoznaczny przepis obliczania wartości odwzorowania φ (x) składający się ze
skończonej liczby kroków.
Jeśli zadanie rozwiązujemy numerycznie to zamiast dokładnymi wartościami x
1
, x
2
, ... x
n
posługujemy się reprezentacjami rd (x
1
), rd (x
2
), ..., rd (x
n
). Operacje elementarne nie są więc
realizowane dokładnie tylko przez odwzorowanie zastępcze fl.
Ta niewielka zmiana danych
powoduje, że różne algorytmy rozwiązania tego samego zadania dają na ogół niejednakowe
wyniki. Dużą rolę odgrywa tu przenoszenie błędów zaokrągleń. Często niewielkie zmiany
danych powodują duże względne zmiany rozwiązania zadania. Zadanie takie nazywamy źle
uwarunkowanym. Wielkości charakteryzujące wpływ zaburzeń danych na zaburzenia
rozwiązania nazywamy wskaźnikami uwarunkowania zadania.
Przykład:
20
1
k
0
1
19
19
20
20
k
x
a
x
a
...
x
a
x
a
)
x
(
w
7
Zerami tego wielomianu będą liczby naturalne 1, 2, ..., 20. Gdy zakłócimy np. a
19
= - 210 i
jego wartość wynosi a
19
(ε) = - (210 + 2
-23
), czyli a
19
(ε) = a
19
(1+ε). Otrzymujemy wtedy
wielomian w
ε
(x) = w (x) – 2
-23
, który posiada już pierwiastki zespolone i np. najbliższe
pierwiastkowi 15 wielomianu w(x) są pierwiastki 13,992358137
2,518830070j.
Metodę badania przenoszenia błędów można rozbudować do analizy różniczkowej błędów
algorytmu.
Niech δ
x1
= rd (x
i
) – x
i
dla i = 1, 2, ..., n
i
n
1
i
i
j
yi
x
x
)
x
(
x
x
rd
dla j = 1, 2, ..., m
We wzorze tym, czynnikiem określającym wrażliwość y
j
na bezwzględną zmianę δx
i
jest
i
j
x
)
x
(
Analogiczny wzór można wyprowadzić dla przenoszenia się błędów względnych. Jeśli
0
y
j
dla j = 0,1, 2, ..., m i
0
x
i
dla i = 1, 2, ..., n to
i
j
x
i
j
n
1
i
j
j
y
x
)
x
(
)
x
(
x
, gdzie
i
j
j
j
x
)
x
(
)
x
(
x
nazywamy wskaźnikiem uwarunkowania.
Jeśli wskaźniki uwarunkowania co do wartości bezwzględnej są duże to zadanie jest źle
uwarunkowane.
Przykład:
Badamy uwarunkowanie sumy y = a + b + c
R
c
b
a
,
,
c
b
a
c
b
a
c
b
a
3
2
1
1
1
1
=
c
b
a
c
c
b
a
b
c
b
a
a
c
b
a
c
b
a
3
2
1
, gdzie ε jest
największą z wartości ε
1
, ε
2
i ε
3
, gdzie
c
b
a
c
c
b
a
b
c
b
a
a
jest wskaźnikiem
uwarunkowania zadania.
Stabilność numeryczna algorytmów
Algorytm jest stabilny, jeżeli posiada tę własność, że małe błędy powstałe w pewnym
kroku algorytmu nie są powiększane w innych krokach oraz nie powodują poważnego
ograniczenia dokładności wszystkich obliczeń. Oznacza to, że algorytm jest numerycznie
stabilny, gdy zwiększając dokładność obliczeń można wyznaczyć (z dowolną dokładnością)
istniejące rozwiązanie zadania. Numeryczną stabilność zadania łatwo sprawdzić rozwiązując
zadanie raz z dokładnością np. 10
-6
, a potem z dokładnością 10
-12
.
8
INTERPOLACJA:
Dana jest funkcja y = f (x),
n
x
x
x
,
0
, dla której znana jest tablica wartości w punktach
zwanych węzłami interpolacji. Należy wyznaczyć taką funkcję W(x), aby:
W(x
0
) = Y
0
, W(x
1
) = Y
1
, ..., W(x
n
) = Y
n
interpolacja funkcji f(x)
Zadaniem interpolacji jest wyznaczenie przybliżonych wartości funkcji zwanej funkcją
interpolową w punktach nie będących węzłami interpolacji. Przybliżoną wartość funkcji
obliczamy za pomocą funkcji zwanej funkcją interpolującą, która w węzłach ma te same
wartości co funkcja interpolowana.
Funkcja interpolująca jest funkcją pewnej klasy. Najczęściej będzie to wielomian
algebraiczny, wielomian trygonometryczny, funkcja wymierna i funkcja sklejana.
Interpolację stosuje się najczęściej, gdy nie znamy analitycznej postaci funkcji (jest ona tylko
stablicowana) lub gdy jej postać analityczna jest zbyt skomplikowana.
Najczęściej stosowaną metodą wyznaczania funkcji W(x) jest jej dobór w postaci kombinacji
liniowej n+1 funkcji bazowych
)
x
(
),...,
x
(
),
x
(
),
x
(
n
2
1
0
n
i
i
i
x
a
x
W
0
)
(
)
(
wyrażenie to nazywamy wielomianem uogólnionym.
Wprowadzając macierz bazową
)
x
(
),...,
x
(
),
x
(
n
1
0
i macierz współczynników
n
1
0
a
a
a
A
mamy
A
x
x
W
)
(
)
(
warunek W(x
0
) = Y
0
, W(x
1
) = Y
1
, ..., W(x
n
) = Y
n
można zapisać w postaci układu równań liniowych X · A = Y, gdzie
9
)
x
(
)
x
(
)
x
(
)
x
(
)
x
(
)
x
(
)
x
(
)
x
(
)
x
(
X
n
n
n
1
n
0
1
n
1
1
1
0
0
n
0
1
0
0
n
1
0
Y
Y
Y
Y
Jeśli macierz X nie jest osobliwa (da się odwrócić), to:
A = X
-1
·Y
co ostatecznie daje
W(x) = Φ(x) · X
-1
· Y
Interpolacja wielomianowa
W praktyce często używa się bazy złożonej z jednomianów
,
)
(
,
,
)
(
,
1
)
(
1
0
n
n
x
x
x
x
x
Baza dla funkcji ciągłych na odcinku skończonym [x
o
, x
n
] jest bazą zamkniętą, tzn., że każda
funkcja tej klasy może być przedstawiona w postaci szeregu złożonego z funkcji bazowych.
Wielomian interpolacyjny ma w tym przypadku postać:
n
n
x
a
x
a
x
a
a
x
W
2
2
1
0
)
(
dodatkowo musi być spełniony warunek:
n
n
n
n
n
1
0
1
n
1
n
1
1
0
0
n
0
n
0
1
0
y
x
a
x
a
a
y
x
a
x
a
a
y
x
a
x
a
a
Powyższy układ równań ma jedyne rozwiązanie względem a
1
, jeżeli wartości x
0
, x
1
, ..., x
n
są
od siebie różne
n
n
n
n
n
x
x
x
x
x
x
X
1
1
1
det
1
1
0
0
Macierz X
-1
dla bazy wielomianowej
,
x
)
x
(
,
,
x
)
x
(
,
1
)
x
(
n
n
1
0
nazywana jest
macierzą Lagrange’a
Zauważyć należy, że każdy zbiór węzłów równoodległych x
i+1
– x
i
= h = const. można
sprowadzić do zbioru podstawowego podstawiając
h
x
x
q
0
, wówczas
n
q
,
0
, a macierz
Φ przyjmuje postać
n
2
q
,
,
q
,
q
,
1
10
Interpolacja Lagrange’a
Przedstawiony powyżej sposób podejścia do interpolacji nie jest zbyt efektywny, ponieważ
macierz X jest macierzą pełną i nie zawsze dobrze uwarunkowaną, co oznacza, że
numeryczna procedura jej odwracania może być obarczona dużym błędem.
W interpolacji wielomianowej Lagrange’a dla n+1 węzłów interpolacji
)
y
,
x
(
,
),
y
,
x
(
,
),
y
,
x
(
),
y
,
x
(
n
n
i
i
1
1
0
0
przyjmuje się funkcje bazowe w postaci
)
x
x
(
)
x
x
)(
x
x
)(
x
x
)(
x
x
(
)
x
(
)
x
x
(
)
x
x
)(
x
x
(
)
x
x
)(
x
x
)(
x
x
)(
x
x
(
)
x
(
)
x
x
(
)
x
x
)(
x
x
)(
x
x
(
)
x
(
)
x
x
(
)
x
x
)(
x
x
)(
x
x
(
)
x
(
1
n
3
2
1
0
n
n
1
i
1
i
3
2
1
0
i
n
3
2
0
1
n
3
2
1
0
Funkcje te są wielomianami stopnia n zbudowanymi w ten sposób, że w funkcji bazowej φ
1
brakuje czynnika (x-x
i
). Zatem wielomian interpolacji wyraża się następującym wzorem:
)
x
x
(
)
x
x
)(
x
x
(
a
...
)
x
x
(
)
x
x
)(
x
x
(
a
)
x
x
(
)
x
x
)(
x
x
(
a
)
x
(
W
1
n
1
0
n
n
2
0
1
n
2
1
0
współczynniki a
0
... a
n
tego wielomianu wyznaczamy z równania:
X · A = Y, przy czym macierz X ma postać:
1
n
n
0
n
n
1
0
1
n
0
1
0
x
x
...
x
x
0
0
0
x
x
...
x
x
0
0
0
x
x
...
x
x
X
Macierz posiada tylko główną przekątną niezerową w związku z tym układ równań X · A = Y
rozwiązuje się natychmiastowo
)
x
x
(
)
x
x
)(
x
x
(
y
a
)
x
x
(
)
x
x
)(
x
x
(
y
a
)
x
x
(
)
x
x
)(
x
x
(
y
a
1
n
n
1
n
0
n
n
n
n
1
2
1
0
1
1
1
n
0
2
0
1
0
0
0
11
Można więc wielomian interpolacyjny Lagrange’a zapisać w postaci ułamka:
n
1
i
n
i
1
i
i
1
i
i
1
i
0
i
n
1
i
1
i
1
0
i
)
x
x
(
)
x
x
)(
x
x
(
)
x
x
)(
x
x
(
)
x
x
(
)
x
x
)(
x
x
(
)
x
x
)(
x
x
(
y
)
x
(
W
lub krócej
n
1
i
i
j
j
i
j
i
)
x
x
(
)
x
x
(
y
)
x
(
W
, j = 0, 1, ..., n
Przykład:
Wyznaczyć wielomian interpolacyjny Lagrange’a funkcji f (x) = e
x
w przedziale [0,2 ; 0,5]
mając dane:
f (0,2) = 1,2214, f (0,4) = 1,4918, f (0,5) = 1,6487
009
,
1
91
,
0
724
,
0
1
,
0
3
,
0
)
4
,
0
)(
2
,
0
(
6487
,
1
)
1
,
0
(
2
,
0
)
5
,
0
)(
2
,
0
(
4918
,
1
)
3
,
0
(
2
,
0
)
5
,
0
)(
4
,
0
(
2214
,
1
)
(
2
x
x
x
x
x
x
x
x
x
W
algorytm do wyznaczania wielomianu Lagrange’a dla podanego punktu:
begin
fx:=0;
for i:=0 to n do
begin
p:=1;
for k:=0 to n do
if k<>i then p:=p*(xx-x[k])/(x[i]-x[k]);
fx:=fx+f[i]*p
end;
end;
Należy pamiętać, ze przed przystąpieniem do obliczeń należy sprawdzić czy w wektorze x
wartości się nie powtarzają, ponieważ grozi to wystąpieniem błędu dzielenia przez zero.
Różnice skończone
Dla funkcji stabelaryzowanej przy stałym kroku h = x
i+1
– x
i
można wyprowadzić pojęcie
różnicy skończonej rzędu k
k
0
j
1
k
i
k
j
j
i
1
k
1
i
1
k
i
1
k
i
k
i
1
i
2
i
i
1
i
i
i
2
i
1
i
i
y
1
y
y
y
y
...
y
y
2
y
y
y
y
y
y
y
y
12
Przykład: Dla funkcji danej w postaci tablicy zbudować tablicę różnic skończonych.
NR
x
y
Δy
Δ
2
y
Δ
3
y
Δ
4
y
Δ
5
y
0
1,2
1,728
0,469
0,078
0,006
0
0
1
1,3
2,197
0,547
0,084
0,006
0
2
1,4
2,744
0,631
0,090
0,006
3
1,5
3,375
0,721
0,096
4
1,6
4,096
0,817
5
1,7
4,913
Wzory interpolacyjne dla argumentów równoodległych
Dla zbioru węzłów tworzących ciąg arytmetyczny x
0
, x
1
= x
0
+ h, x
2
= x
0
+ 2h, ..., x
n
= x
0
+ nh
dane są wartości funkcji f(x
0
), f(x
1
), ..., f(x
n
). Szukany wielomian interpolacyjny zapisać
możemy w następujący sposób:
W(x) = a
0
+ a
1
q + a
2
q(q – 1) + a
3
q(q – 1)(q – 2) + ... + a
n
q(q – 1)(q – 2)...(q – n +1),
gdzie
h
x
x
q
0
Dla:
x = x
0
: q = 0
x = x
1
: q = 1
x = x
2
: q = 2
...
x = x
n
: q = n
Funkcje bazowe dla wielomianu W(x) przyjęto w postaci:
Φ
0
(x) = 1
Φ
1
(x) = q
Φ
2
(x) = q(q – 1)
Φ
3
(x) = q(q – 1)(q – 2)
Φ
n
(x) = q(q – 1)(q – 2)...(q – n + 1)
Otrzymujemy następujący układ równań:
n
3
2
1
0
n
3
2
1
0
y
y
y
y
y
a
a
a
a
a
!
n
)
2
n
)(
1
n
(
n
)
1
n
(
n
n
1
0
6
6
3
1
0
0
2
2
1
0
0
0
1
1
0
0
0
0
1
z którego wyznacza się wartości nieznanych współczynników a
0
, a
1
, a
2
, ..., a
n
13
!
n
y
a
y
a
!
n
...
a
)
1
n
(
n
na
a
...
...
!
3
y
a
y
a
6
a
6
a
3
a
!
2
y
a
y
a
2
a
2
a
y
a
y
a
a
y
a
0
n
n
n
n
2
1
0
0
3
3
3
3
2
1
0
0
2
2
2
2
1
0
0
1
1
1
0
0
0
Ostatecznie otrzymujemy I wzór interpolacyjny Newtona w postaci
0
n
0
2
0
0
y
!
n
)
1
n
q
)...(
1
q
(
q
...
y
!
2
)
1
q
(
q
y
q
y
)
x
(
W
Przykład: Dla zależności f(T)=T log p znaleźć wielomian interpolacyjny stopnia 3 i obliczyć
odchyłki w węzłach interpolacji
T
0,8
1,0
1,2
1,4
1,6
1,8
T log p
0,3566
0,9383
1,5598
2,2169
2,9059
3,6229
Tablica:
T
T log p
Δ
Δ
2
Δ
3
błąd [%]
0,8
0,3566
0,5817
0,0398
-0,0042
0
1,0
0,9383
0,6215
0,0356
-0,0037
0
1,2
1,5598
0,6571
0,0319
-0,0039
0
1,4
2,2169
0,6890
0,0280
0
1,6
2,9059
0,7170
1,9
1,8
3,6229
4,2
przyjmijmy
4
T
5
2
,
0
8
,
0
T
q
8
,
0
T
0
1,5002
1,791T
0,7225T
0,075T
W(T)
2)
1)(q
0,0007q(q
1)
0,0189q(q
0,5817q
0,3566
W(T)
2
3
Zauważmy, że I wzór Newtona daje dobrą dokładność w pobliżu punktu T
0
, natomiast dla
punktów leżących niżej w tabeli błąd wzrasta. Jeśli chcielibyśmy wyzerować odchyłki we
wszystkich węzłach tabeli trzeba by podwyższyć rząd interpolacji i w efekcie wzór
empiryczny staje się wtedy praktycznie mało przydatny. Wzory interpolacyjne stosuje się też
do obliczania wartości funkcji w punktach pośrednich tabeli (do zagęszczenia). Aby obliczyć
T log p dla T = 1,3 z dokładnością do Δ
2
wystarczy przyjąć T
0
= 1,2, a wówczas q = 0,5 i
znajdujemy
88437
,
1
0319
,
0
25
,
0
6571
,
0
5
,
0
5598
,
1
)
3
,
1
(
W
.
14
Zadanie takich wartości T = 1,65 byłoby już niewykonalne, ponieważ brakuje różnic. Widać
tu ograniczenie tzw. interpolacji w przód, określonej I wzorem Newtona. Do interpolacji
wstecz służy II wzór Newtona. Szukamy wielomian interpolacyjny
W(x) = a
0
+ a
1
q + a
2
q(q + 1) + a
3
q(q + 1)(q + 2) + ... + a
n
q(q + 1)(q + 2)...(q + n – 1),
gdzie
h
x
x
q
n
Współczynniki a
0
, a
1
, a
2
, ..., a
n
wyznaczamy jak poprzednio i dochodzimy do wzoru:
0
2
2
1
0
!
)
1
)...(
1
(
...
!
2
)
1
(
)
(
y
n
n
q
q
q
y
q
q
y
q
y
x
W
n
n
n
Przykład: Obliczyć T log p dla T = 1,65 z dokładnością do Δ
2
75
,
0
2
,
0
8
,
1
T
q
8
,
1
T
n
)
65
,
1
(
W
3,6229-0,75*0,717-0,75*0,25*0,028=3,0825
for j:=1 to 2 do
for i:=1 to n-j do
z[i,j]:=z[i+1,j-1]-z[i,j-1];
q0:=(x-x0)/h;
k:=int(q,0);
q:=q0-k;
if k>=n then write (‘Brak danych’)
else wart:=z[k,0]+q*z[k,1]+0.5*q*(q-1)*z[k,2];
Interpolacja trygonometryczna
Załóżmy, że znamy wartości pewnej ciągłej i okresowej funkcji f(x) o okresie 2π w 2n+1
węzłach. Jako bazę interpolacji przyjmujemy zbiór funkcji trygonometrycznych
)
nx
cos(
),
nx
sin(
),...,
x
2
cos(
),
x
2
sin(
),
x
cos(
),
x
sin(
,
2
1
Otrzymujemy zatem wielomian interpolacyjny w postaci
)
nx
cos(
a
)
nx
sin(
b
...
)
x
2
cos(
a
)
x
2
sin(
b
)
x
cos(
a
)
x
sin(
b
2
a
)
x
(
W
n
n
2
2
1
1
0
zawierający 2n+1 nieznanych parametrów.
Ze względu na uproszczenie obliczeń najistotniejszy jest przypadek interpolacji funkcji
określonej na zbiorze równoodległych węzłów x
i
[0,2π] dobranych według następującej
zależności
15
1
n
2
i
2
x
i
, gdzie i = 0, 1, ..., 2n
Czyli
1
n
2
n
4
x
...,
,
1
n
2
2
x
,
0
x
n
1
0
Warunek interpolacji prowadzi do układu równań liniowych w postaci
n
1
0
n
1
0
1
1
2n
2n
1
1
1
1
y
...
y
y
a
...
b
a
nx
cos
nx
sin
...
x
cos
x
sin
2
1
...
...
...
...
...
...
nx
cos
nx
sin
...
x
cos
x
sin
2
1
1
0
...
1
0
2
1
Współczynniki pierwszego wiersza macierzy X wynikają z wartości funkcji sin(hx) i cos (hx)
dla x
0
= 0. Przedstawiony układ równań rozwiązuje się natychmiastowo, ponieważ macierz X
-
1
można wyznaczyć z zależności
T
1
X
1
n
2
2
X
Przykład: Daną funkcję f(x) = 7x – x
2
przybliżyć wielomianem trygonometrycznym
przyjmując n = 3. Współrzędne węzłów interpolacji obliczamy ze wzoru:
1
2n
2iπ
x
i
x
0
0,898
1,795
2,693
3,590
4,488
5,386
y
0
5,478
9,344 11,598 12,242 11,274 8,695
Baza interpolacji – zbiór funkcji:
)
x
3
cos(
),
x
3
sin(
),
x
2
cos(
),
x
2
sin(
),
x
cos(
),
x
sin(
,
2
1
16
Tworzymy macierz X:
901
,
0
434
,
0
223
,
0
975
,
0
623
,
0
782
,
0
707
,
0
623
,
0
782
,
0
901
,
0
434
,
0
223
,
0
975
,
0
707
,
0
223
,
0
975
,
0
623
,
0
782
,
0
901
,
0
434
,
0
707
,
0
223
,
0
975
,
0
623
,
0
782
,
0
901
,
0
434
,
0
707
,
0
623
,
0
782
,
0
901
,
0
434
,
0
223
,
0
975
,
0
707
,
0
901
,
0
434
,
0
223
,
0
975
,
0
623
,
0
782
,
0
707
,
0
1
0
1
0
1
0
707
,
0
Elementy macierzy X mnożymy przez 2/7, otrzymaną macierz transponujemy i obliczamy
współczynniki wzoru interpolacyjnego
T
491
,
1
,
147
,
0
,
961
,
1
,
513
,
0
,
923
,
4
,
336
,
1
,
845
,
11
A
17
Aproksymacja
jest to przybliżanie funkcji f(x) zwanej aproksymowaną inną funkcją Q(x) zwaną funkcją
aproksymującą. Aproksymacja bardzo często występuje w dwóch przypadkach:
gdy funkcja aproksymowana jest przedstawiona w postaci tablicy wartości i
poszukujemy dla niej odpowiedniej funkcji ciągłej
gdy funkcję o dosyć skomplikowanym zapisie analitycznym chcemy przedstawić w
„prostszej” postaci
Dokonując aproksymacji funkcji f(x) musimy rozwiązać dwa ważne problemy:
dobór odpowiedniej funkcji aproksymującej Q(x)
określenie dokładności dokonanej aproksymacji
Dobór odpowiedniej funkcji aproksymującej Q(x)
Najczęściej stosowane funkcje aproksymujące są dobierane w postaci wielomianów
uogólnionych będących kombinacją liniową funkcji q(x)
m
0
j
j
j
m
)
x
(
a
)
x
(
Q
Jako funkcje bazowe stosowane są:
jednomiany
funkcje trygonometryczne
wielomiany ortogonalne
Przyjęcie odpowiednich funkcji bazowych powoduje, że aby wyznaczyć funkcję
aproksymującą należy wyznaczyć wartości współczynników, a
0
, a
1
, ..., a
m
.
Określenie dokładności aproksymacji
Aproksymacja funkcji powoduje powstanie błędów i sposób ich oszacowania wpływa na
wybór metody aproksymacji. Jeśli błąd będzie mierzony na dyskretnym zbiorze punktów x
0
,
x
1
, ..., x
n
to jest oto aproksymacja punktowa, jeśli będzie mierzony w przedziale (a,b) –
aproksymacja integralna lub przedziałowa.
Najczęściej stosowanymi miarami błędów aproksymacji są:
dla aproksymacji średniokwadratowej punktowej
n
0
i
2
)
x
(
Q
)
x
(
f
S
dla aproksymacji średniokwadratowej integralnej lub przedziałowej
b
a
2
dx
)
x
(
Q
)
x
(
f
S
dla aproksymacji jednostajnej
)
x
(
Q
)
x
(
f
sup
S
b
,
a
x
We wszystkich tych przypadkach zadanie aproksymacji sprowadza się do takiego
optymalnego doboru funkcji aproksymującej (dla wielomianu uogólnionych zaś do
optymalnego doboru współczynników a
0
, a
1
, ..., a
m
) aby zdefiniowane wyżej błędy były
minimalne.
18
Aproksymacja średniokwadratowa
Niech dana będzie funkcja y=f(x), która w pewnym zbiorze X punktów x
0
, x
1
, ..., x
n
przyjmuje wartości y
0
, y
1
, ..., y
n
. Wartości te znane tylko w przybliżeniu z pewnym błędem
(np. jako wyniki pomiarów). Poszukujemy takiej funkcji Q(x) przybliżającej daną funkcję
f(x), która umożliwi wygładzenie funkcji f(x), czyli pozwoli z zakłóconych błędami danymi
wartości funkcji przybliżonej otrzymać gładką funkcję przybliżającą.
Niech
j
(x), j=0, 1,...,n będzie układem funkcji bazowych. Poszukujemy wielomianu
uogólnionego Q(x) będącego najlepszym przybliżeniem średniokwadratowym funkcji f(x) na
zbiorze X=(x
j
). Funkcja przybliżająca ma postać
m
0
i
i
i
(x)
a
Q(x)
Przy czym współczynniki a
i
są tak określone, aby wyrażenie
n
0
i
2
i
i
i
)
x
(
Q
)
x
(
f
)
x
(
w
)
x
(
Q
)
x
(
f
0
)
x
(
w
i
dla i=0, 1, ..., n
było minimalne. Funkcja w(x) jest z góry ustaloną funkcją wagową.
Aby wyznaczyć współczynniki a
i
oznaczamy odchylenie
n
0
j
2
j
j
2
m
0
i
i
i
i
j
n
0
j
j
n
1
0
R
)
x
(
w
)
x
(
a
)
x
(
f
)
x
(
w
a
,...,
a
,
a
H
n
0
j
2
j
j
2
m
0
i
j
i
i
j
n
0
j
j
n
1
0
)R
w(x
)
(x
a
)
f(x
)
w(x
a
,...,
a
,
a
H
gdzie R
j
jest odchyleniem w punkcie x
j
.
Następnie obliczamy pochodne cząstkowe funkcji H względem a
i
. Z warunku
0
a
H
k
,
gdzie k = 0, 1, ..., n
otrzymujemy układ m+1 równań o niewiadomych a
i
zwany układem normalnym
0
)
x
(
)
x
(
a
)
x
(
f
)
x
(
w
2
a
H
j
k
m
0
i
i
i
i
j
n
0
j
j
k
,
gdzie k = 0, 1, ..., n
0
)
(x
)
(x
a
)
f(x
)
w(x
2
a
H
j
k
m
0
i
j
i
i
j
n
0
j
j
k
Jeśli wyznacznik tego układu jest różny od zera to rozwiązaniem układu jest minimum funkcji
H. W zapisie macierzowym układ przyjmuje postać
F
D
DA
D
T
T
19
gdzie
n
m
n
0
1
m
1
0
0
m
0
0
x
...
x
...
...
...
x
...
x
x
...
x
D
m
1
0
a
...
a
a
A
)
a
(
f
...
)
x
(
f
)
x
(
f
f
n
1
0
Aproksymacja wielomianowa
Jeżeli za funkcje bazowe przyjmiemy ciąg jednomianów (x
i
), i = 0, 1, ..., n to układ normalny
przyjmuje postać
0
x
)
x
(
a
)
x
(
f
n
0
j
k
j
m
0
i
i
j
i
j
,
k = 0, 1, ..., n
co po zmianie kolejności sumowania daje
n
j
m
i
n
j
k
i
j
i
k
j
j
x
a
x
x
f
0
0
0
)
(
przykład: Odnaleźć zależność między x i y w postaci ax+by=1
x 1 3 4 6 8 9 11 14
y 1 2 4 4 5 7 8
9
Będziemy rozpatrywać odchylenia zarówno wartości x jak i y, ponieważ wiodą one do
różnych wyników. Zakładamy, że dane są obarczone błędami, wówczas minimalizujemy
wyrażenia
2
n
1
i
i
i
y
y
2
n
1
i
i
i
x
x
przy czym
1
y
b
x
a
i
i
n
1
i
i
1
n
1
i
i
0
y
a
x
na
,
gdzie:
b
a
a
1
,
b
1
a
o
n
1
i
i
i
1
n
1
i
2
i
0
n
1
i
i
y
x
a
x
a
x
2
n
1
i
i
n
1
i
2
i
n
1
i
n
1
i
i
i
i
n
1
i
n
1
i
2
i
i
0
x
x
n
y
x
x
x
y
a
2
n
1
i
i
n
1
i
2
i
n
1
i
n
1
i
i
i
n
1
i
i
i
1
x
x
n
x
y
y
x
n
a
20
Układ normalny dla danych z tablicy ma rozwiązanie
a
0
= 6/11
i
a
1
= 7/11,
więc ma postać
11y – 7x = 6
Minimalizując w ten sam sposób wyrażenie
2
n
1
i
i
i
x
x
otrzymamy równanie
2x - 3y = -1
Dla obu otrzymanych równań wykresy się pokrywają.
Algorytm
tablice: tab[1..4,1..5], a, p[1..4], s, suma[1..6]
zmienne rzeczywiste: x, y
zmienne indeksowe: i, j, n, k
for j:=1 to n do
begin
write (‘Podaj x i y dla węzła:’);
readln (x:10:4, y:10:4);
s[1]:=x;
for i:=2 to 6 do s[i]:=x*s[i-1];
for i:=1 to 6 do suma[i]:=suma[i]+s[i];
p[1]:=p[1]+y;
for i:=2 to 4 do p[i]:=p[i]+y*s[i-1];
end;
tab[1,1]:=n;
for i:=1 to 4 do
begin
tab[i,5]:=p[i];
for j:=1 to 4 do
begin
k:=i+j;
if not (k=2 then tab[i,j]:=suma[k-2]
end;
end;
RUKL(tab,a,4,5);
writeln (‘Poszukiwane współczynniki to’);
for i:=1 to 4 do writeln (a[i]:12:6);
Aproksymacja za pomocą wielomianów ortogonalnych ze wzrostem stopnia wielomianu,
obliczenia stają się coraz bardziej pracochłonne a ich wyniki niepewne. Problem ten można
usunąć stosując do aproksymacji wielomiany ortogonalne.
Funkcje f(x) i g(x) nazywa się ortogonalnymi na dyskretnym zbiorze punktów x
0
, x
1
, ..., x
n
jeśli
n
0
i
i
i
)
x
(
g
)
x
(
f
21
przy czym
n
0
i
2
0
)
x
(
f
n
0
i
2
0
)
x
(
g
Analogicznie ciąg funkcyjny
),...
x
(
),...,
x
(
),
x
(
)
x
(
m
1
0
m
nazywamy ortogonalnym na zbiorze punktów x
0
, x
1
, ..., x
n
jeśli
n
0
i
i
k
i
j
0
)
x
(
)
x
(
dla j≠k
Zastosowanie tej metody powoduje, że znika jedna z trudności obliczeniowych przy
aproksymacji wielomianowej, mianowicie złe uwarunkowanie macierzy układu normalnego.
Przy aproksymacji wielomianami ortogonalnymi macierz układu normalnego jest macierzą
diagonalną, a jej elementy położone na głównej przekątnej dane są wzorem
n
0
i
i
2
j
jj
)
x
(
a
Załóżmy, że znamy n+1 równoodległych punktów x
i
(x
i
= x
0
+ih, i = 0, 1, ..., n). Za pomocą
przekształcenia liniowego
h
x
x
q
0
przeprowadzimy te punkty w kolejne liczby całkowite od 0 do n poszukujemy ciągu
wielomianów
)
q
(
F
),...,
q
(
F
,
F
)
q
(
F
)
n
(
m
)
n
(
1
)
n
(
0
)
n
(
i
(dolny indeks oznacza stopień i m ≤ n) spełniających warunek ortogonalności
0
)
i
(
F
)
i
(
F
n
0
i
)
n
(
k
)
n
(
j
dla j≠k
przy czym
k
k
2
2
1
1
0
)
n
(
k
q
a
...
q
a
q
a
a
)
q
(
F
gdzie
)
1
k
q
)...(
1
q
(
q
q
k
22
Często używa się unormowanego ciągu wielomianów spełniających warunek
1
)
0
(
F
)
n
(
k
,
gdzie k = 0, 1, ..., m
]
k
[
k
]
1
[
1
)
n
(
k
q
b
...
q
b
1
)
q
(
F
Co po przekształceniu daje nam wzór na wielomiany Grama, zwane też wielomianami
Czebyszewa stopni k = 0, 1, ..., m w postaci
]
n
[
]
s
[
s
k
s
k
s
k
0
s
s
)
n
(
k
n
q
1
)
q
(
F
,
gdzie k = 0, 1, ..., m
Wzór aproksymacyjny oparty na wielomianach Grama ma postać:
m
0
k
m
0
k
0
)
n
(
k
k
)
n
(
k
k
m
h
x
x
F
s
c
)
q
(
F
s
c
y
gdzie,
n
m
oraz
n
0
q
2
)
n
(
k
k
)
q
(
F
s
)
x
(
F
y
c
i
)
n
(
n
0
i
i
k
Przykład:
Na podstawie tablicy wartości f(x) znaleźć najlepszy w sensie aproksymacji
średniokwadratowej wielomian przybliżający tę funkcję.
x
1
1,5
2
2,5
3
y
3
4,75
7
9,75
13
Ponieważ węzły są równoodległe posłużymy się wielomianem Grama. Konstruujemy
wielomian aproksymacyjny dla n = 4. Wartości F
k
(q), q=1, 2,…, 4.
Obliczamy ze wzoru:
]
[
]
[
0
)
(
)
1
(
)
(
n
S
k
S
s
n
k
h
q
s
s
k
s
k
q
F
1
)
(
0
q
F
n
q
q
F
2
1
)
(
1
1
n
23
)
1
n
(
n
)
1
q
(
q
6
n
q
6
1
)
q
(
F
2
2
n
)
2
n
)(
1
n
(
n
)
2
q
)(
1
q
(
q
20
)
1
n
(
n
)
1
q
(
q
30
n
q
12
1
)
q
(
F
2
3
n
)
3
n
)(
2
n
)(
1
n
(
n
)
3
q
)(
2
q
)(
1
q
(
q
70
)
2
n
)(
1
n
(
n
)
2
q
)(
1
q
(
q
140
)
1
n
(
n
)
1
q
(
q
90
n
q
20
1
)
q
(
F
2
4
n
Zestawienie wyników przedstawia tabela:
q
x
i
y
i
F
0
(q)
F
1
(q)
F
2
(q)
F
3
(q)
F
4
(q)
0
1
3
1
1
1
1
1
1
1,5
4,75
1
0,5
-0,5
-2
-4
2
2
7
1
0
-1
0
6
3
2,5
9,75
1
-0,5
-0,5
2
-4
4
3
13
1
-1
1
-1
1
Następnie korzystając ze wzorów:
n
0
q
2
)
n
(
k
k
)
q
(
F
s
)
x
(
F
y
c
i
)
n
(
n
0
i
i
k
Obliczamy c
i
i s
i
oraz b
i
= c
i
/ s
i
i ze wzoru:
m
j
j
j
m
x
a
x
Q
0
)
(
)
(
otrzymujemy:
c
i
37,5
-12,5
1,75
0
0
s
i
5
2,5
3,5
10
70
b
i
7,5
-5
0,5
0
2
dla
5
,
0
1
x
q
24
Aproksymacja trygonometryczna
Często spotykamy się z przypadkiem, gdy funkcja f(x) jest okresowa. Taką funkcję
wygodniej jest aproksymować, wielomianem trygonometrycznym o bazie:
1, sin x, cos x, sin 2x, cos 2x, …, sin kx, cos kx
Jeżeli f(x) jest funkcją dyskretną określoną w dyskretnym zbiorze równoodległych punktów i
ich liczba jest parzysta i wynosi 2L (dla nieparzystej liczby punktów rozumowanie jest
analogiczne), niech:
L
i
x
i
dla i = 0, 1, …, 2L-1
Baza jest ortogonalna nie tylko na przedziale <0, 2π>, ale też na zbiorze punktów x
i
, przy
czym warunki ortogonalności mają postać:
0
k
m
dla
0
0
k
m
dla
L
k
m
dla
0
)
kx
sin(
)
mx
sin(
1
L
2
0
i
i
i
0
k
m
dla
L
2
0
k
m
dla
L
k
m
dla
0
)
kx
cos(
)
mx
cos(
1
L
2
0
i
i
i
0
)
kx
sin(
)
mx
cos(
1
L
2
0
i
i
i
dla m, k dowolnych, przy czym m i k zmieniają się od 0 do L
Przybliżeniem funkcji f(x) na zbiorze punktów x
i
jest wielomian trygonometryczny:
n
1
j
j
j
0
n
))
jx
sin(
b
)
jx
cos(
a
(
a
2
1
)
x
(
y
, n<L
Współczynniki a
j
i b
j
wielomianu wyznaczamy tak, aby suma kwadratów różnic
1
L
2
0
i
2
i
n
i
)
x
(
y
)
x
(
f
była minimalna.
25
Korzystając z warunku ortogonalności otrzymujemy rozwiązanie układu normalnego
0
)
x
(
)
x
(
a
)
x
(
f
)
x
(
w
2
a
H
j
k
n
0
j
j
i
m
0
i
i
j
j
k
w postaci:
1
L
2
0
i
1
L
2
0
i
i
i
i
j
L
ij
cos
)
x
(
f
L
1
jx
cos
)
x
(
f
L
1
a
dla j=1, 2, …, n
1
L
2
0
i
1
L
2
0
i
i
i
i
j
L
ij
sin
)
x
(
f
L
1
jx
sin
)
x
(
f
L
1
b
Algorytm:
const max: … {zdefiniowanie stałej}
type zakres: 1…max;
wektor: array[zakres] of real;
procedure APR (N, K: zakres; var y, a, b: wektor);
var i, j: integer;
PI, A0, S1, S2: real;
begin
write(‘Podaj wartość y’);
for i:=1 to N do readln(Y[I]);
A0:=0.0;
for j:=1 to N do A0:=A0+Y[I];
A0:=A0/N;
PI:=4.0*arctan(1.0)/N;
for i:=1 to K do
begin
S1:=0.0;
S2:=0.0;
for j:=1 to N do
begin
S1:=S1+Y[J]*cos(PI*J*I);
S2:=S2+Y[J]*sin(PI*J*I);
end;
A[I]:=S1*2/N;
B[I]:=S2*2/N;
end;
end;
26
Szybka transformata Fouriera
Każdą funkcję w tablicy tablicą wartości w n punktach
β
n
2π
x
β
, gdzie
= 0, 1, ..., n-1
można przybliżać wielomianem trygonometrycznym w postaci
1
n
0
k
k
jkx
exp
c
)
x
(
f
,
j - czynnik urojony
lub
m
0
k
1
m
k
k
0
x
)
1
m
(
cos
a
2
1
)
kx
sin(
b
)
kx
cos(
a
2
a
)
x
(
f
gdzie
= 1
2
2
n
m
dla n parzystych
= 1
2
1
n
m
dla n nieparzystych
1
0
1
0
1
0
)
exp(
)
(
1
1
-
n
...,
1,
0,
k
dla
)
sin(
)
(
2
1
)
cos(
)
(
2
1
n
k
n
k
n
k
jkx
x
f
n
c
kx
x
f
n
b
kx
x
f
n
a
Aproksymacja za pomocą funkcji sklejonych stopnia 3 (funkcji giętych)
Każdą funkcję
)
(
)
(
3
n
S
x
s
można przedstawić w postaci kombinacyjnej liniowej
1
n
1
i
3
i
i
)
x
(
c
)
x
(
s
,
a ≤ x ≤ b oraz gdzie
3
i
są funkcjami określonymi wzorem:
pozostaych
dla
0
x
,
x
x
dla
)
x
x
(
x
,
x
x
dla
)
x
x
(
3
)
x
x
(
h
3
)
x
x
(
h
3
h
x
,
x
x
dla
)
x
x
(
3
)
x
x
(
h
3
)
x
x
(
h
3
h
x
,
x
x
dla
)
x
x
(
n
1
2
i
1
i
3
2
i
1
i
i
3
1
i
2
1
i
1
i
2
3
i
1
i
3
1
i
2
1
i
1
i
2
3
1
i
2
i
3
2
i
3
3
i
gdzie
ih
a
x
i
,
n
a
b
h
Aproksymacja na przedziale <a,b>
Współczynniki c
i
dobieramy tak, aby odchylenie kwadratowe dane poniższym wzorem było
jak najmniejsze:
27
dx
)
x
(
C
)
x
(
f
I
2
b
a
1
n
1
i
3
i
i
w celu wyznaczenia minimum funkcji I=I(c
-1
, c
0
, … , c
n+1
) obliczamy pochodne:
j
c
I
j = -1 , 0, … , n+1
i przyrównujemy je do zera. Otrzymane równania tworzą układ n+3 równań z n+3
niewiadomymi c
i
1
n
1
i
b
a
b
a
3
j
3
j
3
i
i
dx
)
x
(
)
x
(
f
dx
)
x
(
)
x
(
C
j= -1, 0, 1, … ,n+1
Funkcje
)
x
(
3
i
są liniowo niezależne na przedziale <a,b>, więc układ równań ma
rozwiązanie w postaci punktu. Funkcji I osiąga minimum. Układ równań możemy zapisać w
postaci:
1
n
1
i
b
a
3
j
i
ij
dx
)
x
(
)
x
(
f
h
1
c
a
j = -1, 0, 1, … , n+1,
gdzie
b
a
3
j
3
i
ij
dx
)
x
(
)
x
(
h
1
a
Aproksymacja funkcji określonej na dyskretnym zbiorze punktów
Niech (x
i
’) dla i= 0, 1, … ,n
1
; n
1
> n + 3 będzie zbiorem punktów, w których dane są wartości
funkcji f(x). Szukane jest minimum wyrażenia:
2
n
0
k
1
n
1
i
k
3
i
i
1
)
'
x
(
c
)
x
(
f
J
W celu wyznaczenia c
i
obliczamy pochodne cząstkowe funkcji J = J (c
-1
, c
0
,… ,c
n+1
)
względem zmiennych c
i
i przyrównujemy je do zera.
Otrzymane równania tworzą układ:
1
n
1
i
n
0
k
k
3
i
k
i
ij
1
)
'
x
(
)
'
x
(
f
c
b
, gdzie
1
n
0
k
k
3
j
k
3
i
ij
)
'
x
(
)
'
x
(
b
28
Jeżeli wyznacznik tego układu jest różny od zera to rozwiązując go możemy wyznaczyć
współczynniki c
j
takie, że funkcja J osiągnie minimum.
Metody numeryczne rozwiązywania układów algebraicznych równań
liniowych
.
Rozpatrujemy układ n równań liniowych zawierających n niewiadomych
n
n
nn
2
2
n
1
1
n
2
n
n
2
2
22
1
21
1
n
n
1
2
12
1
11
b
x
a
...
x
a
x
a
...
...
...
...
...
...
...
b
x
a
...
x
a
x
a
b
x
a
...
x
a
x
a
Układ ten można zapisać także w postaci macierzowej A · X = B, gdzie
n
2
1
n
2
1
nn
2
n
1
n
n
2
22
21
n
1
12
11
x
...
x
x
X
b
...
b
b
B
a
...
a
a
...
...
...
...
a
...
a
a
a
...
a
a
A
metody dokładne: metoda Cramera, metoda eliminacji Gaussa, metoda Crouta (LU)
metody niedokładne: iteracja prosta, Gaussa-Siedla, metoda sukcesywnej nadrelaksacji (SOR)
O wyborze metody decyduje postać macierzy współczynników A. Jeśli macierz A jest
macierzą pełną to na ogół stosujemy metody dokładne. Jeśli macierz współczynników A jest
macierzą niepełną (znacząca ilość współczynników jest równa zero) wtedy stosujemy metody
iteracyjne.
Macierz A nazywana jest macierzą główną układu, X – wektorem niewiadomych, B –
wektorem wyrazów wolnych. Jeśli macierz główna nie jest osobliwa (det A ≠ 0), to układ
równań jest oznaczony (posiada dokładnie jedno rozwiązanie).
Gdy macierz A posiada niezerowe elementy tylko na głównej przekątnej, to układ równań
rozwiązuje się natychmiastowo.
n
n
nn
2
2
22
1
1
11
b
x
a
b
x
a
b
x
a
Niewiadome można wtedy obliczyć
ii
i
i
a
b
x
, a
ii
≠ 0, i = 0, 1, ..., n.
29
Łatwo rozwiązuje się trójkątne układy równań
n
n
nn
2
n
n
2
2
22
1
n
n
1
2
12
1
11
b
x
a
...
...
...
b
x
a
...
x
a
b
x
a
...
x
a
x
a
Z ostatniego równania możemy wyznaczyć x
n
, z przedostatniego x
n-1
nn
n
n
a
b
x
a ogólnie
ii
n
1
i
s
s
is
i
i
a
x
a
b
x
, i = n-1, n-2, ..., 1 przy założeniu, że a
ii
≠ 0, i = 1, 2,
..., n
PRZYKŁAD:
6
1
6
x
x
x
2
0
0
1
2
0
1
1
1
3
2
1
zmienną x
3
mamy daną wprost - x
3
= 3
1
a
x
a
x
a
b
x
2,
a
x
a
b
x
11
3
13
2
12
1
1
22
3
23
2
2
ALGORYTM:
zmienne indeksowe
n, i, j, s
zmienne rzeczywiste
suma
tablice
a[1..n,1..n], b[1..n], c[1..n]
for i:=n downto 1 do
begin
suma:=0;
for s:=i+1 to n do
suma:=suma+a[i,s]*x[s];
x[i]:=(b[i]-suma)/a[i,i];
end;
Wzory Cramera
Jesli oznaczymy symbolem W wyznacznik główny układu równań
n
n
nn
2
2
n
1
1
n
2
n
n
2
2
22
1
21
1
n
n
1
2
12
1
11
b
x
a
...
x
a
x
a
...
...
...
...
...
b
x
a
...
x
a
x
a
b
x
a
...
x
a
x
a
30
nn
2
n
1
n
n
2
22
21
n
1
12
11
a
...
a
a
...
...
...
...
a
...
a
a
a
...
a
a
A
det
W
n
2
n
1
n
2
22
21
1
12
11
j
b
...
a
a
...
...
...
...
b
...
a
a
b
...
a
a
W
to można wykazać, że
W
W
x
j
j
, W ≠ 0, j = 1, 2, ..., n
Metoda ta należy do metod dokładnych. Ze względu na dużą złożoność obliczeniową
praktycznie stosowana do numerycznego rozwiązywania równań.
Metoda eliminacji Gaussa
Metoda polega na zapisaniu układu równań w postaci macierzy C, której n pierwszych
kolumn zawiera elementy a
ij
macierzy głównej A, natomiast kolumnę n+1 tworzą dowolne
wyrazy b
i
.
Zakładając, że C
11
0, odejmujemy pierwsze równanie pomnożone przez
11
1
C
C
i
od i-tego
równania (i=2,3,...,n) a obliczonymi współczynnikami zastępujemy przednie wartości. W
wyniku tej operacji otrzymujemy układ równań i odpowiadającą mu macierz C
1
.
)
1
(
1
,
)
1
(
3
)
1
(
3
2
)
1
(
2
)
1
(
1
,
3
)
1
(
3
3
)
1
(
33
2
)
1
(
32
1
,
2
2
3
23
2
22
1
,
1
1
3
13
2
12
1
11
...
...
...
...
n
n
n
nn
n
n
n
n
n
n
n
n
n
n
n
C
x
C
x
C
x
C
C
x
C
x
C
x
C
C
x
C
x
C
x
C
C
x
C
x
C
x
C
x
C
Za pomocą wzorów określających nowe współczynniki
j
i
ij
ij
C
C
C
C
C
1
11
1
)
1
(
i=2,3,...,n j=2,3,...,n+1
Jeśli
0
)
1
(
22
C
, to odejmiemy drugie równanie układu pomnożone przez
)
1
(
22
)
1
(
2
C
C
i
od i-tego
równania układu i=3,4,...,n. W wyniku otrzymujemy układ:
)
2
(
1
,
)
2
(
)
2
(
3
)
2
(
1
,
3
)
2
(
3
)
2
(
33
)
1
(
1
,
2
)
1
(
2
)
1
(
23
)
1
(
22
1
,
1
1
13
12
11
2
...
0
0
:
:
...
:
:
:
...
0
0
...
0
...
n
n
nn
n
n
n
n
n
n
n
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
A nowe współczynniki wyrażone są wzorami
)
1
(
2
)
1
(
22
)
1
(
2
)
1
(
)
2
(
j
i
ij
ij
C
C
C
C
C
i=3, 4, ..., n j=3, 4, ..., n
Kontynuując takie postępowanie po wykonaniu n kroków dochodzimy do układu trójkątnego
31
)
1
(
1
,
)
1
(
)
2
(
1
,
1
)
2
(
3
3
)
2
(
33
)
1
(
1
,
1
)
1
(
2
3
)
1
(
23
2
)
1
(
22
1
,
1
1
3
13
2
12
1
11
...
...
...
n
n
n
n
n
nn
n
n
n
n
n
n
n
n
n
C
x
C
C
x
C
x
C
C
x
C
x
C
x
C
C
x
C
x
C
x
C
x
C
)
1
(
1
,
)
1
(
)
2
(
1
,
3
)
2
(
3
)
2
(
33
)
1
(
1
,
2
)
1
(
2
)
1
(
23
)
1
(
22
1
,
1
1
13
12
11
2
0
0
0
0
0
0
n
n
n
n
nn
n
n
n
n
n
n
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
Algorytm dochodzenia od układu trójkątnego do układu równań liniowych
)
1
(
)
1
(
)
1
(
)
1
(
)
(
,...,
2
,
1
1
,...,
2
,
1
s
sj
s
ss
s
is
s
ij
s
ij
C
C
C
C
C
n
s
s
i
n
s
Przykład
4
2
3
5
2
3
1
4
3
2
3
2
1
3
2
1
3
2
1
x
x
x
x
x
x
x
x
x
4
2
3
1
5
1
2
3
1
4
3
2
C
Obliczamy elementy macierzy C
1
2
9
0
2
9
2
13
7
2
13
14
11
31
34
)
1
(
34
13
11
31
33
)
1
(
33
12
11
31
32
)
1
(
32
14
11
21
24
)
1
(
24
13
11
21
23
)
1
(
23
12
11
21
22
)
1
(
22
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
32
druki krok s=2
0
13
63
)
1
(
24
)
1
(
22
)
1
(
32
)
1
(
34
)
2
(
34
)
1
(
23
)
1
(
22
)
1
(
32
)
1
(
33
)
2
(
33
C
C
C
C
C
C
C
C
C
C
0
13
63
0
0
2
13
7
2
13
0
1
4
3
2
2
C
co daje nam macierz C
1
:
2
9
0
2
9
0
2
13
7
2
13
0
1
4
3
2
1
C
i dalej ze wzorów dla układów trójkątnych mamy x
3
=0, x
2
=-1, x
1
=1
ALGORYTM:
zmienne indeksowe
n, i, j, s
zmienne rzeczywiste
suma
tablice
c[1..n,1..n+1], x[1..n]
for s:=1 to n-1 do
for i:=s+1 to n do
for j:=s+1 to n+1 do
c[i,j]:=c[i,j]-(c[i,s]*c[s,j])/c[s,s];
x[n]:=c[n,n+1]/c[n,n]
for i:=n-1 downto 1 do
begin
suma:=0;
for s:=i+1 to n do
suma:=suma+c[i,s]*x[s];
x[i]:=(c[i,n+1]-suma)/c[i,i];
end;
Metody iteracyjne
Proces iteracyjnego obliczania wrtości x polega na przyjęciu pewnej wartości początkowej
X
o
, zwanej przybliżeniem początkowym, a następnie wykonaniu określonych operacji
arytmetycznych dających w wyniku X
1
(zwanym pierwszym przybliżeniem). Kolejne
przybliżenia oblicza się wykonując te same operacje na ostatnio obliczonym przybliżeniu.
33
Jeżeli ciąg przybliżeń {X
n
} zmierza do X, to proces iteracyjny jest zbieżny i tylko wtedy
metoda iteracyjna jest efektywna. Metody iteracyjnego rozwiązywania układów równań
liniowych można stosować dla układów typu
...
nn
n2
n1
2n
22
21
1n
12
11
n
2
1
w
w
w
w
w
w
w
w
w
X
X
X
Przykład: Sprowadzić układ do postaci AX=B
2x
1
+x
2
+x
3
=4
x
1
=-0,5x
2
-0,5x
3
+2
x
1
+x
2
-x
3
=1
x
2
=-x
1
+x
3
+1
x
1
+x
2
+x
3
=3
x
3
=-x
1
-2x
2
+4
W pierwszej kolejności macierz A zakładamy na dwie macierze D i R przy czym D zawiera
tylko główną przekątną macierzy A, natomiast R pozostałe jej elementy.
0
2
1
1
0
1
1
1
0
1
0
0
0
1
0
0
0
2
1
2
1
1
1
1
1
1
2
A
=
D
+
R
Równanie A*X=B sprawdziło się dla
D*X+R*X=B D*X=-R*X+B
Ostatnie równanie pomnożymy lewostronnie przez D
-1
i otrzymujemy:
E*X=D
-1
*R*X+D
-1
*B
Oznaczając –D
-1
*R=W i D
-1
*B=Z dochodzimy dorównania: X=W*X+Z
Wracając do zadania:
4
1
4
1
0
0
0
1
0
0
0
2
0
2
1
1
0
1
1
1
0
1
0
0
0
1
0
0
0
2
1
3
2
1
1
3
2
1
x
x
x
x
x
x
ponieważ
1
0
0
0
1
0
0
0
5
,
0
1
0
0
0
1
0
0
0
2
1
więc jak łatwo sprawdzić mamy taki sam układ jaki znaleziono sposobem „naturalnym”
34
Metoda iteracji prostej (Jakobiego):
Metoda ta dla równania X=W*X+Z polega na przyjęciu zerowego przybliżenia wektora
X=X
o
, a następnie przeprowadzenia obliczeń iteracyjnych za pomocą zależności:
X
i+1
=W*X
i
+Z i=0,1,...
Liczba kroków, które należy wykonać, aby uzyskać wymaganą dokładność rozwiązania, jest z
reguły dość duża i istotnie zależy od przyjęcia punktu startowego X
o
. Punkt ten najczęściej się
dobiera na podstawie dodatkowych informacji o fizycznych aspektach problemu.
Przykład: Metodę iteracji prostej wykorzystano do rozwiązania układu równań będącego
różnicową aproksymacją równania Laplace`a z warunkami brzegowymi I rodzaju. Równanie
to opisuje stacjonarne pole temperatury w pewnym dwuwymiarowym obszarze, na którego
brzegach temperatury są znane. Jako przybliżenie zerowe przyjęto wartość temperatur w
wnętrzu obszaru T=50
C, wymiary prostokątne a=b=0,8 cm.
50
0
0
0
0
0
0
0
0
100
50
38
38
38
38
38
25
0
100
63
50
50
50
50
50
38
0
100
63
50
50
50
50
50
38
0
100
63
50
50
50
50
50
38
0
100
63
50
50
50
50
50
38
0
100
63
50
50
50
50
50
38
0
100
75
63
63
63
63
63
50
0
100
100
100
100
100
100
100
100
50
0
0
0
0
0
18
14
10
6
0
32
26
19
10
0
42
35
26
14
0
50
42
32
18
0
58
50
39
22
0
68
61
50
31
0
100
94
90
86
82
78
69
50
0
100
100
100
100
100
100
100
100
50
35
Metoda Gaussa-Seidla
Stanowi ulepszenie metody iteracyjnej prostej polegające na wykorzystaniu obliczonych k
pierwszych składowych wektora X
i+1
do obliczenia składowej k+1. Modyfikacja ta znacznie
przyspiesza proces obliczania.
Przykład: Rozwiązując układ równań:
2
1
5
2
0
3
,
0
1
,
0
2
,
0
1
,
0
0
2
,
0
3
,
0
2
,
0
3
,
0
0
1
,
0
2
,
0
1
,
0
2
,
0
0
4
3
2
1
4
3
2
1
x
x
x
x
x
x
x
x
dla uproszczenia zapisu oznaczamy składowe z wektora iteracji „i” symbolem „x”, a
składowe wektora i iteracji i+1 symbolem „u”. W metodzie iteracji prostej układ ten
zapisujemy:
2
3
,
0
1
,
0
2
,
0
1
1
,
0
2
,
0
3
,
0
5
2
,
0
3
,
0
1
,
0
2
2
,
0
1
,
0
2
,
0
3
2
1
4
4
2
1
3
4
3
1
2
4
3
2
1
x
x
x
u
x
x
x
u
x
x
x
u
x
x
x
u
Przyjmujemy teraz punkt startowy np. X
0
= [3, 6, 4, 3]
T
, składowe tego wektora podstawiamy
do prawych stron równań wyliczając u
1
, ..., u
4
, które w następnej iteracji przejmą rolę punktu
startowego. Kolejne przybliżenia otrzymane tą metodą
78
,
3
13
,
4
33
,
7
53
,
4
37
,
3
05
,
4
27
,
7
45
,
4
15
,
3
4
08
,
7
4
,
4
2
,
3
4
,
3
1
,
7
2
,
4
3
4
6
3
rozwiązaniem dokładnym jest X = [4,59, 7,49, 4,20, 3,44]
Ten sam układ w metodzie Gaussa – Siedla ma postać:
2
3
,
0
1
,
0
2
,
0
1
1
,
0
2
,
0
3
,
0
5
2
,
0
3
,
0
1
,
0
2
2
,
0
1
,
0
2
,
0
3
2
1
4
4
2
1
3
4
3
1
2
4
3
2
1
u
u
u
u
x
u
u
u
x
x
u
u
x
x
x
u
36
a kolejne przybliżenia:
44
,
3
2
,
4
4
,
7
58
,
4
43
,
3
19
,
4
38
,
7
56
,
4
41
,
3
15
,
4
32
,
7
51
,
4
32
,
3
4
22
,
7
2
,
4
3
4
6
3
Formalny zapis algorytmu
Macierz W w równaniu X = W · X + Z przedstawiamy jako sumę dwóch macierzy: górnej
trójkątnej W
u
i dolnej trójkątnej W
l
. Macierz górna trójkątna składa się z elementów macierzy
W leżących nad główną przekątną (pozostałe są zerami), natomiast macierz dolna trójkątna
składa się z elementów macierzy W leżących pod główną przekątną. Macierz W w tym
przykładzie ma zerową główną przekątną
Algorytm metody Gaussa – Siedla realizowany jest następującym wzorem:
Z
X
W
X
W
X
i
i
i
u
i
1
1
czyli mamy:
2
1
5
2
0
3
,
0
1
,
0
2
,
0
0
0
2
,
0
3
,
0
0
0
0
1
,
0
0
0
0
0
0
0
0
0
1
,
0
0
0
0
2
,
0
3
,
0
0
0
2
,
0
1
,
0
2
,
0
0
4
3
2
1
4
3
2
1
x
x
x
x
u
u
u
u
Metoda sukcesywnej nadrelaksacji (SOR)
Jest ulepszeniem metody Gaussa-Seidla mającym poprawić szybkość zbieżności procesu
iteracyjnego. Istota polega na sukcesywnym wprowadzaniu w miejsce składowych u
i
po
prawej stronie układu wartości
x
i
+ α (u
i
– x
i
)
α – współczynnik nadrelaksacji, α
<1,2>
i – liczba niewiadomych a nie iteracja
Sposób dobory optymalnego współczynnika
nazywanego czynnikiem nadrelaksacyjnym
jest trudny do określenia, chociaż wiadomo, że jest on liczbą z przedziału [1 , 2]. Można
zauważyć, że dla
=1 otrzymuje się metodę Gaussa – Siedla. Według źródeł literaturowych
zaleca się przyjmowanie
rzędu 1,8.
Algorytm metody nadrelaksacyjnej sprowadza się do wzoru:
Z
X
X
X
W
X
W
X
i
i
i
i
i
u
i
)]
(
[
1
1
1
37
Zbieżność metod iteracyjnych
Jeżeli macierz A kwadratowa oraz x jest wektorem niezerowym, wówczas istnieje pewien
wektor równoległy do Ax, co oznacza, że istnieje pewna stała
taka, że:
Ax =
x
Możemy wtedy zapisać równanie charakterystyczne
(A -
I)x = 0,
gdzie I – macierz jednostkowa
oraz wielomian charakterystyczny ma postać
p(
) = det (A -
I)
Jeżeli wielomian macierzy A jest
n
n
, wektora x jest n, wówczas p jest wielomianem m-
tego stopnia, którego pierwiastki są pojedyncze i zespolone. Jeżeli wektor
jest
pierwiastkiem wielomianu p wówczas:
0
det
I
A
Jeżeli powyższy warunek jest spełniony, oznacza to, że równanie charakterystyczne ma
rozwiązanie x
0.
Zbieżność metod iteracyjnych (2)
Promień spektralny
(A) macierzy A jest zdefiniowany jako:
max
)
A
(
,
gdzie
- oznacza moduł wektora
wówczas równanie charakterystyczne jest zbieżne gdy:
1
)
A
(
Można wykazać, że:
A
macierzy
normę
oznacza
A
gdzie
A
A
)
(
ostatecznie dla metod iteracyjnych można zapisać twierdzenie:
Warunkiem koniecznym i wystarczającym zbieżności procesu iteracyjnego danego wzorem
Z
X
W
X
i
i
1
przy dowolnym wektorze początkowym
)
0
(
X
i dowolnym wektorze Z jest,
aby wszystkie wartości własne macierzy W były co do modułu mniejsze od jedności.
38
Zbieżność metod iteracyjnych
W praktyce jest na ogół wygodniej korzystać z twierdzenia:
Jeżeli norma macierzy W jest mniejsza od jedności to proces iteracyjny
określany wzorem X
i+1
= W · X
i
+ Z jest zbieżny.
Z przytoczonych twierdzeń wynika, że zbieżność procesu iteracyjnego nie zależy ani od
wartości początkowej
)
0
(
X
, ani też od wartości danego wektora Z, a wystarczy jedynie
spełnienie jednego z następujących warunków:
n
j
ij
i
w
1
1
max
norma wierszowa
n
i
ij
j
w
1
1
max
norma kolumnowa
n
i
n
j
ij
w
1
1
2
1
norma energetyczna
Kryterium „stopu” metod iteracyjnych
Obliczenia iteracyjne trwają do momentu kiedy zostanie spełniony warunek
i
1
i
X
X
gdzie: X
(i+1)
oraz X
(i)
są kolejnymi przybliżeniami rozwiązania X
i+1
= W · X
i
+ Z oraz δ jest
dokładnością obliczeń – dowolna liczba większa od 0
j
1
i
j
i
1
i
x
x
x
,
gdzie ε
<10
-1
, 10
-256
> – przyjęta dokładność obliczeń
39
Numeryczne rozwiązywanie równań nieliniowych
Metoda Newtona
Dany jest układ równań nieliniowych z n niewiadomymi
0
)
x
...
x
(
f
..
..........
..........
0
)
x
...
x
(
f
0
)
x
...
x
(
f
n
1
n
n
1
2
n
1
1
który możemy zapisać w postaci wektorowej
f(x)=0
gdzie:
n
x
x
x
:
:
1
)
(
:
:
)
(
)
(
1
x
f
x
f
x
f
n
O funkcjach fi(x
1
...x
n
) dla i = 1, 2, ..., n zakładamy, że mają ciągłe pochodne pierwszego
rzędu w pewnym obszarze zawierającym odosobniony pierwiastek układu równań. Niech
x
(k)
= {x
1
(k)
, ..., x
n
(k)
}
będzie k-tym przybliżeniem pierwiastka a = {a
1
, ..., a
n
} równania.
Dokładną wartość pierwiastka wyraża wzór a = x
(k)
+ ε
(k)
Gdzie
}
,...,
{
)
(
)
(
1
k
n
k
jest błędem pierwiastka przybliżonego
)
(k
x
. Skoro istnieje ciągła
pochodna funkcji f(x) możemy zapisać:
)
(
)
(
)
(
)
(
'
)
(
)
(
k
k
k
x
f
x
f
a
f
Zastępując błąd
)
(k
przyrostem
)
(
)
1
(
)
(
k
k
k
x
x
h
i porównując prawą stronę powyższego
wyrażenia do zera otrzymujemy równanie liniowe:
0
)
)(
(
'
)
(
)
(
)
1
(
)
(
)
(
k
k
k
k
x
x
x
f
x
f
Wzór ten stanowi zapis rekurencyjny dla metody Newtona w postaci macierzowej.
40
Pochodna f’(x) jest macierzą Jakobiego
n
n
2
n
1
n
n
2
2
2
1
2
n
1
2
1
1
1
x
f
...
x
f
x
f
...
...
...
...
x
f
...
x
f
x
f
x
f
...
x
f
x
f
)
x
(
W
)
x
(
'
f
Przyjmując, że
)
(
)
(k
x
W
jest macierzą nieosobliwą możemy równanie liniowe przekształcić do
postaci:
)
(
)
(
)
(
)
(
1
)
(
)
1
(
k
k
k
k
x
f
x
W
x
x
stąd przyjmując dowolną wartość
)
0
(
x
otrzymujemy ciąg kolejnych przybliżeń pierwiastka
równania
)
(
)
2
(
)
1
(
)
0
(
...,
,
,
,
k
x
x
x
x
uzyskanych metodą Newtona.
Metoda iteracji
Dany jest układ równań nieliniowych z n niewiadomymi
)
x
...
x
(
g
x
..
..........
..........
)
x
...
x
(
g
x
)
x
...
x
(
g
x
n
1
n
3
n
1
2
2
n
1
1
1
który możemy zapisać w postaci wektorowej
x=g(x)
gdzie:
n
x
x
x
:
:
1
)
(
:
:
)
(
)
(
1
x
g
x
g
x
f
n
Zakładamy, że funkcja g
1
, g
2
, ..., g
n
są rzeczywiste i ciągłe w pewnym otoczeniu
odosobnionego pierwiastka a = {a
1
, ..., a
n
} układu równań.
41
Metoda iteracji polega na stosowaniu następującego wzoru
x
k+1
= g(x
(k)
)
Po przyjęciu przybliżenia
)
0
(
x
, w rezultacie otrzymujemy ciąg kolejnych przybliżeń
...
,
,
,
)
2
(
)
1
(
)
0
(
x
x
x
pierwiastka a równania.
Metoda siecznych
Rozpatrujemy układ równań nieliniowych w postaci
0
)
x
...,
,
x
,
x
(
f
.
..........
..........
..........
0
)
x
...,
,
x
,
x
(
f
0
)
x
...,
,
x
,
x
(
f
n
2
1
n
n
2
1
2
n
2
1
1
Który możemy ogólnie zapisać w postaci:
f(x)=0
gdzie x jest wektorem n zmiennych.
Metoda iteracji polega na stosowaniu wzoru wyliczającego k-te przybliżenie i-tej zmiennej
tych:
)
(
)
(
)
(
)
(
1
1
1
k
i
k
i
k
i
k
i
k
i
k
i
k
i
x
f
x
f
x
f
x
x
x
x
Metoda ma dwie wady:
- potrzebne są 2 zestawy punktów startowych: x
0
, f(x
0
), x
1
, f(x
1
)
- wykrywamy jeden wektor x-ów (globalny) (jeden zestaw x-ów), a pierwiastki mogą być
wielokrotne
Przykład: Dany jest układ równań w postaci
0
2
y
x
0
4
x
2
Napisać program wyznaczający rozwiązanie tego układu.
Niezbędne deklaracje typów:
42
wektor array[1..2] of Extended;
zmienne:
x, X1, Xn, F, F1, Fn: Extended;
iteracja: LongInt;
Funkcja reprezentująca układ równań
function funkcja (z:word):wektor;
begin
funkcja[1]:=sqrt(z[1])-4;
funkcja[2]:=z[1]*z[2]-2;
end;
Realizacja obliczeń:
Procedure oblicz
Begin
x[1]:=4;
x[2]:=7;
X1[1]:=9;
X1[2]:=11;
Repeat
Inc(iteracje);
F:=funkcja(x);
F1:=funkcja(X1);
xn[1]:=x[1]
-
((x[1]-x1[1])*F[1])/F[1]-F1[1];
Xn[2]:=x[2]
-
((x[2]-X1[2])*F[2])/F[2]-F1[2];
Fn:=funkcja(xn);
X1:=x;
x:=Xn;
until (abs(Fn[1]+abs(Fn[2])<=0.0000000001;
Po zakończeniu obliczeń wektor x przechowuje wartości zmiennych, natomiast zmienna
iteracje zawiera liczbę wykonanych iteracji.
43
Rozważania dotyczą numerycznego obliczania wartości:
- całek pojedynczych funkcji ciągłych
- całek pojedynczych funkcji niewłaściwych
- całek funkcji osobliwych
- całek wielokrotnych
- całek Monte Carlo
Całkowanie pojedyncze – metody obliczeń
a) kwadratury z ustalonymi węzłami
kwadratury Newtona-Cotesa
złożone kwadratury Newtona-Cotesa
metody Romberga
metoda adaptacyjna
b) kwadratury Gaussa
- kwadratury Gaussa
- złożone kwadratury Gaussa
Całki pojedyncze właściwe
Całka oznaczona Riemanna
Niech funkcja f(x) będzie określona i ograniczona na przedziale [a,b]. Dokonajmy podziału
przedziału [a,b] takiego, że:
a = x
0
< x
1
< ... < x
n+1
= b
i utwórzmy sumę
1
c
f
x
x
S
n
0
i
i
i
1
i
n
,
gdzie
1
i
i
i
x
,
x
c
są dowolnymi punktami pośrednimi
Jeżeli ciąg {S
N
} dla
n
jest zbieżny do tej samej granicy S przy każdym przedziale
odcinka [a ; b] takim, że
0
max
1
0
i
i
n
i
x
x
niezależnie od wyboru punktów C
i
to funkcję f(x)
nazywamy całkowalną w przedziale [a ; b], a granicę S ciągu (1) nazywamy:
Całką oznaczoną Riemanna funkcji f o oznaczamy:
b
a
dx
x
f
)
(
Można wykazać, że funkcja ograniczona w przedziale domkniętym jest całkowalna wtedy i
tylko wtedy, gdy jej punkt nieciągłości w przedziale [a ; b] tworzą zbiór miary zero.
Jeżeli przez F(x) oznaczymy funkcję pierwotną funkcji f(x) w przedziale [a,b] (jeżeli F’(x) =
f(x)) to ma miejsce wzór
44
2
a
F
b
F
dx
x
f
f
I
b
a
3
f
E
S
f
I
n
, gdzie E(f) – reszta kwadratury
Interpretacją graficzną całki oznaczonej jest pole pomiędzy osią OX, krzywą.
Wzór (1) jest inaczej nazywany kwadraturą. Nazwa pochodzi od sumowania wyrazów ciągu,
które w najprostszym przypadku graficzną interpretację mają w postaci czworokątów.
45
Całka oznaczona Riemanna
Interpretacją graficzną całki oznaczonej jest pole powierzchni pomiędzy osią OX a krzywą.
Rysunek obok to ilustruje
Wyznaczmy S sumę ciągu n-wyrazów, które są określone jako iloczyn wartości funkcji w
punktach x
i
<a,b> oraz współczynników niezależnych od rodzaju funkcji.
Wówczas
n
i
i
i
x
f
a
f
S
0
)
(
)
(
(1)
b
a
dx
x
f
f
I
)
(
)
(
(2)
Można wykazać, że suma ciągu S(f) jest zbieżna do całki oznaczonej I(f) dla
n
46
Uwagi ogólne o całkowaniu numerycznym
Całkowanie numeryczne stosujemy wtedy gdy trudno obliczyć całkę w sposób analityczny
lub w przypadku gdy znane są tylko wartości funkcji podcałkowej w punktach zwanych
węzłami
W celu przybliżonego obliczenia całki oznaczonej funkcji f(x) na skończonym przedziale
]
;
[ b
a
x
, funkcję podcałkową f(x) zastępujemy interpolującą funkcją φ(x), którą łatwo
można całkować. Funkcja φ(x) będzie wielomianem interpolującym z węzłami interpolacji x
0
,
x1, …, x
n
. Wówczas całkę możemy zastąpić sumą, współczynniki a
i
zależą od metody
obliczeń
b
a
b
a
n
i
i
i
x
f
a
dx
x
dx
x
f
0
)
(
)
(
)
(
(4)
Kwadratury z ustalonymi węzłami
Kwadratury Newtona-Cotesa
Całkowanie z zastosowaniem kwadratur o ustalonych węzłach polega na zastąpieniu funkcji
podcałkowej wielomianem interpolacyjnym Lagrange’a Li(x).
Jeżeli dla skończonego przedziału [a,b] wybierzemy zbiór punktów węzłowych {x
0
, ..., x
n
}
takich, że:
x
i
= x
0
+ i · h
47
gdzie:
n
a
b
h
)
(
dla wówczas=(0, 1, …, n)
n - oznacza również stopień wielomianu interpolacyjnego
n
i
i
i
x
L
x
f
x
0
)
(
)
(
)
(
(5), gdzie:
n
k
i
k
k
i
k
i
x
x
x
x
x
L
0
)
(
)
(
)
(
(6)
wówczas
7
dx
!
1
n
x
f
x
x
dx
x
L
x
f
dx
x
f
b
a
n
0
i
b
a
n
0
i
1
n
i
i
i
b
a
8
dx
x
f
x
x
!
1
n
1
dx
x
f
a
b
a
n
0
i
b
a
n
0
i
1
n
i
i
i
ξ – dowolny punkt pośredni z przedziału (x,x
i
)
gdzie:
]
;
[
)
(
b
a
x
oraz
b
a
n
i
n
i
i
i
t
n
t
t
t
i
n
i
h
dx
x
L
a
0
)
)...(
1
(
)!
(
!
)
1
(
)
(
(9)
ostatecznie wzór kwadratury Newtona – Cotesa oraz reszta kwadratury
b
a
n
i
i
i
f
E
x
f
a
dx
x
f
0
)
(
)
(
)
(
(10)
b
a
n
i
n
i
dx
x
f
x
x
n
f
E
0
)
1
(
))
(
(
)
(
)!
1
(
1
)
(
(11)
Twierdzenie 1
Kwadratury Newtona-Cotesa oparte na n+1 węzłach są rzędu
ych
nieparzyst
n
dla
1
n
parzystych
n
dla
2
n
wówczas wzory kwadratur dla n parzystych mają postać
48
12
b
,
a
gdzie
,
dt
n
t
...
1
t
t
!
2
n
f
h
dx
x
f
a
dx
x
f
b
a
n
0
i
n
0
2
n
2
n
3
n
i
i
b
a
oraz dla nieparzystych mają postać:
b
a
n
n
i
n
n
n
i
i
dt
n
t
t
t
n
f
h
x
f
a
dx
x
f
0
0
)
1
(
)
2
(
)
)...(
1
(
)!
1
(
)
(
)
(
)
(
(13)
Kwadratury z ustalonymi węzłami – wzór trapezów
Wyznaczmy wzór kwadratury dla n = 1
14
f
E
dx
x
f
x
x
x
x
x
f
x
x
x
x
dx
x
f
1
0
x
x
1
0
1
0
0
1
0
1
b
a
gdzie:
1
1
3
1
1
)
)(
(
)
(
2
1
)
)(
))(
(
(
2
1
)
(
x
x
x
x
o
n
o
n
o
o
dx
x
x
x
x
f
dx
x
x
x
x
x
f
f
E
)
(
12
2
)
(
3
)
(
2
1
3
0
1
1
0
2
0
1
3
n
n
f
h
x
x
x
x
x
x
x
x
x
f
wówczas:
b
a
n
f
h
x
x
x
f
x
x
x
x
x
f
x
x
x
x
dx
x
f
)
(
12
)
(
)
(
2
)
(
)
(
)
(
2
)
(
)
(
3
0
1
1
0
1
2
0
0
1
0
2
1
49
)
(
12
)
(
)
(
2
)
(
3
1
0
0
1
n
f
h
x
f
x
f
x
x
(16)
Ostatecznie otrzymujemy
17
b
a
,
''
f
12
h
b
f
a
f
2
h
dx
x
f
3
b
a
powyższe równanie jest znane jako wzór trapezów – widać, że wzór pozwala dokładnie
obliczyć całkę dla dowolnej funkcji, dla której druga pochodna jest zero (wtedy reszta jest 0).
Zauważmy, że dla powyższych rozwiązań końce przedziału całkowania są węzłami
kwadratury x
0
=a oraz x
n
=b. Wzory kwadratur oparte na tych węzłach nazywamy wzorami
zamkniętymi Newtona – Cotesa.
Podobnie wyznaczyć można wzory dla innych n.
Kwadratury Newtona-Cotesa – wzory zamknięte
Poniżej znajdują się wzory zamknięte dla n = (1..3):
wzór trapezów
18
,
''
12
2
1
1
0
3
1
0
x
x
f
h
x
f
x
f
h
dx
x
f
n
b
a
wzór parabol (Simpsona)
19
,
90
4
3
2
2
0
)
3
(
5
2
1
0
x
x
f
h
x
f
x
f
x
f
h
dx
x
f
n
b
a
wzór Bessela
20
,
80
3
3
3
8
3
3
3
0
)
4
(
7
3
2
1
0
x
x
f
h
x
f
x
f
x
f
x
f
h
dx
x
f
n
b
a
50
Interpretacja graficzna
wzór parabol (Simpsona)
wzór Bessela
Kwadratury Newtona - Cotesa - wzory otwarte
Wzory kwadratur nie opartych na węzłach będących końcami przedziału całkowania
nazywamy wzorami otwartymi Newtona - Cotesa
51
Jeżeli dla skończonego [a ; b] przedziału wybierzemy zbiór punktów węzłowych {x
0
, …, x
n
}
takich, że :
h
i
x
x
i
0
gdzie:
2
)
(
n
a
b
h
dla i=(0, …, n)
oraz
h
a
x
0
W oparciu o Twierdzenie 1 można również wyznaczyć wzory kwadratur otwartych.
Wówczas
wzory kwadratur dla n parzystych mają postać
21
]
b
,
a
[
gdzie
,
dt
n
t
...
1
t
t
!
2
n
f
h
dx
x
f
a
dx
x
f
n
0
i
1
n
1
2
2
n
3
n
i
i
b
a
oraz dla n nieparzystych mają postać:
b
a
n
i
n
n
n
n
i
i
b
a
gdzie
dt
n
t
t
t
n
f
h
x
f
a
dx
x
f
0
1
1
)
1
(
)
1
(
)
2
(
)
22
(
]
;
[
,
)
)...(
1
(
)!
1
(
)
(
)
(
)
(
wzór prostokątów
23
x
x
,
''
f
3
h
x
f
2
h
dx
x
f
0
n
1
1
3
0
b
a
24
x
x
,
''
f
3
h
3
x
f
x
f
2
h
3
dx
x
f
1
n
2
1
3
1
0
b
a
25
,
45
14
2
2
3
4
2
3
1
)
4
(
5
2
1
0
x
x
f
h
x
f
x
f
x
f
h
dx
x
f
n
b
a
26
,
144
95
11
11
24
5
3
4
1
)
4
(
5
3
2
1
0
x
x
f
h
x
f
x
f
x
f
x
f
h
dx
x
f
n
b
a
52
Interpretacja geometryczna
wzór prostokątów
wzór trapezów
Kwadratury Newtona-Cotesa – porównanie
Przykładowe wyniki obliczeń całki
29289922
,
0
2
2
1
x
cos
dx
x
sin
4
0
4
0
Wyniki obliczeń dla wzorów otwartych i zamkniętych:
n
0
1
2
3
wzory zamknięte
0,27768018
0,29293264
0,29291070
błąd
0,01521303
0,00003924
0,00001748
wzory otwarte
0,30055887
0,29798754
0,29285866
0,29286923
błąd
0,007666565
0,00509432
0,00003456
0,00002399
53
Kwadratury Newtona-Cotesa – podsumowanie
Kwadratury Newtona-Cotesa dają coraz lepszą dokładność wraz ze wzrostem n. Jednak wraz
ze wzrostem n wzór na sumę posiada również rosnącą liczbę czynników. Dlatego kwadratur
Newtona-Cotesa nie stosuje się dla dużych n oraz z uwagi na pojawianie się ujemnych
współczynników (dla n = 7). Co więcej dla dużych n istnieją funkcje ciągłe, dla których ciąg
kwadratur Newtona – Cotesa nie jest zbieżny do całki
b
a
dx
x
f
)
(
W praktyce nie stosuje się kwadratur Newtona – Cotesa wysokiego rzędu.
Jak widać aby obliczyć wartość całki stosując kwadratury Newtona – Cotesa wystarczy
wyznaczyć wartość funkcji w punktach węzłowych, a następnie podstawić je do wzoru.
Większą trudność sprawia oszacowanie reszty kwadratury, gdyż należy wyznaczyć
maksimum pochodnej k rzędu funkcji dla przedziału całkowania.
Obliczyć poniższą całkę stosując metodę Simpsona
4
4
5
0
d
)
sin(
d
2
90
1
sin
2
sin
4
0
sin
6
dx
)
x
sin(
cos
d
)
sin(
d
4
4
, wtedy
1
)
cos(
max
,
0
wówczas wyniki:
09439510
,
2
3
2
dx
)
x
sin(
0
oraz oszacowanie błędu obliczeń:
106256835
,
0
2880
1
2
90
1
5
2
rozwiązanie dokładne:
2
)
0
cos(
)
cos(
)
x
cos(
dx
)
x
sin(
0
0
Kody źródłowe:
przykładowa całkowana funkcja:
function f(x:extender):extender,
begin
result:=sin (2*x*x)+2*cos(3x);
end;
54
wzór prostokątów:
function calka_n_0_wzor_otwarty(a,b:extended):extended;
var h:extended;
begin
h:=(b-a)/2;
result:=2*h*f(a+h);
end;
wzór trapezów:
function calka_n_1_wzor_zamkniety(a,b:extended):extended;
var h:extended;
begin
h:=(b-a);
result:=(h/2)*(f(a)+f(b));
end;
wzór parabol:
function calka_n_2_wzor_zamkniety(a,b:extended):extended;
var h:extended;
begin
h:=(b-a)/2;
result:=(h/3)*(f(a)+f(a+h)=f(b));
end;
55
Kwadratury z ustalonymi węzłami
Złożone kwadratury Newtona – Cotesa:
W poprzednim punkcie zostały przedstawione kwadratury Newtona-Cotesa. Wiemy, że błąd
kwadratury zależy od długości przedziału całkowania oraz rzędu kwadratury. Dlatego w
oparciu właśnie o kwadratury Newtona-Cotesa niskich rzędów zdefiniujemy tzw. wzory
złożone Newtona-Cotesa. W celu zwiększenia dokładności obliczeń należy zmniejszyć
długość przedziału. Wynika z tego, że dzieląc przedział całkowania na m podprzedziałów i
obliczając dla każdego przedziału kwadraturę Newtona-Cotesa niskiego rzędu, następnie, gdy
zsumujemy wyniki podprzedziałów, uzyskamy wynik zbieżny do rozwiązania dokładnego,
natychmiast zwiększy się dokładność obliczeń.
m=1
Otrzymane w ten sposób kwadratury określono na całym przedziale całkowania nazywamy
złożonymi kwadraturami Newtona-Cotesa
Twierdzenie 2:
Podzielmy przedział całkowania [a,b] na m części przy pomocy punktów:
h
i
a
x
i
, dla
)
1
m
,...,
1
,
0
(
i
,
m
)
a
b
(
h
mamy wtedy:
1
m
0
i
x
x
b
a
1
i
i
dx
)
x
(
f
dx
)
x
(
f
56
stosując do każdej z całek po prawej stronie równania dowolny wzór kwadratury Newtona-
Cotesa
...
Na podstawie poprzedniej zależności możemy wyprowadzić wzory złożonych kwadratur:
złożony wzór prostokątów:
E
2
h
x
f
h
E
2
h
x
f
h
dx
)
x
(
f
1
m
0
i
i
1
m
0
i
i
1
m
0
i
i
b
a
(28)
gdzie błąd przybliżenia
1
m
0
i
)
2
(
3
f
3
h
E
(29)
oraz
1
i
i
i
x
,
x
złożony wzór trapezów:
E
x
f
2
)
x
(
f
)
x
(
f
h
E
)
x
(
f
)
x
(
f
2
h
dx
)
x
(
f
1
m
1
i
i
m
0
1
m
0
i
i
1
m
0
i
1
i
i
b
a
(30)
gdzie błąd przybliżenia
i
1
m
0
i
)
2
(
3
f
12
h
E
(31)
oraz
1
i
i
i
x
,
x
57
złożony wzór parabol – dla m parzystych:
1
m
0
i
i
1
m
0
i
1
i
1
i
i
b
a
E
)
x
(
f
)
2
h
x
(
f
4
)
x
(
f
3
h
dx
)
x
(
f
(32)
E
)
x
(
f
4
)
x
(
f
2
)
x
(
f
)
x
(
f
3
h
dx
)
x
(
f
1
2
m
1
i
1
i
2
1
2
m
1
i
i
2
m
0
b
a
(33)
gdzie błąd przybliżenia
)
(
f
h
180
)
a
b
(
E
)
4
(
4
(34)
oraz
1
i
i
i
x
,
x
Podsumowanie:
Z powyższych wzorów wynika, że jeśli funkcja f jest wystarczająco regularna, to całka tej
funkcji może być z dowolnie dużą dokładnością przybliżona przy pomocy złożonych
kwadratur. Dla złożonych kwadratur Newtona-Cotesa, również największą trudność sprawia
oszacowanie błędu obliczeń.
Twierdzenie 3:
Jeżeli
]
b
;
a
[
c
f
(jest jednokrotnie całkowalna w przedziale [a,b]), to dla każdego ciągu m
złożonych kwadratur Newtona-Cotesa n-tego rzędu dla funkcji f jest zbieżny do całki:
b
a
dx
)
x
(
f
, przy
m
inaczej
1
m
1
i
x
x
m
b
a
1
i
i
dx
)
x
(
f
lim
dx
)
x
(
f
h
i
a
x
i
,
)
1
m
,...,
1
,
0
(
i
,
m
)
a
b
(
h
58
Przykład:
Obliczyć poniższą całkę:
Rozwiązanie dokładne
59815
,
53
0
4
4
0
e
e
dx
e
x
Metoda Simpsona m=1, h=2
76958
,
56
4
3
2
4
2
0
4
0
e
e
e
dx
e
x
błąd = -3,17143
Metoda Simpsona m=1, h=1
86385
,
53
4
3
1
4
3
1
4
3
2
2
1
0
4
2
2
0
4
0
e
e
e
e
e
e
dx
e
dx
e
dx
e
x
x
x
błąd = -0,26570
Metoda Simpsona m=4, h=0,5
61622
,
53
4
3
3
2
2
1
1
0
4
0
dx
e
dx
e
dx
e
dx
e
dx
e
x
x
x
x
x
błąd = -0,01807
Przykładowe wyniki obliczeń całki:
2
)
cos(
)
sin(
0
0
x
dx
x
Wyniki obliczeń dla wzorów złożonych w zależności od m:
m
2
4
8
16
32
64
wzór
prostokątów
wynik
1,110720735 1,751785441 1,936297295 1,983970758 1,995986709 1,998996147
błąd
0,889279265 0,248214559 0,066102705 0,016029242 0,004013791 0,001003853
wzór
trapezów
wynik
1,570796327 1,896118898 1,1974231602 1,993570344 1,998393361 1,999598389
błąd
0,429203673 0,103881102 0,025768398 0,006429656 0,001606639 0,000401611
wzór parabol
wynik
2,094395102 2,004559755 2,000269170 2,000016591 2,000001033 2,000000065
błąd
-0,094395102 -0,004559755 -0,000269170 -0,000016591 0,000001033 0,000000065
Kody źródłowe:
Złożony wzór prostokątów:
function
Calka_Prostokat_Zlozony(a,b:extended;m:integer):extended;
var h,x,I_0:extended;
begin
h:=(b-a)/m;
59
I_0:=0.0;
for i:=0 to m-1 do
begin
x:=a+(i*h);
x:=x+(h/2);
I_0:=I_0+f(x)
end;
result:=I_0*h
end;
Złożony wzór trapezów:
function
Calka_Trapez_Zlozony(a,b:extended;m:integer):extended;
var h,x,I_0,I_1:extended;
begin
h:=(b-a)/m;
I_0:=f(a)+f(b);
I_1:=0.0;
for i:=0 to m-1 do
begin
x:=a+(i*h);
I_1:=I_1+f(x);
end;
result:=(0.5*I_0+I_1)*h;
end;
Złożony wzór parabol:
function
Calka_Simpson_Zlozony(a,b:extended;m:integer):extended;
var h,x,I_0,I_1,I_2:extended;
begin
h:=(b-a)/m;
I_0:=f(a)+f(b);
I_1:=0.0;
I_2:=0.0;
for i:=0 to m-1 do
begin
x:=a+(i*h);
if (i mod 2)=0
then I_2:=I_2+f(x);
else I_1:=I_1+f(x);
end;
result:=(I_0+2.0*I_2+4.0*I_1)*h/3;
end;
60
Kwadratury z ustalonymi węzłami
Metoda ekstrapolacyjna – Romberga
Metoda adaptacyjna
Kwadratury Gaussa
W poprzednim rozdziale omawiane były metody Newtona-Cotesa, które bazowały na
całkowaniu wielomianu interpolacyjnego stopnia n. Pozwalał on na dokładne przybliżenie
całki dla funkcji stopnia (n+1) lub (n+2). Wszystkie metody Newtona-Cotesa
wykorzystywały wartości funkcji dla punktów węzłów równoległych. Ten warunek jest
dogodny dla konstrukcji złożonych wzorów, ale może znacząco zmniejszyć dokładność
przybliżenia całek.
Poniżej widać kwadraturę Newtona-Cotesa opartą na końcach przedziału – przybliżenie
funkcji nie jest dokładne.
Kwadratury Gaussa dobierają takie węzły, aby osiągnąć optymalne przybliżenie całki:
n
1
i
i
i
b
a
)
x
(
f
c
dx
)
x
(
f
Węzły kwadratury x
1
, x
2
,..., x
n
z przedziału całkowania [a,b] oraz współczynniki c
1
, c
2
,...,c
n
są
tak dobrane, aby zminimalizować błąd przybliżenia. Nie zakładamy jednak żadnych
ograniczeń na węzły x
i
o współczynnikach c
i
natomiast chcemy zmaksymalizować rząd
kwadratury.
61
Mamy znaleźć 2n stałe we wzorze, przypuszczamy więc, że rząd kwadratury powinien
wynosić 2n, a więc wzór powinien być dokładny dla wielomianów stopnia (2n-1). Dla
wielomianu stopnia równego lub wyższego od (2n-1) należy uwzględnić resztę kwadratury E:
E
)
x
(
f
c
dx
)
x
(
f
n
1
i
i
i
b
a
Aby pokazać proces wyboru odpowiednich węzłów oraz współczynników zilustrujemy to na
przykładzie:
n=2 – liczba węzłów oraz przedział [-1,1]
Należy wyznaczyć x
1
, x
2
oraz c
1
, c
2
, ponieważ:
)
x
(
f
c
)
x
(
f
c
dx
)
x
(
f
2
2
1
1
1
1
(64)
Zgodnie z założeniami, powyższe przybliżenie będzie dokładne tylko gdy stopień funkcji f(x)
mniejszy lub równy
3
1
)
2
(
2
1
2
n
wówczas
3
3
2
2
1
0
x
a
x
a
x
a
a
)
x
(
f
(65)
gdzie a
0
, a
1
, a
2
, a
3
– dowolne stałe
a więc możemy zapisać:
dx
x
a
dx
x
a
xdx
a
dx
1
a
dx
x
a
x
a
x
a
a
3
3
2
2
1
0
3
3
2
2
1
0
Poszukujemy x
1
, x
2
oraz c
1
, c
2
, więc można zapisać:
2
1
c
1
c
2
1
(67)
0
dx
x
x
c
x
c
1
1
2
2
1
1
(68)
3
2
dx
x
x
c
x
c
1
1
2
2
2
2
2
1
1
(69)
0
dx
x
c
x
c
1
1
3
2
2
3
1
1
(70)
otrzymaliśmy układ liniowych równań algebraicznych, którego rozwiązanie jest:
1
c
1
1
c
2
3
3
x
1
3
3
x
2
62
ostatecznie otrzymaliśmy wzór na przybliżenie całki dla funkcji danej wzorem (64):
3
3
f
3
3
f
dx
)
x
(
f
1
1
Metoda poszukiwania węzłów oraz współczynników może zostać zastosowana dla funkcji
wyższych rzędów, lecz istnieje łatwiejsza metoda. Bazuje ona na wykorzystaniu jednego z
wielomianów ortogonalnych. Mają one taką właściwość, że całka z iloczynu dwóch
dowolnych wielomianów jest równa 0:
b
a
dx
x
g
x
f
g
f
0
)
(
)
(
,
Dla naszych dalszych rozważań, zajmiemy się tylko ciągiem wielomianów Legendre’a
)
(
)
(
),
(
1
0
x
P
x
P
x
P
n
, który posiada właściwości:
1. Dla każdego n, P
n
(x) jest wielomianem stopnia n
2.
1
1
0
)
(
)
(
dx
x
P
x
P
n
, gdzie P(x) jest wielomianem stopnia wyższego
(niższego???)
od n
Wzory kilku pierwszych wielomianów Legendre’a:
1
)
(
0
x
P
x
x
P
)
(
1
3
1
2
3
)
(
2
2
x
x
P
x
x
x
P
3
5
2
5
)
(
3
3
oraz wzór rekurencyjny do wyznaczania wielomianów krzywych
???
oraz ich pochodnych:
)
(
1
)
(
1
2
)
(
2
1
x
P
n
n
x
P
x
n
n
x
P
n
n
n
)
(
1
)
(
1
)
(
'
1
2
2
x
P
x
n
x
P
x
x
n
x
P
n
n
n
Pierwiastki tego wielomianu są tylko pojedyncze, leżą w przedziale <-1,1> oraz są
symetryczne względem początku układu współrzędnych. Najważniejszą właściwością jest to,
że bardzo dobrze nadają się do wyznaczania punktów węzłowych oraz współczynników.
Punkty węzłowe potrzebne w celu dokładnego przybliżenia całki funkcji stopnia mniejszego
od 2n są pierwiastkami wielomianu Legendre’a
n-tego stopnia.
Jeżeli x
1
, x
2
, ... x
n
są pierwiastkami wielomianu Lagrange’a, P
n
(x) stopnia n dla każdego i =
(1,2,...n) wtedy współczynniki c
i
można obliczyć ze wzoru:
2
1
2
1
1
1
1
1
)
(
)
1
(
1
2
)
(
'
)
(
)
1
(
2
i
n
i
i
n
i
n
n
i
k
k
k
i
k
x
P
n
x
x
P
x
P
n
dx
x
x
x
x
c
wówczas jeżeli P(x) jest dowolnym wielomianem stopnia mniejszego od 2n wtedy:
1
1
1
)
(
)
(
n
i
i
i
x
P
c
dx
x
P
(75)
63
reszta dla stopnia większego od (2n-1):
)
2
2
(
1
1
2
!
1
2
)
(
n
f
n
dx
x
P
E
Pierwiastki wielomianu oraz współczynniki c
i
są tabelaryzowane – nie obliczamy ich przy
każdym całkowaniu.
Prawa strona wzoru (75) jest nazwana kwadraturą Gaussa-Lagrange’a.
Wzór ten pozwala na obliczenie całki tylko w przedziale
1
,
1
. W celu obliczenia
przybliżenia całki w dowolnym przedziale
b
a,
dla którego f(x) jest ciągła należy wykonać
przekształcenie poprzez podstawienie.
Dokonujemy przekształcenia dowolnego przedziału
b
a,
do przedziału
1
,
1
.
2
2
a
b
t
a
b
x
dt
a
b
dx
2
Wówczas możemy zapisać:
b
a
dt
a
b
t
a
b
f
a
b
dx
x
f
1
1
2
2
2
)
(
po zastosowaniu kwadratury Gaussa-Lagrange’a otrzymujemy
2
2
2
2
2
2
)
(
1
1
1
a
b
x
a
b
f
c
a
b
dt
a
b
t
a
b
f
a
b
dx
x
f
i
n
i
i
b
a
gdzie x
i
– pierwiastek wielomianu stopnia n
c
i
– współczynnik ze wzoru (74)
przykład:
Obliczyć całkę wykorzystując kwadraturę Gaussa:
5
,
1
1
1093643
,
0
2
dx
e
x
dokończyć ???
najpierw przekształcenia przedziału całkowania:
dt
t
f
dt
t
f
dx
e
x
1
1
5
,
1
1
1
1
2
16
5
4
1
2
1
5
,
1
2
1
5
,
1
2
1
5
,
1
2
po zastosowaniu kwadratury Gaussa-Lagrange’a, n=2 otrzymujemy:
c
1
=1 c
2
=1
3
3
1
x
3
3
2
x
64
1094003
,
0
1
1
4
1
4
1
16
5
3
3
16
5
3
3
16
5
2
16
5
1
5
,
1
1
2
1
2
e
e
e
c
e
c
dx
e
x
x
x
Po zastosowaniu kwadratury Gaussa-Lagrange’a, n=3 otrzymujemy:
9
5
1
c
9
8
2
c
9
5
3
c
3
15
1
x
0
2
x
3
15
3
x
1093642
,
0
9
5
9
8
9
5
4
1
4
1
16
5
3
15
16
5
0
16
5
3
15
16
5
3
16
5
2
16
5
1
5
,
1
1
3
2
1
2
e
e
e
e
c
e
c
e
c
dx
e
x
x
x
x
Tablica pierwiastków wielomianu Legendre’a:
n
x
0
x
1
x
2
x
3
x
4
1
0
-
-
-
-
2
-0,5773509
0,5773509
-
-
-
3
0,7745966
0
-0,77459666
-
-
4
0,86113631
0,33998104
-0,33998104
-0,86113631
-
5
-0,53846931
-0,90617984
0
0,90617984
0,53846931
Kod źródłowy kwadratury Gaussa:
Przykład f. całką wykorzystując kwadraturę G w przedziale
b
a,
oraz n- oznacz. stopnia
met. L.
function Gauss_Legendre (n:Integer; var x:vector); Extended
var k,l,m:Integer;
p0,p1,p2,s,x1,ci:Extended
begin
if n>0 then
begin
s:=0
l:=n div2;
for k:=0 to n do
begin
p0:=1;
x1:=x[k];
65
p1:=x1
for m:=1 to n-1 do
begin
p2:=((2*m+1)*x1*p1-m*p0)/(m+1);
p0:=p1;
p1:=p2;
end;
ci:=-2*(x1*x1-1)/((n+1)*(n+1)*p1*p1);
s:=s+f(x1)*ci;
end;
Result:=s;
end;
end;
gdzie x – tabl. zawiera pierwiastki wielomianu Legendre’a
oraz n – stopień wielomianu Legendre’a
Obliczanie całek niewłaściwych:
Przy numerycznym całkowaniu w przedziale obustronnie lub jednostronnie nieskończonym,
należy stosować takie funkcje wagowe, które zapewniają zbieżność całki, W zależności od
przedziału stosujemy odpowiedni wielomian interpolacyjny:
n
i
i
i
b
a
x
P
c
dx
x
f
x
w
I
1
)
(
)
(
)
(
Całkowanie funkcji osobliwych:
W celu obliczenia przybliżonej całki funkcji posiadającej punkty nieciągłości, należy
zmodyfikować poznanie metody. Aby dokonać takiego całkowania, należy zlokalizować
wszystkie punkty nieciągłości należące do przedziału całkowania.
Przykłady funkcji nieciągłych:
b
a
dx
x
f
I
)
(
,
gdzie f(x):
przypadek 1*
przypadek 2*
0
66
b
a
c
a
b
a
dx
x
f
c
dx
x
f
dx
x
f
)
(
)
(
)
(
Powyższe metody wymagają wielu przekształceń matematycznych, które nie dają wyznaczyć
się informatycznie.
Ma to charakter historyczny.
Alternatywną metodą wyliczenia całki jest metoda polegająca na oddalaniu się od punktu
nieciągłości
o małe ε ???
Metoda zwykłego szeregu Taylora:
9253141
,
2
0017691
,
0
9235450
,
2
1
0
dx
x
e
x
Metoda adaptacyjna:
925324091
,
2
1
0000000001
,
0
1
0
dx
x
e
dx
x
e
x
x
zał. dokładne 0,0001 liczba przedzi. 235 il. wywoł. funkcji 941
???
Całki wielokrotne:
rozpatrzmy całkę podwójną:
b
a
x
d
x
c
dx
dy
y
x
f
)
(
)
(
)
,
(
x – stała i stosunek kwadratury do obliczonej całki
)
(
)
(
)
,
(
)
(
x
d
x
c
dy
y
x
f
x
s
67
otrzymujemy:
1
1
1
1
2
2
2
0
)
,
(
4
)
,
(
2
)
,
(
)
,
(
6
)
(
)
(
)
(
m
j
m
j
j
j
m
y
x
f
y
x
f
y
x
f
y
x
f
m
x
c
x
d
x
s
m – liczba przedziałów przedziału [c(x) ; d(x)]
wykorzystanie z kolei kwadratury Simpsona do obliczenia całki:
b
a
dx
x
S
s
)
(
ostatecznie:
1
1
1
1
2
2
2
0
)
(
4
)
(
2
)
(
)
(
6
n
i
n
i
i
i
n
x
S
x
S
x
S
x
S
n
a
b
s
np.:
5
,
1
0
,
1
0
,
2
4
,
1
4295549269
,
0
)
2
ln(
dx
dy
y
x
stos. złoż. kwadratury Simpsona:
dla n=4 i m=2
5
,
1
0
,
1
0
,
2
4
,
1
7
4295522438
,
0
)
2
ln(
dx
dy
y
x
1800
6
,
0
5
,
0
E
6
10
72
,
4
E
Całkowanie wielokrotne – całka podwójna
kod źródłowy:
function f (x,y:real):real;
begin
Result:=
explg (x)???
end;
function c(x:real):real;
begin
Result:=x*x*x
end;
function d(x:real):real;
begin
Result:=x*x
end;
68
function simpson-Podwojna
(a,b:real; m,n:integer):real;
var i,j:integer
h,hx,j1,j2,j3
k1,k2,k3,x,y,x:real;
begin
h:=(b-a)/(2*n)
j1:=0; j2:=0;
j3:=0
for i:=0 to 2*n do
begin
x:=a+i*h;
nx:=(d(x)-c(x))/(2*m);
k1:=f(x,c(x))+f(x,d(x));
k2:=0
k3:=0
for j:=1 to 2*m-1 do
begin
y:=c(x)+j*hx
z:=f(x,y);
if Odd(j) ???
then k3 :=k3+z
else k2:=k2+z
end
i:=(k1+2*k2+4*k3)+hx/3
if(i=0) or (i=2*n)
then j1:=j1+1
else if Odd (i)
then j3:=j3+l
else j2:=j2+l
end;
Result:=(j1+2*j2+4*h/3)
end;
Kody źródłowe:
Function f(x,y:real):real;
Begin
Result:=exp(y/x);
End;
Function c (x:real):real;
Begin
Result:=x*x*x;
End;
Function d (x:real):real;
Begin
Result:=x*x;
69
End;
function Simpson_Podwojna (a,b:real; m,n:integer):real;
var i,j:integer;
h,hx,j1,j2,j3,k1,k2,k3,I,x,y,z:real;
begin
h:=(b-a)/(2*n);
j1:=0;
j2:=0;
j3:=0;
for i:=0 to 2*n do
begin
x:=a+i*h;
tu tomek napisał ... j*h a ja mam i*h sprawdźcie!!!
hx:=(d(x)-c(x))/(2*m)
k1:=f(x,c(x)+f(x,d(x));
k2:=0;
k3:=0;
for j:=1 to 2*m-1 do
begin
y:=c(x) +j*hx;
z:=f(x,y);
if odd(j)
then k3:=k3+z
else k2:=k2+z
end;
I:=(k1+2*k2+4*k3)*hx/3;
If (i=0) or (i=2*n)
then j1:=j1+I
else if Odd(i)
tu też mam troche inaczej
then j3:=j3+l
else j2:=j2+l
end;
Result:=(j1+2*j2+4*j3)*h/3;
End;
Obliczanie całek metodą Monte-Carlo
Metoda Monte Carlo – całka pojedyncza
Mamy obliczyć przybliżoną wartość In całki oznaczonej
a
b
f(x)dx
I
(129)
przy założenie, że f(x) jest funkcją ciągłą w przedziale domkniętym [a;b]
Następnie n-krotnie generujemy realizację x zmiennej losowej X o rozkładzie równomiernym
w przedziale [a,b], otrzymujemy w rezultacie ciąg liczb x
1
, x
2
, ..., x
n
,
70
Przybliżoną wartość całki określa wzór
n
i
i
n
x
f
n
a
b
I
1
)
(
)
(
(130)
Zbieżność przybliżonej wartości całki In do jej dokładnej wartości jest gwarantowana w
sensie probalistycznym, tzn., że dla każdego ε>0 zachodzi relacja:
1
)
1
n
(
P
lim
n
(131)
gdzie P jest prawdopodobieństwem, przy czym błąd metody monte Carlo można określić
wzorem
n
I
I
n
1
(132)
Maszyny typowe zazwyczaj dysponują generatorem liczb losowych o rozkładzie
równomiernym w przedziale [0,1], zachodzi zatem konieczność przeliczania wylosowanej
realizacji Y z przedziału [0,1], na realizację zmiennej losowej X o rozkładzie równomiernym
w przedziale [a, b]. Wzory poniżej, to ilustrują:
Przeliczanie zmiennych losowych Y na X
X=(b-a)Y+a
(133)
Przeliczanie realizacji zmiennych losowych Y na X
x=(b-a)y+a
(134)
Kod źródłowy:
Function Calka_MonteCarlo (a,b:extended; n:lingint):extended
Var i:integer;
h,x:extended;
begin
Result:=0;
h:=(b-a)/n;
for :=1 to n do
begin
x:=a+((b-a)*Random);
Result:=Result+f(x);
End;
Result:=Result*h;
end;
71
Całkowanie – Metody Monte Carlo. Całki wielokrotne właściwe
Dana jest całka wielokrotna:
3
2
1
2
1
...
)
,...,
,
(
dx
dx
dx
x
x
x
f
I
n
(135)
Każda ze zmiennych x
i
zawiera się w przedziale
i
i
i
b
x
a
(136)
Losujemy n-krotnie punkty Pk kostki, przyjmując że współrzędnymi tych punktów są
zmienne losowe X
1
(k)
, X
2
(k)
,..., X
m
(k)
o rozkładzie równomiernym w odpowiednich
przedziałach [a
i
,b
i
]. Każdemu punktowi Pk odpowiada zatem ciąg m-elementowy
wylosowanych realizacji x
1
(k)
, x
2
(k)
,..., x
m
(k)
. Po wylosowaniu sprawdzamy czy każdy punkt Pk
należy do obszaru Ω. Możemy zatem po zakończeniu losowań uporządkować wylosowane
punkty w następujący sposób:
k
P
dla k=(1, 2, ..., n),
k
P
dla k=(n+1, n+2, ..., N)
(137)
Przy dostatecznej liczbie losowań możemy obliczyć przybliżoną wartość V
Ω
objętości m-
wymiarowanego obszaru całkowania:
V
N
n
V
(138),
gdzie
)
(
1
i
m
i
i
a
b
V
(139)
jest objętością m-wymiarowej kostki. Przybliżona wartość I
n
całki (135) określana wzorem:
n
k
k
n
P
f
N
V
I
1
)
(
(140)
Aby uzyskać ciąg realizacji X1, X2, ..., Xm losujemy m-krotnie zmienną losową Y o
rozkładzie równomiernym w przedziale [0,1], a następnie przeliczamy realizacje y na
realizację x
i
(k)
w przedziale [a
i
; b
i
] korzystając ze zależności:
x
i
(k)
=(b
i
-a
i
)y
i
(k)
+a
i
Rozwiązywanie równań różniczkowych
Rozważania dotyczą:
- równań różniczkowych zwyczajnych z warunkami początkowymi i brzegowymi,
- układów równań różniczkowych zwyczajnych z warunkami początkowymi i
brzegowymi,
- równań różniczkowych cząstkowych,
- układów równań różniczkowych cząstkowych
72
Numeryczne obliczanie pochodnej
Aby obliczyć pochodną funkcji f w punkcie x
0
, mamy
h
x
f
h
x
f
x
f
h
)
(
)
(
lim
)
(
'
0
0
0
(1)
W celu wyznaczenia przybliżonej wartości pochodnej f’(x) dla małych wartości h możemy
użyć wyrażenia
h
x
f
h
x
f
)
(
)
(
0
0
(2)
Możemy to udowodnić. Dokonujemy interpolacji f(x); zdefiniujemy wielomian Lagrange`a
stopnia 1, który będzie oparty na węzłach x
0
oraz x
0
+h pod warunkiem, że funkcja f jest
dwukrotnie różniczkowalna w przedziale [a; b]
!
2
)
)(
(
)
(
)
(
0
0
1
h
x
x
x
x
x
P
x
f
ktoś nie zdązył do końca, a ja nie wiem czy w
mianowniku jest 2!, czy co??
Wówczas obliczając pochodną funkcji f, mamy:
h
)
x
(
f
)
h
x
(
f
)
x
(
'
f
0
0
0
(3)
oraz po uproszczeniu
)
(
"
2
f
h
(4)
73
Jeżeli dokonamy interpolacji funkcji f wielomianem stopnia 2, wówczas
h
h
x
f
h
x
f
x
f
2
)
(
)
(
)
(
'
0
0
0
(5)
oraz po uproszczeniu
)
(
6
)
3
(
2
f
h
I(6)
Za pomocą wielomianów wyższych stopni możemy wyznaczyć wzory na przybliżenie
pochodnej.
Przykład: Obliczyć pochodną funkcji f(x)=ln(x) dla x
0
=1,8 stosując metodę w przód i wstecz
stopnia 1.
2
)
(
)
(
)
(
'
0
0
0
h
h
x
f
h
x
f
x
f
M
h
h
f
h
f
f
2
)
8
,
1
(
)
8
,
1
(
)
(
"
,
gdzie:
)
(
"
max
)
8
,
1
,
8
,
1
(
f
M
h
brak textu
Metoda w przód
Metoda wstecz
h
Wynik
Błąd
Wynik
Błąd
0,1
0,5406772
0,0154321
0,5715841
0,0173010
0,01
0,5540180
0,0015432
0,5571045
0,0015605
0,001
0,5554013
0,0001543
0,5557099
0,0001545
74
Numeryczne obliczanie drugiej pochodnej.
Jeżeli dokonamy rozwinięcia funkcji f w szereg Taylora 3 rzędu w otoczeniu punktu x
0
dla
punktów x
0
+h oraz x
0
-h, zakładając, że f jest 4-krotnie różniczkowalna na przedziale [x
0
-h,
x
0
+h], wówczas dla x
0
-h:
Wzor (8)
gdzie
]
x
h;
[x
ξ
0
0
dodając oba równania stronami otrzymujemy:
wzór (9)
Rozwiązując równanie ze względu na f”(x
0
), otrzymujemy:
)]
(
)
(
2
)
(
[
1
)
(
"
0
0
0
2
0
h
x
f
x
f
h
x
f
h
x
f
(11)
oraz
)
(
[
12
)
(
n
n
f
h
, gdzie
]
;
[
0
0
h
x
h
x
Powyższy wzór otrzymamy różniczkując obustronnie wzór na pierwszą pochodną (5) kroku
h/2.
Przykład: Obliczyć drugą pochodną funkcji f(x)=ln(x
3
) dla x
0
=2 stosując powyższą metodę
różnych wartości kroku h.
Wzór
Rozwiązanie dokładne
2
3
)
(
"
x
x
f
f”(2)=-0,1875
h
Wynik
Błąd
Oszacowanie ε
0,1
-0,1875586182
0,0000586182
0,01
-0,1875005859
0,001
-0,1875000049
Równania różniczkowe zwyczajne z warunkami początkowymi
Równania różniczkowe są wykorzystywane do odwzorowania wielu problemów z zakresu
różnych nauk, które dotyczą zmian jednej zmiennej względem innej zmiennej. Wiele
problemów świata rzeczywistego opisywanych przez równania różniczkowe jest za bardzo
złożonych, aby uzyskać rozwiązanie dokładne (analityczne), należy znaleźć rozwiązanie
przybliżone. Są dwa sposoby aby uzyskać rozwiązanie przybliżone. Po pierwsze można
dokonać uproszczenia zadanego równania różniczkowego do takiej postaci, dla której znamy
rozwiązanie dokładne, wówczas uproszczone równanie będzie przybliżało rozwiązanie
równania zadanego. Drugim sposobem jest zastosowanie jednej z metod numerycznych w
celu przybliżenia rozwiązania zadanego równania.
Metody różnicowe jednokrokowe:
- metody Eulera;
- metody Rundego-Kutty
- metody Rundego-Kutty-Fehenberga
75
Metody różnicowe wielokrokowe:
- metody bezpośrednie
- metody pośrednie
- metoda Predictor-Corrector
- metoda trapezów
Metoda Eulera
Rozpatrzmy zadanie znalezienia funkcji y=y(t), które dla
]
,
[
b
a
t
spełnia równanie
różniczkowe zwyczajne pierwszego rzędu:
)
,
( y
t
f
dx
dy
(7)
oraz warunek początkowy
y(a)=y
0
(8),
gdzie f jest daną funkcją dwu zmiennych. Warunek wystarczający istnienia i jednoznaczności
rozwiązania zadania (7) oraz (8) jest warunek ciągłości funkcji f(t,y)
Metoda ta rzadko jest używana w praktyce, natomiast ze względu na jej prostotę, zostanie
przedstawiona w celu zobrazowania i zrozumienia technik stosowanych w bardziej
zaawansowanych metodach, pomijając przy tym złożone działania matematyczne.
Przedmiotem rozważań będzie poniższe równanie.
)
,
( y
t
f
dx
dy
, gdzie
]
,
[
b
a
t
oraz y(a)=a (10)
Na rozwiązanie powyższego zagadnienia będziemy obliczać przybliżone wartości funkcji
yi=y(ti), gdzie ti=a+ih, h=(b-a)/N, dla którego i=(0,1, ..., N), gdzie ti nazywane są punktami
węzłowymi, natomiast h odległością między nimi.
Rozwiniemy y(t) w szereg Taylora w celu wyprowadzenia metody Eulera. Zakładając, że y(t)
jest jedynym rozwiązaniem (10) oraz posiada drugą pochodną, która jest ciągła w przedziale
[a, b] wówczas dla każdego i=1, 2, ..., N.
Wykorzystując powyższe założenia możemy zapisać:
)
(
"
2
)
(
)
(
'
)
(
)
(
)
(
2
1
1
1
i
i
i
i
i
i
i
i
y
t
t
t
y
t
t
t
y
t
y
(11), gdzie
)
,
(
1
i
i
i
t
t
Oznaczamy h=t
i+1
-t
i
, wówczas otrzymujemy:
)
(
"
2
)
(
'
)
(
)
(
2
1
i
i
i
i
y
h
t
hy
t
y
t
y
(12)
Użyjemy podstawienia y’(t)=f(t, y), wówczas otrzymujemy:
)
(
"
2
))
(
,
(
)
(
)
(
2
1
i
i
i
i
i
y
h
t
y
t
hf
t
y
t
y
(13)
Zapisując ω
i
≈y(t
i
) oraz pomijając błąd przybliżenia otrzymujemy ω
i+1
=ω
i
+hf(t
i
, ω
i
) (14)
76
Powyższy wzór nazywamy metodą Eulera – wzór ten nazywany jest inaczej równaniem
różniczkowym, gdyż można zapisać:
h
t
f
i
i
i
i
1
)
,
(
(15)
Aby wyznaczyć wartość szukanej funkcji y(x) w następnym kroku h, wykorzystujemy
poprzednią wartość funkcji oraz wielkości zmian funkcji – dzięki pochodnej. Natomiast
uwzględniając błąd przybliżenia wzór (13) przyjmuje postać:
i
i
i
i
i
h
t
y
t
f
1
))
(
,
(
, gdzie
)
(
"
2
2
i
i
y
h
(16)
Lokalny błąd dyskretyzacji τ
i+1
(h)
)
(
1
)
,
(
(
1
)
(
1
1
1
1
i
i
i
i
i
i
i
y
h
t
f
y
y
h
h
(17)
dodatkowo możemy określić krok h, dla którego błąd lokalny jest mniejszy od zadanej
dokładności δ
M
2
h
, gdzie
)
(
"
y
max
M
)
t
,
t
(
1
i
i
Globalny błąd dyskretyzacji g(x)
g(t)= ω(t)-y(t) (19)
Dla wybranego punktu t
i
możemy zapisać:
g
i
=ω
i
-y(t
i
) (20), wówczas
1
2
)
(
)
(
0
t
t
L
i
i
i
i
e
L
hM
t
y
g
, gdzie L- liczba Lipschnitz`a
77
Lokalny błąd dyskretyzacji
Globalny błąd dyskretyzacji
nie będę wskazywał na tego, kto nie narysował go ;)))
Stosując metodę Eulera rozwiązać następujące równanie różniczkowe
y’=y-t
2
+1 0<t<2 y(0)=0,5
dla N=10, wtedy h=0,2; t
i
=ih, ω=0,5 oraz ω
i+1
=ω
i
+hf(t
i
, ω
i
), gdzie f(t,y)=y-t
2
+1
ω
i+1
=ω
i
+h(ω
i
-t
2
i
+1)= ω
i
+0,2[ω
i
-0,04i
2
+1]=1,2ω
i
-0,008i
2
+0,02 dla i=0, 1, ..., 9
Rozwiązanie dokładne ma postać
y(t)=(t
2
-1)-0,5e
t
78
Metoda Eulera – wstecz
Powyższe rozważania dotyczyły metody Eulera w przód, ponieważ krok spełniał warunek
h>0. W analogiczny sposób można wyprowadzić metodą Eulera wstecz przyjmując h<0,
wówczas otrzymujemy:
w
i+1
= w
i
+ hf (t
i+1
, w
i+1
)
(22)
w
i+1
(k)
= w
i
+ hf (t
i+1
, w
i+1
(k-1)
)
(23)
Metoda wstecz różni się w stosunku do metody w przód argumentami funkcji f.
Metoda w przód wykorzystuje do obliczenia przybliżenia wartości z poprzedniego kroku,
natomiast metoda wstecz jest równaniem uwikłanym, ponieważ aby obliczyć kolejne
przybliżenie w
i+1
wykorzystujemy wartości z poprzedniego kroku oraz poszukiwaną wartość
w
i+1
. Takiego równania nie można rozwiązać w sposób bezpośredni. Aby rozwiązać takie
równanie (23) należy zastosować proces iteracyjny, czyli poszukujemy kilkakrotnie w
i+1
(k)
,
stojącej po prawej stronie równania podstawiając jako w
i+1
(k-1)
- lewa strona równania, wynik
przybliżenia z poprzedniej iteracji (k-1). Proces trwa do momentu, kiedy spełniony zostanie
warunek:
|w
i+1
(k)
- w
i+1
(k-1)
| ≤ ε (24)
gdzie ε - tolerancja obliczeń
Kody źródłowe:
Metoda Eulera w przód:
Function Euler_order_One_Forward (t0, tn, y0:extended;
n:integer): extended;
var i:integer;
h, wi, wi_1,t :extended;
Begin
h:= (tn-t0)/n;
yi_1:=y0;
for i:=0 to n-1 do
begin
t:=t0+i*h;
wi:=wi_1+h*f(t,wi_1);
wi_1:=wi;
end;
Result:=wi;
end;
79
Przykładowa funkcja definiująca równanie różniczkowe
)
,
( y
t
f
dt
dy
function f(t,y:extended):extended;
begin
result:=y-sqr(t)+1;
end;
Metoda Eulera “wstecz”
function Euler_Order_One_Backward
t0, tn, y0, eps : extended; n:integer):extended;
var i:integer;
h,wi,wi_p,wi_1,t :extended
begin
h:=(tn-t0)/n
wi_1:=y0;
for i:=0 to n-1 do
begin
t:=t0 + i*h;
wi:=wi_1 + h*f (t, wi_1);
repeat
wi_p:=wi;
wi:=wi_1 + h*f (t, wi_p);
until abs(wi-wi_p)<eps;
wi_1:=wi;
end;
Result:=wi;
end;
Przykładowa funkcja definiująca równanie różniczkowe
)
,
( y
t
f
dt
dy
function f(t,y:extended):extended;
begin
result:=y-sqr(t)+1;
end;
Metody Rungego-Kutty
Powszechnie na całym świecie stosuje się metody Rungego–Kutty czwartego rzędu. Polegają
one na rozwiązaniu zagadnienia:
)
,
( y
t
f
dt
dy
gdzie t
a,b
oraz y(a)=
(25)
i przedstawieniu różnicy funkcji y(t) w punktach t
i+1
oraz t
i
w postaci:
w
i+1
- w
i
=
m
i
j
j
k
c
1
lub inaczej w
i+1
– w
i
= h0 (t
i
,w
i
,h) (26)
80
gdzie m oznacza rząd metody, c
j
są stałymi, a
1
1
1
,
,
,
j
j
i
j
i
j
i
l
li
i
j
i
j
k
w
t
hf
k
w
h
t
hf
k
gdzie α
j
, β
j
, γ
j
,
δlj
– stałe
h - krok całkowania
Określenie stałych c
j
, α
j
, β
j
otrzymujemy poprzez rozwinięcie funkcji f(t,y) w szereg Taylora
w otoczeniu punktu t
i
Metody R-K – metoda 2 rzędu
Metoda wyprowadzona przez rozwinięcie f(t,y) w szereg Taylora 2 rzędu pozwala określić
stałe c
1
,
1
,
1
, c
2
,
2
,
2
:
Metoda punktu środkowego:
C
1
= 0
C
2
= 1
α
1
= 0
α
2
= h/2
1
= 0
2
= ½ (28)
wówczas możemy zapisać:
k
1
=hf (t
i
,w
i
) k
2
=hf (t
i
+ ½ h, w
i
+ ½ k
1
)
ostatecznie:
w
i+1
=w
i
+k
2
lub
w
i+1
=w
i
+ hf(t
i
+ h/2, w
i
+ h/2 * f(t
i
, w
i
))
81
Interpretacja graficzna:
f
1
=f(t
i
,w
i
) f
2
=f(t
i
+ h/2, w
i
+ h/2 * f
1
)
Stosując metody Rungego–Kutty 2 rzędu , rozwiązać następujące równanie różniczkowe:
y’=y-t
2
+ 1 0
t
2 y(0)=0,5
dla N = 10 wtedy h = 0,2; t
i
= i
h ; w
0
= 0,5 oraz
metoda punktu środkowego w
i+1
= 1,22 w
i
– 0,0088 i
2
+ 0,218
Metoda zmodyfikowana Eulera
w
i+1
= 1,22 w
i
- 0,0088 i
2
+ 0,216
Metoda Heana dla i = 0, 1, ..., 9
w
i+1
= 1,22 w
i
- 0,0088 i
2
+ 0,2173
Rozwiązanie dokładne ma postać:
y(t) = (t
2
+1) - 0,5 e
t
Metody R-K rzędu 4
Metoda wyprowadzona przez rozwinięcie f(t,y) w szereg Taylora 4 rzędu, pozwala określić
stałe we wzorze (wzór ogólny R-K). Poniżej przedstawiono najczęściej stosowaną metodę 4
rzędu
k
1
= hf (t
i
, w
i
)
k
2
= hf (t
i
+ ½ *h, w
i
+ ½ *k
1
)
k
3
= hf (t
i
+ ½ *h, w
i
+ ½ *k
2
) (37)
82
k
4
= hf (t
i
+h, w
i
+k
3
)
ostatecznie:
w
i+1
=w
i
+ 1/6 (k
1
+ 2k
2
+ 2k
3
+ k
4
) (38)
Interpretacja graficzna:
f
1
= f (t
i
, w
i
)
f
2
= f (t
i
+ ½ * h, w
i+1
+ h f
1
)
f
3
= f (t
i
+ ½ * h, w
i
+ ½ * h f
2
)
f
4
= f (t
i
+ h, w
i
+ h f
3
)
linia 3-4 nie jest łamana – jest prosta !!! (na 99%) :)
Metoda R-K jest najczęściej stosowaną metodą do rozwiązywania układów równań
różniczkowych
Przykład:
Stosując metodę R-K 4-tego rzędu rozwiązać następujące równanie różniczkowe:
y’ = y-t
2
+ 1 0
t
2 y(0) = 0,5
dla N = 10 wtedy h = 0,2 t
1
= i*h w
0
= 0,5 oraz:
Metoda Rungego-Kutty 4-rzedu:
w
i+1
= 1,2214 w
i
- 0,008856 i
2
+ 0,00856 i + 0,218593
dla i = 0, 1, ..., 9
Rozwiązanie dokładne ma postać y(t) = (t
2
+1) - 0,5e
t
83
Kody źródłowe:
function RK_Rzedu_4 (t0, tn, y0: extended;
n:integer):extended;
var i:integer;
h, yi, y_i, t, k1, k2, k3, k4: extended;
begin
h:= (tn-t0)/n
y_i = y0;
for i:=0 to n-1 do
begin
t:=t0 + i*h;
k
1 :
= h*ft, y_i);
k
2 :
= h*f(t + 0.5 * h, y_i+0.5 * k1);
k
3 :
= h*f(t+0.5*h, y_i + 0.5 * k
2
);
k
4 :
= h*f(t+h, y_i + k3);
yi:= y_i+(k1 + 2*k2 + 2*k3 + k4)/6;
y_i:=yi;
end;
Result := yi;
end;
Metody wykorzystujące więcej niż jedną przybliżoną wartość dla poprzednich punktów
węzłowych w celu przybliżenia wartości w następnym punkcie węzłowym nazywamy
metodami wielokrokowymi.
Ogólny wzór na m-krotną metodę wielokrokową dla rozwiązania równania
)
,
( y
t
f
dt
dy
gdzie t
a,b
oraz
y(a) =
(55)
jest równaniem różniczkowym dla znalezienia przybliżenia w
i+1
, dla punktu węzłowego t
i+1
,
gdzie m>1, które przyjmuje postać:
w
i+1
= a
m-1
* w
i
+ a
m-2
* w
i-1
+ ... + a
0
* w
i+1_m
+ h [b
n
f(t
i+1
,w
i+1
) + b
n-1
f(t
i
, w
i
) + ... + t + b
0
f(t
i+1_m
, w
i+1_m
)]
dla i = (m-1, m, ..., N-1),
N - liczba węzłów w przedziale, gdzie h = (b-a) / N
84
Równania różniczkowe cząstkowe z warunkami brzegowymi – rodzaje:
a)
równania różniczkowe eliptyczne
•
równanie Laplace’a,
•
równanie Poisson’a,
b)
równanie różniczkowe paraboliczne
•
równanie dyfuzji
c) równanie różniczkowe hiperboliczne
•
równanie falowe
Równania różniczkowe cząstkowe z warunkiem brzegowym – wstęp
Równania różniczkowe prezentowane w poprzednich rozdziałach pozwalały na wyznaczenie
poszukiwanej funkcji jednej zmiennej dla podanych warunków początkowych lub
brzegowych. Problemy przedstawione w tym rozdziale będą opisywane funkcjami wielu
zmiennych.
W zależności od specyfiki problemu można opisać go za pomocą odpowiedniego
równania różniczkowego: Laplace’a, Poisson’a, Dyfuzji, Falowego, itp… . Aby rozwiązać
problem opisany odpowiednim równaniem różniczkowym należy określić warunki brzegowe.
Warunki brzegowe mogą być opisane w formie warunków Dirichleta lub Neumanna.
Równania różniczkowe cząstkowe – warunek brzegowy Dirichleta
Jeżeli funkcja u(x,y) jest poszukiwana na obszarze R, wówczas musi mieć zdefiniowane
warunki w postaci funkcji g(x,y) na brzegu S.
Warunek brzegowy Dirichleta
y
,
x
g
y
,
x
u
dla
S
y
,
x
gdzie: u(x,y) – poszukiwana funkcja w punktach wewnątrz obszaru R, g(x,y) – zadana
funkcja dla punktów (x,y) należących do brzegu S obszaru R.
Wynika z tego, że w trakcie obliczeń wartości funkcji g w poszczególnych punktach (x,y)
S
nie mogą się zmienić.
85
Równania różniczkowe cząstkowe – warunek brzegowy Neumanna
Warunek brzegowy Neumanna
)
y
,
x
(
g
)
y
,
x
(
n
u
dla
S
y
,
x
gdzie:
)
y
,
x
(
n
u
– pochodna normalna poszukiwanej funkcji w punktach należących do
brzegu S obszaru R, g(x,y) – zadana funkcja dla punktów (x,y) należących do brzegu S
obszaru R.
Oznacza to, że na brzegu obszaru zmienność funkcji u(x,y) w kierunku normalnym dla
punktów (x,y)
S jest równa funkcji g(x,y).
Dla zmiennej x
)
y
,
x
(
g
h
)
y
,
h
x
(
u
)
y
,
x
(
u
Dla zmiennej y
)
y
,
x
(
g
h
)
h
y
,
x
(
u
)
y
,
x
(
u
Równania różniczkowe cząstkowe z warunkiem brzegowym - Równania eliptyczne -
Poissona i Laplace’a
Równania różniczkowe cząstkowe eliptyczne znane jako równanie Poissona , dla dwóch
wymiarów i prostokątnego układu współrzędnych przyjmuje postać:
)
y
,
x
(
f
)
y
,
x
(
dy
u
d
)
y
,
x
(
dx
u
d
)
y
,
x
(
u
2
2
2
2
2
(1)
na
d
y
c
b
x
a
y
x
R
;
)
,
(
z warunkami brzegowymi
)
,
(
)
,
(
y
x
g
y
x
u
dla
S
y
x
)
,
(
W powyższym równaniu funkcja f opisuje dane wejściowe na płaszczyźnie dla obszaru R
ograniczonego brzegiem S. Równania tego typu są stosowane dla różnych problemów
fizycznych, które są niezależne od czasu. Najczęściej stosowane są do obliczenia rozkładu
potencjału (np. temperatury) w stanie ustalonym.
Szczególnym przypadkiem równania Poissona gdy f(x,y)=0 jest równanie Laplace’a.
Aby rozwiązać równanie cząstkowe eliptyczne (1) zastosujemy metodę różnic skończonych
MRS.
W tym celu weźmiemy dwie liczby całkowite m i n, i zdefiniujemy kroki całkowania h = (b-
a)/n oraz k = (c-d)/m. Dzięki temu możemy podzielić: przedział [a,b] na n równych części o
szerokości h oraz przedział [c,d] na m równych części o szerokości k.
Umieśćmy siatkę na obszarze R poprzez narysowanie pionowych i poziomych linii
przechodzących przez punkty (xi, yj), takich że:
86
ih
a
x
i
dla i = 0,1, .. n
oraz
jk
c
y
j
dla j = 0,1, .. m
Linie x = xi oraz y = yj nazywamy liniami siatki, natomiast ich punkty przecięcia (xi,yj)
nazywamy punkami węzłowymi siatki. Dla wewnętrznych punktów węzłowych (xi,yj),
które nie należą do brzegu obszaru R, zastosujemy metodę różnic skończonych – MRS.
Metoda ta opiera się na zastąpieniu pochodnych cząstkowych rozwinięciem funkcji w szereg
Taylora w otoczeniu punktu (xi,yj) – (patrz druga pochodna numeryczna) . Wówczas
otrzymujemy:
Dla zmiennej x
)
,
(
12
)
,
(
)
,
(
2
)
,
(
)
,
(
4
4
2
2
1
1
2
2
j
i
j
i
j
i
j
i
j
i
y
x
u
h
h
y
x
u
y
x
u
y
x
u
y
x
x
u
(2)
gdzie:
)
,
(
1
1
i
i
i
x
x
Dla zmiennej y
)
,
(
12
)
,
(
)
,
(
2
)
,
(
)
,
(
4
4
2
2
1
1
2
2
j
i
j
i
j
i
j
i
j
i
x
y
u
k
k
y
x
u
y
x
u
y
x
u
y
x
y
u
(3)
gdzie:
)
,
(
1
1
j
j
j
y
y
87
Podstawiając wzory (2) oraz (3) do równania Poissona (1) otrzymujemy:
2
1
1
2
1
1
)
,
(
)
,
(
2
)
,
(
)
,
(
)
,
(
2
)
,
(
k
y
x
u
y
x
u
y
x
u
h
y
x
u
y
x
u
y
x
u
j
i
j
i
j
i
j
i
j
i
j
i
)
,
(
12
)
,
(
12
)
,
(
4
4
2
4
4
2
j
i
j
i
x
y
u
k
y
x
u
h
y
x
f
(4)
dla i = 1,2 … n-1 oraz j = 1,2, … m-1
natomiast warunki brzegowe mają postać:
)
,
(
)
,
(
0
0
j
j
y
x
g
y
x
u
oraz
)
,
(
)
,
(
j
n
j
n
y
x
g
y
x
u
dla j = 0,1, .. m (5)
)
,
(
)
,
(
m
i
m
i
y
x
g
y
x
u
oraz
)
,
(
)
,
(
0
0
y
x
g
y
x
u
i
i
dla i = 1,2, .. n (6)
Ostatecznie wzór na MRS, podstawiając za (xi,yj) = wij oraz wyłączając oszacowanie błędu
zapisujemy:
)
,
(
)
(
)
(
1
2
2
1
,
1
,
2
,
1
,
1
2
j
i
j
i
j
i
j
i
j
i
ij
y
x
f
h
w
w
k
h
w
w
w
k
h
(7)
dla i = 1,2 … n-1 oraz j = 1,2, … m-1
natomiast warunki brzegowe mają postać:
)
,
(
0
0
j
j
y
x
g
w
oraz
)
,
(
j
n
nj
y
x
g
w
dla j = 0,1, .. m (8)
)
,
(
0
0
y
x
g
w
i
i
oraz
)
,
(
m
i
im
y
x
g
w
dla i = 1,2, .. n (9)
gdzie w
ij
jest przybliżeniem u(x
i
,y
j
).
88
Analizując równanie (7) można zauważyć, że w celu wyznaczenia przybliżenie rozwiązania w
punkcie (xi,yj), potrzebne są wartości przybliżenia rozwiązania w czterech sąsiednich
punktach – patrz rysunek obok.
Postać macierzowa
Jak widać poszukujemy rozwiązania równania Poissona dla wewnętrznych węzłów siatki, w
tym celu dla każdego punktu siatki układamy równanie (7). Po uwzględnieniu warunków
brzegowych otrzymujemy układ (n-1)
(m-1) równań liniowych. Układ równań możemy
rozwiązać metodami deterministycznymi bądź iteracyjnymi.
Zakładając, że h = k możemy wyprowadzić uproszczoną postać równania (7) :
)
,
(
)
(
)
(
4
2
1
,
1
,
,
1
,
1
j
i
j
i
j
i
j
i
j
i
ij
y
x
f
h
w
w
w
w
w
(10)
dla i = 1,2 … n-1 oraz j = 1,2, … m-1
natomiast warunki brzegowe są takie same jak (8) i (9) :
Postać iteracyjna
Równanie (7) i (10) można przedstawić również w postaci iteracyjnej, wyprowadzając wij z
każdego równania, gdzie l oznacza numer iteracji :
z równania (7)
1
2
)
(
)
(
)
,
(
2
)
(
1
,
)
(
1
,
2
)
(
,
1
)
(
,
1
2
)
1
(
k
h
l
j
i
l
j
i
k
h
l
j
i
l
j
i
j
i
l
ij
w
w
w
w
y
x
f
h
w
(11)
z równania (10)
4
)
(
)
,
(
)
(
1
,
)
(
1
,
)
(
,
1
)
(
,
1
2
)
1
(
l
j
i
l
j
i
l
j
i
l
j
i
j
i
l
ij
w
w
w
w
y
x
f
h
w
(12)
89
Równania eliptyczne - przykład 1
Wyznaczyć rozkład temperatury w stanie ustalonym dla cienkiej kwadratowej
metalowej płytki o wymiarach 0,5 m na 0,5 m. Na brzegu płytki znajdują się źródła
ciepła utrzymujące temperaturę na poziomie 0°C dla boku dolnego i prawego, natomiast
temperatura boku górnego i lewego zmienia się liniowo od 0°C do 100°C Problem
rozwiązać układając układ równań liniowych (postać macierzowa) dla wewnętrznych
węzłów siatki 5 x 5 – układ równań rozwiązać metodą Gaussa -Siedla.
Siatka dyskretyzacyjna 5 x 5
Problem ten opisuje równanie Laplace’a
0
)
,
(
)
,
(
)
,
(
2
2
2
2
2
y
x
dy
T
d
y
x
dx
T
d
y
x
T
z warunkami brzegowymi:
1) T(0,y) = 0 [°C ]
2) T(x,0) = 0 [°C ]
3) T(0,5,y) = 200y [°C ]
4) T(x,0,5) = 200x [°C ]
Układ równań tworzymy w oparciu o wzór (10)
25
0
0
50
0
0
100
50
25
4
1
0
1
0
0
0
0
0
1
4
1
0
1
0
0
0
0
0
1
4
0
0
1
0
0
0
1
0
0
4
1
0
1
0
0
0
1
0
1
4
1
0
1
0
0
0
1
0
1
4
0
0
1
0
0
0
1
0
0
4
1
0
0
0
0
0
1
0
1
4
1
0
0
0
0
0
1
0
1
4
3
,
1
2
,
1
1
,
1
3
,
2
2
,
2
1
,
2
3
,
3
2
,
3
1
,
3
T
T
T
T
T
T
T
T
T
90
Rozwiązanie układu równań metodą Gaussa-Siedla, daje wyniki
75
,
18
50
,
12
25
,
6
50
,
37
00
,
25
50
,
12
25
,
56
50
,
38
75
,
18
3
,
1
2
,
1
1
,
1
3
,
2
2
,
2
1
,
2
3
,
3
2
,
3
1
,
3
T
T
T
T
T
T
T
T
T
Ostatecznie wyniki dla założonej siatki mają postać:
00
,
0
00
,
0
00
,
0
00
,
0
00
,
0
0
00
,
25
75
,
18
50
,
12
25
,
6
00
,
0
125
,
0
00
,
50
50
,
37
00
,
25
50
,
12
00
,
0
25
,
0
00
,
75
25
,
56
50
,
37
75
,
18
00
,
0
375
,
0
00
,
100
00
,
75
00
,
50
00
,
25
00
,
0
5
,
0
5
,
0
375
,
0
25
,
0
125
,
0
0
y
y
y
y
y
x
x
x
x
x
Rozwiązanie dokładne problemu ma postać:
xy
y
x
T
400
)
,
(
Natomiast oszacowanie błędu obliczeń :
0
)
,
(
)
,
(
12
4
4
4
4
2
j
i
j
i
x
y
T
y
x
T
h
, bo
0
)
400
(
)
,
(
4
4
4
4
xy
x
y
x
x
T
i
0
)
400
(
)
,
(
4
4
4
4
xy
y
y
x
y
T
Wynika z tego, że rozwiązanie numeryczne jest równe z rozwiązaniu dokładnemu.
Równania eliptyczne - przykład 1 - kody źródłowe
Funkcje charakteryzujące dane równanie różniczkowe eliptyczne
function F( X,Y : real ) : real;
begin
F := 0;
end;
function G( X,Y : real ) : real;
begin
if X = 0 then G := 0;
if X = 0,5 then G := 200*Y;
if Y = 0 then G := 0;
if Y = 0,5 then G := 200*X;
end;
Definicja wymaganych stałych i typów
91
const max = 4;
type
Macierz = array [ 0..max, 0..max ] of real;
Vector = array [ 0..max ] of real;
Funkcja obliczająca przybliżenie równania różniczkowego eliptycznego
function PoissonEq(A,B,C,D,TOL : real;N,M,NMAX : integer; var W : Macierz)
:boolean;
var
X,Y : Vector; H,K,V,VV,Z,E : real;
M1,M2,N1,N2,I,J,L,LL: integer;
Ok :boolean;
begin
M1 := M - 1; M2 := M - 2;
N1 := N - 1; N2 := N - 2;
H := ( B - A ) / N; K := ( D - C ) / M;
for I := 0 to N do X[I] := A + I * H;
for J := 0 to M DO Y[J] := C + J * K;
for I := 1 to N1 do begin
W[I,0] := g(X[I],Y[0]);
W[I,M] := g(X[I],Y[M])
end;
for J := 0 to M do begin
W[0,J] := g(X[0],Y[J]);
W[N,J] := g(X[N],Y[J])
end;
for I := 1 to N1 do for J := 1 to M1 do W[I,J] := 0.0;
V := H * H / ( K * K ); VV := 2.0 * ( 1.0 + V );
L := 1; OK := false;
while ( ( L <= NMAX ) and ( not OK ) ) do begin
Z := (-H*H*F(X[1],Y[M1])+G(A,Y[M1])+V*G(X[1],D)+V*W[1,M2]+W[2,M1])/VV;
E := abs( W[1,M1] - Z ); W[1,M1] := Z;
for I := 2 to N2 do begin
Z := (-H*H*F(X[I],Y[M1])+V*G(X[I],D)+W[I-1,M1]+W[I+1,M1]+V*W[I,M2])/VV;
if ( abs( W[I,M1] - Z ) > E ) then E := abs( W[I,M1] - Z );
W[I,M1] := Z
end;
Z := (-
H*H*F(X[N1],Y[M1])+G(B,Y[M1])+V*G(X[N1],D)+W[N2,M1]+V*W[N1,M2])/VV;
if ( abs( W[N1,M1] - Z ) > E ) then E := abs( W[N1,M1] - Z ); W[N1,M1] := Z;
for LL := 2 to M2 do begin
J := M2 - LL + 2; Z := (-H*H*F(X[1],Y[J])+G(A,Y[J])+V*W[1,J+1]+V*W[1,J-
1]+W[2,J])/VV;
if ( abs( W[1,J] - Z ) > E ) then E := abs( W[1,J] - Z );
W[1,J] := Z;
for I := 2 to N2 do begin
Z := (-H*H*F(X[I],Y[J])+W[I-1,J]+V*W[I,J+1]+V*W[I,J-1]+W[I+1,J])/VV;
if ( abs( W[I,J] - Z ) > E ) then E := abs( W[I,J] - Z );
92
W[I,J] := Z
end;
Z := (-H*H*F(X[N1],Y[J])+G(B,Y[J])+W[N2,J]+V*W[N1,J+1]+V*W[N1,J-1])/VV;
if ( abs( W[N1,J] - Z ) > E ) then E := abs( W[N1,J] - Z ); W[N1,J] := Z
end;
Z := ( -H * H * F( X[1],Y[1] ) + V * G( X[1], C ) + G( A, Y[1] ) + V * W[1,2] + W[2,1] ) /
VV;
if ( abs( W[1,1] - Z ) > E ) then E := abs( W[1,1] - Z );
W[1,1] := Z;
for I := 2 to N2 do begin
Z := (-H*H*F(X[I],Y[1])+V*G(X[I],C)+ W[I+1,1]+W[I-1,1]+V*W[I,2])/VV;
if ( abs( W[I,1] - Z ) > E ) then E := abs( W[I,1] - Z );
W[I,1] := Z
end;
Z := (-H*H*F(X[N1],Y[1])+V*G(X[N1],C)+G(B,Y[1])+W[N2,1]+V*W[N1,2])/VV;
if ( abs( W[N1,1] - Z ) > E ) then E := ABS( W[N1,1] - Z );
W[N1,1] := Z;
if ( E <= TOL ) then OK := true
else L := L + 1;
end;
Result := OK;
end;
93
Równania eliptyczne - przykład 2
Wyznaczyć rozkład temperatury w stanie ustalonym dla cienkiej kwadratowej metalowej
płytki o wymiarach 1 m na 1 m. Na trzech brzegach płytki znajdują się źródła ciepła
utrzymujące temperaturę na poziomie 100°C dla boku górnego, 50°C dla boku lewego, 0°C
dla boku dolnego (warunki Dirichleta), natomiast przy czwartym lewym boku jest
umieszczony izolator termiczny (warunek Neumanna – pochodna normalna równa jest 0).
Problem rozwiązać wg postaci iteracyjnej dla wewnętrznych węzłów siatki 101 x 101.
Tolerancję obliczeń przyjąć na poziomie 10
-6
.
Problem ten opisuje równanie Laplace’a
0
)
,
(
)
,
(
)
,
(
2
2
2
2
2
y
x
dy
T
d
y
x
dx
T
d
y
x
T
z warunkami brzegowymi:
1) T(0,y) = 0 [°C ]
2) T(x,0) = 50 [°C ]
3) T(1,y) = T(1-h,y) [°C ]
4) T(x,1) = 100 [°C ]
94
Przybliżenie rozwiązania równania w postaci mapy barwnej po 13484 iteracjach
Tabela wyników dla wybranych punków prawego boku płytki - izolator
Numer węzła
Numer iteracji
Dokładne
0
1000
3000
7000
13000
T100,0
100,00
100,00
100,00
100,00
100,00
100,00
T100,10
0,00
75,17
86,03
89,59
89,92
90,00
T100,20
0,00
52,58
72,46
79,24
79,87
80,00
T100,30
0,00
34,03
59,66
69,00
69,86
70,00
T100,40
0,00
20,27
47,91
58,89
59,90
60,00
T100,50
0,00
11,08
37,37
48,90
49,97
50,00
T100,60
0,00
5,55
28,07
39,03
40,05
40,00
T100,70
0,00
2,55
19,91
29,22
30,09
30,00
T100,80
0,00
1,07
12,70
19,47
20,10
20,00
T100,90
0,00
0,37
6,17
9,73
10,06
10,00
T100,100
0,00
0,00
0,00
0,00
0,00
0,00
95
Równania eliptyczne - przykład 2 - kody źródłowe
Funkcje charakteryzujące dane równanie różniczkowe eliptyczne
function F( X,Y : real ) : real;
begin
F := 0;
end;
Definicja wymaganych stałych i typów
const max = 100;
type
Macierz = array [ 0..max, 0..max ] of real;
Vector = array [ 0..max ] of real;
function PoissonEquationIteration (a,b,c,d, TOL : real ; N,M ,nmax : integer; var V :
Macierz) : boolean;
var EPS,H,K,Vs : real; X,Y : Vector; I,J,L : integer;
begin
H := ( B - A ) / N; K := ( D - C ) / M;
for I := 0 to N do X[I] := A + I * H; for J := 0 to M do Y[J] := C + J * K;
for i:=0 to N do for j:=0 to M do V[i,j]:=0;
L := 0; Result := true;
repeat
eps:=0; L := L + 1;
for i:=0 to N do for j:=0 to M do begin
Vs:=V[i,j];
if j=0 then V[i,j]:=0 else // bok dolny
if j=M then V[i,j]:=100 else // bok gora
if i=0 then V[i,j]:= 50 else // bok lewo
if i=N then V[i,j]:= V[i-1,j] else // bok prawo
V[i,j]:=(- h*h* F(X[i],Y[j] ) + (V[i+1,j]+V[i-1,j])+sqr(h/k)*(V[i,j+1]+V[i,j-1])
)/(2*(sqr(h/k)+1));
If abs(V[i,j]-Vs) > eps then eps := abs(V[i,j]-Vs);
end;
until (eps< TOL ) or ( L > nmax);
if (L > nmax) then Result := false;
end;
96
Równania różniczkowe cząstkowe z warunkiem brzegowym
Równania paraboliczne - Dyfuzja
Równania różniczkowe cząstkowe paraboliczne znane jako równanie dyfuzji lub
przewodnictwa, w zależności od jednego wymiaru oraz czasu przyjmuje postać:
dla oraz
(13)
z warunkami brzegowymi dla
i początkowymi dla
Równania tego typu są stosowane dla różnych problemów fizycznych, które są zależne od
czasu. Najczęściej stosowane są do obliczenia przepływu ciepła wzdłuż pręta o długości l,
przy założeniu że, powierzchnia boczna pręta jest odizolowana od otoczenia. Stała
jest
niezależna od położenia oraz czasu i najczęściej określa przewodność cieplną materiału z
którego zrobiony jest pręt. Funkcja f określa początkowy rozkład temperatury w pręcie.
Aby rozwiązać równanie cząstkowe paraboliczne (13) zastosujemy metodę różnic
skończonych MRS.
W tym celu weźmiemy liczbę
całkowitą m i zdefiniujemy
krok całkowania h = (b-a)/m.
Dodatkowo ustalimy wartość
kroku czasowego k. Dzięki
temu możemy wyznaczyć
węzły siatki (xi , tj ), gdzie:
dla i = 0,1, .. m oraz dla j = 0,1, ..
Dla wewnętrznych punktów węzłowych (xi,tj), które nie należą do brzegu, zastosujemy
metodę różnic skończonych – MRS. Metoda ta opiera się na zastąpieniu pochodnych
cząstkowych rozwinięciem funkcji w szereg Taylora w otoczeniu punktu (xi,tj) – (patrz
pierwsza i druga pochodna numeryczna) . Wówczas otrzymujemy:
Dla zmiennej x (14)
gdzie:
Dla zmiennej t (15)
gdzie:
0
)
,
(
)
,
(
2
2
2
t
x
dx
u
d
t
x
dt
du
l
x
0
0
t
0
)
,
(
)
,
0
(
t
l
u
t
u
0
t
)
(
)
0
,
(
x
f
x
u
l
x
0
ih
x
i
jk
t
j
)
,
(
12
)
,
(
)
,
(
2
)
,
(
)
,
(
4
4
2
2
1
1
2
2
j
i
j
i
j
i
j
i
j
i
t
x
u
h
h
t
x
u
t
x
u
t
x
u
y
x
x
u
)
,
(
1
1
i
i
i
x
x
)
,
(
2
)
,
(
)
,
(
)
,
(
2
2
1
j
i
j
i
j
i
j
i
x
t
u
k
k
t
x
u
t
x
u
t
x
t
u
)
,
(
1
j
j
j
t
t
97
j
i
j
i
j
i
j
i
w
w
h
k
w
h
k
w
,
1
,
1
2
2
,
2
2
1
,
2
1
Podstawiając wzory (14) oraz (15) do równania dyfuzji (13) oraz wyłączając oszacowanie
błędu zapisujemy:
(16)
warunek brzegowy : (17)
dla i = 1,2 … m-1 oraz j = 1,2, …
natomiast lokalny błąd odcięcia jest równy :
(18)
gdzie: oraz
Ostatecznie wzór na MRS, podstawiając za ( xi , tj ) = wi,j oraz wyłączając oszacowanie
błędu zapisujemy:
(19)
przekształcając powyższy wzór ze względu na wi,j+1 ,wówczas otrzymujemy:
(20)
lub
gdzie: (21)
dla i = 1,2 … m-1 oraz j = 1,2, …
0
)
,
(
)
,
(
2
)
,
(
)
,
(
)
,
(
2
1
1
2
1
h
t
x
u
t
x
u
t
x
u
k
t
x
u
t
x
u
j
i
j
i
j
i
j
i
j
i
0
)
,
(
)
,
0
(
j
j
t
l
u
t
u
)
,
(
12
)
,
(
2
4
4
2
2
2
2
j
i
j
i
ij
t
x
u
h
x
t
u
k
)
,
(
1
1
i
i
i
x
x
)
,
(
1
j
j
j
t
t
0
2
2
,
1
,
,
1
2
,
1
,
h
w
w
w
k
w
w
j
i
j
i
j
i
j
i
j
i
j
i
j
i
j
i
j
i
w
w
w
w
,
1
,
1
,
1
,
2
1
2
h
k
98
Równania paraboliczne – Schemat jawny
Analizując równanie (21) można zauważyć,
że w celu wyznaczenia przybliżenia
rozwiązania w punkcie (xi,tj+1),
potrzebne są wartości trzech punktów z
poprzedniego kroku czasowego tj - patrz
rysunek obok. Metoda ta nazywana jest
rozwiązywaniem równania parabolicznego
wg schematu jawnego
Jak
widać
poszukujemy
rozwiązania
równania dyfuzji dla wewnętrznych węzłów siatki, w tym celu dla każdego punktu siatki
układamy równanie (21). Obliczenia rozpoczynamy od kroku czasowego j = 1, gdyż najpierw
musimy wygenerować wartości początkowe dla t = 0 dla wszystkich punktów u( xi , 0).
Obliczmy je zgodnie z warunkami początkowymi u( xi , 0) = f(xi).
A więc mamy do rozwiązania układ (m-1) równań, który powstał poprzez ułożenie równania
(21) dla każdego punktu wi,j siatki dla i = 1,2, .. m-1 :
(22)
co rozpisując daje nam:
(23)
Czyli w ten sposób możemy wyznaczyć poszukiwane wartości w(j+1). Wystarczy tylko
przemnożyć macierz A przez wektor w(j) - wartości dla poprzedniego kroku.
)
(
)
1
(
j
j
Aw
w
j
m
j
j
j
m
j
j
w
w
w
w
w
w
,
1
,
2
,
1
1
,
1
1
,
2
1
,
1
)
2
1
(
0
0
0
0
0
0
)
2
1
(
99
Równania paraboliczne - Dyfuzja - przykład 1
Rozwiązać poniższe równanie przewodnictwa:
dla oraz
z warunkami brzegowymi dla
i początkowymi dla
Obliczenia wykonać dla czasu t = 0.5 używając schematu jawnego najpierw dla h = 0.1 oraz k
= 0.0005 , następnie dla h = 0.1 oraz k = 0.01 ,przyjąć
2 = 1
Rozwiązanie dokładne ma postać:
Tabela wyników – schemat jawny
Mały krok czasu
Duży krok czasu
xi
u( xi , 0.5 )
wi,1000
k = 0.0005
| u( xi , 0.5 ) –
wi,1000|
wi,50
k = 0.01
| u( xi , 0.5 ) –
wi,50|
0,0
0
0
0
0,1
0,00222241 0,00228652 6,411
10-5
8,19876
107 8,199
107
0,2
0,00422728 0,00434922 1,219
10-4
-1,55719
108 1,557
108
0,3
0,00581836 0,00598619 1,678
10-4
2,13833
108 2,138
108
0,4
0,00683989 0,00703719 1,973
10-4
-2,50642
108 2,506
108
0,5
0,00719788 0,00739934 2,075
10-4
2,62685
107 2,627
108
0,6
0,00683989 0,00703719 1,973
10-4
-2,49015
108 2,490
108
0,7
0,00581836 0,00598619 1,678
10-4
2,11200
107 2,112
108
0,8
0,00422728 0,00434922 1,219
10-4
-1,53086
108 1,531
108
0,9
0,00222241 0,00228652 6,411
10-5
8,03604
107 8,036
107
1,0
0
0
0
0
)
,
(
)
,
(
2
2
2
t
x
dx
u
d
t
x
dt
du
1
0
x
0
t
0
)
,
1
(
)
,
0
(
t
u
t
u
0
t
)
sin(
)
0
,
(
x
x
u
1
0
x
)
sin(
)
,
(
2
x
e
t
x
u
t
100
Interpretacja graficzna
Równania paraboliczne – Schemat jawny
Analizując wyniki z przykładu 1 można zauważyć, że metoda wykorzystująca schemat jawny
jest warunkowo stabilna, gdyż dla zbyt dużych wartości kroku wyniki są rozbieżne w
stosunku do rozwiązania dokładnego. Wynika to wartości własnej macierzy A , dla której
powinien być spełniony warunek
(23)
z powyższego wyrażenia wynika, że warunek stabilności tej metody jest:
(24)
Wykorzystując powyższy wzór można wyznaczyć dla przykładu 1 graniczną wartość kroku
czasowego, dla której schemat jawny jest stabilny:
wtedy
1
sin
4
1
max
)
(
2
2
1
1
m
i
m
i
A
2
1
2
2
h
k
2
2
1
h
k
0,005
1
1
,
0
2
1
2
k
101
0
)
,
(
)
,
(
2
)
,
(
)
,
(
)
,
(
2
1
1
2
1
h
t
x
u
t
x
u
t
x
u
k
t
x
u
t
x
u
j
i
j
i
j
i
j
i
j
i
Równania paraboliczne – Schemat niejawny
Poprzednia metoda wykorzystująca schemat jawny, była warunkowo stabilna, gdyż należy
zwracać uwagę na wielkości kroków całkowania. Metoda opisana poniżej nie ma tej wady, i
jest bezwarunkowo stabilna dla różnych wielkości kroków całkowania. Metoda ta
wykorzystuje schemat niejawny gdyż zalicz się do metod pośrednich wymagających
wielokrotnych iteracji dla danej chwili czasu. Różnica polega na innym określeniu pochodnej
cząstkowej po czasie:
Dla zmiennej t (25)
gdzie:
Wówczas podstawiając wzór (25) oraz (14) do równania parabolicznego (13), pomijając błąd
obcięcia otrzymamy :
(26)
Modyfikując powyższy wzór, podstawiając za ( xi , tj ) = wi,j zapisujemy:
(27)
lub:
gdzie: (28)
dla i = 1,2 … m-1 oraz j = 1,2, …
Analizując równanie (28) można
zauważyć, że w celu wyznaczenia
przybliżenia rozwiązania w punkcie
(xi,tj), potrzebne są wartości trzech
punktów:
dwóch
sąsiednich
z
aktualnego kroku czasowego oraz
wartość dla tego samego węzła, ale
dla poprzedniej iteracji tj-1 - patrz
rysunek obok. Metoda ta nazywana
jest
rozwiązywaniem
równania
parabolicznego
wg
schematu
niejawnego
Jak widać poszukujemy rozwiązania równania dyfuzji dla wewnętrznych węzłów siatki, w
tym celu dla każdego punktu siatki układamy równanie (28). Obliczenia rozpoczynamy od
kroku czasowego j = 1, gdyż najpierw musimy wygenerować wartości początkowe dla t = 0
dla wszystkich punktów u( xi , 0). Obliczmy je zgodnie z warunkami początkowymi u( xi , 0)
= f(xi).
)
,
(
2
)
,
(
)
,
(
)
,
(
2
2
1
j
i
j
i
j
i
j
i
x
t
u
k
k
t
x
u
t
x
u
t
x
t
u
)
,
(
1
j
j
j
t
t
0
2
2
,
1
,
,
1
2
1
,
,
h
w
w
w
k
w
w
j
i
j
i
j
i
j
i
j
i
1
,
,
1
,
1
,
)
(
)
2
1
(
j
i
j
i
j
i
j
i
w
w
w
w
2
h
k
102
1
,
1
1
,
2
1
,
1
,
1
,
2
,
1
)
2
1
(
0
0
0
0
0
0
)
2
1
(
j
m
j
j
j
m
j
j
w
w
w
w
w
w
A więc mamy do rozwiązania układ (m-1) równań, który powstał poprzez ułożenie równania
(28) dla każdego punktu wi,j siatki dla i = 1,2, .. m-1 :
(29)
co rozpisując daje nam:
(30)
Czyli w ten sposób możemy wyznaczyć poszukiwane wartości w(j). Najlepiej do rozwiązania
tego typu układów użyć rozkładu Crouta dla macierzy trójdiagonalnych
Równania paraboliczne - Dyfuzja - przykład 2 (1)
Rozwiązać poniższe równanie przewodnictwa:
dla oraz
z warunkami brzegowymi dla
i początkowymi dla
Obliczenia wykonać dla czasu t = 0.5 używając schematu niejawnego dla h = 0.1 oraz k
= 0.01 ,przyjąć
2 = 1
Rozwiązanie dokładne ma postać:
)
1
(
)
(
j
j
w
Aw
0
)
,
(
)
,
(
2
2
2
t
x
dx
u
d
t
x
dt
du
1
0
x
0
t
0
)
,
1
(
)
,
0
(
t
u
t
u
0
t
)
sin(
)
0
,
(
x
x
u
1
0
x
)
sin(
)
,
(
2
x
e
t
x
u
t
103
Tabela wyników – schemat niejawny
xi
u( xi , 0.5 )
wi,50
k = 0.01
| u( xi , 0.5 )
– wi,50|
0,0
0
0
0,1
0,00222241 0,00289802 6,756
10-4
0,2
0,00422728 0,00551236 1,285
10-3
0,3
0,00581836 0,00758711 1,769
10-3
0,4
0,00683989 0,00891918 2,079
10-3
0,5
0,00719788 0,00937818 2,186
10-3
0,6
0,00683989 0,00891918 2,079
10-3
0,7
0,00581836 0,00758711 1,769
10-3
0,8
0,00422728 0,00551236 1,285
10-3
0,9
0,00222241 0,00289802 6,756
10-4
1,0
0
0
Równania paraboliczne – Schemat niejawny
Analizując wyniki z przykładu 2 można zauważyć, że metoda wykorzystująca schemat
niejawny jest bezwarunkowo stabilna, gdyż jest ona zawsze zbieżna nie zależnie od wielkości
kroków całkowania. Wynika to ze spełniania poniższego warunku dla promienia spektralnego
macierzy A:
(31)
Równania paraboliczne – kody źródłowe
Funkcja określająca warunki początkowe równania paraboliczne
function F( X,Y : real ) : real;
begin
F := sin( pi * x );
end;
Definicja wymaganych stałych i typów
const
max = 100;
type
Vector = array [ 0..max ] of real;
1
sin
4
1
max
)
(
2
2
1
1
m
i
m
i
A
104
Funkcji „HeatingForward” - schemat jawny
function HeatingForward (FX,ALPHA ,FT : real; M,N : integer; var W : Vector):
boolean ;
var
H,K,VV : real; I,J : integer; W_1 : Vector;
begin
H := FX / M; K := FT / N;
VV := ALPHA * ALPHA * K / ( H * H );
if K <= (0.5*SQR(H/ALPHA)) then begin
W_1[0] := 0; W_1[M] := 0; Result := True;
W[0] := 0; W[M] := 0;
for I := 1 to M-1 do W_1[I] := F( I * H );
for J := 1 to N do begin
for I := 1 to M-1 do W[I] := (1.0 - 2.0 * VV)*W_1[I] + VV*( W_1[I-1] + W_1[I+1]
);
W_1 := W;
end;
end
else Result := False;
end;
Funkcji „HeatingBackward” - schemat niejawny
function HeatingBackward (FX,alpha ,FT : real; m,n : integer; var W : Vector):
boolean ;
var
L,U,Z : Vector; H,K,VV : real; I1,I,J : integer;
begin
Result := true;
H := FX / M; K := FT / N;
VV := ALPHA * ALPHA * K / ( H * H );
W[0] := 0; W[M] := 0;
for I := 1 to M-1 do W[I] := F( I * H );
L[1] := 1.0 + 2.0 * VV; U[1] := -VV / L[1];
for I := 2 to M-2 do begin
L[I] := 1.0 + 2.0 * VV + VV * U[I-1]; U[I] := -VV / L[I]
end;
L[M-1] := 1.0 + 2.0 * VV + VV * U[M-2];
for J := 1 to N do begin
Z[1] := W[1] / L[1];
for I := 2 to M-1 do Z[I] := ( W[I] + VV * Z[I-1] ) / L[I];
W[M-1] := Z[M-1];
105
for I1 := 1 to M-2 do begin
I := M-2 - I1 + 1; W[I] := Z[I] - U[I] * W[I+1]
end;
end;
end;
106
Równanie hiperboliczne – Falowe
Równania różniczkowe cząstkowe hiperboliczne znane jako równanie falowe, w zależności
od jednego wymiaru oraz czasu przyjmuje postać:
(32)
dla oraz
z warunkami brzegowymi
dla
i początkowymi
oraz
dla
Równania tego typu są stosowane dla różnych problemów fizycznych, które są zależne od
czasu. Najczęściej stosowane są do obliczenia wartości fali płaskiej w dowolnym punkcie
pomiędzy 0 a l dla dowolnej chwili czasu większej od zera. Funkcje f(x) i g(x) określają
warunki początkowe położenia oraz prędkości rozchodzenia się fali.
Aby rozwiązać równanie cząstkowe paraboliczne (32) zastosujemy metodę różnic
skończonych MRS.
W tym celu weźmiemy liczbę całkowitą m i zdefiniujemy krok całkowania h = (b-a)/m.
Dodatkowo ustalimy wartość kroku czasowego k. Dzięki temu możemy wyznaczyć węzły
siatki (x
i
, t
j
), gdzie:
dla i = 0,1, .. m oraz dla j = 0,1, ..
Dla wewnętrznych punktów węzłowych (x
i
,t
j
), które nie należą do brzegu, zastosujemy
metodę różnic skończonych – MRS. Metoda ta opiera się na zastąpieniu pochodnych
cząstkowych rozwinięciem funkcji w szereg Taylora w otoczeniu punktu (x
i
,t
j
) – (patrz druga
pochodna numeryczna) . Wówczas otrzymujemy:
Dla zmiennej x
(33)
0
)
,
(
)
,
(
2
2
2
2
2
t
x
dx
u
d
t
x
dt
u
d
l
x
0
0
t
0
)
,
(
)
,
0
(
t
l
u
t
u
0
t
)
(
)
0
,
(
x
f
x
u
)
(
)
0
,
(
x
g
x
t
u
l
x
0
ih
x
i
jk
t
j
)
,
(
12
)
,
(
)
,
(
2
)
,
(
)
,
(
4
4
2
2
1
1
2
2
j
i
j
i
j
i
j
i
j
i
t
x
u
h
h
t
x
u
t
x
u
t
x
u
t
x
x
u
107
gdzie
Dla zmiennej t
(34)
gdzie
Podstawiając wzory (33) oraz (34) do równania falowego (31) oraz wyłączając oszacowanie
błędu zapisujemy:
(35)
warunek brzegowy :
(36)
dla
i = 1,2 … m-1 oraz j = 1,2, …
natomiast lokalny błąd odcięcia jest równy :
(37)
gdzie:
oraz
Ostatecznie wzór na MRS, podstawiając za ( x
i
, t
j
) = w
i,j
zapisujemy:
(38)
przekształcając powyższy wzór ze względu na w
i,j+1
,wówczas otrzymujemy:
gdzie:
(39)
dla i = 1,2 … m-1 oraz j = 1,2, …
z warunkami brzegowymi
dla j = 1,2, …
i początkowymi
dla i = 1,2, … m-1
)
,
(
1
1
i
i
i
x
x
)
,
(
12
)
,
(
)
,
(
2
)
,
(
)
,
(
4
4
2
2
1
1
2
2
j
i
j
i
j
i
j
i
j
i
x
t
u
k
k
t
x
u
t
x
u
t
x
u
t
x
y
u
)
,
(
1
1
j
j
j
t
t
0
)
,
(
)
,
(
2
)
,
(
)
,
(
)
,
(
2
)
,
(
2
1
1
2
2
1
1
h
t
x
u
t
x
u
t
x
u
k
t
x
u
t
x
u
t
x
u
j
i
j
i
j
i
j
i
j
i
j
i
0
)
,
(
)
,
0
(
j
j
t
l
u
t
u
)
,
(
)
,
(
12
1
4
4
2
2
4
4
2
j
i
j
i
ij
t
x
u
h
x
t
u
k
)
,
(
1
1
i
i
i
x
x
)
,
(
1
1
j
j
j
t
t
0
2
2
2
,
1
,
,
1
2
2
1
,
,
1
,
h
w
w
w
k
w
w
w
j
i
j
i
j
i
j
i
j
i
j
i
1
,
,
1
,
1
2
,
2
1
,
1
2
j
i
j
i
j
i
j
i
j
i
w
w
w
w
w
h
k
0
,
,
0
j
m
j
w
w
)
(
0
,
i
i
x
f
w
108
Analizując równanie (39) można zauważyć, że w celu wyznaczenia przybliżenia rozwiązania
w punkcie (x
i
,t
j+1
), potrzebne są wartości trzech punktów z poprzedniego kroku
czasowego t
j
oraz jednej wartości dla czasu t
j-1
– patrz rysunek poniżej.
Jak widać poszukujemy rozwiązania równania falowego dla wewnętrznych węzłów siatki, w
tym celu dla każdego punktu siatki układamy równanie (39). Obliczenia rozpoczynamy od
kroku czasowego j = 2, gdyż najpierw musimy wygenerować wartości początkowe dla t = 0
dla wszystkich punktów u( x
i
, 0)=f(x
i
), następnie wygenerujemy z drugiego warunku
początkowego wartości punktów dla t = t
1
jako du( x
i
, t
1
) / dt = g(x
i
).
A więc mamy do rozwiązania układ (m-1) równań, który powstał poprzez ułożenie równania
(39) dla każdego punktu w
i,j
siatki dla i = 1,2, .. m-1 :
(40)
co rozpisując daje nam:
(41)
Czyli w ten sposób możemy wyznaczyć poszukiwane wartości w
(j+1)
. Wystarczy tylko
przemnożyć macierz A przez wektor w
(j)
- wartości dla poprzedniego kroku, a następnie odjąć
od tego wektor wynikowy dla t
j-1
czyli w
(j-1)
.
)
1
(
)
(
)
1
(
j
j
j
w
Aw
w
1
,
1
1
,
2
1
,
1
,
1
,
2
,
1
2
2
2
2
2
2
1
,
1
1
,
2
1
,
1
)
1
(
2
0
0
0
0
0
0
)
1
(
2
j
m
j
j
j
m
j
j
j
m
j
j
w
w
w
w
w
w
w
w
w
109
Ponieważ obliczania rozpoczynamy od kroku czasowego j = 2 musimy oprócz warunków
początkowych dla j = 0 określić wartości dla j = 1, skorzystamy z drugiego warunku
początkowej prędkości fali:
dla
(42)
Zamienimy pochodną cząstkową po czasie wykorzystując wzór schematu jawnego dla
równań parabolicznych:
(43)
Wyprowadźmy u(x
i
,t
1
) z powyższego równania, wtedy otrzymujemy:
(44)
Podstawmy do równania (44) wartość pochodnej z równania (42) :
(45)
Podstawmy za u( x
i
,t
j
) = w
i,j
, pomijając błąd odcięcia, ostatecznie mamy:
dla i = 1,2, .. m-1
:
(46)
Powyższy wzór pozwala na określenie warunków początkowych dla j = 1, wadą tego jest
duży błąd oszacowania rzędu O(k). Alternatywą do tego jest wzór wyprowadzony jako
rozwiniecie funkcji u(x
i
,t
1
) w szereg Maclurina 2 rzędu:
dla i = 1,2, .. m-1
(47)
gdzie:
:
l
x
0
)
(
)
0
,
(
x
g
x
t
u
)
~
,
(
2
)
0
,
(
)
,
(
)
0
,
(
2
2
1
j
i
i
i
x
t
u
k
k
x
u
t
x
u
x
t
u
)
,
0
(
~
1
t
j
)
~
,
(
2
)
0
,
(
)
0
,
(
)
,
(
2
2
1
j
i
i
i
x
t
u
k
x
t
u
k
x
u
t
x
u
)
~
,
(
2
)
(
)
0
,
(
)
,
(
2
2
1
j
i
i
i
i
x
t
u
k
x
g
k
x
u
t
x
u
)
(
0
,
1
,
i
i
i
x
g
k
w
w
)
(
2
)
1
(
0
,
1
0
,
1
2
0
,
2
1
,
i
i
i
i
i
x
g
k
w
w
w
w
h
k
110
Równania hiperboliczne - Falowe - przykład 1
Rozwiązać poniższe równanie falowe:
dla
oraz
z warunkami brzegowymi
dla
i początkowymi
oraz
dla
Obliczenia wykonać dla czasu t = 1 wykorzystując metodę MRS dla równań falowych,
przyjąć: h = 0.1 oraz k = 0.05 ,przyjąć
2
= 4
Rozwiązanie dokładne ma postać:
Tabela wyników
0
)
,
(
)
,
(
2
2
2
2
2
t
x
dx
u
d
t
x
dt
u
d
0
)
,
1
(
)
,
0
(
t
u
t
u
1
0
x
0
t
0
t
)
sin(
)
0
,
(
x
x
u
1
0
x
0
)
0
,
(
x
t
u
)
2
cos(
)
sin(
)
,
(
t
x
t
x
u
0
0,8090169944
0,8090169944
0,3
0
0,3090169944
0,5877852523
0,8090169944
0,9510565163
1,0000000000
0,9510565163
0,5877852523
0,3090169944
0
u( x
i
, 1
)
5,421
10
-20
0
0
1,1102
10
-16
1,1102
10
-16
0
0
0
| u( x
i
, 1
)
– w
i,20
|
0
0,3090169944
0,5877852523
0,8090169944
0,9510565163
1,0000000000
0,9510565163
0,5877852523
0,3090169944
0
w
i,20
k = 0.05
0,5
0,6
0,7
0,8
0,9
1,0
0,4
0,2
0,1
0,0
x
i
111
Interpretacja graficzna
112
Równania hiperboliczne – kody źródłowe
Funkcja określająca warunki początkowe równania paraboliczne
function F( X,Y : real ) : real;
begin
F := sin( pi * x );
end;
function G( X,Y : real ) : real;
begin
G := 0;
end;
Definicja wymaganych stałych i typów
const
max = 100;
type
Vector = array [ 0..max ] of real;
Macierz = array [ 0..max ] of Vector;
Funkcja „WaveEquation” – oblicza przybliżenie rozwiązania równania falowego
function WaveEquation(FX,ALPHA ,FT : real; M,N : integer; var W : Macierz):
boolean ;
var H,K,V : real; M1,N1,I,J : integer;
begin
M1 := M + 1; N1 := N + 1; Result := True;
H := FX / M; K := FT / N;
V := ALPHA * K / H;
for J := 1 to N do begin
W[0,J] := 0; W[M,J] := 0;
end;
W[0,0] := F( 0 ); W[M,0] := F ( FX );
for I := 1 to M-1 do begin
W[I,0] := F( H * ( I ) );
W[I,1] := (1-V*V)*F(H*(I))+V*V*(F((I+1)*H) + F(H*(I-1)))/2+K*G(H*(I));
end;
for J := 1 to N-1 do for I := 1 to M-1 do
W[I,J+1] := 2*(1-V*V)*W[I,J]+V*V*(W[I+1,J]+W[I-1,J])-W[I,J-1];
end;