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):
, gdzie ei = 0 lub 1
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:
, może być zatem przechowywana w słowie o długości D + 1 > 5 bitów jako
Liczby rzeczywiste (zmiennopozycyjne):
, gdzie ½ < m < 1, m-mantysa, i-cecha
W maszynie cyfrowej mantysa zapisywana jest w t-bitach mt
w wyniku zaokrąglenia do t-bitów mantysy
Zmiennopozycyjna reprezentacja liczby rzeczywistej x oznaczana jest symbolem rd(x) i jest równa
! |
Gdy elementy macierzy są rzędu μ lub n musimy dokonać przekształcenia całego układu. Uciekamy się od operacji na małych liczbach. |
Nie prowadzi się operacji na małych liczbach, co oznacza, że:
, gdzie
liczbę 2-t nazywamy dokładnością maszynową
przykład: Liczba 18,5 daje się przedstawić w postaci:
może być więc przechowywana w słowie o długości d + 1 > 11 bitów o t > 6 bitach 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. Źródłem 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 równego szeregowi:
zastępuje je sumą częściową o odpowiednio dobranej wartości 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:
gdzie εi są co do modułu niewiększe niż dokładność maszyny ε. Przedrostek fl oznacza, że są to wyniki działań wyznaczone przy użyciu maszyny cyfrowej.
Oszacowanie błędów w obliczeniach iteracyjnych:
Wiele algorytmów obliczeniowych polega na wyznaczeniu ciągu liczb zbieżnego do poszukiwanej wartości.
Przykład:
wyznaczyć wartość funkcji
, 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:
, dla i = 1, 2, 3, ...
(niewiadoma musi być po obu stronach równania)
, współczynnik przed funkcją (tutaj 1) musi być < 1 dlatego:
- iteracja i-ta
np.
takiego równania nie da się rozwiązać bo 7/2
ale można zrobić tak:
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ż
, to
, a ciąg po > 0
, przy i = 1, 2, 3, ... jest zbieżny do 1
Przykład:
rozwiązanie iteracyjne
wszystkie składniki przy niewiadomych są < 1 (7/8, -7/8, 3/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:
W przypadku dokładnych działań, gdy
i
, więc qi < 0 dla i = 1, 2, 3, ... jeśli pi>1/3. Oznacza to, że maleje wartość bezwzględna różnicy
odpowiednio dużego numeru iteracji i.
Dla obliczeń wykonywanych na maszynie cyfrowej, mamy:
, gdzie
, k = 1, 2, 3, ...
Dla liczby
uzyskamy
Jeżeli yi jest już dostatecznie dobrym przybliżeniem liczby
, 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
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 ϕ: 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 źle uwarunkowanym. Wielkości charakteryzujące wpływ zaburzeń danych na zaburzenia rozwiązania nazywamy wskaźnikami uwarunkowania zadania.
Przykład:
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
dla j = 1, 2, ..., m
We wzorze tym, czynnikiem określającym wrażliwość yj na bezwzględną zmianę δxi jest
Analogiczny wzór można wyprowadzić dla przenoszenia się błędów względnych. Jeśli
dla j = 0,1, 2, ..., m i
dla i = 1, 2, ..., n to
, gdzie
nazywamy wskaźnikiem uwarunkowania. Jeśli wskaźniki uwarunkowania co do wartości bezwzględnej są duże to zadanie jest źle uwarunkowane.
Przykład:
Badamy uwarunkowanie sumy y = a + b + c
=
, gdzie ε jest
największą z wartości ε1, ε2 i ε3, gdzie
jest wskaźnikiem uwarunkowania zadania.
Stabilność numeryczna algorytmów
Algorytm jest stabilny, jeżeli posiada tę własność, że małe błędy powstałe w pewnym kroku algorytmu nie są powiększane w innych krokach oraz nie powodują poważnego ograniczenia dokładności wszystkich obliczeń. Oznacza to, że algorytm jest numerycznie stabilny, gdy zwiększając dokładność obliczeń można wyznaczyć (z dowolną dokładnością) istniejące rozwiązanie zadania. Numeryczną stabilność zadania łatwo sprawdzić rozwiązując zadanie raz z dokładnością np. 10-6, a potem z dokładnością 10-12.
INTERPOLACJA:
Dana jest funkcja y = f (x),
, 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 liniowej n+1 funkcji bazowych
wyrażenie to nazywamy wielomianem uogólnionym.
Wprowadzając macierz bazową
i macierz współczynników
mamy
warunek W(x0) = Y0, W(x1) = Y1, ..., W(xn) = Yn
można zapisać w postaci układu równań liniowych X · A = Y, gdzie
Jeśli macierz X nie jest osobliwa (da się odwrócić), to:
A = X-1·Y
co ostatecznie daje
W(x) = Φ(x) · X-1· Y
Interpolacja wielomianowa
W praktyce często używa się bazy złożonej z jednomianów
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ć:
dodatkowo musi być spełniony warunek:
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
Macierz X-1 dla bazy wielomianowej
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 sprowadzić do zbioru podstawowego podstawiając
, wówczas
, a macierz Φ przyjmuje postać
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
przyjmuje się funkcje bazowe w postaci
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:
współczynniki a0 ... an tego wielomianu wyznaczamy z równania:
X · A = Y, przy czym macierz X ma postać:
Macierz posiada tylko główną przekątną niezerową w związku z tym układ równań X · A = Y rozwiązuje się natychmiastowo
Można więc wielomian interpolacyjny Lagrange'a zapisać w postaci ułamka:
lub krócej
, j = 0, 1, ..., n
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
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
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),
gdzie
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ń:
z którego wyznacza się wartości nieznanych współczynników a0, a1, a2, ..., an
Ostatecznie otrzymujemy I wzór interpolacyjny Newtona w postaci
Przykład: Dla zależności f(T)=T log p znaleźć wielomian interpolacyjny stopnia 3 i obliczyć odchyłki w węzłach interpolacji
T |
0,8 |
1,0 |
1,2 |
1,4 |
1,6 |
1,8 |
T log p |
0,3566 |
0,9383 |
1,5598 |
2,2169 |
2,9059 |
3,6229 |
Tablica:
T |
T log p |
Δ |
Δ2 |
Δ3 |
błąd [%] |
0,8 |
0,3566 |
0,5817 |
0,0398 |
-0,0042 |
0 |
1,0 |
0,9383 |
0,6215 |
0,0356 |
-0,0037 |
0 |
1,2 |
1,5598 |
0,6571 |
0,0319 |
-0,0039 |
0 |
1,4 |
2,2169 |
0,6890 |
0,0280 |
|
0 |
1,6 |
2,9059 |
0,7170 |
|
|
1,9 |
1,8 |
3,6229 |
|
|
|
4,2 |
przyjmijmy
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
.
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),
gdzie
Współczynniki a0, a1, a2, ..., an wyznaczamy jak poprzednio i dochodzimy do wzoru:
Przykład: Obliczyć T log p dla T = 1,65 z dokładnością do Δ2
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
Otrzymujemy zatem wielomian interpolacyjny w postaci
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
, gdzie i = 0, 1, ..., 2n
Czyli
Warunek interpolacji prowadzi do układu równań liniowych w postaci
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
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:
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:
Tworzymy macierz X:
Elementy macierzy X mnożymy przez 2/7, otrzymaną macierz transponujemy i obliczamy współczynniki wzoru interpolacyjnego
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)
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
dla aproksymacji średniokwadratowej integralnej lub przedziałowej
dla aproksymacji jednostajnej
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.
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 ϕj(x), j=0, 1,...,n będzie układem funkcji bazowych. Poszukujemy wielomianu uogólnionego Q(x) będącego najlepszym przybliżeniem średniokwadratowym funkcji f(x) na zbiorze X=(xj). Funkcja przybliżająca ma postać
Przy czym współczynniki ai są tak określone, aby wyrażenie
dla i=0, 1, ..., n
było minimalne. Funkcja w(x) jest z góry ustaloną funkcją wagową.
Aby wyznaczyć współczynniki ai oznaczamy odchylenie
gdzie Rj jest odchyleniem w punkcie xj.
Następnie obliczamy pochodne cząstkowe funkcji H względem ai. Z warunku
, gdzie k = 0, 1, ..., n
otrzymujemy układ m+1 równań o niewiadomych ai zwany układem normalnym
, gdzie k = 0, 1, ..., n
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ć
gdzie
Aproksymacja wielomianowa
Jeżeli za funkcje bazowe przyjmiemy ciąg jednomianów (xi), i = 0, 1, ..., n to układ normalny przyjmuje postać
, k = 0, 1, ..., n
co po zmianie kolejności sumowania daje
przykład: Odnaleźć zależność między x i y w postaci ax+by=1
x |
1 |
3 |
4 |
6 |
8 |
9 |
11 |
14 |
y |
1 |
2 |
4 |
4 |
5 |
7 |
8 |
9 |
Będziemy rozpatrywać odchylenia zarówno wartości x jak i y, ponieważ wiodą one do różnych wyników. Zakładamy, że dane są obarczone błędami, wówczas minimalizujemy wyrażenia
przy czym
, gdzie:
,
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
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
przy czym
Analogicznie ciąg funkcyjny
nazywamy ortogonalnym na zbiorze punktów x0, x1, ..., xn jeśli
dla j≠k
Zastosowanie tej metody powoduje, że znika jedna z trudności obliczeniowych przy aproksymacji wielomianowej, mianowicie złe uwarunkowanie macierzy układu normalnego. Przy aproksymacji wielomianami ortogonalnymi macierz układu normalnego jest macierzą diagonalną, a jej elementy położone na głównej przekątnej dane są wzorem
Załóżmy, że znamy n+1 równoodległych punktów xi (xi = x0+ih, i = 0, 1, ..., n). Za pomocą przekształcenia liniowego
przeprowadzimy te punkty w kolejne liczby całkowite od 0 do n poszukujemy ciągu wielomianów
(dolny indeks oznacza stopień i m ≤ n) spełniających warunek ortogonalności
dla j≠k
przy czym
gdzie
Często używa się unormowanego ciągu wielomianów spełniających warunek
, gdzie k = 0, 1, ..., m
Co po przekształceniu daje nam wzór na wielomiany Grama, zwane też wielomianami Czebyszewa stopni k = 0, 1, ..., m w postaci
, gdzie k = 0, 1, ..., m
Wzór aproksymacyjny oparty na wielomianach Grama ma postać:
gdzie,
oraz
Przykład:
Na podstawie tablicy wartości f(x) znaleźć najlepszy w sensie aproksymacji średniokwadratowej wielomian przybliżający tę funkcję.
x |
1 |
1,5 |
2 |
2,5 |
3 |
y |
3 |
4,75 |
7 |
9,75 |
13 |
Ponieważ węzły są równoodległe posłużymy się wielomianem Grama. Konstruujemy wielomian aproksymacyjny dla n = 4. Wartości Fk(q), q=1, 2,…, 4.
Obliczamy ze wzoru:
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:
Obliczamy ci i si oraz bi= ci / si i ze wzoru:
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 |
dla
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:
dla i = 0, 1, …, 2L-1
Baza jest ortogonalna nie tylko na przedziale <0, 2π>, ale też na zbiorze punktów xi, przy czym warunki ortogonalności mają postać:
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<L
Współczynniki aj i bj wielomianu wyznaczamy tak, aby suma kwadratów różnic
była minimalna.
Korzystając z warunku ortogonalności otrzymujemy rozwiązanie układu normalnego
w postaci:
dla j=1, 2, …, n
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;
Szybka transformata Fouriera
Każdą funkcję w tablicy tablicą wartości w n punktach
, gdzie β = 0, 1, ..., n-1 można przybliżać wielomianem trygonometrycznym w postaci
, j - czynnik urojony
lub
gdzie
Θ = 1
dla n parzystych
Θ = 1
dla n nieparzystych
Aproksymacja za pomocą funkcji sklejonych stopnia 3 (funkcji giętych)
Każdą funkcję
można przedstawić w postaci kombinacyjnej liniowej
, a ≤ x ≤ b oraz gdzie
są funkcjami określonymi wzorem:
gdzie
,
Aproksymacja na przedziale <a,b>
Współczynniki ci dobieramy tak, aby odchylenie kwadratowe dane poniższym wzorem było jak najmniejsze:
w celu wyznaczenia minimum funkcji I=I(c-1 , c0 , … , cn+1) obliczamy pochodne:
j = -1 , 0, … , n+1
i przyrównujemy je do zera. Otrzymane równania tworzą układ n+3 równań z n+3 niewiadomymi ci
j= -1, 0, 1, … ,n+1
Funkcje
są liniowo niezależne na przedziale <a,b>, więc układ równań ma rozwiązanie w postaci punktu. Funkcji I osiąga minimum. Układ równań możemy zapisać w postaci:
j = -1, 0, 1, … , n+1,
gdzie
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:
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:
, gdzie
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
Układ ten można zapisać także w postaci macierzowej A · X = B, gdzie
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.
Niewiadome można wtedy obliczyć
, aii ≠ 0, i = 0, 1, ..., n.
Łatwo rozwiązuje się trójkątne układy równań
Z ostatniego równania możemy wyznaczyć xn, z przedostatniego xn-1
a ogólnie
, i = n-1, n-2, ..., 1 przy założeniu, że aii ≠ 0, i = 1, 2, ..., n
PRZYKŁAD:
zmienną x3 mamy daną wprost - x3 = 3
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ń
to można wykazać, że
, W ≠ 0, j = 1, 2, ..., n
Metoda ta należy do metod dokładnych. Ze względu na dużą złożoność obliczeniową praktycznie stosowana do numerycznego rozwiązywania równań.
Metoda eliminacji Gaussa
Metoda polega na zapisaniu układu równań w postaci macierzy C, której n pierwszych kolumn zawiera elementy aij macierzy głównej A, natomiast kolumnę n+1 tworzą dowolne wyrazy bi.
Zakładając, że C11≠0, odejmujemy pierwsze równanie pomnożone przez
od i-tego równania (i=2,3,...,n) a obliczonymi współczynnikami zastępujemy przednie wartości. W wyniku tej operacji otrzymujemy układ równań i odpowiadającą mu macierz C1.
Za pomocą wzorów określających nowe współczynniki
i=2,3,...,n j=2,3,...,n+1
Jeśli
, to odejmiemy drugie równanie układu pomnożone przez
od i-tego równania układu i=3,4,...,n. W wyniku otrzymujemy układ:
A nowe współczynniki wyrażone są wzorami
i=3, 4, ..., n j=3, 4, ..., n
Kontynuując takie postępowanie po wykonaniu n kroków dochodzimy do układu trójkątnego
Algorytm dochodzenia od układu trójkątnego do układu równań liniowych
Przykład
Obliczamy elementy macierzy C1
druki krok s=2
co daje nam macierz C1:
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. 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...
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.
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:
ponieważ
więc jak łatwo sprawdzić mamy taki sam układ jaki znaleziono sposobem „naturalnym”
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=50°C, wymiary prostokątne a=b=0,8 cm.
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ń:
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:
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ą
rozwiązaniem dokładnym jest X = [4,59, 7,49, 4,20, 3,44]
Ten sam układ w metodzie Gaussa - Siedla ma postać:
a kolejne przybliżenia:
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:
czyli mamy:
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 α nazywanego czynnikiem nadrelaksacyjnym jest trudny do określenia, chociaż wiadomo, że jest on liczbą z przedziału [1 , 2]. Można zauważyć, że dla α=1 otrzymuje się metodę Gaussa - Siedla. Według źródeł literaturowych zaleca się przyjmowanie α rzędu 1,8.
Algorytm metody nadrelaksacyjnej sprowadza się do wzoru:
Zbieżność metod iteracyjnych
Jeżeli macierz A kwadratowa oraz x jest wektorem niezerowym, wówczas istnieje pewien wektor równoległy do Ax, co oznacza, że istnieje pewna stała λ taka, że:
Ax = λx
Możemy wtedy zapisać równanie charakterystyczne
(A - λI)x = 0, gdzie I - macierz jednostkowa
oraz wielomian charakterystyczny ma postać
p(λ) = det (A - λI)
Jeżeli wielomian macierzy A jest
, wektora x jest n, wówczas p jest wielomianem m-tego stopnia, którego pierwiastki są pojedyncze i zespolone. Jeżeli wektor λ jest pierwiastkiem wielomianu p wówczas:
Jeżeli powyższy warunek jest spełniony, oznacza to, że równanie charakterystyczne ma rozwiązanie x≠0.
Zbieżność metod iteracyjnych (2)
Promień spektralny ρ(A) macierzy A jest zdefiniowany jako:
, gdzie
- oznacza moduł wektora λ
wówczas równanie charakterystyczne jest zbieżne gdy:
Można wykazać, że:
ostatecznie dla metod iteracyjnych można zapisać twierdzenie:
Warunkiem koniecznym i wystarczającym zbieżności procesu iteracyjnego danego wzorem
|
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 wartości początkowej
, ani też od wartości danego wektora Z, a wystarczy jedynie spełnienie jednego z następujących warunków:
norma wierszowa
norma kolumnowa
norma energetyczna
Kryterium „stopu” metod iteracyjnych
Obliczenia iteracyjne trwają do momentu kiedy zostanie spełniony warunek
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
, gdzie ε ∈ <10-1, 10-256> - przyjęta dokładność obliczeń
Numeryczne rozwiązywanie równań nieliniowych
Metoda Newtona
Dany jest układ równań nieliniowych z n niewiadomymi
który możemy zapisać w postaci wektorowej
f(x)=0
gdzie:
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)
Gdzie
jest błędem pierwiastka przybliżonego
. Skoro istnieje ciągła pochodna funkcji f(x) możemy zapisać:
Zastępując błąd
przyrostem
i porównując prawą stronę powyższego wyrażenia do zera otrzymujemy równanie liniowe:
Wzór ten stanowi zapis rekurencyjny dla metody Newtona w postaci macierzowej.
Pochodna f'(x) jest macierzą Jakobiego
Przyjmując, że
jest macierzą nieosobliwą możemy równanie liniowe przekształcić do postaci:
stąd przyjmując dowolną wartość
otrzymujemy ciąg kolejnych przybliżeń pierwiastka równania
uzyskanych metodą Newtona.
Metoda iteracji
Dany jest układ równań nieliniowych z n niewiadomymi
który możemy zapisać w postaci wektorowej
x=g(x)
gdzie:
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ń.
Metoda iteracji polega na stosowaniu następującego wzoru
xk+1 = g(x(k))
Po przyjęciu przybliżenia
, w rezultacie otrzymujemy ciąg kolejnych przybliżeń
pierwiastka a równania.
Metoda siecznych
Rozpatrujemy układ równań nieliniowych w postaci
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:
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
Napisać program wyznaczający rozwiązanie tego układu.
Niezbędne deklaracje typów:
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.
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ń
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ę
, gdzie
są dowolnymi punktami pośrednimi
Jeżeli ciąg {SN} dla
jest zbieżny do tej samej granicy S przy każdym przedziale odcinka [a ; b] takim, że
niezależnie od wyboru punktów Ci to funkcję f(x) nazywamy całkowalną w przedziale [a ; b], a granicę S ciągu (1) nazywamy:
Całką oznaczoną Riemanna funkcji f o oznaczamy:
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
, 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.
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∈<a,b> oraz współczynników niezależnych od rodzaju funkcji.
Wówczas
(1)
(2)
Można wykazać, że suma ciągu S(f) jest zbieżna do całki oznaczonej I(f) dla
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
, 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ń
(4)
Kwadratury z ustalonymi węzłami
Kwadratury Newtona-Cotesa
Całkowanie z zastosowaniem kwadratur o ustalonych węzłach polega na zastąpieniu funkcji podcałkowej wielomianem interpolacyjnym Lagrange'a Li(x).
Jeżeli dla skończonego przedziału [a,b] wybierzemy zbiór punktów węzłowych {x0, ..., xn} takich, że:
xi = x0 + i · h
gdzie:
dla wówczas=(0, 1, …, n)
n - oznacza również stopień wielomianu interpolacyjnego
(5), gdzie:
(6)
wówczas
ξ - dowolny punkt pośredni z przedziału (x,xi)
gdzie:
oraz
(9)
ostatecznie wzór kwadratury Newtona - Cotesa oraz reszta kwadratury
(10)
(11)
Twierdzenie 1
Kwadratury Newtona-Cotesa oparte na n+1 węzłach są rzędu
|
wówczas wzory kwadratur dla n parzystych mają postać
oraz dla nieparzystych mają postać:
(13)
Kwadratury z ustalonymi węzłami - wzór trapezów
Wyznaczmy wzór kwadratury dla n = 1
gdzie:
wówczas:
(16)
Ostatecznie otrzymujemy
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
wzór parabol (Simpsona)
wzór Bessela
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
Jeżeli dla skończonego [a ; b] przedziału wybierzemy zbiór punktów węzłowych {x0, …, xn} takich, że :
gdzie:
dla i=(0, …, n)
oraz
W oparciu o Twierdzenie 1 można również wyznaczyć wzory kwadratur otwartych. Wówczas wzory kwadratur dla n parzystych mają postać
oraz dla n nieparzystych mają postać:
wzór prostokątów
Interpretacja geometryczna
wzór prostokątów
wzór trapezów
Kwadratury Newtona-Cotesa - porównanie
Przykładowe wyniki obliczeń całki
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 |
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
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
, wtedy
wówczas wyniki:
oraz oszacowanie błędu obliczeń:
rozwiązanie dokładne:
Kody źródłowe:
przykładowa całkowana funkcja:
function f(x:extender):extender,
begin
result:=sin (2*x*x)+2*cos(3x);
end;
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;
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:
, dla
,
mamy wtedy:
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:
(28)
gdzie błąd przybliżenia
(29)
oraz
złożony wzór trapezów:
(30)
gdzie błąd przybliżenia
(31)
oraz
złożony wzór parabol - dla m parzystych:
(32)
(33)
gdzie błąd przybliżenia
(34)
oraz
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
(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:
, przy
inaczej
,
,
Przykład:
Obliczyć poniższą całkę:
Rozwiązanie dokładne
Metoda Simpsona m=1, h=2
błąd = -3,17143
Metoda Simpsona m=1, h=1
błąd = -0,26570
Metoda Simpsona m=4, h=0,5
błąd = -0,01807
Przykładowe wyniki obliczeń całki:
Wyniki obliczeń dla wzorów złożonych w zależności od m:
m |
|
2 |
4 |
8 |
16 |
32 |
64 |
wzór prostokątów |
wynik |
1,110720735 |
1,751785441 |
1,936297295 |
1,983970758 |
1,995986709 |
1,998996147 |
|
błąd |
0,889279265 |
0,248214559 |
0,066102705 |
0,016029242 |
0,004013791 |
0,001003853 |
wzór trapezów |
wynik |
1,570796327 |
1,896118898 |
1,1974231602 |
1,993570344 |
1,998393361 |
1,999598389 |
|
błąd |
0,429203673 |
0,103881102 |
0,025768398 |
0,006429656 |
0,001606639 |
0,000401611 |
wzór parabol |
wynik |
2,094395102 |
2,004559755 |
2,000269170 |
2,000016591 |
2,000001033 |
2,000000065 |
|
błąd |
-0,094395102 |
-0,004559755 |
-0,000269170 |
-0,000016591 |
0,000001033 |
0,000000065 |
Kody źródłowe:
Złożony wzór prostokątów:
function Calka_Prostokat_Zlozony(a,b:extended;m:integer):extended;
var h,x,I_0:extended;
begin
h:=(b-a)/m;
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;
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:
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.
Mamy znaleźć 2n stałe we wzorze, przypuszczamy więc, że rząd kwadratury powinien wynosić 2n, a więc wzór powinien być dokładny dla wielomianów stopnia (2n-1). Dla wielomianu stopnia równego lub wyższego od (2n-1) należy uwzględnić resztę kwadratury E:
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ż:
(64)
Zgodnie z założeniami, powyższe przybliżenie będzie dokładne tylko gdy stopień funkcji f(x) mniejszy lub równy
wówczas
(65)
gdzie a0, a1, a2, a3 - dowolne stałe
a więc możemy zapisać:
Poszukujemy x1, x2 oraz c1, c2, więc można zapisać:
(67)
(68)
(69)
(70)
otrzymaliśmy układ liniowych równań algebraicznych, którego rozwiązanie jest:
ostatecznie otrzymaliśmy wzór na przybliżenie całki dla funkcji danej wzorem (64):
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:
Dla naszych dalszych rozważań, zajmiemy się tylko ciągiem wielomianów Legendre'a
, który posiada właściwości:
1. Dla każdego n, Pn(x) jest wielomianem stopnia n
2.
, gdzie P(x) jest wielomianem stopnia wyższego (niższego???) od n
Wzory kilku pierwszych wielomianów Legendre'a:
oraz wzór rekurencyjny do wyznaczania wielomianów krzywych ??? oraz ich pochodnych:
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:
wówczas jeżeli P(x) jest dowolnym wielomianem stopnia mniejszego od 2n wtedy:
(75)
reszta dla stopnia większego od (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
. W celu obliczenia przybliżenia całki w dowolnym przedziale
dla którego f(x) jest ciągła należy wykonać przekształcenie poprzez podstawienie.
Dokonujemy przekształcenia dowolnego przedziału
do przedziału
.
Wówczas możemy zapisać:
po zastosowaniu kwadratury Gaussa-Lagrange'a otrzymujemy
gdzie xi - pierwiastek wielomianu stopnia n
ci - współczynnik ze wzoru (74)
przykład:
Obliczyć całkę wykorzystując kwadraturę Gaussa:
dokończyć ??? najpierw przekształcenia przedziału całkowania:
po zastosowaniu kwadratury Gaussa-Lagrange'a, n=2 otrzymujemy:
c1=1 c2=1
Po zastosowaniu kwadratury Gaussa-Lagrange'a, n=3 otrzymujemy:
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 źródłowy kwadratury Gaussa:
Przykład f. całką wykorzystując kwadraturę G w przedziale
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];
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:
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:
, gdzie f(x):
przypadek 1*
|
przypadek 2*
|
|
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:
Metoda adaptacyjna:
zał. dokładne 0,0001 liczba przedzi. 235 il. wywoł. funkcji 941 ???
Całki wielokrotne:
rozpatrzmy całkę podwójną:
x - stała i stosunek kwadratury do obliczonej całki
otrzymujemy:
m - liczba przedziałów przedziału [c(x) ; d(x)]
wykorzystanie z kolei kwadratury Simpsona do obliczenia całki:
ostatecznie:
np.:
stos. złoż. kwadratury Simpsona:
dla n=4 i m=2
Całkowanie wielokrotne - całka podwójna
kod źródłowy:
function f (x,y:real):real;
begin
Result:=explg (x)???
end;
function c(x:real):real;
begin
Result:=x*x*x
end;
function d(x:real):real;
begin
Result:=x*x
end;
function simpson-Podwojna
(a,b:real; m,n:integer):real;
var i,j:integer
h,hx,j1,j2,j3
k1,k2,k3,x,y,x:real;
begin
h:=(b-a)/(2*n)
j1:=0; j2:=0;
j3:=0
for i:=0 to 2*n do
begin
x:=a+i*h;
nx:=(d(x)-c(x))/(2*m);
k1:=f(x,c(x))+f(x,d(x));
k2:=0
k3:=0
for j:=1 to 2*m-1 do
begin
y:=c(x)+j*hx
z:=f(x,y);
if Odd(j) ???
then k3 :=k3+z
else k2:=k2+z
end
i:=(k1+2*k2+4*k3)+hx/3
if(i=0) or (i=2*n)
then j1:=j1+1
else if Odd (i)
then j3:=j3+l
else j2:=j2+l
end;
Result:=(j1+2*j2+4*h/3)
end;
Kody źródłowe:
Function f(x,y:real):real;
Begin
Result:=exp(y/x);
End;
Function c (x:real):real;
Begin
Result:=x*x*x;
End;
Function d (x:real):real;
Begin
Result:=x*x;
End;
function Simpson_Podwojna (a,b:real; m,n:integer):real;
var i,j:integer;
h,hx,j1,j2,j3,k1,k2,k3,I,x,y,z:real;
begin
h:=(b-a)/(2*n);
j1:=0;
j2:=0;
j3:=0;
for i:=0 to 2*n do
begin
x:=a+i*h; tu tomek napisał ... j*h a ja mam i*h sprawdźcie!!!
hx:=(d(x)-c(x))/(2*m)
k1:=f(x,c(x)+f(x,d(x));
k2:=0;
k3:=0;
for j:=1 to 2*m-1 do
begin
y:=c(x) +j*hx;
z:=f(x,y);
if odd(j)
then k3:=k3+z
else k2:=k2+z
end;
I:=(k1+2*k2+4*k3)*hx/3;
If (i=0) or (i=2*n)
then j1:=j1+I
else if Odd(i) tu też mam troche inaczej
then j3:=j3+l
else j2:=j2+l
end;
Result:=(j1+2*j2+4*j3)*h/3;
End;
Obliczanie całek metodą Monte-Carlo
Metoda Monte Carlo - całka pojedyncza
Mamy obliczyć przybliżoną wartość In całki oznaczonej
(129)
przy założenie, że f(x) jest funkcją ciągłą w przedziale domkniętym [a;b]
Następnie n-krotnie generujemy realizację x zmiennej losowej X o rozkładzie równomiernym w przedziale [a,b], otrzymujemy w rezultacie ciąg liczb x1, x2, ..., xn,
Przybliżoną wartość całki określa wzór
(130)
Zbieżność przybliżonej wartości całki In do jej dokładnej wartości jest gwarantowana w sensie probalistycznym, tzn., że dla każdego ε>0 zachodzi relacja:
(131)
gdzie P jest prawdopodobieństwem, przy czym błąd metody monte Carlo można określić wzorem
(132)
Maszyny typowe zazwyczaj dysponują generatorem liczb losowych o rozkładzie równomiernym w przedziale [0,1], zachodzi zatem konieczność przeliczania wylosowanej realizacji Y z przedziału [0,1], na realizację zmiennej losowej X o rozkładzie równomiernym w przedziale [a, b]. Wzory poniżej, to ilustrują:
Przeliczanie zmiennych losowych Y na X
X=(b-a)Y+a (133)
Przeliczanie realizacji zmiennych losowych Y na X
x=(b-a)y+a (134)
Kod źródłowy:
Function Calka_MonteCarlo (a,b:extended; n:lingint):extended
Var i:integer;
h,x:extended;
begin
Result:=0;
h:=(b-a)/n;
for :=1 to n do
begin
x:=a+((b-a)*Random);
Result:=Result+f(x);
End;
Result:=Result*h;
end;
Całkowanie - Metody Monte Carlo. Całki wielokrotne właściwe
Dana jest całka wielokrotna:
(135)
Każda ze zmiennych xi zawiera się w przedziale
(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:
dla k=(1, 2, ..., n),
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:
(138),
gdzie
(139)
jest objętością m-wymiarowej kostki. Przybliżona wartość In całki (135) określana wzorem:
(140)
Aby uzyskać ciąg realizacji X1, X2, ..., Xm losujemy m-krotnie zmienną losową Y o rozkładzie równomiernym w przedziale [0,1], a następnie przeliczamy realizacje y na realizację 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
Numeryczne obliczanie pochodnej
Aby obliczyć pochodną funkcji f w punkcie x0, mamy
(1)
W celu wyznaczenia przybliżonej wartości pochodnej f'(x) dla małych wartości h możemy użyć wyrażenia
(2)
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]
ktoś nie zdązył do końca, a ja nie wiem czy w mianowniku jest 2!, czy co??
Wówczas obliczając pochodną funkcji f, mamy:
(3)
oraz po uproszczeniu
(4)
Jeżeli dokonamy interpolacji funkcji f wielomianem stopnia 2, wówczas
(5)
oraz po uproszczeniu
I(6)
Za pomocą wielomianów wyższych stopni możemy wyznaczyć wzory na przybliżenie pochodnej.
Przykład: Obliczyć pochodną funkcji f(x)=ln(x) dla x0=1,8 stosując metodę w przód i wstecz stopnia 1.
,
gdzie:
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 |
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
dodając oba równania stronami otrzymujemy:
wzór (9)
Rozwiązując równanie ze względu na f”(x0), otrzymujemy:
(11)
oraz
, gdzie
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
f”(2)=-0,1875
h |
Wynik |
Błąd |
Oszacowanie ε |
0,1 |
-0,1875586182 |
0,0000586182 |
|
0,01 |
-0,1875005859 |
|
|
0,001 |
-0,1875000049 |
|
|
Równania różniczkowe zwyczajne z warunkami początkowymi
Równania różniczkowe są wykorzystywane do odwzorowania wielu problemów z zakresu różnych nauk, które dotyczą zmian jednej zmiennej względem innej zmiennej. Wiele problemów świata rzeczywistego opisywanych przez równania różniczkowe jest za bardzo złożonych, aby uzyskać rozwiązanie dokładne (analityczne), należy znaleźć rozwiązanie przybliżone. Są dwa sposoby aby uzyskać rozwiązanie przybliżone. Po pierwsze można dokonać uproszczenia zadanego równania różniczkowego do takiej postaci, dla której znamy rozwiązanie dokładne, wówczas uproszczone równanie będzie przybliżało rozwiązanie równania zadanego. Drugim sposobem jest zastosowanie jednej z metod numerycznych w celu przybliżenia rozwiązania zadanego równania.
Metody różnicowe jednokrokowe:
metody Eulera;
metody Rundego-Kutty
metody Rundego-Kutty-Fehenberga
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
spełnia równanie różniczkowe zwyczajne pierwszego rzędu:
(7)
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.
, gdzie
oraz y(a)=a (10)
Na rozwiązanie powyższego zagadnienia będziemy obliczać przybliżone wartości funkcji yi=y(ti), gdzie ti=a+ih, h=(b-a)/N, dla którego i=(0,1, ..., N), gdzie ti nazywane są punktami węzłowymi, natomiast h odległością między nimi.
Rozwiniemy y(t) w szereg Taylora w celu wyprowadzenia metody Eulera. Zakładając, że y(t) jest jedynym rozwiązaniem (10) oraz posiada drugą pochodną, która jest ciągła w przedziale [a, b] wówczas dla każdego i=1, 2, ..., N.
Wykorzystując powyższe założenia możemy zapisać:
(11), gdzie
Oznaczamy h=ti+1-ti, wówczas otrzymujemy:
(12)
Użyjemy podstawienia y'(t)=f(t, y), wówczas otrzymujemy:
(13)
Zapisując ωi≈y(ti) oraz pomijając błąd przybliżenia otrzymujemy ωi+1=ωi+hf(ti, ωi) (14)
Powyższy wzór nazywamy metodą Eulera - wzór ten nazywany jest inaczej równaniem różniczkowym, gdyż można zapisać:
(15)
Aby wyznaczyć wartość szukanej funkcji y(x) w następnym kroku h, wykorzystujemy poprzednią wartość funkcji oraz wielkości zmian funkcji - dzięki pochodnej. Natomiast uwzględniając błąd przybliżenia wzór (13) przyjmuje postać:
, gdzie
(16)
Lokalny błąd dyskretyzacji τi+1(h)
(17)
dodatkowo możemy określić krok h, dla którego błąd lokalny jest mniejszy od zadanej dokładności δ
, gdzie
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
, gdzie L- liczba Lipschnitz`a
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 0<t<2 y(0)=0,5
dla 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,2ωi-0,008i2+0,02 dla i=0, 1, ..., 9
Rozwiązanie dokładne ma postać
y(t)=(t2-1)-0,5et
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)| ≤ ε (24)
gdzie ε - tolerancja obliczeń
Kody źródłowe:
Metoda Eulera w przód:
Function Euler_order_One_Forward (t0, tn, y0:extended; n:integer): extended;
var i:integer;
h, wi, wi_1,t :extended;
Begin
h:= (tn-t0)/n;
yi_1:=y0;
for i:=0 to n-1 do
begin
t:=t0+i*h;
wi:=wi_1+h*f(t,wi_1);
wi_1:=wi;
end;
Result:=wi;
end;
Przykładowa funkcja definiująca równanie różniczkowe
function f(t,y:extended):extended;
begin
result:=y-sqr(t)+1;
end;
Metoda Eulera “wstecz”
function Euler_Order_One_Backward
t0, tn, y0, eps : extended; n:integer):extended;
var i:integer;
h,wi,wi_p,wi_1,t :extended
begin
h:=(tn-t0)/n
wi_1:=y0;
for i:=0 to n-1 do
begin
t:=t0 + i*h;
wi:=wi_1 + h*f (t, wi_1);
repeat
wi_p:=wi;
wi:=wi_1 + h*f (t, wi_p);
until abs(wi-wi_p)<eps;
wi_1:=wi;
end;
Result:=wi;
end;
Przykładowa funkcja definiująca równanie różniczkowe
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:
gdzie t∈ [a,b] oraz y(a)=α (25)
i przedstawieniu różnicy funkcji y(t) w punktach ti+1 oraz ti w postaci:
wi+1 - wi =
lub inaczej wi+1 - wi = h0 (ti,wi,h) (26)
gdzie m oznacza rząd metody, cj są stałymi, a
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, α1,β1, c2, α2, β2 :
Metoda punktu środkowego:
C1 = 0 C2 = 1 α1 = 0 α2 = h/2 β1 = 0 β2 = ½ (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))
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 = i⋅h ; 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)
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
Kody źródłowe:
function RK_Rzedu_4 (t0, tn, y0: extended; n:integer):extended;
var i:integer;
h, yi, y_i, t, k1, k2, k3, k4: extended;
begin
h:= (tn-t0)/n
y_i = y0;
for i:=0 to n-1 do
begin
t:=t0 + i*h;
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
gdzie t∈ [a,b] oraz y(a) = α (55)
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
Równania różniczkowe cząstkowe z warunkami brzegowymi - rodzaje:
równania różniczkowe eliptyczne
równanie Laplace'a,
równanie Poisson'a,
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
dla
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ć.
Równania różniczkowe cząstkowe - warunek brzegowy Neumanna
Warunek brzegowy Neumanna
dla
gdzie:
- pochodna normalna poszukiwanej funkcji w punktach należących do brzegu S obszaru R, g(x,y) - zadana funkcja dla punktów (x,y) należących do brzegu S obszaru R.
Oznacza to, że na brzegu obszaru zmienność funkcji u(x,y) w kierunku normalnym dla punktów (x,y) S jest równa funkcji g(x,y).
Dla zmiennej x
Dla zmiennej y
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ć:
(1)
na
z warunkami brzegowymi
dla
W powyższym równaniu funkcja f opisuje dane wejściowe na płaszczyźnie dla obszaru R ograniczonego brzegiem S. Równania tego typu są stosowane dla różnych problemów fizycznych, które są niezależne od czasu. Najczęściej stosowane są do obliczenia rozkładu potencjału (np. temperatury) w stanie ustalonym.
Szczególnym przypadkiem równania Poissona gdy f(x,y)=0 jest równanie Laplace'a.
Aby rozwiązać równanie cząstkowe eliptyczne (1) zastosujemy metodę różnic skończonych MRS.
W tym celu weźmiemy dwie liczby całkowite m i n, i zdefiniujemy kroki całkowania h = (b-a)/n oraz k = (c-d)/m. Dzięki temu możemy podzielić: przedział [a,b] na n równych części o szerokości h oraz przedział [c,d] na m równych części o szerokości k.
Umieśćmy siatkę na obszarze R poprzez narysowanie pionowych i poziomych linii przechodzących przez punkty (xi, yj), takich że:
dla i = 0,1, .. n
oraz
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
(2)
gdzie:
Dla zmiennej y
(3)
gdzie:
Podstawiając wzory (2) oraz (3) do równania Poissona (1) otrzymujemy:
(4)
dla i = 1,2 … n-1 oraz j = 1,2, … m-1
natomiast warunki brzegowe mają postać:
oraz
dla j = 0,1, .. m (5)
oraz
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:
(7)
dla i = 1,2 … n-1 oraz j = 1,2, … m-1
natomiast warunki brzegowe mają postać:
oraz
dla j = 0,1, .. m (8)
oraz
dla i = 1,2, .. n (9)
gdzie wij jest przybliżeniem u(xi,yj).
Analizując równanie (7) można zauważyć, że w celu wyznaczenia przybliżenie rozwiązania w punkcie (xi,yj), potrzebne są wartości przybliżenia rozwiązania w czterech sąsiednich punktach - patrz rysunek obok.
Postać macierzowa
Jak widać poszukujemy rozwiązania równania Poissona dla wewnętrznych węzłów siatki, w tym celu dla każdego punktu siatki układamy równanie (7). Po uwzględnieniu warunków brzegowych otrzymujemy układ (n-1)Ⴔ(m-1) równań liniowych. Układ równań możemy rozwiązać metodami deterministycznymi bądź iteracyjnymi.
Zakładając, że h = k możemy wyprowadzić uproszczoną postać równania (7) :
(10)
dla i = 1,2 … n-1 oraz j = 1,2, … m-1
natomiast warunki brzegowe są takie same jak (8) i (9) :
Postać iteracyjna
Równanie (7) i (10) można przedstawić również w postaci iteracyjnej, wyprowadzając wij z każdego równania, gdzie l oznacza numer iteracji :
z równania (7)
(11)
z równania (10)
(12)
Równania eliptyczne - przykład 1
Wyznaczyć rozkład temperatury w stanie ustalonym dla cienkiej kwadratowej metalowej płytki o wymiarach 0,5 m na 0,5 m. Na brzegu płytki znajdują się źródła ciepła utrzymujące temperaturę na poziomie 0°C dla boku dolnego i prawego, natomiast temperatura boku górnego i lewego zmienia się liniowo od 0°C do 100°C Problem rozwiązać układając układ równań liniowych (postać macierzowa) dla wewnętrznych węzłów siatki 5 x 5 - układ równań rozwiązać metodą Gaussa -Siedla.
Siatka dyskretyzacyjna 5 x 5
Problem ten opisuje równanie Laplace'a
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)
Rozwiązanie układu równań metodą Gaussa-Siedla, daje wyniki
Ostatecznie wyniki dla założonej siatki mają postać:
Rozwiązanie dokładne problemu ma postać:
Natomiast oszacowanie błędu obliczeń :
, bo
i
Wynika z tego, że rozwiązanie numeryczne jest równe z rozwiązaniu dokładnemu.
Równania eliptyczne - przykład 1 - kody źródłowe
Funkcje charakteryzujące dane równanie różniczkowe eliptyczne
function F( X,Y : real ) : real;
begin
F := 0;
end;
function G( X,Y : real ) : real;
begin
if X = 0 then G := 0;
if X = 0,5 then G := 200*Y;
if Y = 0 then G := 0;
if Y = 0,5 then G := 200*X;
end;
Definicja wymaganych stałych i typów
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 );
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;
Równania eliptyczne - przykład 2
Wyznaczyć rozkład temperatury w stanie ustalonym dla cienkiej kwadratowej metalowej płytki o wymiarach 1 m na 1 m. Na trzech brzegach płytki znajdują się źródła ciepła utrzymujące temperaturę na poziomie 100°C dla boku górnego, 50°C dla boku lewego, 0°C dla boku dolnego (warunki Dirichleta), natomiast przy czwartym lewym boku jest umieszczony izolator termiczny (warunek Neumanna - pochodna normalna równa jest 0). Problem rozwiązać wg postaci iteracyjnej dla wewnętrznych węzłów siatki 101 x 101. Tolerancję obliczeń przyjąć na poziomie 10-6.
Problem ten opisuje równanie Laplace'a
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 ]
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 |
Równania eliptyczne - przykład 2 - kody źródłowe
Funkcje charakteryzujące dane równanie różniczkowe eliptyczne
function F( X,Y : real ) : real;
begin
F := 0;
end;
Definicja wymaganych stałych i typów
const max = 100;
type
Macierz = array [ 0..max, 0..max ] of real;
Vector = array [ 0..max ] of real;
function PoissonEquationIteration (a,b,c,d, TOL : real ; N,M ,nmax : integer; var V : Macierz) : boolean;
var EPS,H,K,Vs : real; X,Y : Vector; I,J,L : integer;
begin
H := ( B - A ) / N; K := ( D - C ) / M;
for I := 0 to N do X[I] := A + I * H; for J := 0 to M do Y[J] := C + J * K;
for i:=0 to N do for j:=0 to M do V[i,j]:=0;
L := 0; Result := true;
repeat
eps:=0; L := L + 1;
for i:=0 to N do for j:=0 to M do begin
Vs:=V[i,j];
if j=0 then V[i,j]:=0 else // bok dolny
if j=M then V[i,j]:=100 else // bok gora
if i=0 then V[i,j]:= 50 else // bok lewo
if i=N then V[i,j]:= V[i-1,j] else // bok prawo
V[i,j]:=(- h*h* F(X[i],Y[j] ) + (V[i+1,j]+V[i-1,j])+sqr(h/k)*(V[i,j+1]+V[i,j-1]) )/(2*(sqr(h/k)+1));
If abs(V[i,j]-Vs) > eps then eps := abs(V[i,j]-Vs);
end;
until (eps< TOL ) or ( L > nmax);
if (L > nmax) then Result := false;
end;
Równania różniczkowe cząstkowe z warunkiem brzegowym
Równania paraboliczne - Dyfuzja
Równania różniczkowe cząstkowe paraboliczne znane jako równanie dyfuzji lub przewodnictwa, w zależności od jednego wymiaru oraz czasu przyjmuje postać:
dla oraz (13)
z warunkami brzegowymi dla
i początkowymi dla
Równania tego typu są stosowane dla różnych problemów fizycznych, które są zależne od czasu. Najczęściej stosowane są do obliczenia przepływu ciepła wzdłuż pręta o długości l, przy założeniu że, powierzchnia boczna pręta jest odizolowana od otoczenia. Stała ၡ jest niezależna od położenia oraz czasu i najczęściej określa przewodność cieplną materiału z którego zrobiony jest pręt. Funkcja f określa początkowy rozkład temperatury w pręcie.
Aby rozwiązać równanie cząstkowe paraboliczne (13) zastosujemy metodę różnic skończonych MRS.
W tym celu weźmiemy liczbę całkowitą m i zdefiniujemy krok całkowania h = (b-a)/m. Dodatkowo ustalimy wartość kroku czasowego k. Dzięki temu możemy wyznaczyć węzły siatki (xi , tj ), gdzie:
dla i = 0,1, .. m oraz dla j = 0,1, ..
Dla wewnętrznych punktów węzłowych (xi,tj), które nie należą do brzegu, zastosujemy metodę różnic skończonych - MRS. Metoda ta opiera się na zastąpieniu pochodnych cząstkowych rozwinięciem funkcji w szereg Taylora w otoczeniu punktu (xi,tj) - (patrz pierwsza i druga pochodna numeryczna) . Wówczas otrzymujemy:
Dla zmiennej x (14)
gdzie:
Dla zmiennej t (15)
gdzie:
Podstawiając wzory (14) oraz (15) do równania dyfuzji (13) oraz wyłączając oszacowanie błędu zapisujemy:
(16)
warunek brzegowy : (17)
dla i = 1,2 … m-1 oraz j = 1,2, …
natomiast lokalny błąd odcięcia jest równy :
(18)
gdzie: oraz
Ostatecznie wzór na MRS, podstawiając za ( xi , tj ) = wi,j oraz wyłączając oszacowanie błędu zapisujemy:
(19)
przekształcając powyższy wzór ze względu na wi,j+1 ,wówczas otrzymujemy:
(20)
lub
gdzie: (21)
dla i = 1,2 … m-1 oraz j = 1,2, …
Równania paraboliczne - Schemat jawny
Analizując równanie (21) można zauważyć, że w celu wyznaczenia przybliżenia rozwiązania w punkcie (xi,tj+1), potrzebne są wartości trzech punktów z poprzedniego kroku czasowego tj - patrz rysunek obok. Metoda ta nazywana jest rozwiązywaniem równania parabolicznego wg schematu jawnego
Jak widać poszukujemy rozwiązania równania dyfuzji dla wewnętrznych węzłów siatki, w tym celu dla każdego punktu siatki układamy równanie (21). Obliczenia rozpoczynamy od kroku czasowego j = 1, gdyż najpierw musimy wygenerować wartości początkowe dla t = 0 dla wszystkich punktów u( xi , 0). Obliczmy je zgodnie z warunkami początkowymi u( xi , 0) = f(xi).
A więc mamy do rozwiązania układ (m-1) równań, który powstał poprzez ułożenie równania (21) dla każdego punktu wi,j siatki dla i = 1,2, .. m-1 :
(22)
co rozpisując daje nam:
(23)
Czyli w ten sposób możemy wyznaczyć poszukiwane wartości w(j+1). Wystarczy tylko przemnożyć macierz A przez wektor w(j) - wartości dla poprzedniego kroku.
Równania paraboliczne - Dyfuzja - przykład 1
Rozwiązać poniższe równanie przewodnictwa:
dla oraz
z warunkami brzegowymi dla
i początkowymi dla
Obliczenia wykonać dla czasu t = 0.5 używając schematu jawnego najpierw dla h = 0.1 oraz k = 0.0005 , następnie dla h = 0.1 oraz k = 0.01 ,przyjąć ၡ 2 = 1
Rozwiązanie dokładne ma postać:
Tabela wyników - schemat jawny
|
|
Mały krok czasu |
Duży krok czasu |
||
xi |
u( xi , 0.5 ) |
wi,1000 k = 0.0005 |
| u( xi , 0.5 ) - wi,1000| |
wi,50 k = 0.01 |
| u( xi , 0.5 ) - wi,50| |
0,0 |
0 |
0 |
|
0 |
|
0,1 |
0,00222241 |
0,00228652 |
6,411 Ⴔ 10-5 |
8,19876 Ⴔ 107 |
8,199 Ⴔ 107 |
0,2 |
0,00422728 |
0,00434922 |
1,219 Ⴔ 10-4 |
-1,55719 Ⴔ 108 |
1,557 Ⴔ 108 |
0,3 |
0,00581836 |
0,00598619 |
1,678 Ⴔ 10-4 |
2,13833 Ⴔ 108 |
2,138 Ⴔ 108 |
0,4 |
0,00683989 |
0,00703719 |
1,973 Ⴔ 10-4 |
-2,50642 Ⴔ 108 |
2,506 Ⴔ 108 |
0,5 |
0,00719788 |
0,00739934 |
2,075 Ⴔ 10-4 |
2,62685 Ⴔ 107 |
2,627 Ⴔ 108 |
0,6 |
0,00683989 |
0,00703719 |
1,973 Ⴔ 10-4 |
-2,49015 Ⴔ 108 |
2,490 Ⴔ 108 |
0,7 |
0,00581836 |
0,00598619 |
1,678 Ⴔ 10-4 |
2,11200 Ⴔ 107 |
2,112 Ⴔ 108 |
0,8 |
0,00422728 |
0,00434922 |
1,219 Ⴔ 10-4 |
-1,53086 Ⴔ 108 |
1,531 Ⴔ 108 |
0,9 |
0,00222241 |
0,00228652 |
6,411 Ⴔ 10-5 |
8,03604 Ⴔ 107 |
8,036 Ⴔ 107 |
1,0 |
0 |
0 |
|
0 |
|
Interpretacja graficzna
Równania paraboliczne - Schemat jawny
Analizując wyniki z przykładu 1 można zauważyć, że metoda wykorzystująca schemat jawny jest warunkowo stabilna, gdyż dla zbyt dużych wartości kroku wyniki są rozbieżne w stosunku do rozwiązania dokładnego. Wynika to wartości własnej macierzy A , dla której powinien być spełniony warunek
(23)
z powyższego wyrażenia wynika, że warunek stabilności tej metody jest:
(24)
Wykorzystując powyższy wzór można wyznaczyć dla przykładu 1 graniczną wartość kroku czasowego, dla której schemat jawny jest stabilny:
wtedy
Równania paraboliczne - Schemat niejawny
Poprzednia metoda wykorzystująca schemat jawny, była warunkowo stabilna, gdyż należy zwracać uwagę na wielkości kroków całkowania. Metoda opisana poniżej nie ma tej wady, i jest bezwarunkowo stabilna dla różnych wielkości kroków całkowania. Metoda ta wykorzystuje schemat niejawny gdyż zalicz się do metod pośrednich wymagających wielokrotnych iteracji dla danej chwili czasu. Różnica polega na innym określeniu pochodnej cząstkowej po czasie:
Dla zmiennej t (25)
gdzie:
Wówczas podstawiając wzór (25) oraz (14) do równania parabolicznego (13), pomijając błąd obcięcia otrzymamy :
(26)
Modyfikując powyższy wzór, podstawiając za ( xi , tj ) = wi,j zapisujemy:
(27)
lub:
gdzie: (28)
dla i = 1,2 … m-1 oraz j = 1,2, …
Analizując równanie (28) można zauważyć, że w celu wyznaczenia przybliżenia rozwiązania w punkcie (xi,tj), potrzebne są wartości trzech punktów: dwóch sąsiednich z aktualnego kroku czasowego oraz wartość dla tego samego węzła, ale dla poprzedniej iteracji tj-1 - patrz rysunek obok. Metoda ta nazywana jest rozwiązywaniem równania parabolicznego wg schematu niejawnego
Jak widać poszukujemy rozwiązania równania dyfuzji dla wewnętrznych węzłów siatki, w tym celu dla każdego punktu siatki układamy równanie (28). Obliczenia rozpoczynamy od kroku czasowego j = 1, gdyż najpierw musimy wygenerować wartości początkowe dla t = 0 dla wszystkich punktów u( xi , 0). Obliczmy je zgodnie z warunkami początkowymi u( xi , 0) = f(xi).
A więc mamy do rozwiązania układ (m-1) równań, który powstał poprzez ułożenie równania (28) dla każdego punktu wi,j siatki dla i = 1,2, .. m-1 :
(29)
co rozpisując daje nam:
(30)
Czyli w ten sposób możemy wyznaczyć poszukiwane wartości w(j). Najlepiej do rozwiązania tego typu układów użyć rozkładu Crouta dla macierzy trójdiagonalnych
Równania paraboliczne - Dyfuzja - przykład 2 (1)
Rozwiązać poniższe równanie przewodnictwa:
dla oraz
z warunkami brzegowymi dla
i początkowymi dla
Obliczenia wykonać dla czasu t = 0.5 używając schematu niejawnego dla h = 0.1 oraz k = 0.01 ,przyjąć ၡ 2 = 1
Rozwiązanie dokładne ma postać:
Tabela wyników - schemat niejawny
xi |
u( xi , 0.5 ) |
wi,50 k = 0.01 |
| u( xi , 0.5 ) - wi,50| |
0,0 |
0 |
0 |
|
0,1 |
0,00222241 |
0,00289802 |
6,756 Ⴔ 10-4 |
0,2 |
0,00422728 |
0,00551236 |
1,285 Ⴔ 10-3 |
0,3 |
0,00581836 |
0,00758711 |
1,769 Ⴔ 10-3 |
0,4 |
0,00683989 |
0,00891918 |
2,079 Ⴔ 10-3 |
0,5 |
0,00719788 |
0,00937818 |
2,186 Ⴔ 10-3 |
0,6 |
0,00683989 |
0,00891918 |
2,079 Ⴔ 10-3 |
0,7 |
0,00581836 |
0,00758711 |
1,769 Ⴔ 10-3 |
0,8 |
0,00422728 |
0,00551236 |
1,285 Ⴔ 10-3 |
0,9 |
0,00222241 |
0,00289802 |
6,756 Ⴔ 10-4 |
1,0 |
0 |
0 |
|
Równania paraboliczne - Schemat niejawny
Analizując wyniki z przykładu 2 można zauważyć, że metoda wykorzystująca schemat niejawny jest bezwarunkowo stabilna, gdyż jest ona zawsze zbieżna nie zależnie od wielkości kroków całkowania. Wynika to ze spełniania poniższego warunku dla promienia spektralnego macierzy A:
(31)
Równania paraboliczne - kody źródłowe
Funkcja określająca warunki początkowe równania paraboliczne
function F( X,Y : real ) : real;
begin
F := sin( pi * x );
end;
Definicja wymaganych stałych i typów
const
max = 100;
type
Vector = array [ 0..max ] of real;
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];
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;
Równanie hiperboliczne - Falowe
Równania różniczkowe cząstkowe hiperboliczne znane jako równanie falowe, w zależności od jednego wymiaru oraz czasu przyjmuje postać:
(32)
dla oraz
z warunkami brzegowymi dla
i początkowymi oraz dla
Równania tego typu są stosowane dla różnych problemów fizycznych, które są zależne od czasu. Najczęściej stosowane są do obliczenia wartości fali płaskiej w dowolnym punkcie pomiędzy 0 a l dla dowolnej chwili czasu większej od zera. Funkcje f(x) i g(x) określają warunki początkowe położenia oraz prędkości rozchodzenia się fali.
Aby rozwiązać równanie cząstkowe paraboliczne (32) zastosujemy metodę różnic skończonych MRS.
W tym celu weźmiemy liczbę całkowitą m i zdefiniujemy krok całkowania h = (b-a)/m. Dodatkowo ustalimy wartość kroku czasowego k. Dzięki temu możemy wyznaczyć węzły siatki (xi , tj ), gdzie:
dla i = 0,1, .. m oraz dla j = 0,1, ..
Dla wewnętrznych punktów węzłowych (xi,tj), które nie należą do brzegu, zastosujemy metodę różnic skończonych - MRS. Metoda ta opiera się na zastąpieniu pochodnych cząstkowych rozwinięciem funkcji w szereg Taylora w otoczeniu punktu (xi,tj) - (patrz druga pochodna numeryczna) . Wówczas otrzymujemy:
Dla zmiennej x (33)
gdzie
Dla zmiennej t (34)
gdzie
Podstawiając wzory (33) oraz (34) do równania falowego (31) oraz wyłączając oszacowanie błędu zapisujemy:
(35)
warunek brzegowy : (36)
dla i = 1,2 … m-1 oraz j = 1,2, …
natomiast lokalny błąd odcięcia jest równy :
(37)
gdzie: oraz
Ostatecznie wzór na MRS, podstawiając za ( xi , tj ) = wi,j zapisujemy:
(38)
przekształcając powyższy wzór ze względu na wi,j+1 ,wówczas otrzymujemy:
gdzie: (39)
dla i = 1,2 … m-1 oraz j = 1,2, …
z warunkami brzegowymi dla j = 1,2, …
i początkowymi dla i = 1,2, … m-1
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 :
(40)
co rozpisując daje nam:
(41)
Czyli w ten sposób możemy wyznaczyć poszukiwane wartości w(j+1). Wystarczy tylko przemnożyć macierz A przez wektor w(j) - wartości dla poprzedniego kroku, a następnie odjąć od tego wektor wynikowy dla tj-1 czyli w(j-1).
Ponieważ obliczania rozpoczynamy od kroku czasowego j = 2 musimy oprócz warunków początkowych dla j = 0 określić wartości dla j = 1, skorzystamy z drugiego warunku początkowej prędkości fali:
dla (42)
Zamienimy pochodną cząstkową po czasie wykorzystując wzór schematu jawnego dla równań parabolicznych:
(43)
Wyprowadźmy u(xi,t1) z powyższego równania, wtedy otrzymujemy:
(44)
Podstawmy do równania (44) wartość pochodnej z równania (42) :
(45)
Podstawmy za u( xi,tj ) = wi,j , pomijając błąd odcięcia, ostatecznie mamy:
dla i = 1,2, .. m-1 : (46)
Powyższy wzór pozwala na określenie warunków początkowych dla j = 1, wadą tego jest duży błąd oszacowania rzędu O(k). Alternatywą do tego jest wzór wyprowadzony jako rozwiniecie funkcji u(xi,t1) w szereg Maclurina 2 rzędu:
dla i = 1,2, .. m-1 (47)
gdzie:
:
Równania hiperboliczne - Falowe - przykład 1
Rozwiązać poniższe równanie falowe:
dla oraz
z warunkami brzegowymi dla
i początkowymi oraz dla
Obliczenia wykonać dla czasu t = 1 wykorzystując metodę MRS dla równań falowych, przyjąć: h = 0.1 oraz k = 0.05 ,przyjąć α 2 = 4
Rozwiązanie dokładne ma postać:
Tabela wyników
Interpretacja graficzna
Równania hiperboliczne - kody źródłowe
Funkcja określająca warunki początkowe równania paraboliczne
function F( X,Y : real ) : real;
begin
F := sin( pi * x );
end;
function G( X,Y : real ) : real;
begin
G := 0;
end;
Definicja wymaganych stałych i typów
const
max = 100;
type
Vector = array [ 0..max ] of real;
Macierz = array [ 0..max ] of Vector;
Funkcja „WaveEquation” - oblicza przybliżenie rozwiązania równania falowego
function WaveEquation(FX,ALPHA ,FT : real; M,N : integer; var W : Macierz): boolean ;
var H,K,V : real; M1,N1,I,J : integer;
begin
M1 := M + 1; N1 := N + 1; Result := True;
H := FX / M; K := FT / N;
V := ALPHA * K / H;
for J := 1 to N do begin
W[0,J] := 0; W[M,J] := 0;
end;
W[0,0] := F( 0 ); W[M,0] := F ( FX );
for I := 1 to M-1 do begin
W[I,0] := F( H * ( I ) );
W[I,1] := (1-V*V)*F(H*(I))+V*V*(F((I+1)*H) + F(H*(I-1)))/2+K*G(H*(I));
end;
for J := 1 to N-1 do for I := 1 to M-1 do
W[I,J+1] := 2*(1-V*V)*W[I,J]+V*V*(W[I+1,J]+W[I-1,J])-W[I,J-1];
end;
28
0
0,8090169944
0,8090169944
0,3
0
0,3090169944
0,5877852523
0,8090169944
0,9510565163
1,0000000000
0,9510565163
0,5877852523
0,3090169944
0
u( xi , 1 )
5,421 × 10-20
0
0
1,1102 × 10-16
1,1102 × 10-16
0
0
0
| u( xi , 1 ) - wi,20|
0
0,3090169944
0,5877852523
0,8090169944
0,9510565163
1,0000000000
0,9510565163
0,5877852523
0,3090169944
0
wi,20
k = 0.05
0,5
0,6
0,7
0,8
0,9
1,0
0,4
0,2
0,1
0,0
xi