METODY NUMERYCZNE
W ELEKTROTECHNICE
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):
n
l = znak 2i , gdzie ei = 0 lub 1
ei
i=0
Jeżeli rejestr ma d-bitów, wówczas liczba całkowita n może zawierać się w przedziale
-2d-1 < n < 2d-1-1
przykład: Zapisać liczbę 18 w systemie dwójkowym:
l =1 24 + 0 23 + 0 22 +1 21 + 0 20, może być zatem przechowywana w słowie o długości D
+ 1 > 5 bitów jako
0 004243
{
1...0010010
4 4
znak
d cyfr rozwiniecia
Liczby rzeczywiste (zmiennopozycyjne):
x = zn m 2c , gdzie < m < 1, m-mantysa, i-cecha
Ą
m = 2-i
e-i
i=1
2
W maszynie cyfrowej mantysa zapisywana jest w t-bitach mt
t
mt = e-(t+2) 2-t + 2-i
e-i
i=1
w wyniku zaokrąglenia do t-bitów mantysy
1
m - mt Ł 2-t
2
Zmiennopozycyjna reprezentacja liczby rzeczywistej x oznaczana jest symbolem rd(x) i jest
równa
rd(x) = zn 2c mt
rd(x) - x
Ł 2-t
x
Gdy elementy macierzy są rzędu m 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:
rd(x) = x(1+ e), gdzie e Ł 2-t
liczbę 2-t nazywamy dokładnością maszynową
przykład: Liczba 18,5 daje się przedstawić w postaci:
2
x = 2(12 +021+120 )(1 2-1 + 0 2-2 + 0 2-3 +1 2-4 + 0 2-5 +1 2-6)
może być więc przechowywana w słowie o długości d + 1 > 11 bitów o t > 6 bitach mantysy:
0 0 000...001011001010...0
{ {
14243 14243
znak znakc d-t bitów t bitów
cechy mantysy
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. yródłem
3
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 ex
1 1
równego szeregowi: 1+ x + x2 + ... + xN + ... zastępuje je sumą częściową o
2 N!
1 1
odpowiednio dobranej wartości N: 1+ x + x2 +...+ xN
2 N!
Jeżeli liczba N będzie niedostatecznie duża, to uzyskana w ten sposób wartość liczby ex
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:
fl(x1 ą x2) = x1(1+ e1) ą x2(1+ e2)
fl(x1 x2) = x1 x2(1+ e3) e3 = e1 + e2
ć
x1 x1(1+ e4) x1 e1 e2
fl = = e4 = e5 =
x2 x2 x2(1+ e5) e2 e1
Ł ł
gdzie są co do modułu niewiększe niż dokładność maszyny . Przedrostek fl oznacza, że są
i
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 y = x , przekształcając do równania iteracyjnego. Można to
uczynić obliczając ciąg y0, y1, y2..., gdzie y0 jest dowolnie wybraną liczbą, natomiast każdy
element ciągu yi jest dany wzorem:
ć
1 x
yi = yi-1 + , dla i = 1, 2, 3, ...
2 yi-1
Ł ł
(niewiadoma musi być po obu stronach równania)
2
y = x
y2 = x
y y = x
x
y =1 , współczynnik przed funkcją (tutaj 1) musi być < 1 dlatego:
y
x
2y - y =
y
4
x
2y = + y : 2
y
ć
1 x
yi = + yi-1 - iteracja i-ta
2 yi-1
Ł ł
ć
7x 1 7x
np. y2 = 7x (2y - y) = yi = + yi-1 takiego równania nie da się rozwiązać
y 2 yi-1
Ł ł
bo 7/2
ale można zrobić tak:
ć
7x 1 7x
(8y - 7y) = yi = + 7yi-1
y 8 yi-1
Ł ł
jeśli założyć, że działania wykonywane są dokładnie, to ciąg y0, y1, y2, ... będzie zbieżny do
liczby x, ponieważ yi = pi x , to
ć
1 1
yi+1 = + x , a ciąg po > 0
pi
2 pi
Ł ł
ć
1 1
pi+1 = + , przy i = 1, 2, 3, ... jest zbieżny do 1
pi
2 pi
Ł ł
Przykład:
rozwiązanie iteracyjne
x + 7y = 3 8x - 7x + 7y = 3
4x - 2y = 5
- 2(5y - 4y) + 4x = 5
1
x = - 7yi + 3)
(7xi
i+1
8
wszystkie składniki przy niewiadomych są < 1 (7/8, -7/8, 3/5)
yi+1 = ć2xi + 4yi - 5 1
2 5
Ł ł
przyjmujemy dowolnie xi = 1, yi = 2
Działania potrzebne do wyznaczania kolejnych wartości ciągu y0, y1, y2, ... nazywamy
iteracją. Informację o tym czy w kolejnej iteracji nastąpiło przybliżenie do rozwiązania
uzyskujemy obliczając wyrażenie:
yi+1 - x
qi =
yi x
pi -1
W przypadku dokładnych działań, gdy yi = pi x i qi = , więc qi < 0 dla i = 1, 2, 3, ...
2pi
jeśli pi>1/3. Oznacza to, że maleje wartość bezwzględna różnicy y1 - x odpowiednio
dużego numeru iteracji i.
5
Dla obliczeń wykonywanych na maszynie cyfrowej, mamy:
ł
1 x(1+ ei )(1+ ei )ś , gdzie ei Ł 3 , k = 1, 2, 3, ...
i
2 3
fl(yi+1) = (1+ e1)+
i
ęy k
2 yi
y' - x
' '
i+1
Dla liczby yi+1 = fl(yi+1) uzyskamy qi =
yi - x
Jeżeli yi jest już dostatecznie dobrym przybliżeniem liczby x , a więc pi 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 pi.
Przykład: dla liczby 4 x=4
1 4 1 ć 4 1 ć 4
ć4
y0 = 1 y1 = + = 2,5 y2 = 2,5 + = 2,05 y3 = 2,05 + = 2,000609
2 4 2 2,5 2 2,05
Ł ł
Ł ł Ł ł
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 = ( x1, x2, ... xn ), na ich
podstawie chcemy obliczyć skończenie wiele wyników rzeczywistych y = ( y1, y2, ... ym ).
Będziemy więc chcieli określić wartość y według odwzorowania y = Ć (x), gdzie j: D Rm
jest ciągłym odwzorowaniem i D Rm.
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 x1, x2, ... xn
posługujemy się reprezentacjami rd (x1), rd (x2), ..., rd (xn). 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 zle
uwarunkowanym. Wielkości charakteryzujące wpływ zaburzeń danych na zaburzenia
rozwiązania nazywamy wskaznikami uwarunkowania zadania.
Przykład:
20
w(x) = a x20 + a19x19 + ... + a1x + a0 = - k)
20 (x
k=1
6
Zerami tego wielomianu będą liczby naturalne 1, 2, ..., 20. Gdy zakłócimy np. a19 = - 210 i
jego wartość wynosi a19() = - (210 + 2-23), czyli a19 () = a19 (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 (xi) xi dla i = 1, 2, ..., n
n
śjj(x)
dyi = j(rd(x))- j(x)- dxi dla j = 1, 2, ..., m
śxi
i=1
We wzorze tym, czynnikiem określającym wrażliwość yj na bezwzględną zmianę xi jest
śjj(x)
śxi
Analogiczny wzór można wyprowadzić dla przenoszenia się błędów względnych. Jeśli
yj ą 0dla j = 0,1, 2, ..., m i xi ą 0 dla i = 1, 2, ..., n to
n
xj śjj(x) xj śjj(x)
ey = ex , gdzie nazywamy wskaznikiem uwarunkowania.
j (x)
j i
śxi jj(x) śxi
i=1
j
Jeśli wskazniki uwarunkowania co do wartości bezwzględnej są duże to zadanie jest zle
uwarunkowane.
Przykład:
Badamy uwarunkowanie sumy y = a + b + c (a,b,c R)
a(1+ e1)+ b(1+ e2)+ c(1+ e3)- a - b - c
=
a + b + c
ae1 + be2 + ce3 ć a b c
Ł e + + , gdzie jest
a + b + c a + b + c a + b + c a + b + c
Ł ł
a b c
największą z wartości 1, 2 i 3, gdzie + + jest wskaznikiem
a + b + c a + b + c a + b + c
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.
7
INTERPOLACJA:
Dana jest funkcja y = f (x), x[x0, xn], dla której znana jest tablica wartości w punktach
zwanych węzłami interpolacji. Należy wyznaczyć taką funkcję W(x), aby:
W(x0) = Y0, W(x1) = Y1, ..., W(xn) = Yn
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
n
liniowej n+1 funkcji bazowych j0(x),j1(x),j2(x),...,jn (x) W (x) = ji (x)
ai
i=0
wyrażenie to nazywamy wielomianem uogólnionym.
Wprowadzając macierz bazową
a0
ł
ęa ś
1
ę ś
F = [j0(x),j1(x),...,jn (x)] i macierz współczynników A =
ę ś
M
ęa ś
n
mamy W(x) = F(x) A
warunek W(x0) = Y0, W(x1) = Y1, ..., W(xn) = Yn
można zapisać w postaci układu równań liniowych X A = Y, gdzie
8
j0(x0) j1(x0) L jn (x0) Y0
ł ł
ęj (x1) j1(x1) L jn (x1)ś ęY ś
0 1
ę ś ę ś
X = Y =
ę ś ę ś
M M M M M
ę ęY ś
(xn ) j1(xn ) L jn (xn )ś
j0 n
Jeśli macierz X nie jest osobliwa (da się odwrócić), to:
A = X-1Y
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
j0(x) =1,j1(x) = x,K,jn(x) = xn,
Baza dla funkcji ciągłych na odcinku skończonym [xo, xn] 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ć:
W(x) = a0 + a1x + a2x2 +K+ anxn dodatkowo musi być spełniony warunek:
a0 + a1x0 +K+ anx0n = y0
a0 + a1x1 +K+ anx1n = y1
K K K K K K
a0 + a1xn +K+ anxn n = yn
Powyższy układ równań ma jedyne rozwiązanie względem a1, jeżeli wartości x0, x1, ..., xn są
od siebie różne
1 x0 K x0n
1 x1 K x1n
det X =
M M M M
1 xn K xn n
Macierz X-1 dla bazy wielomianowej j0(x) =1, j1(x) = x, K, jn(x) = xn, nazywana jest
macierzą Lagrange a
Zauważyć należy, że każdy zbiór węzłów równoodległych xi+1 xi = h = const. można
x - x0
sprowadzić do zbioru podstawowego podstawiając q = , wówczas q[0,n], a macierz
h
Ś przyjmuje postać
F =[1, q, q2, K, qn]
9
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
(x0, y0), (x1, y1),K, (xi, yi),K, (xn, yn)
przyjmuje się funkcje bazowe w postaci
j0(x) = (x - x1)(x - x2)(x - x3)K(x - xn )
j1(x) = (x - x0)(x - x2)(x - x3)K(x - xn )
K
ji(x) = (x - x0)(x - x1)(x - x2)(x - x3)K(x - xi-1)(x - xi+1)K(x - xn )
K
jn (x) = (x - x0)(x - x1)(x - x2)(x - x3)K(x - xn -1)
Funkcje te są wielomianami stopnia n zbudowanymi w ten sposób, że w funkcji bazowej Ć1
brakuje czynnika (x-xi). Zatem wielomian interpolacji wyraża się następującym wzorem:
W(x) = a0(x - x1)(x - x2)K(x - xn) + a1(x - x0)(x - x2)K(x - xn) + ... + an(x - x0)(x - x1)K(x - xn-1)
współczynniki a0 ... an tego wielomianu wyznaczamy z równania:
X A = Y, przy czym macierz X ma postać:
(x0 - x1)...(x0 - xn) 0 L 0
ł
ę ś
0 (x1 - x0)...(x1 - xn) L 0
ę ś
X =
ę ś
M M M M
ę
0 0 L (xn - x0)...(xn - xn -1)ś
Macierz posiada tylko główną przekątną niezerową w związku z tym układ równań X A = Y
rozwiązuje się natychmiastowo
y0
a0 =
(x0 - x1)(x0 - x2)K(x0 - xn )
y1
a1 =
(x1 - x0)(x1 - x2)K(x1 - xn )
L
yn
an =
(xn - x0)(xn - x1)K(xn - xn -1)
10
Można więc wielomian interpolacyjny Lagrange a zapisać w postaci ułamka:
n
(x - x0)(x - x1)K(x - xi-1)(x - xi+1)K(x - xn )
W(x) =
yi
(xi - x0)(xi - x1)K(xi - xi-1)(xi - xi+1)K(xi - xn )
i=1
lub krócej
n
(x - xj)
W(x) =
yi(x - xj) , j = 0, 1, ..., n
i=1 jąi
i
Przykład:
Wyznaczyć wielomian interpolacyjny Lagrange a funkcji f (x) = ex w przedziale [0,2 ; 0,5]
mając dane:
f (0,2) = 1,2214, f (0,4) = 1,4918, f (0,5) = 1,6487
(x - 0,4)(x - 0,5) (x - 0,2)(x - 0,5) (x - 0,2)(x - 0,4)
W (x) = 1,2214 +1,4918 +1,6487 =
- 0,2 (-0,3) 0,2 (-0,1) 0,3 0,1
= 0,724x2 + 0,91x +1,009
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 = xi+1 xi można wyprowadzić pojęcie
różnicy skończonej rzędu k
Dyi = yi+1 - yi
D2yi = D[Dyi]= Dyi+1 - Dyi = yi+2 - 2yi+1 + yi
...
k
j
k
Dkyi = D[Dk-1yi]= Dk-1yi+1 - Dk-1yi = ( )yi+k-1
(-1) j
j=0
11
Przykład: Dla funkcji danej w postaci tablicy zbudować tablicę różnic skończonych.
NR x y "y "2y "3y "4y "5y
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 x0, x1 = x0 + h, x2 = x0 + 2h, ..., xn = x0 + nh
dane są wartości funkcji f(x0), f(x1), ..., f(xn). Szukany wielomian interpolacyjny zapisać
możemy w następujący sposób:
W(x) = a0 + a1q + a2q(q 1) + a3q(q 1)(q 2) + ... + anq(q 1)(q 2)...(q n +1),
x - x0
gdzie q =
h
Dla:
x = x0: q = 0
x = x1: q = 1
x = x2: q = 2
...
x = xn: 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ń:
1 0 0 0 L 0 a0 y0
ł ł ł
ę1 1 0
0 L 0ś ęa1 ś ęy1 ś
ę ś ę ś ę ś
ę ś ę ś ę ś
1 2 2 0 L 0 a2 y2
=
ę1 3 6
6 L 0ś ęa3 ś ęy3 ś
ę ś ę ś ę ś
ęM M M ś ęLś ęLś
M L M
ę ś ę ś ę ś
ę ś
1 n n(n -1) n(n -1)(n - 2) L n! ę ś ę ś
an yn
z którego wyznacza się wartości nieznanych współczynników a0, a1, a2, ..., an
12
a0 = y0
a0 + a1 = y1 a1 = Dy0
D2y0
a0 + 2a1 + 2a2 = y2 a2 =
2!
D3y0
a0 + 3a1 + 6a2 + 6a3 = y3 a3 =
3!
... ...
Dny0
a0 + na1 + n(n -1)a2 + ... + n!an = yn an =
n!
Ostatecznie otrzymujemy I wzór interpolacyjny Newtona w postaci
q(q -1) q(q -1)...(q - n +1)
W(x) = y0 + qDy0 + D2y0 + ... + Dny0
2! n!
Przykład: Dla zależności f(T)=T log p znalezć 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
T - 0,8
przyjmijmy T0 = 0,8 -q = = 5T - 4
0,2
W(T) = 0,3566 + 0,5817q + 0,0189q(q -1) - 0,0007q(q -1)(q - 2)
W(T) = -0,075T3 + 0,7225T2 +1,791T -1,5002
Zauważmy, że I wzór Newtona daje dobrą dokładność w pobliżu punktu T0, 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ąć T0 = 1,2, a wówczas q = 0,5 i
znajdujemy W(1,3) =1,5598 + 0,5 0,6571- 0,25 0,0319 =1,88437 .
13
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) = a0 + a1q + a2q(q + 1) + a3q(q + 1)(q + 2) + ... + anq(q + 1)(q + 2)...(q + n 1),
x - xn
gdzie q =
h
Współczynniki a0, a1, a2, ..., an wyznaczamy jak poprzednio i dochodzimy do wzoru:
q(q +1) q(q +1)...(q + n -1)
W (x) = y0 + qDyn-1 + D2 yn-2 + ... + Dn y0
2! n!
Przykład: Obliczyć T log p dla T = 1,65 z dokładnością do "2
T -1,8
Tn =1,8 q = = -0,75
0,2
W(1,65) = 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
1
F = ,sin(x),cos(x),sin(2x),cos(2x),...,sin(nx),cos(nx)ł
ę ś
2
Otrzymujemy zatem wielomian interpolacyjny w postaci
a0
W(x) = + b1 sin(x) + a1 cos(x) + b2 sin(2x) + a2 cos(2x) +...+ bn sin(nx) + an cos(nx)
2
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 xi[0,2Ą] dobranych według następującej
zależności
14
2ip
xi = , gdzie i = 0, 1, ..., 2n
2n +1
Czyli
2p 4np
x0 = 0, x1 = , ..., xn =
2n +1 2n +1
Warunek interpolacji prowadzi do układu równań liniowych w postaci
1
ł ł ł
0 1 ... 0 1 a0 y0
ę ś ę ś ę ś
2
ę ś ę ś ę ś
1
ę
sin(x1) cos(x1) ... sin(nx1) cos(nx1)ś ę b1 ś ę y1 ś
=
ę ś ę ś ę ś
2
ę ś ę ś ę ś
... ... ... ... ... ... ... ...
ę ś ę ś
1
sin(x2n) cos(x2n) ... sin(nx1) cos(nx1)ś ę an ś ę yn ś
ę ś ę
2
Współczynniki pierwszego wiersza macierzy X wynikają z wartości funkcji sin(hx) i cos (hx)
dla x0 = 0. Przedstawiony układ równań rozwiązuje się natychmiastowo, ponieważ macierz X-
1
można wyznaczyć z zależności
2
X-1 = XT
2n +1
Przykład: Daną funkcję f(x) = 7x x2 przybliżyć wielomianem trygonometrycznym
przyjmując n = 3. Współrzędne węzłów interpolacji obliczamy ze wzoru:
2iĄ
xi =
2n +1
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:
1
F = ,sin(x),cos(x),sin(2x),cos(2x),sin(3x),cos(3x)ł
ę ś
2
15
Tworzymy macierz X:
0,707 0 1 0 1 0 1
ł
ę0,707 0,782 0,623 0,975 - 0,223 0,434 - 0,901ś
ę ś
ę ś
0,707 0,975 - 0,223 - 0,434 - 0,901 - 0,782 0,623
ę0,707 0,434 - 0,901 - 0,782 0,623 0,975 - 0,223ś
ę ś
ę0,707 - 0,434 - 0,901 0,782 0,623 - 0,975 - 0,223ś
ę ś
ę0,707 - 0,975 - 0,223 0,434 - 0,901 0,782 0,623 ś
ę0,707 - 0,782 0,623 - 0,975 - 0,223 - 0,434 - 0,901ś
Elementy macierzy X mnożymy przez 2/7, otrzymaną macierz transponujemy i obliczamy
współczynniki wzoru interpolacyjnego
T
A = [11,845, -1,336, - 4,923, - 0,513, -1,961, - 0,147, -1,491]
16
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
Qm (x) = jj(x)
a j
j=0
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, a0, a1, ..., am.
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 x0,
x1, ..., xn 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
2
S = (x) - Q(x)}
{f
i=0
- dla aproksymacji średniokwadratowej integralnej lub przedziałowej
b
2
S = (x) - Q(x)} dx
{f
a
- dla aproksymacji jednostajnej
S = f (x) - Q(x)
sup
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 a0, a1, ..., am) aby zdefiniowane wyżej błędy były
minimalne.
17
Aproksymacja średniokwadratowa
Niech dana będzie funkcja y=f(x), która w pewnym zbiorze X punktów x0, x1, ..., xn
przyjmuje wartości y0, y1, ..., yn. 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 jj(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=(xj). Funkcja przybliżająca ma postać
m
Q(x) = ji (x)
ai
i=0
Przy czym współczynniki ai są tak określone, aby wyrażenie
n
2
f (x) - Q(x) = )[f (xi ) - Q(xi )] w(xi ) ł 0 dla i=0, 1, ..., n
w(xi
i=0
było minimalne. Funkcja w(x) jest z góry ustaloną funkcją wagową.
Aby wyznaczyć współczynniki ai oznaczamy odchylenie
2
n m n
ł
2
H(a0,a1,...,an ) = )ęf (x ) - ji (xi )ś = )R
w(x j j ai w(x j j
j=0 i=0 j=0
2
n m n
ł
H(a0,a1,...,an )= )ęf(x ) - ji (x )ś = )R2
w(xj j ai j w(xj j
j=0 i=0 j=0
gdzie Rj jest odchyleniem w punkcie xj.
Następnie obliczamy pochodne cząstkowe funkcji H względem ai. Z warunku
śH
= 0 , gdzie k = 0, 1, ..., n
śak
otrzymujemy układ m+1 równań o niewiadomych ai zwany układem normalnym
n m
śH ł
= -2 ji (xi )śjk (x ) = 0 , gdzie k = 0, 1, ..., n
w(x )ęf (x ) -ai j
śak j=0 j j i=0
n m
śH ł
= -2 )ęf(x ) - ji (x )śjk (x ) = 0
w(xj ai
śak j=0 j i=0 j j
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ć
DTDA = DTF
18
gdzie
j0x0 ... jmx0 a0 f (x0 )
ł ł ł
ęj x1 ... jmx1 ś ę ś ęf (x1)ś
a1
0
ę ś ę ś ę ś
D = A = f =
ę ś ę ś ę ś
... ... ... ... ...
ę ęa ś ęf (a )ś
xn ... jmxn ś
j0 m n
Aproksymacja wielomianowa
Jeżeli za funkcje bazowe przyjmiemy ciąg jednomianów (xi), i = 0, 1, ..., n to układ normalny
przyjmuje postać
n m
ł
) - (xij)śxk = 0 , k = 0, 1, ..., n
f (x j ai j
ę
j=0 i=0
co po zmianie kolejności sumowania daje
n m n
ć
i+k
f (xj )xk =
j ai xj
j=0 i=0 j=0
Ł ł
przykład: Odnalezć 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 2
n n
(yi - yi) (xi - x )
i
i=1 i=1
przy czym ax + byi = 1
i
n n
ć a 1
na0 + a1 = , gdzie: a1 = - , ao =
xi yi
b b
Ł i=1 ł i=1
n n n
ć ć
2
a0 + a1 = yi
xi xi xi
Ł i=1 ł Ł i=1 ł i=1
n n n n
n n n
yi 2 - yi
n yi -
xi xixi
xi yixi
i=1 i=1 i=1 i=1
i=1 i=1 i=1
a0 =
a1 =
2
2
n n
n n
ć
2 ć
2
n -
n -
xi xi
xi xi
i=1 Ł i=1 ł i=1 Ł i=1 ł
19
Układ normalny dla danych z tablicy ma rozwiązanie
a0 = 6/11 i a1 = 7/11,
więc ma postać
11y 7x = 6
Minimalizując w ten sam sposób wyrażenie
2
n
(xi - x )
i
i=1
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 x0, x1, ..., xn
jeśli
n
f (xi ) - g(xi )
i=0
20
przy czym
n n
2 2
[f (x)] > 0 [g(x)] > 0
i=0 i=0
Analogicznie ciąg funkcyjny
(jm(x)) j0(x),j1(x),...,jm(x),...
nazywamy ortogonalnym na zbiorze punktów x0, x1, ..., xn jeśli
n
(xi )jk (xi ) = 0 dla j`"k
jj
i=0
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
2
a = (xi )
jj jj
i=0
Załóżmy, że znamy n+1 równoodległych punktów xi (xi = x0+ih, i = 0, 1, ..., n). Za pomocą
przekształcenia liniowego
x - x0
q =
h
przeprowadzimy te punkty w kolejne liczby całkowite od 0 do n poszukujemy ciągu
wielomianów
(Fi(n)(q)) F0(n),F1(n)(q),...,Fm(n)(q)
(dolny indeks oznacza stopień i m d" n) spełniających warunek ortogonalności
n
(n) (n)
(i)Fk (i) = 0 dla j`"k
Fj
i=0
przy czym
Fk(n)(q) = a0 + a1q[1] + a2q[2] + ... + akq[k]
gdzie
q[k] = q(q -1)...(q - k +1)
21
Często używa się unormowanego ciągu wielomianów spełniających warunek
Ł
(n)
F (0) = 1, gdzie k = 0, 1, ..., m
k
Ł
(n)
F (q) =1+ b1q[1] + ... + bkq[k]
k
Co po przekształceniu daje nam wzór na wielomiany Grama, zwane też wielomianami
Czebyszewa stopni k = 0, 1, ..., m w postaci
k
Ł
q[s]
(n)
s k k+s
F (q) = ( )( )n , gdzie k = 0, 1, ..., m
k
-1 s s
[n]
s=0
Wzór aproksymacyjny oparty na wielomianach Grama ma postać:
m
ck Ł (n) m k Ł (n) x - x0
ym = F (q) = F
c ć h
sk sk Ł
ł
k=0 k=0
gdzie,
m Ł n
oraz
(n)
n n
Ł
2
sk = [Fk(n)(q)] ck = F (xi )
yi
q=0 i=0
Przykład:
Na podstawie tablicy wartości f(x) znalezć 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 Fk(q), q=1, 2,& , 4.
Obliczamy ze wzoru:
(n) ]
k
Ł
k k + s
ć ć
s
F (q) =
k
(-1) s s q[S
h[n]
S =0
Ł łŁ ł
F0 (q) = 1
q
F1(q) = 1- 2 n ł 1
n
22
q q(q -1)
F2(q) =1- 6 + 6 n ł 2
n n(n -1)
q q(q -1) q(q -1)(q - 2)
F2(q) =1-12 + 30 - 20 n ł 3
n n(n -1) n(n -1)(n - 2)
q q(q -1) q(q -1)(q - 2) q(q -1)(q - 2)(q - 3)
F2(q) =1- 20 + 90 -140 + 70 n ł 4
n n(n -1) n(n -1)(n - 2) n(n -1)(n - 2)(n - 3)
Zestawienie wyników przedstawia tabela:
q xi yi F0(q) F1(q) F2(q) F3(q) F4(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)
n n
Ł
2
sk = [Fk (n) (q)] ck = F (xi)
y
i
q=0 i=0
Obliczamy ci i si oraz bi= ci / si i ze wzoru:
m
Qm (x) = j (x)
a j j
j=0
otrzymujemy:
ci 37,5 -12,5 1,75 0 0
si 5 2,5 3,5 10 70
bi 7,5 -5 0,5 0 2
x -1
q =
dla
0,5
23
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:
pi
xi =
dla i = 0, 1, & , 2L-1
L
Baza jest ortogonalna nie tylko na przedziale <0, 2Ą>, ale też na zbiorze punktów xi, przy
czym warunki ortogonalności mają postać:
0 dla m ą k
2L-1
sin(mx )sin(kxi ) = L dla m = k ą 0
i
i=0
0 dla m = k = 0
0 dla m ą k
2L-1
cos(mx )cos(kxi ) = L dla m = k ą 0
i
i=0
2L dla m = k = 0
2L-1
)sin(kxi ) = 0
cos(mxi
i=0
dla m, k dowolnych, przy czym m i k zmieniają się od 0 do L
Przybliżeniem funkcji f(x) na zbiorze punktów xi jest wielomian trygonometryczny:
n
1
yn (x) = a0 + cos(jx) + bj sin( jx)) , n(aj
2
j=1
Współczynniki aj i bj wielomianu wyznaczamy tak, aby suma kwadratów różnic
2L-1
2
była minimalna.
[f (xi) - yn(xi)]
i=0
24
Korzystając z warunku ortogonalności otrzymujemy rozwiązanie układu normalnego
n m
śH ł
= -2 )ęf (xj) - ji(xj)śjk (xj) = 0
w(xj ai
śak j=0
i=0
w postaci:
2L-1 2L-1
1 1 pij
a = (xi)cos jxi = (xi)cos
j f f
L L L
i=0 i=0
dla j=1, 2, & , n
2L-1 2L-1
1 1 pij
bj = (xi)sin jxi = (xi)sin
f f
L L L
i=0 i=0
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;
25
Szybka transformata Fouriera
2Ą
Każdą funkcję w tablicy tablicą wartości w n punktach x = , gdzie b = 0, 1, ..., n-1
n
można przybliżać wielomianem trygonometrycznym w postaci
n-1
f (x) = exp(jkx), j - czynnik urojony
ck
k=0
lub
a0 m 1
f (x) = + cos(kx) + bk sin(kx))+ Qam+1 cos[(m +1)x]
(ak
2 2
k=0
gdzie
n - 2
Q = 1 m = dla n parzystych
2
n -1
Q = 1 m = dla n nieparzystych
2
n-1
1
ak = f (xb ) cos(kxb )
2n
b =0
n-1
1
bk = f (xb )sin(kxb ) dla k = 0, 1, ..., n -1
2n
b =0
n-1
1
ck = f (xb ) exp(- jkxb )
n
b =0
Aproksymacja za pomocą funkcji sklejonych stopnia 3 (funkcji giętych)
Każdą funkcję s(x) S3(Dn ) można przedstawić w postaci kombinacyjnej liniowej
n +1
s(x) = F3(x) , a d" x d" b oraz gdzie F3 są funkcjami określonymi wzorem:
ci i i
i=-1
(x - xi-2)3 dla x xi-2,xi-1
h3 + 3h2(x - xi-1) + 3h(x - xi-1)2 - 3(x - xi-1)3 dla x xi-1,xi
1
F3 =
h3 + 3h2(xi+1 - x) + 3h(xi+1 - x)2 - 3(xi+1 - x)3 dla x xi,xi+1
i
n3
(x - x)3
dla x xi+1,xi+2
i+2
dla pozostaych
0
b - a
gdzie xi = a + ih , h =
n
Aproksymacja na przedziale
Współczynniki ci dobieramy tak, aby odchylenie kwadratowe dane poniższym wzorem było
jak najmniejsze:
26
2
b
n +1
ł
I = (x) - f3(x)ś dx
Ci i
ęf
i=-1
a
w celu wyznaczenia minimum funkcji I=I(c-1 , c0 , & , cn+1) obliczamy pochodne:
śI
j = -1 , 0, & , n+1
ścj
i przyrównujemy je do zera. Otrzymane równania tworzą układ n+3 równań z n+3
niewiadomymi ci
b b
n +1
3
Ci i j j
F (x)F3(x)dx = f (x)F3(x)dx j= -1, 0, 1, & ,n+1
i=-1
a a
Funkcje F3(x) są liniowo niezależne na przedziale , więc układ równań ma
i
rozwiązanie w postaci punktu. Funkcji I osiąga minimum. Układ równań możemy zapisać w
postaci:
b
n+1
1
ci = (x)F3(x)dx j = -1, 0, 1, & , n+1,
aij j
f
h
i=-1
a
b
1
3
gdzie aij = (x)F3(x)dx
i j
F
h
a
Aproksymacja funkcji określonej na dyskretnym zbiorze punktów
Niech (xi ) dla i= 0, 1, & ,n1 ; n1 > n + 3 będzie zbiorem punktów, w których dane są wartości
funkcji f(x). Szukane jest minimum wyrażenia:
2
n1
n =1
ł
J = (x) - F3(xk ')ś
f ci i
ę
k=0 i=-1
W celu wyznaczenia ci obliczamy pochodne cząstkowe funkcji J = J (c-1, c0,& ,cn+1)
względem zmiennych ci i przyrównujemy je do zera.
Otrzymane równania tworzą układ:
n1
n +1
ci =
bij f (xk ')F3(xk ') , gdzie
i
i=-1 k =0
n1
3
bij = (xk ')F3(xk ')
Fi j
k =0
27
Jeżeli wyznacznik tego układu jest różny od zera to rozwiązując go możemy wyznaczyć
współczynniki cj 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
a11x1 + a12x2 + ... + a1nxn = b1
a21x1 + a22x2 + ... + a xn = b2
2n
... ... ... ... ... ... ...
an1x1 + a x2 + ... + a xn = bn
n2 nn
Układ ten można zapisać także w postaci macierzowej A X = B, gdzie
a11 a12 ... a1n b1 x1
ł ł ł
ęa a22 ... a ś ęb ś ęx ś
21 2n 2 2
ę ś ę ś ę ś
A = B = X =
ę ś ę ś ę ś
... ... ... ... ... ...
ęa a ... a ś ęb ś ęx ś
n1 n2 nn n n
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.
a11x1 = b1
a22x2 = b2
L L
annxn = bn
bi
Niewiadome można wtedy obliczyć xi = , aii `" 0, i = 0, 1, ..., n.
aii
28
Aatwo rozwiązuje się trójkątne układy równań
a11x1 + a12x2 + ... + a1n xn = b1
a x2 + ... + a xn = b2
22 2n
+ ... + ... = ...
a xn = bn
nn
Z ostatniego równania możemy wyznaczyć xn, z przedostatniego xn-1
n
bi - xs
ais
bn
s=i+1
xn = a ogólnie xi = , i = n-1, n-2, ..., 1 przy założeniu, że aii `" 0, i = 1, 2,
a aii
nn
..., n
PRZYKAAD:
1 1 1 x1 6
ł ł ł
ę0 2 -1ś ęx ś ę1ś
=
2
ę ś ę ś ę ś
ę
0 0 2 ś ę ś ę
x3 6ś
zmienną x3 mamy daną wprost - x3 = 3
b2 - a23x3 b1 - a12x2 - a13x3
x2 = = 2, x1 = = 1
a22 a11
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ń
a11x1 + a12x2 + ... + a1n xn = b1
a x1 + a22x2 + ... + a xn = b2
21 2n
... + ... + ... + ... = ...
an1x1 + an2x2 + ... + annxn = bn
29
a11 a12 ... a1n a11 a12 ... b1
a a ... a a a ... b2
21 22 2n 21 22
W = det A Wj =
... ... ... ... ... ... ... ...
a a ... a a a ... bn
n1 n2 nn n1 n 2
Wj
to można wykazać, że x = , W `" 0, j = 1, 2, ..., n
j
W
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 aij macierzy głównej A, natomiast kolumnę n+1 tworzą dowolne
wyrazy bi.
Ci1
Zakładając, że C11ą0, odejmujemy pierwsze równanie pomnożone przez od i-tego
C11
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 C1.
C11x1 + C12x2 + C13x3 + ... + C1n xn = C1,n+1
C22x2 + C23x3 + ... + C2n xn = C2,n+1
(1 (1 ( ( )
C32) x2 + C33) x3 + ... + C31) xn = C31n+1
n ,
( ( (1 ( )
Cn1) x2 + Cn1) x3 + ... + Cnn) xn = Cn1n+1
2 3 ,
Za pomocą wzorów określających nowe współczynniki
Ci1
(
Cij1) = Cij - C1 j i=2,3,...,n j=2,3,...,n+1
C11
1
Ci(2)
(1
Jeśli C22) ą 0, to odejmiemy drugie równanie układu pomnożone przez od i-tego
(1
C22)
równania układu i=3,4,...,n. W wyniku otrzymujemy układ:
C11 C12 C13 ... C1n C1,n+1
ł
ę (1 (1 ( ( )
0 C22) C23) ... C21) C21n+1ś
n ,
ę ś
(2 (2 (
ę ś
C2 = 0 0 C33) ... C3n) C3,2)+1
n
ę ś
: : : ... : :
ę ś
( (2 ( )
ę
0 0 Cn2) ... Cnn) Cn2n+1ś
3 ,
A nowe współczynniki wyrażone są wzorami
1
Ci(2) (
( (
Cij2) = Cij1) - C21j) i=3, 4, ..., n j=3, 4, ..., n
(1
C22)
Kontynuując takie postępowanie po wykonaniu n kroków dochodzimy do układu trójkątnego
30
C11x1 + C12x2 + C13x3 + ... + C1n xn = C1,n+1
(1 (1 ( (
C22) x2 + C23) x3 + ... + C21) xn = C1,1)+1
n n
(2 (2 (
C33) x3 + ... + C3n) xn = C1,2)+1
n
M
(n ( )
Cnn-1) xn = Cnnn-11
, +
C11 C12 C13 L C1n C1,n+1
ł
ę (1 (1 ( ( )
0 C22) C23) L C21) C21n+1 ś
n ,
ę ś
(2 (2 (
ę ś
C2 = 0 0 C33) L C3n) C3,2)+1
n
ę ś
M M M M M M
ę ś
(n ( )
ę
0 0 0 L Cnn-1) Cnnn-11 ś
, +
Algorytm dochodzenia od układu trójkątnego do układu równań liniowych
s = 1,2,...,n -1
i = s +1, s + 2,...,n
(
( ( Ciss-1) (
Cijs) = Cijs-1) - Csjs-1)
(s-1)
Css
Przykład
2x1 + 3x2 - 4x3 = -1 2 3 - 4 -1
ł
ę ś
3x1 - 2x2 + x3 = 5 ą C = 3 - 2 1 5
ę ś
ę ś
- x1 + 3x2 + 2x3 = -4 -1 3 2 - 4
Obliczamy elementy macierzy C1
C21 13
(1
C22) = C22 - C12 = -
C11 2
C21
(1
C23) = C23 - C13 = 7
C11
C21 13
(1
C24) = C24 - C14 =
C11 2
C31 9
(1
C32) = C32 - C12 =
C11 2
C31
(1
C33) = C33 - C13 = 0
C11
C31 9
(1
C34) = C34 - C14 = -
C11 2
31
druki krok s=2
ł
(1
C32) (1 63
(2 (1
ę2 3 - 4 -1ś
C33) = C33) - C23) = -
(1
ę
C22) 13
13 13ś
ą C2 = 0 - 7
ę ś
(1
2 2
C32) (1
ę ś
(2 (1
C34) = C34) - C24) = 0
(1
ę0 0 63 0 ś
C22)
ę ś
13
co daje nam macierz C1:
ł
ę2 3 - 4 -1ś
ę ś
13 13
C1 = 0 - 7
ę ś
2 2
ę ś
ę0 9 0 - 9ś
ę
2 2ś
i dalej ze wzorów dla układów trójkątnych mamy x3=0, x2=-1, x1=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
Xo, zwanej przybliżeniem początkowym, a następnie wykonaniu określonych operacji
arytmetycznych dających w wyniku X1 (zwanym pierwszym przybliżeniem). Kolejne
przybliżenia oblicza się wykonując te same operacje na ostatnio obliczonym przybliżeniu.
32
Jeżeli ciąg przybliżeń {Xn} 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...
X1 w11 w12 L w1n
ł ł ł
ęX ś ęw w22 L w2n ś ę ś
2 21
ę ś ę ś ę ś
=
ę ś ę ś ę ś
M M M L M
ęX ś ęw wn2 L wnn ś ę ś
n n1
Przykład: Sprowadzić układ do postaci AX=B
2x1+x2+x3=4 x1=-0,5x2-0,5x3+2
x1+x2-x3=1 ą x2=-x1+x3+1
x1+x2+x3=3 x3=-x1-2x2+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.
2 1 1 2 0 0 0 1 1
ł ł ł
ę1 1 -1ś ę0 1 0ś + ę1 0 -1ś
=
ę ś ę ś ę ś
ę
1 2 1 ś ę 0 1 ę 2 0 ś
0 ś 1
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:
-1 -1
x1 2 0 0 0 1 1 x1 - 2 0 0 4
ł ł ł ł ł ł
ęx ś ę1 ęx ś ę ę1ś
= -ę0 1 0ś 0 -1ś + 0 1 0ś
2 2
ę ś ę ś ę ś ę ś ę ś ę ś
ę ś ę ś ś
x3 0 0 1 ę 2 0 ś ę ś ę 0 0 1 ę
1 x3 4ś
ponieważ
-1
2 0 0 0,5 0 0
ł ł
ę0 1 0ś ę
= 0 1 0ś
ę ś ę ś
ę ś ś
0 0 1 ę 0 0 1
więc jak łatwo sprawdzić mamy taki sam układ jaki znaleziono sposobem naturalnym
33
Metoda iteracji prostej (Jakobiego):
Metoda ta dla równania X=W*X+Z polega na przyjęciu zerowego przybliżenia wektora
X=Xo, a następnie przeprowadzenia obliczeń iteracyjnych za pomocą zależności:
Xi+1=W*Xi+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 Xo. 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=50C, wymiary prostokątne a=b=0,8 cm.
50 100 100 100 100 100 100 100 100
ł
ę
0 50 63 63 63 63 63 75 100ś
ę ś
ę ś
0 38 50 50 50 50 50 63 100
ę ś
0 38 50 50 50 50 50 63 100ś
ę
ę
0 38 50 50 50 50 50 63 100ś
ę ś
0 38 50 50 50 50 50 63 100ś
ę
ę
0 38 50 50 50 50 50 63 100ś
ę ś
0 25 38 38 38 38 38 50 100ś
ę
ę ś
0 0 0 0 0 0 0 0 50
50 100 100 100 100 100 100 100 100
ł
ę
0 50 69 78 82 86 90 94 100ś
ę ś
ę ś
0 31 50 61 68
ę ś
0 22 39 50 58
ę ś
ę ś
0 18 32 42 50
ę ś
0 14 26 35 42
ę ś
ę ś
0 10 19 26 32
ę ś
0 6 10 14 18
ę ś
ę ś
0 0 0 0 0
34
Metoda Gaussa-Seidla
Stanowi ulepszenie metody iteracyjnej prostej polegające na wykorzystaniu obliczonych k
pierwszych składowych wektora Xi+1 do obliczenia składowej k+1. Modyfikacja ta znacznie
przyspiesza proces obliczania.
Przykład: Rozwiązując układ równań:
x1 0 0,2 0,1 0,2 x1 2
ł ł ł ł
ęx ś ę0,1 0 0,3 0,2ś ęx ś ę5ś
2 2
ę ś ę ś ę ś ę ś
= +
ę ś ę ś ę ś ę ś
x3 0,3 0,2 0 0,1 x3 1
ęx ś ę0,2 - 0,1 0,3 0 ś ęx ś ę2ś
4 4
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:
u1 = 0,2x2 + 0,1x3 + 0,2x4 + 2
u2 = 0,1x1 + 0,3x3 + 0,2x4 + 5
u3 = 0,3x1 + 0,2x2 + 0,1x4 +1
u4 = 0,2x1 - 0,1x2 + 0,3x3 + 2
Przyjmujemy teraz punkt startowy np. X0 = [3, 6, 4, 3]T, składowe tego wektora podstawiamy
do prawych stron równań wyliczając u1, ..., u4, które w następnej iteracji przejmą rolę punktu
startowego. Kolejne przybliżenia otrzymane tą metodą
3 4,2 4,4 4,45 4,53
ł ł ł ł ł
ę6ś ę7,1ś ę7,08ś ę7,27ś ę7,33ś
ę ś ę ś ę ś ę ś ę ś
ę ś ę ś ę ś ę ś ę ś
4 3,4 4 4,05 4,13
ę3ś ę3,2ś ę3,15ś ę3,37ś ę3,78ś
rozwiązaniem dokładnym jest X = [4,59, 7,49, 4,20, 3,44]
Ten sam układ w metodzie Gaussa Siedla ma postać:
u1 = 0,2x2 + 0,1x3 + 0,2x4 + 2
u2 = 0,1u1 + 0,3x3 + 0,2x4 + 5
u3 = 0,3u1 + 0,2u2 + 0,1x4 +1
u4 = 0,2u1 - 0,1u2 + 0,3u3 + 2
35
a kolejne przybliżenia:
3 4,2 4,51 4,56 4,58
ł ł ł ł ł
ę6ś ę7,22ś ę7,32ś ę7,38ś ę ś
7,4
ę ś ę ś ę ś ę ś ę ś
ę ś ę ś ę ś ę ś ę ś
4 4 4,15 4,19 4,2
ę3ś ę3,32ś ę3,41ś ę3,43ś ę3,44ś
Formalny zapis algorytmu
Macierz W w równaniu X = W X + Z przedstawiamy jako sumę dwóch macierzy: górnej
trójkątnej Wu i dolnej trójkątnej Wl. 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:
i+1 i i+1
X = Wu X +Wi X + Z
czyli mamy:
u1 0 0,2 0,1 0,2 x1 0 0 0 0 2
ł ł ł ł ł
ęu ś ę0 0 0,3 0,2ś ęx ś ę0,1 0 0 0ś ę5ś
2 2
ę ś ę ś ę ś ę ś ę ś
= + +
ę ś ę ś ę ś ę ś ę ś
u3 0 0 0 0,1 x3 0,3 0,2 0 0 1
ęu ś ę0 0 0 0 ś ęx ś ę0,2 - 0,1 0,3 0ś ę2ś
4 4
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 ui po
prawej stronie układu wartości
xi + ą (ui xi) ą współczynnik nadrelaksacji, ą<1,2>
i liczba niewiadomych a nie iteracja
Sposób dobory optymalnego współczynnika a 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 a=1 otrzymuje się metodę Gaussa Siedla. Według zródeł literaturowych
zaleca się przyjmowanie a rzędu 1,8.
Algorytm metody nadrelaksacyjnej sprowadza się do wzoru:
i+1 i i+1 i+1 i
X = Wu X +Wi[X +a(X - X )] + Z
36
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 l taka, że:
Ax = lx
Możemy wtedy zapisać równanie charakterystyczne
(A - lI)x = 0, gdzie I macierz jednostkowa
oraz wielomian charakterystyczny ma postać
p(l) = det (A - lI)
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 l jest
pierwiastkiem wielomianu p wówczas:
det(A - lI) = 0
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 r(A) macierzy A jest zdefiniowany jako:
r(A) = max l , gdzie l - oznacza moduł wektora l
wówczas równanie charakterystyczne jest zbieżne gdy:
r(A) < 1
Można wykazać, że:
r(A) < A gdzie A oznacza normę macierzy A
ostatecznie dla metod iteracyjnych można zapisać twierdzenie:
Warunkiem koniecznym i wystarczającym zbieżności procesu iteracyjnego danego wzorem
i+1 i (0)
X = W X + Z przy dowolnym wektorze początkowym X i dowolnym wektorze Z jest,
aby wszystkie wartości własne macierzy W były co do modułu mniejsze od jedności.
37
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 Xi+1 = W Xi + Z jest zbieżny.
Z przytoczonych twierdzeń wynika, że zbieżność procesu iteracyjnego nie zależy ani od
(0)
wartości początkowej X , ani też od wartości danego wektora Z, a wystarczy jedynie
spełnienie jednego z następujących warunków:
n
max wij < 1 norma wierszowa
i
j=1
n
max wij < 1 norma kolumnowa
j
i=1
n n
2
wij < 1 norma energetyczna
i=1 j=1
Kryterium stopu metod iteracyjnych
Obliczenia iteracyjne trwają do momentu kiedy zostanie spełniony warunek
X(i+1) - X(i) < d
gdzie: X(i+1) oraz X(i) są kolejnymi przybliżeniami rozwiązania Xi+1 = W Xi + Z oraz jest
dokładnością obliczeń dowolna liczba większa od 0
xi+1 - xi
j
Ł e , gdzie <10-1, 10-256> przyjęta dokładność obliczeń
xi+1
j
38
Numeryczne rozwiązywanie równań nieliniowych
Metoda Newtona
Dany jest układ równań nieliniowych z n niewiadomymi
f1(x1...xn ) = 0
f2 (x1...xn ) = 0
......................
fn (x1...xn ) = 0
który możemy zapisać w postaci wektorowej
f(x)=0
gdzie:
x1 f1(x)
ł ł
ę ś ę ś
: :
ę ś ę ś
x = f (x) =
ę ś ę ś
: :
ęx ś ę
fn (x)ś
n
O funkcjach fi(x1...xn) 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) = {x1(k), ..., xn(k)}
będzie k-tym przybliżeniem pierwiastka a = {a1, ..., an} równania.
Dokładną wartość pierwiastka wyraża wzór a = x(k) + (k)
(k (
Gdzie e = {e1 ) ,...,enk )} jest błędem pierwiastka przybliżonego x(k ) . Skoro istnieje ciągła
pochodna funkcji f(x) możemy zapisać:
(k)
f (a) @ f (x(k ) ) + f '(x(k ) ) e
(k )
Zastępując błąd e przyrostem h(k ) = x(k+1) - x(k) i porównując prawą stronę powyższego
wyrażenia do zera otrzymujemy równanie liniowe:
f (x(k ) ) + f '(x(k) )(x(k+1) - x(k) ) = 0
Wzór ten stanowi zapis rekurencyjny dla metody Newtona w postaci macierzowej.
39
Pochodna f (x) jest macierzą Jakobiego
śf1 śf1 śf1
ł
...
ę
śx1 śx2 śxn ś
ę ś
śf2 śf2 śf2 ś
ę
...
ę
śx1 śx2 śxn ś
f '(x) = W(x) =
ę ś
ę ś
... ... ... ...
ę ś
ę ś
śfn śfn śfn
...
ę
śx1 śx2 śxn ś
Przyjmując, że W (x(k ) ) jest macierzą nieosobliwą możemy równanie liniowe przekształcić do
postaci:
-1
x(k+1) = x(k ) -W (x(k) ) f (x(k ) )
stąd przyjmując dowolną wartość x(0) otrzymujemy ciąg kolejnych przybliżeń pierwiastka
równania x(0) , x(1) , x(2) , ..., x(k ) uzyskanych metodą Newtona.
Metoda iteracji
Dany jest układ równań nieliniowych z n niewiadomymi
x1 = g1(x1...xn )
x2 = g2 (x1...xn )
......................
x3 = gn (x1...xn )
który możemy zapisać w postaci wektorowej
x=g(x)
gdzie:
x1 g1(x)
ł ł
ę ś ę ś
: :
ę ś ę ś
x = f (x) =
ę ś ę ś
: :
ęx ś ęg (x)ś
n n
Zakładamy, że funkcja g1, g2, ..., gn są rzeczywiste i ciągłe w pewnym otoczeniu
odosobnionego pierwiastka a = {a1, ..., an} układu równań.
40
Metoda iteracji polega na stosowaniu następującego wzoru
xk+1 = g(x(k))
Po przyjęciu przybliżenia x(0) , w rezultacie otrzymujemy ciąg kolejnych przybliżeń
x(0), x(1), x(2), ... pierwiastka a równania.
Metoda siecznych
Rozpatrujemy układ równań nieliniowych w postaci
f1(x1, x2 , ..., xn ) = 0
f (x1, x2 , ..., xn ) = 0
2
...............................
fn (x1, x2 , ..., xn ) = 0
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:
(xik - xik+1) fi (xk )
xik+1 = xik -
fi (xk ) - fi (xk-1)
Metoda ma dwie wady:
- potrzebne są 2 zestawy punktów startowych: x0, f(x0), x1, f(x1)
- 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
x2 - 4 = 0
x y - 2 = 0
Napisać program wyznaczający rozwiązanie tego układu.
Niezbędne deklaracje typów:
41
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.
42
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 = x0 < x1 < ... < xn+1 = b
i utwórzmy sumę
n
Sn = - xi ) f(ci ) (1), gdzie ci xi , xi+1 są dowolnymi punktami pośrednimi
(xi+1
i=0
Jeżeli ciąg {SN} dla n Ąjest zbieżny do tej samej granicy S przy każdym przedziale
odcinka [a ; b] takim, że max xi+1 - xi 0 niezależnie od wyboru punktów Ci to funkcję f(x)
0ŁiŁn
nazywamy całkowalną w przedziale [a ; b], a granicę S ciągu (1) nazywamy:
Całką oznaczoną Riemanna funkcji f o oznaczamy:
b
f (x)dx
a
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
43
b
I(f ) =
f(x)dx = F(b)- F(a) (2)
a
I(f ) = Sn - E(f ) (3), 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.
44
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 xi oraz współczynników niezależnych od rodzaju funkcji.
Wówczas
b
n
S( f ) = f (xi ) (1) I( f ) = f (x)dx (2)
ai
i=0
a
Można wykazać, że suma ciągu S(f) jest zbieżna do całki oznaczonej I(f) dla n Ą
45
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
x [a;b] , 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 x0,
x1, & , xn. Wówczas całkę możemy zastąpić sumą, współczynniki ai zależą od metody
obliczeń
b b
n
f (x)dx =
ai
j(x)dx f (xi ) (4)
i=0
a a
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 {x0, ..., xn}
takich, że:
xi = x0 + i h
46
gdzie:
(b - a)
h = dla wówczas=(0, 1, & , n)
n
n - oznacza również stopień wielomianu interpolacyjnego
n
j(x) = f (xi )Li (x) (5), gdzie:
i=0
n
(x - xk )
Li (x) =
(x - xk ) (6)
k =0
i
iąk
wówczas
b b b
(n+1)
n n
(x(x))dx (7)
)Li(x)dx +
f(xi (x - xi )f (n +1)!
f(x)dx =
i=0 i=0
a a a
b b
n n
1
= f(xi )dx +
ai (x - xi )f (n+1)(x(x))dx (8)
(n +1)!
i=0 i=0
a a
dowolny punkt pośredni z przedziału (x,xi)
gdzie:
b
(-1)n-i n t(t -1)...(t - n)
x(x) [a;b] oraz ai = (x)dx = h (9)
i
L
i!(n - i)! t - i
a 0
ostatecznie wzór kwadratury Newtona Cotesa oraz reszta kwadratury
b
n
f (x)dx = f (xi ) + E( f ) (10)
ai
i=0
a
b
n
1
(n+1)
E( f ) =
(x - xi ) f (x (x))dx (11)
(n +1)!
i=0
a
Twierdzenie 1
Kwadratury Newtona-Cotesa oparte na n+1 węzłach są rzędu
n + 2 dla n parzystych
n +1 dla n nieparzyst
ych
wówczas wzory kwadratur dla n parzystych mają postać
47
b b n n
(n+2)
n
h(n+3)f (x)
2
f(xi )dx +
ai
f(x)dx = t (t -1)...(t - n)dt, gdzie x [a,b] (12)
(n + 2)!
i=0
a a 0
oraz dla nieparzystych mają postać:
b
(n+1)
n
h(n+2) f (x)n n
f (x)dx = f (xi ) +
ai
t(t -1)...(t - n)dt (13)
(n +1)!
i=0
a 0
Kwadratury z ustalonymi węzłami wzór trapezów
Wyznaczmy wzór kwadratury dla n = 1
x1
b
(x - x1)f(x (x - x0 )
ę 0 ś
f(x)dx = x0 - x1 )+ (x1 - x0 )f(x1)łdx + E(f ) (14)
a x0
gdzie:
x1 x1
1 1
n n
E( f ) = f (x(x))(x - xo )(x - x1)dx = f (x) - xo )(x - x1)dx =
(x
2 2
xo xo 3
x1 h3 n
(x1 + x0 ) ł
1 x3
n
= f (x)ę - x2 + x0x1xś = - f (x)
2 3 2 x0 12
wówczas:
b
ł x1 h3 n
(x - x1)2 (x - x0 )2
f (x)dx =
ę2(x - x1) f (x0 ) + 2(x1 - x0 ) f (x1)ś x0 - 12 f (x) =
a 0
48
(x1 - x0 ) h3 n
= [f (x0 ) + f (x1)]- f (x) (16)
2 12
Ostatecznie otrzymujemy
b
h h3
f(x)dx = 2 [f(a)+ f(b)]- 12 f ''(x), a < x < b (17)
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 x0=a oraz xn=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
b
h h3
n = 1 f (x)dx = [f (x0 )+ f (x1)]- f ''(x), x0 < x < x1 (18)
2 12
a
wzór parabol (Simpsona)
b
h h5 (3)
n = 2 f (x)dx = [f (x0 )+ 4 f (x1)+ f (x2 )]- f (x), x0 < x < x2 (19)
3 90
a
wzór Bessela
b
3h 3h7 (4)
n = 3 f (x)dx = [f (x0 )+ 3 f (x1)+ 3 f (x2 )+ f (x3)]- f (x), x0 < x < x3 (20)
8 80
a
49
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
50
Jeżeli dla skończonego [a ; b] przedziału wybierzemy zbiór punktów węzłowych {x0, & , xn}
takich, że :
xi = x0 + i h
gdzie:
(b - a)
h = dla i=(0, & , n)
n + 2
oraz
x0 = a + h
W oparciu o Twierdzenie 1 można również wyznaczyć wzory kwadratur otwartych. Wówczas
wzory kwadratur dla n parzystych mają postać
b n+1
(n+2)
n
hn+3f (x)
2
ai
f(x)dx = f(xi )dx + (n + 2)! t (t -1)...(t - n)dt, gdzie x [a,b] (21)
i=0
a -1
oraz dla n nieparzystych mają postać:
b
(n+1)
n
h(n+2) f (x)(n+1) n+1
f (x)dx = f (xi ) +
ai
t(t -1)...(t - n)dt, gdzie x [a;b] (22)
(n +1)!
i=0
a -1
wzór prostokątów
b
h h3
n = 0
f(x)dx = 2 f(x0 )+ 3 f ''(x), x-1 < x < x1 (23)
a
b
3
3h
n = 1
f(x)dx = 2 [f(x0 )+ f(x1)]3h f ''(x), x-1 < x < x2 (24)
3
a
b
4h 14h5 (4)
n = 2 f (x)dx = [2 f (x0 )- f (x1)+ 2 f (x2 )]+ f (x), x-1 < x < x3 (25)
3 45
a
b
5h 95h5 (4)
n = 3 f (x)dx = [11f (x0 )+ f (x1)+ f (x2 )+11f (x3)]+ f (x), x-1 < x < x4 (26)
24 144
a
51
Interpretacja geometryczna
wzór prostokątów
wzór trapezów
Kwadratury Newtona-Cotesa porównanie
Przykładowe wyniki obliczeń całki
p
p
4
4
2
sin(x)dx = [- cos(x)] = 1- 2 0,29289922
0
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
52
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
f (x)dx
a
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
5
p
p p 1 p d4(sin(x))
ć
-
ęsin ś
sin(x)dx 6 (0)+ 4sinć 2 + sin(p)ł 90 2 dx4
Ł ł Ł ł
0
d4(sin(x))
= - cos(x), wtedy max(- cos(x))=1 x 0, p
dx4
wówczas wyniki:
p
2p
sin(x)dx 3 = 2,09439510
0
oraz oszacowanie błędu obliczeń:
2
1 p p5
ć
- 1 = - = -0,106256835
90 2 2880
Ł ł
rozwiązanie dokładne:
p
p
sin(x)dx = [- cos(x)] = -cos(p) + cos(0) = 2
0
0
Kody zródłowe:
przykładowa całkowana funkcja:
function f(x:extender):extender,
begin
result:=sin (2*x*x)+2*cos(3x);
end;
53
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;
54
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:
(b - a)
xi = a + i h , dla i = (0,1,..., m -1) , h =
m
mamy wtedy:
xi+1
b
m-1
f (x)dx = f (x)dx
i=0
a xi
55
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:
b
m-1 m-1 m-1
h
h ćxi Ei f ć h
f (x)dx = f + 2 + = h xi + 2 + E (28)
Ł ł Ł ł
i=0 i=0 i=0
a
h3 m-1 (2)
gdzie błąd przybliżenia E =
f (x) (29)
3
i=0
oraz xi (xi , xi+1)
złożony wzór trapezów:
b
m-1 m-1 m-1
f (x0 ) + f (xm )
h
)ł + E (30)
Ei f(xi
f (x)dx = 2 (f (xi ) + f (xi+1))+ = hę 2 + ś
i=0 i=0 i=1
a
h3 m-1 (2)
gdzie błąd przybliżenia E = -
f (xi ) (31)
12
i=0
oraz xi (xi , xi+1)
56
złożony wzór parabol dla m parzystych:
b
m-1 m-1
h h
ćf
Ei
f (x)dx = 3 (xi ) + 4f (xi+1 + 2) + f (xi+1) + (32)
Ł ł
i=0 i=0
a
m m
-1 -1 ł
b
2 2
h
f f ś
f (x)dx = 3 ęf (x0 ) + f (xm ) + 2 (x2i ) + 4 (x2i+1)ś + E (33)
ę
i=1 i=1
a
ę ś
(b - a)
(4)
gdzie błąd przybliżenia E = - h4f (x) (34)
180
oraz xi (xi , xi+1)
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 f c[a;b](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
f (x)dx , przy m Ą
a
xi+1
b
m-1
inaczej (x)dx = lim (x)dx
f f
mĄ
i=1
a xi
(b - a)
xi = a + i h , i = (0,1,...,m -1) , h =
m
57
Przykład:
Obliczyć poniższą całkę:
Rozwiązanie dokładne
4
x
e dx =[e4 - e0]= 53,59815
0
Metoda Simpsona m=1, h=2
4
2
x
[e0
e dx 3 + 4e2 + e4]= 56,76958
0
błąd = -3,17143
Metoda Simpsona m=1, h=1
4 2 4
1 1
x x x
[e0 [e2
e dx = e dx + e dx 3 + 4e1 + e2]+ 3 + 4e3 + e4]= 53,86385
0 0 2
błąd = -0,26570
Metoda Simpsona m=4, h=0,5
4 1 2 3 4
x x x x x
e dx = e dx + e dx + e dx + e dx K = 53,61622
0 0 1 2 3
błąd = -0,01807
Przykładowe wyniki obliczeń całki:
p
p
0
sin(x)dx = [- cos(x)] = 2
0
Wyniki obliczeń dla wzorów złożonych w zależności od m:
m 2 4 8 16 32 64
1,110720735 1,751785441 1,936297295 1,983970758 1,995986709 1,998996147
wzór wynik
0,889279265 0,248214559 0,066102705 0,016029242 0,004013791 0,001003853
prostokątów błąd
1,570796327 1,896118898 1,1974231602 1,993570344 1,998393361 1,999598389
wzór wynik
0,429203673 0,103881102 0,025768398 0,006429656 0,001606639 0,000401611
trapezów błąd
2,094395102 2,004559755 2,000269170 2,000016591 2,000001033 2,000000065
wynik
wzór parabol
błąd -0,094395102 -0,004559755 -0,000269170 -0,000016591 0,000001033 0,000000065
Kody zró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;
58
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;
59
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:
b
n
ci
f (x)dx = f (xi )
i=1
a
Węzły kwadratury x1, x2,..., xn z przedziału całkowania [a,b] oraz współczynniki c1, c2,...,cn są
tak dobrane, aby zminimalizować błąd przybliżenia. Nie zakładamy jednak żadnych
ograniczeń na węzły xi o współczynnikach ci natomiast chcemy zmaksymalizować rząd
kwadratury.
60
Mamy znalezć 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:
b
n
ci
f (x)dx = f (xi ) + E
i=1
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ć x1, x2 oraz c1, c2, ponieważ:
1
f (x)dx c1 f (x1) + c2 f (x2 ) (64)
-1
Zgodnie z założeniami, powyższe przybliżenie będzie dokładne tylko gdy stopień funkcji f(x)
mniejszy lub równy 2n -1= 2(2) -1= 3 wówczas
f (x) = a0 + a1x + a2x2 + a3x3 (65)
gdzie a0, a1, a2, a3 dowolne stałe
a więc możemy zapisać:
(a0 + a1x + a2x2 + a3x3)dx = a0 + a1 + a2 2dx + a3 3dx
1dx xdx x x
Poszukujemy x1, x2 oraz c1, c2, więc można zapisać:
c1 1+ c2 1 = 2 (67)
1
c1 x1 + c2 x2 = x dx = 0 (68)
-1
1
2
2
c1 x12 + c2 x2 2 = dx = (69)
x
3
-1
1
c1 x13 + c2 x23 = = 0 (70)
dx
-1
otrzymaliśmy układ liniowych równań algebraicznych, którego rozwiązanie jest:
3 3
c1 = 1 c2 = 1 x1 = - x2 =
3 3
61
ostatecznie otrzymaliśmy wzór na przybliżenie całki dla funkcji danej wzorem (64):
1
ć ć
3 3
f (x)dx f - 3 + f 3
-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
( f , g)= f (x) g(x)dx = 0
a
Dla naszych dalszych rozważań, zajmiemy się tylko ciągiem wielomianów Legendre a
{P0(x), P1(x)KPn(x)K}, który posiada właściwości:
1. Dla każdego n, Pn(x) jest wielomianem stopnia n
1
2. P(x) Pn(x)dx = 0 , gdzie P(x) jest wielomianem stopnia wyższego (niższego???) od n
-1
Wzory kilku pierwszych wielomianów Legendre a:
3 1 5 5
ć ć
P0(x) =1 P1(x) = x P2(x) = x2 - P3(x) = x3 - x
2 3 2 3
Ł ł Ł ł
oraz wzór rekurencyjny do wyznaczania wielomianów krzywych ??? oraz ich pochodnych:
2n -1 n -1 n n
Pn(x) = x Pn-1(x) - Pn-2(x) Pn'(x) = x Pn(x) - Pn-1(x)
n n x2 -1 x2 -1
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 x1, x2, ... xn są pierwiastkami wielomianu Lagrange a, Pn(x) stopnia n dla każdego i =
(1,2,...n) wtedy współczynniki ci można obliczyć ze wzoru:
1
n
x - xk - 2 - 2(xi 2 -1)
c1 = dx = =
2
xi - xk (n +1)Pn+1(xi )Pn'(xi )
[(n +1) Pn-1(xi )]
k =1
-1
k ąi
wówczas jeżeli P(x) jest dowolnym wielomianem stopnia mniejszego od 2n wtedy:
1
n
ci
P(x)dx = P(xi ) (75)
i=1
-1
62
1
2
(P(x)) dx
(2n+2)
-1
reszta dla stopnia większego od (2n-1): E = f (x)
(2n +1)!
Pierwiastki wielomianu oraz współczynniki ci 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 a,b dla którego f(x) jest ciągła należy wykonać
przekształcenie poprzez podstawienie.
Dokonujemy przekształcenia dowolnego przedziału a,b do przedziału -1,1 .
b - a b + a b - a
x = t + dx = dt
2 2 2
Wówczas możemy zapisać:
b 1
b - a b - a b + a
ć dt
f (x)dx = f t +
2 2 2
Ł ł
a -1
po zastosowaniu kwadratury Gaussa-Lagrange a otrzymujemy
b 1
b - a b - a b + a b - a b - a b + a
ć dt n f ć
f (x)dx = f t + xi +
ci
2 2 2 2 2 2
Ł ł Ł ł
i=1
a -1
gdzie xi pierwiastek wielomianu stopnia n
ci współczynnik ze wzoru (74)
przykład:
Obliczyć całkę wykorzystując kwadraturę Gaussa:
1,5
-x2
e dx = 0,1093643
1
dokończyć ??? najpierw przekształcenia przedziału całkowania:
1,5 1 1
ć1,5 -1 1,5 +12ł
1,5 -1 1 - t + 5
ć dt
x2
f t +
ę śdt = f
e dx 2
2 4 16
ł Ł ł
ęŁ 2 ś
1 -1 -1
po zastosowaniu kwadratury Gaussa-Lagrange a, n=2 otrzymujemy:
3 3
c1=1 c2=1 x1 = - x2 =
3 3
63
ć ć
ć ć
3
ł
+5 - 3 +5
- -
3
3
Ł ł Ł ł
ę ś
16 16
-x1+5 -x2 +5
ć ć
1,5 ę ś
ł
1 1
x2 16 16
ł ł ł
+1eŁ ł ś = 0,1094003
e dx 4 ęc1 eŁ + c2 eŁ ś = 4 ę1eŁ
ę ś
ę ś
1
ę ś
ę ś
Po zastosowaniu kwadratury Gaussa-Lagrange a, n=3 otrzymujemy:
5 8 5 15 - 15
c1 = c2 = c3 = x1 = x2 = 0 x3 =
9 9 9 3 3
ć
ć
15
ł
+5
- - ć
15
3 - +5
Ł ł
ę ś
3
-x1+5 -x2 +5 -x3 +5 -(0)+5
ć ć ć ę5 16 8 ć 16 ś
1,5
ł
1 1 5
x2 16 16 16
ł ł ł ł
+ eŁ 16 ł + eŁ ł ś =
e dx 4 ęc1 eŁ + c2 eŁ + c3 eŁ ś = 4 ę9 eŁ
ę ś
9 9
ę ś
1
ę ś
ę ś
= 0,1093642
Tablica pierwiastków wielomianu Legendre a:
n x0 x1 x2 x3 x4
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 zródłowy kwadratury Gaussa:
Przykład f. całką wykorzystując kwadraturę G w przedziale a,b 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];
64
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:
b
n
I =
ci
w(x) f (x) dx = P(xi )
i=1
a
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
I = f (x) dx , gdzie f(x):
a
przypadek 1* przypadek 2*
e 0
65
b c b
f (x) dx = f (x) dx + c f (x) dx
a a a
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:
1
ex
dx 2,9235450 + 0,0017691 = 2,9253141
x
0
Metoda adaptacyjna:
1 1
ex ex
dx dx =2,925324091
x x
0 0,0000000001
zał. dokładne 0,0001 liczba przedzi. 235 il. wywoł. funkcji 941 ???
Całki wielokrotne:
rozpatrzmy całkę podwójną:
d ( x)
b
f (x, y) dy dx
a c( x)
x stała i stosunek kwadratury do obliczonej całki
d ( x)
s(x) = f (x, y) dy
c( x)
66
otrzymujemy:
m-1 m
ł
d(x) - c(x)
s(x) f (x, y0) + f (x, y2m) + 2 f (x, y2 j ) + 4 f (x, y2 j-1)ś
ę
6m
j=1 j=1
m liczba przedziałów przedziału [c(x) ; d(x)]
wykorzystanie z kolei kwadratury Simpsona do obliczenia całki:
b
s dx
S(x)
a
ostatecznie:
n-1 n
b - a ł
s ) + 4 )ś
0 S(x2i S(x2i-1
ęS(x ) + S(x2n) + 2
6n
i=1 i=1
2,0 1,5
np.:
ln(x + 2y) dy dx = 0,4295549269
1,4 1,0
stos. złoż. kwadratury Simpsona:
dla n=4 i m=2
2,0 1,5
ln(x + 2y) dy dx 0,42955224387
1,4 1,0
(- 0,5)(0,6)
E = [ ]
1800
E Ł -4,7210-6
Całkowanie wielokrotne całka podwójna
kod zró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;
67
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 zró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;
68
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 sprawdzcie!!!
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
I =
f(x)dx (129)
b
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 x1, x2, ..., xn,
69
Przybliżoną wartość całki określa wzór
n
(b - a)
In = f (xi ) (130)
n
i =1
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:
lim = P( n -1 < e) 1 (131)
n Ą
gdzie P jest prawdopodobieństwem, przy czym błąd metody monte Carlo można określić
wzorem
1
In - I < (132)
n
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 zró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;
70
Całkowanie Metody Monte Carlo. Całki wielokrotne właściwe
Dana jest całka wielokrotna:
I = f (x1, x2,..., xn )dx1dx2...dx3 (135)
W
Każda ze zmiennych xi zawiera się w przedziale ai Ł xi Ł bi (136)
Losujemy n-krotnie punkty Pk kostki, przyjmując że współrzędnymi tych punktów są
zmienne losowe X1(k), X2(k),..., Xm(k) o rozkładzie równomiernym w odpowiednich
przedziałach [ai,bi]. Każdemu punktowi Pk odpowiada zatem ciąg m-elementowy
wylosowanych realizacji x1(k), x2(k),..., xm(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:
Pk W dla k=(1, 2, ..., n), Pk W 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:
n
VW @ V (138),
N
gdzie
m
V = - a ) (139)
(bi
i=1
i
jest objętością m-wymiarowej kostki. Przybliżona wartość In całki (135) określana wzorem:
n
V
In f (Pk ) (140)
N
k =1
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ę xi(k) w przedziale [ai; bi] korzystając ze zależności:
xi(k)=(bi-ai)yi(k)+ai
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
71
Numeryczne obliczanie pochodnej
Aby obliczyć pochodną funkcji f w punkcie x0, mamy
f (x0 + h) - f (x0 )
f '(x) = lim = (1)
h 0
h
W celu wyznaczenia przybliżonej wartości pochodnej f (x) dla małych wartości h możemy
użyć wyrażenia
f (x0 + h) - f (x0 )
(2)
h
Możemy to udowodnić. Dokonujemy interpolacji f(x); zdefiniujemy wielomian Lagrange`a
stopnia 1, który będzie oparty na węzłach x0 oraz x0+h pod warunkiem, że funkcja f jest
dwukrotnie różniczkowalna w przedziale [a; b]
(x - x0 )(x - x0 - h)
f (x) = P1 (x) + ktoś nie zdązył do końca, a ja nie wiem czy w
2!
mianowniku jest 2!, czy co??
Wówczas obliczając pochodną funkcji f, mamy:
f (x + h) - f (x )
0 0
f '(x ) = + e (3)
0
h
oraz po uproszczeniu
h
e = - f "(x) (4)
2
72
Jeżeli dokonamy interpolacji funkcji f wielomianem stopnia 2, wówczas
f (x0 + h) - f (x0 - h)
f '(x0 ) = + e (5)
2h
oraz po uproszczeniu
h2 (3)
e - f (x ) I(6)
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 x0=1,8 stosując metodę w przód i wstecz
stopnia 1.
f (x0 + h) - f (x0 ) f (1,8 + h) - f (1,8)
h h
f '(x0 ) = - f "(x) = - M ,
h 2 h 2
gdzie:
M = max f "(x)
x (1,8,1,8 + 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
73
Numeryczne obliczanie drugiej pochodnej.
Jeżeli dokonamy rozwinięcia funkcji f w szereg Taylora 3 rzędu w otoczeniu punktu x0 dla
punktów x0+h oraz x0-h, zakładając, że f jest 4-krotnie różniczkowalna na przedziale [x0-h,
x0+h], wówczas dla x0-h:
Wzor (8)
gdzie [x0 - h; x0 ]
dodając oba równania stronami otrzymujemy:
wzór (9)
Rozwiązując równanie ze względu na f (x0), otrzymujemy:
1
f "(x0 ) = [ f (x0 + h) - 2 f (x0 ) + f (x0 - h)] + e (11)
h2
oraz
hn (n)
e = [ f (x ) , gdzie x [x0 - h; x0 + h]
12
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(x3) dla x0=2 stosując powyższą metodę
różnych wartości kroku h.
Wzór
Rozwiązanie dokładne
3
f "(x) = - f (2)=-0,1875
x2
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 znalezć 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
74
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 t [a, b] spełnia równanie
różniczkowe zwyczajne pierwszego rzędu:
dy
= f (t, y) (7)
dx
oraz warunek początkowy
y(a)=y0 (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.
dy
= f (t, y) , gdzie t [a, b] oraz y(a)=a (10)
dx
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ć:
(ti +1 - ti )2
y(ti +1) = y(ti ) + (ti +1 - ti )y' (ti ) + y" (xi ) (11), gdzie xi (ti , ti +1)
2
Oznaczamy h=ti+1-ti, wówczas otrzymujemy:
h2
y(ti +1) = y(ti ) + hy' (ti ) + y" (xi ) (12)
2
Użyjemy podstawienia y (t)=f(t, y), wówczas otrzymujemy:
h2
y(ti +1) = y(ti ) + hf (ti , y(ti )) + y" (xi ) (13)
2
Zapisując iH"y(ti) oraz pomijając błąd przybliżenia otrzymujemy i+1=i+hf(ti, i) (14)
75
Powyższy wzór nazywamy metodą Eulera wzór ten nazywany jest inaczej równaniem
różniczkowym, gdyż można zapisać:
wi +1 - wi
f (ti , wi ) = (15)
h
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ć:
wi +1 - wi
h2
f (ti , y(ti )) = + e , gdzie ei = y" (xi ) (16)
i
h 2
Lokalny błąd dyskretyzacji i+1(h)
1 1
t (h) = (yi +1 - yi + f (ti , wi ) = (yi +1 - wi +1) (17)
i +1
h h
dodatkowo możemy określić krok h, dla którego błąd lokalny jest mniejszy od zadanej
dokładności
2d
h = , gdzie M = max y" (x)
x(t , t )
i i + 1
M
Globalny błąd dyskretyzacji g(x)
g(t)= (t)-y(t) (19)
Dla wybranego punktu ti możemy zapisać:
gi=i-y(ti) (20), wówczas
hM
i 0
gi = y(ti ) - wi Ł [eL(t - t ) -1], gdzie L- liczba Lipschnitz`a
2L
76
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-t2+1 0dla N=10, wtedy h=0,2; ti=ih, =0,5 oraz i+1=i+hf(ti, i), gdzie f(t,y)=y-t2+1
i+1=i+h(i-t2i+1)= i+0,2[i-0,04i2+1]=1,2i-0,008i2+0,02 dla i=0, 1, ..., 9
Rozwiązanie dokładne ma postać
y(t)=(t2-1)-0,5et
77
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:
wi+1 = wi + hf (ti+1, wi+1) (22)
wi+1(k) = wi + hf (ti+1, wi+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 wi+1 wykorzystujemy wartości z poprzedniego kroku oraz poszukiwaną wartość
wi+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 wi+1(k),
stojącej po prawej stronie równania podstawiając jako wi+1(k-1) - lewa strona równania, wynik
przybliżenia z poprzedniej iteracji (k-1). Proces trwa do momentu, kiedy spełniony zostanie
warunek:
|wi+1(k) - wi+1(k-1)| d" (24)
gdzie - tolerancja obliczeń
Kody zró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;
78
Przykładowa funkcja definiująca równanie różniczkowe
dy
= f (t, y)
dt
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)wi_1:=wi;
end;
Result:=wi;
end;
Przykładowa funkcja definiująca równanie różniczkowe
dy
= f (t, y)
dt
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:
dy
= f (t, y) gdzie t [a,b] oraz y(a)=a (25)
dt
i przedstawieniu różnicy funkcji y(t) w punktach ti+1 oraz ti w postaci:
m
wi+1 - wi = k lub inaczej wi+1 wi = h0 (ti,wi,h) (26)
cj j
i=1
79
gdzie m oznacza rząd metody, cj są stałymi, a
j-1
ć
k = hf + g h, wi + kl = hf (ti,a , wi + b k )
j dli j j j-1
ti j
i=1
Ł ł
gdzie ąj, j, łj, lj stałe
h - krok całkowania
Określenie stałych cj, ąj, j otrzymujemy poprzez rozwinięcie funkcji f(t,y) w szereg Taylora
w otoczeniu punktu ti
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 c1, a1,b1, c2, a2, b2 :
Metoda punktu środkowego:
C1 = 0 C2 = 1 ą1 = 0 ą2 = h/2 b1 = 0 b2 = (28)
wówczas możemy zapisać:
k1=hf (ti,wi) k2=hf (ti + h, wi + k1)
ostatecznie:
wi+1=wi+k2
lub
wi+1=wi + hf(ti + h/2, wi + h/2 * f(ti, wi))
80
Interpretacja graficzna:
f1=f(ti,wi) f2=f(ti + h/2, wi + h/2 * f1)
Stosując metody Rungego Kutty 2 rzędu , rozwiązać następujące równanie różniczkowe:
y =y-t2 + 1 0Ł t Ł2 y(0)=0,5
dla N = 10 wtedy h = 0,2; ti = ih ; w0 = 0,5 oraz
metoda punktu środkowego wi+1 = 1,22 wi 0,0088 i2 + 0,218
Metoda zmodyfikowana Eulera
wi+1 = 1,22 wi - 0,0088 i2 + 0,216
Metoda Heana dla i = 0, 1, ..., 9
wi+1 = 1,22 wi - 0,0088 i2 + 0,2173
Rozwiązanie dokładne ma postać:
y(t) = (t2+1) - 0,5 et
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
k1 = hf (ti, wi)
k2 = hf (ti + *h, wi + *k1)
k3 = hf (ti + *h, wi + *k2) (37)
81
k4 = hf (ti+h, wi+k3)
ostatecznie:
wi+1=wi + 1/6 (k1 + 2k2 + 2k3 + k4) (38)
Interpretacja graficzna:
f1 = f (ti, wi)
f2 = f (ti + * h, wi+1 + h f1)
f3 = f (ti + * h, wi + * h f2)
f4 = f (ti + h, wi + h f3)
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-t2 + 1 0Ł t Ł 2 y(0) = 0,5
dla N = 10 wtedy h = 0,2 t1 = i*h w0 = 0,5 oraz:
Metoda Rungego-Kutty 4-rzedu:
wi+1 = 1,2214 wi - 0,008856 i2 + 0,00856 i + 0,218593
dla i = 0, 1, ..., 9
Rozwiązanie dokładne ma postać y(t) = (t2+1) - 0,5et
82
Kody zró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;
k1 := h*ft, y_i);
k2 := h*f(t + 0.5 * h, y_i+0.5 * k1);
k3 := h*f(t+0.5*h, y_i + 0.5 * k2);
k4 := 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
dy
= f (t, y) gdzie t [a,b] oraz y(a) = a (55)
dt
jest równaniem różniczkowym dla znalezienia przybliżenia wi+1, dla punktu węzłowego ti+1,
gdzie m>1, które przyjmuje postać:
wi+1 = am-1 * wi + am-2 * wi-1 + ... + a0 * wi+1_m + h [bn f(ti+1,wi+1) + bn-1 f(ti, wi) + ... + t + b0
f(ti+1_m, wi+1_m)]
dla i = (m-1, m, ..., N-1),
N - liczba węzłów w przedziale, gdzie h = (b-a) / N
83
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
u(x, y) = g(x, y) dla (x, y)S
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ć.
84
Równania różniczkowe cząstkowe warunek brzegowy Neumanna
Warunek brzegowy Neumanna
śu
(x, y) = g(x, y) dla (x, y)S
śn
śu
gdzie: (x, y) pochodna normalna poszukiwanej funkcji w punktach należących do
śn
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
u(x, y) - u(x - h, y)
= g(x, y)
h
Dla zmiennej y
u(x, y) - u(x, y - h)
= g(x, y)
h
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ć:
d2u d2u
Ń2u(x, y) = (x, y) + (x, y) = f (x, y) (1)
dx2 dy2
na R ={(x, y) a < x < b;c < y < d }
z warunkami brzegowymi u(x, y) = g(x, y) dla (x, y)S
W powyższym równaniu funkcja f opisuje dane wejściowe na płaszczyznie 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 wezmiemy 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:
85
xi = a + ih dla i = 0,1, .. n
oraz yj = c + jk 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
u(xi+1, yj ) - 2u(xi, yj ) + u(xi-1, yj )
ś2u h2 ś4u
(xi, yj ) = - (xi, yj ) (2)
śx2 h2 12 śx4
gdzie: xi (xi-1, xi+1)
Dla zmiennej y
u(xi, yj+1) - 2u(xi, yj ) + u(xi, yj-1)
ś2u k2 ś4u
(xi, yj ) = - (xi,h )(3)
śy2 k2 12 śy4 j
gdzie: h (yj-1, yj+1)
j
86
Podstawiając wzory (2) oraz (3) do równania Poissona (1) otrzymujemy:
u(xi+1, yj ) - 2u(xi, yj ) + u(xi-1, yj ) u(xi, yj+1) - 2u(xi, yj ) + u(xi, yj-1)
+
h2 k2
h2 ś4u k2 ś4u
= f (x, y) + (xi, yj ) + (xi,h ) (4)
12 śx4 12 śy4 j
dla i = 1,2 & n-1 oraz j = 1,2, & m-1
natomiast warunki brzegowe mają postać:
u(x0, yj ) = g(x0, yj ) oraz u(xn, yj ) = g(xn, yj ) dla j = 0,1, .. m (5)
u(xi, ym) = g(xi, ym) oraz u(xi, y0) = g(xi, y0) 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:
2 2
ł
h h
ć
2ęć +1świj - (wi+1, j + wi-1, j ) - (wi, j+1 + wi, j-1) = -h2 f (xi, yj ) (7)
k
Ł ł
ęŁ k ł ś
dla i = 1,2 & n-1 oraz j = 1,2, & m-1
natomiast warunki brzegowe mają postać:
w0 j = g(x0, yj ) oraz wnj = g(xn, yj ) dla j = 0,1, .. m (8)
wi0 = g(xi, y0) oraz wim = g(xi, ym) dla i = 1,2, .. n (9)
gdzie wij jest przybliżeniem u(xi,yj).
87
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ądz iteracyjnymi.
Zakładając, że h = k możemy wyprowadzić uproszczoną postać równania (7) :
4wij - (wi+1, j + wi-1, j ) - (wi, j+1 + wi, j-1) = -h2 f (xi, yj ) (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 :
2
l) l) ) )
h
- h2 f (xi, yj ) + (wi(+1, j + wi(-1, j ) +( ) (wi(,lj+1 + wi(,lj-1)
( k
z równania (7) wijl+1) = (11)
2
h
2[( ) +1]
k
l) l) ) )
- h2 f (xi, yj ) + (wi(+1, j + wi(-1, j + wi(,lj+1 + wi(,lj-1)
(
z równania (10) wijl+1) = (12)
4
88
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ę zródła
ciepła utrzymujące temperaturę na poziomie 0C dla boku dolnego i prawego, natomiast
temperatura boku górnego i lewego zmienia się liniowo od 0C do 100C 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
2 2
d T d T
Ń2T(x, y) = (x, y) + (x, y) = 0
dx2 dy2
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)
T3,1 ł
4 -1 0 -1 0 0 0 0 0 ł 25
ł
ęT ś
ę ę
50ś
ę ś
ę-1 4 -1 0 -1 0 0 0 0ś 3,2 ę ś
ś
ę ś
ę ś T3,3 ę ś
0 -1 4 0 0 -1 0 0 0 100
ęT ś
ę ę
0ś
ę ś
ę-1 0 0 4 -1 0 -1 0 0ś 2,1 ę ś
ś
ę ę
0 -1 0 -1 4 -1 0 -1 0ś ęT2,2 ś = 0ś
ę ś
ę ś ę ś
0 0 -1 0 -1 4 0 0 -1ś ęT2,3 ś ę 50ś
ę
ę
0 0 0 -1 0 0 4 -1 0ś ęT1,1 ś ę 0ś
ę ś
ę ś ę ś
0 0 0 0 -1 0 -1 4 -1ś 0ś
ęT1,2 ś
ę ę
ę
0 0 0 0 0 -1 0 -1 4ś ęT1,3 ś ę 25ś
89
Rozwiązanie układu równań metodą Gaussa-Siedla, daje wyniki
T3,1 T3,2 T3,3 T2,1 T2,2 T2,3 T1,1 T1,2 T1,3
18,75 38,50 56,25 12,50 25,00 37,50 6,25 12,50 18,75
Ostatecznie wyniki dla założonej siatki mają postać:
x = 0 x = 0,125 x = 0,25 x = 0,375 x = 0,5
y = 0,5 0,00 25,00 50,00 75,00 100,00
y = 0,375 0,00 18,75 37,50 56,25 75,00
y = 0,25 0,00 12,50 25,00 37,50 50,00
y = 0,125 0,00 6,25 12,50 18,75 25,00
y = 0 0,00 0,00 0,00 0,00 0,00
Rozwiązanie dokładne problemu ma postać:
T(x, y) = 400xy
Natomiast oszacowanie błędu obliczeń :
ł
h2 ś4T ś4T
t = (xi, yj ) + (xi,h )ś = 0
ę
12 śx4 śy4 j , bo
ś4T ś4
ś4T ś4
(x, y) = (400xy) = 0
(x, y) = (400xy) = 0
śy4 śy4
śx4 śx4
i
Wynika z tego, że rozwiązanie numeryczne jest równe z rozwiązaniu dokładnemu.
Równania eliptyczne - przykład 1 - kody zró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
90
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 );
91
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;
92
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ę zródła ciepła
utrzymujące temperaturę na poziomie 100C dla boku górnego, 50C dla boku lewego, 0C
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
2 2
d T d T
Ń2T(x, y) = (x, y) + (x, y) = 0
dx2 dy2
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 ]
93
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
94
Równania eliptyczne - przykład 2 - kody zró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;
95
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ć:
2
du d u
2
(x,t) -a (x,t) = 0
dla oraz (13)
0 < x < l t > 0
dt dx2
z warunkami brzegowymi dla
u(0,t) = u(l,t) = 0 t > 0
u(x,0) = f (x) 0 < x < l
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 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 wezmiemy 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:
xi = ih t = jk
dla i = 0,1, .. m oraz dla j = 0,1, ..
j
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:
u(xi+1,t ) - 2u(xi,t ) + u(xi-1,t )
ś2u h2 ś4u
j j j
Dla zmiennej x (xi, yj ) = - (xi,t ) (14)
śx2 h2 12 śx4 j
gdzie: xi (xi-1, xi+1)
u(xi,t ) - u(xi,t )
śu k ś2u
j+1 j
(xi,t ) = - (xi, m )
Dla zmiennej t (15)
j
śt k 2 śt2 j
m (t ,t )
gdzie:
j j j+1
96
Podstawiając wzory (14) oraz (15) do równania dyfuzji (13) oraz wyłączając oszacowanie
błędu zapisujemy:
u(xi,t ) - u(xi,t ) u(xi+1,t ) - 2u(xi,t ) + u(xi-1,t )
j+1 j j j j
2
(16)
-a = 0
k h2
u(0,t ) = u(l,t ) = 0
warunek brzegowy : (17)
j j
dla i = 1,2 & m-1 oraz j = 1,2, &
natomiast lokalny błąd odcięcia jest równy :
k ś2u h2 ś4u
(18)
2
tij = (xi, m ) -a (xi,t )
2 śt2 j 12 śx4 j
gdzie: oraz
xi (xi-1, xi+1) m (t ,t )
j j j+1
Ostatecznie wzór na MRS, podstawiając za ( xi , tj ) = wi,j oraz wyłączając oszacowanie
błędu zapisujemy:
wi, j+1 - wi, j 2 wi+1, j - 2wi, j + wi-1, j
-a = 0
(19)
k h2
przekształcając powyższy wzór ze względu na wi,j+1 ,wówczas otrzymujemy:
2
ć
2a k k
2
(20)
wi, j+1 = ( )
wi,
1- h2 j +a h2 wi+1, j + wi-1, j
Ł ł
lub
2
a
wi, j+1 = (1- 2l)wi, j + l(wi+1, j + wi-1, j)
l = kć
gdzie: (21)
h
Ł ł
dla i = 1,2 & m-1 oraz j = 1,2, &
97
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 :
w( j+1) = Aw( j)
(22)
co rozpisując daje nam:
(1- 2l) l
ł
0 0
w1, j+1 ę w1,
ł ł
ś j
ę ś ę ś
ś
w2, j+1 ę l w2, j
ę ś ę ś
ę ś
0 (23)
ę ś ę ś
M M
=
ę ś
ę ś 0 ę ś
ę ś
M M
ę ś ę ś
ę ś
l
ęwm-1, j+1ś ęwm-1, ś
ę ś
l (1- 2l)ś j
ę 0 0
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.
98
Równania paraboliczne - Dyfuzja - przykład 1
Rozwiązać poniższe równanie przewodnictwa:
2
du d u
2
0 < x <1 t > 0
(x,t) -a (x,t) = 0
dla oraz
dt dx2
z warunkami brzegowymi u(0,t) = u(1,t) = 0 dla t > 0
u(x,0) = sin(p x) 0 < x <1
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ąć a 2 = 1
Rozwiązanie dokładne ma postać:
2
u(x,t) = e-p t sin(p x)
Tabela wyników schemat jawny
Mały krok czasu Duży krok czasu
| u( xi , 0.5 )
xi u( xi , 0.5 ) wi,1000 | u( xi , 0.5 ) wi,50
k = 0.0005 wi,1000| k = 0.01 wi,50|
0,0 0 0 0
0,1 0,002222410,002286526,411 10-5 8,19876 107 8,199 107
0,2 0,004227280,004349221,219 10-4 -1,55719 108 1,557 108
0,3 0,005818360,005986191,678 10-4 2,13833 108 2,138 108
0,4 0,006839890,007037191,973 10-4 -2,50642 108 2,506 108
0,5 0,007197880,00739934 2,075 10-4 2,62685 107 2,627 108
0,6 0,006839890,007037191,973 10-4 -2,49015 108 2,490 108
0,7 0,005818360,00598619 1,678 10-4 2,11200 107 2,112 108
0,8 0,004227280,004349221,219 10-4 -1,53086 108 1,531 108
0,9 0,002222410,00228652 6,411 10-5 8,03604 107 8,036 107
1,0 0 0 0
99
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
2
ip
r(A) = max 1- 4l(sin( )) Ł1
(23)
2m
1ŁiŁm-1
z powyższego wyrażenia wynika, że warunek stabilności tej metody jest:
k 1
2
a Ł (24)
h2 2
Wykorzystując powyższy wzór można wyznaczyć dla przykładu 1 graniczną wartość kroku
czasowego, dla której schemat jawny jest stabilny:
2 2
1 h 1 0,1
ć ć
k Ł k Ł = 0,005
wtedy
2 2 1
Ła ł Ł ł
100
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:
u(xi,t ) -u(xi,t )
śu k ś2u
j j-1
Dla zmiennej t (25)
(xi,t ) = - (xi, m )
j
śt k 2 śt2 j
gdzie:
m (t ,t )
j j-1 j
Wówczas podstawiając wzór (25) oraz (14) do równania parabolicznego (13), pomijając błąd
obcięcia otrzymamy :
u(xi,t ) -u(xi,t ) u(xi+1,t ) - 2u(xi,t ) + u(xi-1,t )
j j-1 j j j
2
-a = 0
(26)
k h2
Modyfikując powyższy wzór, podstawiając za ( xi , tj ) = wi,j zapisujemy:
wi, j - wi, j-1 2 wi+1, j - 2wi, j + wi-1, j
(27)
-a = 0
k h2
lub:
2
(28)
(1+ 2l)wi, j - l(wi+1, j + wi-1, j ) = wi, j-1 gdzie: l = kća
h
Ł ł
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).
101
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 :
Aw( j) = w( j-1)
(29)
co rozpisując daje nam:
(1+ 2l)
ł
- l 0 0
w1, j w1, j-1
ł ł
ę ś
- l
ę
ę ś
w2, j ś ę w2, j-1 ś
ę ś ę ś
ę ś
0
(30)
ę ś ę ś
M M
=
ę ś
0 ę ś ę ś
ę ś
M M
ę ś ę ś
ę ś
- l
ęwm-1, ś ęwm-1, ś
ę ś
(1+ 2l)ś j j-1
ę 0 0 - l
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:
2
du d u
2
0 < x <1 t > 0
(x,t) -a (x,t) = 0
dla oraz
dt dx2
z warunkami brzegowymi dla
u(0,t) = u(1,t) = 0 t > 0
i początkowymi dla
u(x,0) = sin(p x) 0 < x <1
Obliczenia wykonać dla czasu t = 0.5 używając schematu niejawnego dla h = 0.1 oraz k
= 0.01 ,przyjąć a 2 = 1
2
u(x,t) = e-p t sin(p x)
Rozwiązanie dokładne ma postać:
102
Tabela wyników schemat niejawny
wi,50 | u( xi , 0.5 )
xi u( xi , 0.5 )
k = 0.01 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:
2
ip
r(A) = max 1+ 4l(sin( )) Ł1
(31)
2m
1ŁiŁm-1
Równania paraboliczne kody zró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;
103
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];
104
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;
105
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)
2 2
d u d u
2 dla 0 < x < l oraz t > 0
(x,t) -a (x,t) = 0
2
dt dx2
u(0,t) = u(l,t) = 0
z warunkami brzegowymi dla t > 0
u(x,0) = f (x) śu 0 < x < l
i początkowymi oraz dla
(x, 0) = g(x)
śt
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 wezmiemy 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:
xi = ih t = jk
dla i = 0,1, .. m oraz j 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 druga
pochodna numeryczna) . Wówczas otrzymujemy:
ś2u u(xi+1 ,t j ) - 2 u(xi,t j ) + u(xi-1,t j ) h2 ś4u
Dla zmiennej x (33)
(xi,t ) = - (xi,t )
śx2 j h2 12 śx4 j
106
xi (xi-1, xi+1)
gdzie
ś2u
Dla zmiennej t (xi,t j ) = u(xi,t j+1) - 2u(xi,t j ) + u(xi,t j-1) - k2 ś4u (xi, m j ) (34)
śy2 k2 12 śt4
gdzie m j (t j-1,t j+1)
Podstawiając wzory (33) oraz (34) do równania falowego (31) oraz wyłączając oszacowanie
błędu zapisujemy:
u(xi ,t ) - 2u(xi,t j ) + u(x i,t ) u(x i+1,t ) - 2u(xi,t j ) + u(x i-1,t )
j+1 j-1 j j
2
(35)
-a = 0
k2 h2
warunek brzegowy : (36)
u(0,t ) = u(l,t ) = 0
j j
dla i = 1,2 & m-1 oraz j = 1,2, &
natomiast lokalny błąd odcięcia jest równy :
ł
1 ś4u 2 ś4u
2
tij =
j
ęk (xi, m ) -a h2 śx4 (xi,t )ś (37)
12 śt4 j
gdzie: oraz
xi (xi-1, xi+1) m (t ,t )
j j-1 j+1
Ostatecznie wzór na MRS, podstawiając za ( xi , tj ) = wi,j zapisujemy:
wi, j+1 - 2wi, j + wi, j-1 2 wi+1, j - 2wi, j + wi-1, j (38)
- a = 0
k2 h2
przekształcając powyższy wzór ze względu na wi,j+1 ,wówczas otrzymujemy:
a
l = k
wi, j+1 = 2(1- l2)wi, j + l2(wi+1 , j + wi-1 , j)- wi, j-1 gdzie: (39)
h
dla i = 1,2 & m-1 oraz j = 1,2, &
w0 , j = wm, j = 0
z warunkami brzegowymi dla j = 1,2, &
wi ,0 = f (xi )
i początkowymi dla i = 1,2, & m-1
107
Analizując równanie (39) 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 oraz jednej wartości dla czasu tj-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( xi , 0)=f(xi), następnie wygenerujemy z drugiego warunku
początkowego wartości punktów dla t = t1 jako du( xi , t1) / dt = g(xi).
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 wi,j siatki dla i = 1,2, .. m-1 :
w( j+1) = Aw( j) - w( j-1)
(40)
co rozpisując daje nam:
ł
2(1- l2) l2 0
0
w1, j+1 ę w1, j w1, j-1
ł ł ł
ś
ę ś ę ś ę ś
ś
w2, j+1 ę l2 w2, j w2, j-1
ę ś ę
ę ś
0 ś ę M ś (41)
ę ś ę ś ę ś
M ę ś M
= -
ę ś ę ś ę ś
ę 0 ś
M M M
ę ś ę ś ę ś
ę
l2 ś ęwm-1, ś ęwm-1, ś
ęwm-1, j+1ś
ę ś
j j-1
ę l2 2(1- l2)ś
0 0
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 tj-1 czyli w(j-1).
108
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:
śu
0 < x < l
(x,0) = g(x) dla (42)
śt
Zamienimy pochodną cząstkową po czasie wykorzystując wzór schematu jawnego dla
równań parabolicznych:
śu u( xi,t1) - u(xi,0) k ś2u
~
~
m (0,t1)
(43)
(x,0) = - (xi, m )
j
śt k 2 śt2 j
Wyprowadzmy u(xi,t1) z powyższego równania, wtedy otrzymujemy:
śu k ś 2u
~ (44)
u(xi,t1) = u(xi,0) + k (x,0) + (xi, m j )
śt 2 śt2
Podstawmy do równania (44) wartość pochodnej z równania (42) :
k ś2u
~
u( xi,t1) = u(xi,0) + k g(x i ) + (xi, m )
2 ś t2 j (45)
Podstawmy za u( xi,tj ) = wi,j , pomijając błąd odcięcia, ostatecznie mamy:
wi,1 = wi,0 + k g(xi )
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(xi,t1) w szereg Maclurina 2 rzędu:
l2
wi,1 = (1- l2)wi,0 + (wi+1 ,0 + wi-1, 0)+ k g(xi )
dla i = 1,2, .. m-1 (47)
2
a
l = k
gdzie:
h
:
109
Równania hiperboliczne - Falowe - przykład 1
Rozwiązać poniższe równanie falowe:
d 2u d 2u
2
0 < x <1 t > 0
dla oraz
(x,t ) -a (x,t) = 0
dt2 dx2
z warunkami brzegowymi dla
u( 0,t) = u (1,t) = 0 t > 0
ś u
i początkowymi oraz dla
u( x,0) = sin(p x) 0 < x <1
(x,0) = 0
śt
Obliczenia wykonać dla czasu t = 1 wykorzystując metodę MRS dla równań falowych,
2
przyjąć: h = 0.1 oraz k = 0.05 ,przyjąć a = 4
Rozwiązanie dokładne ma postać:
u(x,t) = sin(p x)cos(2p t)
Tabela wyników
xi u( xi , 1 ) wi,20 | u( xi , 1 ) wi,20|
k = 0.05
0,0 0 0
0,1 0,3090169944 0,3090169944 0
0,2 0,5877852523 0,5877852523 0
0,3 0,8090169944 0,8090169944 0
0,4 0,9510565163 0,9510565163 0
0,5 1,0000000000 1,0000000000
1,1102 10-16
0,6 0,9510565163 0,9510565163
1,1102 10-16
0,7 0,8090169944 0,8090169944 0
0,8 0,5877852523 0,5877852523 0
0,9 0,3090169944 0,3090169944
5,421 10-20
1,0 0 0
110
Interpretacja graficzna
111
Równania hiperboliczne kody zró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;
112
Wyszukiwarka
Podobne podstrony:
komputerowe projektowania nowoczesnych instalacji elektrycznych
Komputerowe projektowanie technologii Bartosz Drzazga
49 przyklad projektu elektryki
Komputerowe projektowanie biżuterii
Narzędzia Komputerowe w Projektowaniu i Symulacji Info
Projekt elektrownia wiatrowa
Katedra Opto Komputerowe Systemy Elektroniczne
Projekt elektrownia
Projekt elektryczny
Komputerowe projektowanie odzieży i wizualizacja 3D(1)
projekt sieci komputerowej
Projekt umowy dzierzawy gruntow rolnych pod elektrownie wiatrowe
więcej podobnych podstron