1. Liczby rzeczywiste, zmiennopozycyjne w maszynie

Reprezentacja liczb w maszynie cyfrowej

Liczby całkowite (stałopozycyjne = stałoprzecinkowe):

, gdzie e_i = 0 lub 1

Jeżeli rejestr ma d-bitów, wówczas liczba całkowita n może zawierać się w przedziale

-2^d-1 < n < 2^d-1-1

Liczby rzeczywiste (zmiennopozycyjne):

, gdzie ½ < m < 1, m-mantysa, i-cecha

W maszynie cyfrowej mantysa zapisywana jest w t-bitach m_t

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ą

2. 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. Błędy tego typu powstają też często przy obliczaniu wielkości będących granicami.

• 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.

3. Uwarunkowanie zadania i stabilność algorytmów.

Uwarunkowanie zadania

Jeśli zadanie rozwiązujemy numerycznie to zamiast dokładnymi wartościami x₁, x₂, ... x_n posługujemy się reprezentacjami rd (x₁), rd (x₂), ..., rd (x_n). Operacje elementarne nie są więc realizowane dokładnie tylko przez odwzorowanie zastępcze fl. Ta niewielka zmiana danych powoduje, że różne algorytmy rozwiązania tego samego zadania dają na ogół niejednakowe wyniki. Dużą rolę odgrywa tu przenoszenie błędów zaokrągleń. Często niewielkie zmiany danych powodują duże względne zmiany rozwiązania zadania. Zadanie takie nazywamy źle uwarunkowanym. Wielkości charakteryzujące wpływ zaburzeń danych na zaburzenia rozwiązania nazywamy wskaźnikami uwarunkowania zadania. Jeśli wskaźniki uwarunkowania co do wartości bezwzględnej są duże to zadanie jest źle uwarunkowane.

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.

4. 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(x₀) = Y₀, W(x₁) = Y₁, ..., W(x_n) = Y_n

interpolacja funkcji f(x)

Zadaniem interpolacji jest wyznaczenie przybliżonych wartości funkcji zwanej funkcją interpolową w punktach nie będących węzłami interpolacji. Przybliżoną wartość funkcji obliczamy za pomocą funkcji zwanej funkcją interpolującą, która w węzłach ma te same wartości co funkcja interpolowana.

Interpolacja wielomianowa

Baza dla funkcji ciągłych na odcinku skończonym [x_o, x_n] jest bazą zamkniętą, tzn., że każda funkcja tej klasy może być przedstawiona w postaci szeregu złożonego z funkcji bazowych. Wielomian interpolacyjny ma w tym przypadku postać:

dodatkowo musi być spełniony warunek:

Powyższy układ równań ma jedyne rozwiązanie względem a₁, jeżeli wartości x₀, x₁, ..., x_n są od siebie różne

Macierz X^-1 dla bazy wielomianowej
nazywana jest macierzą Lagrange'a

Interpolacja Lagrange'a

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 φ₁ brakuje czynnika (x-x_i). Zatem wielomian interpolacji wyraża się następującym wzorem:

współczynniki a₀ ... a_n 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

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 x_i∈[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 x₀ = 0. Przedstawiony układ równań rozwiązuje się natychmiastowo, ponieważ macierz X^-1 można wyznaczyć z zależności

5. Aproksymacja

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

Określenie dokładności aproksymacji

Aproksymacja funkcji powoduje powstanie błędów i sposób ich oszacowania wpływa na wybór metody aproksymacji. Jeśli błąd będzie mierzony na dyskretnym zbiorze punktów x₀, x₁, ..., x_n to jest oto aproksymacja punktowa, jeśli będzie mierzony w przedziale (a,b) - aproksymacja integralna lub przedziałowa.

Najczęściej stosowanymi miarami błędów aproksymacji są:

• dla aproksymacji średniokwadratowej punktowej

• 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 a₀, a₁, ..., a_m) 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 x₀, x₁, ..., x_n przyjmuje wartości y₀, y₁, ..., y_n. Poszukujemy wielomianu uogólnionego Q(x) będącego najlepszym przybliżeniem średniokwadratowym funkcji f(x) na zbiorze X=(x_j). Funkcja przybliżająca ma postać

Przy czym współczynniki a_i 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 a_i oznaczamy odchylenie

gdzie R_j jest odchyleniem w punkcie x_j.

Następnie obliczamy pochodne cząstkowe funkcji H względem a_i. Z warunku

, gdzie k = 0, 1, ..., n

otrzymujemy układ m+1 równań o niewiadomych a_i 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 (xⁱ), i = 0, 1, ..., n to układ normalny przyjmuje postać

, k = 0, 1, ..., n

co po zmianie kolejności sumowania daje

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

Przybliżeniem funkcji f(x) na zbiorze punktów x_i jest wielomian trygonometryczny:

, n<L

Współczynniki a_j i b_j 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

6. 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ć
, a_ii ≠ 0, i = 0, 1, ..., n.

Łatwo rozwiązuje się trójkątne układy równań

Z ostatniego równania możemy wyznaczyć x_n, z przedostatniego x_n-1

a ogólnie
, i = n-1, n-2, ..., 1 przy założeniu, że a_ii ≠ 0, i = 1, 2, ..., n

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ń.

Metody iteracyjne

Proces iteracyjnego obliczania wrtości x polega na przyjęciu pewnej wartości początkowej X^o, zwanej przybliżeniem początkowym, a następnie wykonaniu określonych operacji arytmetycznych dających w wyniku X¹ (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ń {Xⁿ} 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...

Metoda iteracji prostej (Jakobiego):

Metoda ta dla równania X=W*X+Z polega na przyjęciu zerowego przybliżenia wektora X=X^o, a następnie przeprowadzenia obliczeń iteracyjnych za pomocą zależności:

Xⁱ⁺¹=W*Xⁱ+Z i=0,1,...

Liczba kroków, które należy wykonać, aby uzyskać wymaganą dokładność rozwiązania, jest z reguły dość duża i istotnie zależy od przyjęcia punktu startowego X^o. Punkt ten najczęściej się dobiera na podstawie dodatkowych informacji o fizycznych aspektach problemu.

Metoda Gaussa-Seidla

Stanowi ulepszenie metody iteracyjnej prostej polegające na wykorzystaniu obliczonych k pierwszych składowych wektora Xⁱ⁺¹ do obliczenia składowej k+1. Modyfikacja ta znacznie przyspiesza proces obliczania.

Macierz W w równaniu X = W · X + Z przedstawiamy jako sumę dwóch macierzy: górnej trójkątnej W_u i dolnej trójkątnej W_l. Macierz górna trójkątna składa się z elementów macierzy W leżących nad główną przekątną (pozostałe są zerami), natomiast macierz dolna trójkątna składa się z elementów macierzy W leżących pod główną przekątną. Macierz W w tym przykładzie ma zerową główną przekątną

Algorytm metody Gaussa - Siedla realizowany jest następującym wzorem:

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 u_i po prawej stronie układu wartości

x_i + α (u_i - x_i) α - współczynnik nadrelaksacji, α∈<1,2>

i - liczba niewiadomych a nie iteracja

Sposób dobory optymalnego współczynnika α nazywanego czynnikiem nadrelaksacyjnym jest trudny do określenia, chociaż wiadomo, że jest on liczbą z przedziału [1 , 2]. Można zauważyć, że dla α=1 otrzymuje się metodę Gaussa - Siedla. Według źródeł literaturowych zaleca się przyjmowanie α rzędu 1,8.

Algorytm metody nadrelaksacyjnej sprowadza się do wzoru:

Zbieżność metod iteracyjnych

Warunkiem koniecznym i wystarczającym zbieżności procesu iteracyjnego danego wzorem przy dowolnym wektorze początkowym i dowolnym wektorze Z jest, aby wszystkie wartości własne macierzy W były co do modułu mniejsze od jedności.

W praktyce jest na ogół wygodniej korzystać z twierdzenia:

Jeżeli norma macierzy W jest mniejsza od jedności to proces iteracyjny określany wzorem X^i+1 = W · 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

7. 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(x₁...x_n) dla i = 1, 2, ..., n zakładamy, że mają ciągłe pochodne pierwszego rzędu w pewnym obszarze zawierającym odosobniony pierwiastek układu równań. Niech

x^(k) = {x₁^(k), ..., x_n^(k)}

będzie k-tym przybliżeniem pierwiastka a = {a₁, ..., a_n} 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 g₁, g₂, ..., g_n są rzeczywiste i ciągłe w pewnym otoczeniu odosobnionego pierwiastka a = {a₁, ..., a_n} układu równań.

Metoda iteracji polega na stosowaniu następującego wzoru

x^k+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: x₀, f(x₀), x₁, f(x₁)

- wykrywamy jeden wektor x-ów (globalny) (jeden zestaw x-ów), a pierwiastki mogą być wielokrotne

8. Całkowanie numeryczne

Całkowanie pojedyncze - metody obliczeń

1. kwadratury z ustalonymi węzłami

• kwadratury Newtona-Cotesa

• złożone kwadratury Newtona-Cotesa

• metody Romberga

• metoda adaptacyjna

b) kwadratury Gaussa

- kwadratury Gaussa

- złożone kwadratury Gaussa

Całki pojedyncze właściwe

Całka oznaczona Riemanna

Niech funkcja f(x) będzie określona i ograniczona na przedziale [a,b]. Dokonajmy podziału przedziału [a,b] takiego, że:

a = x₀ < x₁ < ... < x_n+1 = b

i utwórzmy sumę

, gdzie
są dowolnymi punktami pośrednimi

Jeżeli ciąg {S_N} dla
jest zbieżny do tej samej granicy S przy każdym przedziale odcinka [a ; b] takim, że
niezależnie od wyboru punktów C_i to funkcję f(x) nazywamy całkowalną w przedziale [a ; b], a granicę S ciągu (1) nazywamy:

Całką oznaczoną Riemanna funkcji f o oznaczamy:

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 x_i∈<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 x₀, x1, …, x_n. Wówczas całkę możemy zastąpić sumą, współczynniki a_i 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 {x₀, ..., x_n} takich, że:

x_i = x₀ + i · h

Twierdzenie 1

Kwadratury Newtona-Cotesa oparte na n+1 węzłach są rzędu

Zauważmy, że dla powyższych rozwiązań końce przedziału całkowania są węzłami kwadratury x₀=a oraz x_n=b. Wzory kwadratur oparte na tych węzłach nazywamy wzorami zamkniętymi Newtona - Cotesa.

Podobnie wyznaczyć można wzory dla innych n.

Kwadratury Newtona-Cotesa - wzory zamknięte

Poniżej znajdują się wzory zamknięte dla n = (1..3):

wzór trapezów

wzór parabol (Simpsona)

Wzór Bessla

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 {x₀, …, x_n} takich, że :

gdzie:

dla i=(0, …, n) oraz

wzór prostokątów

Wzór trapezów:

Wzór parabol:

Wzór Bessla:

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.

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:

Istnieje kilka metod obliczania całek funkcji osobliwych, jedną z nich jest Metoda adaptacyjna w której dokonuje się modyfikacji polegającej na pominięciu małego otoczenia punktu osobliwego np.

zał. dokładne 0,0001 liczba przedzi. 235 il. wywoł. funkcji 941 ???

Całki wielokrotne:

Całkę podwójną obliczamy poprzez dwukrotne zastosowanie kwadratury złożonej Simpsona.

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:

n- liczba przedziałów przedziału [a,b]

Całka potrójna: obliczamy tak samo jak podwójną tylko stosujemy kwadraturę Simpsona 3-krotnie.

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 x₁, x₂, ..., x_n,

Przybliżoną wartość całki określa wzór

(130)

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.

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)

Całkowanie - Metody Monte Carlo. Całki wielokrotne właściwe

Dana jest całka wielokrotna:

(135)

Każda ze zmiennych x_i zawiera się w przedziale
(136)

Losujemy n-krotnie punkty Pk kostki, a następnie 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ść I_n 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ę x_i^(k) w przedziale [a_i; b_i] korzystając ze zależności:

x_i^(k)=(b_i-a_i)y_i^(k)+a_i

9. Równania różniczkowe I rzędu

Numeryczne obliczanie pochodnej

Aby obliczyć pochodną funkcji f w punkcie x₀, 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 x₀ oraz x₀+h pod warunkiem, że funkcja f jest dwukrotnie różniczkowalna w przedziale [a; b]

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.

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)=y₀ (8),

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.

Zapisując ω_i≈y(t_i) oraz pomijając błąd przybliżenia otrzymujemy

ω_i+1=ω_i+hf(t_i, ω_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
- błąd przybliżenia (16)

Lokalny błąd dyskretyzacji τ_i+1(h)

(17)

Globalny błąd dyskretyzacji g(x)

g(t)= ω(t)-y(t) (19)

dodatkowo możemy określić krok h, dla którego błąd lokalny jest mniejszy od zadanej dokładności δ

, gdzie

Metoda Eulera - wstecz

Powyższe rozważania dotyczyły metody Eulera w przód, ponieważ krok spełniał warunek h>0. W analogiczny sposób można wyprowadzić metodą Eulera wstecz przyjmując h<0, wówczas otrzymujemy:

w_i+1 = w_i + hf (t_i+1, w_i+1) (22)

w_i+1^(k) = w_i + hf (t_i+1, w_i+1^(k-1)) (23)

Metoda wstecz różni się w stosunku do metody w przód argumentami funkcji f.

Metoda w przód wykorzystuje do obliczenia przybliżenia wartości z poprzedniego kroku, natomiast metoda wstecz jest równaniem uwikłanym, ponieważ aby obliczyć kolejne przybliżenie w_i+1 wykorzystujemy wartości z poprzedniego kroku oraz poszukiwaną wartość w_i+1. Takiego równania nie można rozwiązać w sposób bezpośredni. Aby rozwiązać takie równanie (23) należy zastosować proces iteracyjny, czyli poszukujemy kilkakrotnie w_i+1^(k), stojącej po prawej stronie równania podstawiając jako w_i+1^(k-1) - lewa strona równania, wynik przybliżenia z poprzedniej iteracji (k-1). Proces trwa do momentu, kiedy spełniony zostanie warunek:

|w_i+1^(k) - w_i+1^(k-1)| ≤ ε (24)

gdzie ε - tolerancja obliczeń

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 t_i+1 oraz t_i w postaci:

w_i+1 - w_i =
lub inaczej w_i+1 - w_i = hϕ (t_i,w_i,h) (26)

gdzie m oznacza rząd metody, c_j są stałymi, a

gdzie α_j, β_j, γ_j, _δlj - stałe

h - krok całkowania

Określenie stałych c_j, α_j, β_j otrzymujemy poprzez rozwinięcie funkcji f(t,y) w szereg Taylora w otoczeniu punktu t_i

Metody R-K rzędu 4

Metoda wyprowadzona przez rozwinięcie f(t,y) w szereg Taylora 4 rzędu, pozwala określić stałe we wzorze (wzór ogólny R-K). Poniżej przedstawiono najczęściej stosowaną metodę 4 rzędu

k₁ = hf (t_i, w_i)

k₂ = hf (t_i + ½ *h, w_i + ½ *k₁)

k₃ = hf (t_i + ½ *h, w_i + ½ *k₂) (37)

k₄ = hf (t_i+h, w_i+k₃)

ostatecznie:

w_i+1=w_i + 1/6 (k₁ + 2k₂ + 2k₃ + k₄) (38)

Interpretacja graficzna:

f₁ = f (t_i, w_i)

f₂ = f (t_i + ½ * h, w_i+1 + h f₁)

f₃ = f (t_i + ½ * h, w_i + ½ * h f₂)

f₄ = f (t_i + h, w_i + h f₃)

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

10. Układy równań różniczkowych 1. Rzędu

Układ równań różniczkowych zwyczajnych pierwszego rzędu z warunkami początkowymi. Mamy układ równań różniczkowych:

dla t należy (a,b)

Z warunkami początkowymi: u₁(a)=α₁ ; u₂(a)=α₂, … u_m(a)=α_m

Należy obliczyć przybliżenie m=funkcji u_1,u₂, … u_m­ które spełniają każde z równań różniczkowych.

Aby obliczyć taki układ należy zastosować dowolną z poznanych wcześniej metod np. Rungego-Kutty 4 rzędu.

Jeżeli podzielimy przedział t należy [a,b] na N przedziałów gdzie
oraz t_j=a+jh dla j=(1,2,…,N).

Użyjemy notacji w_ij dla każdego i=(1,2,…,m) oraz j=(0,1,…,N) w celu oznaczenia przybliżenia u_i(t­_j­).

Warunki początkowe oznaczamy: w₁₀= α w₂₀= α₂ w_m0= α_m

Wówczas słuszne są wzory:

k_1,i=hf_i(t_j,w_1j,w_2j,…,w_mj)

k_2,i=hf_i(t_j+1/2*h,w_ij+1/2*k₁₁;w_2j+1/2*k₁₂,…, w_mj+1/2*k_1m)

k_3,i=hf_i(t_j+1/2*h,w_ij+1/2*k₂₁;w_2j+1/2*k₂₂,…, w_mj+1/2*k_2m)

k_4,i=hf_i(t_j+h,w_ij+k₃₁;w_2j+k₃₂,…, w_mj+k_3m)

w_i,j+1=w_ij+1/6(k_1i+2k_2i+2k_3i+k_4i)

11.Równania różniczkowe wyższych rzędów

Rozwiązujemy poprzez sprowadzenie równania różniczkowego wyższego rzędu do układu równań różniczkowych rzędu I.

Równanie różniczkowe m-rzędu ma postać

dla t należy [a,b]

Warunki początkowe (jest ich m-1):

y(a)=α₁ y'(a)=α₂ , …., y_­^(m-1)(a)=α_m

W celu dokonania konwersji rozwiązania do układu równani oznaczamy:

u₁(t)=y(t), u₂(t)=y'(t), … , u­_m(t)=y^(m-1)(t)

oraz

Z warunkami początkowymi u₁(a)=y(a)=α₁

13. 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 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.

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 w_ij jest przybliżeniem u(x_i,y_j).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.

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.

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, …

26

---