Metody numeryczne w elektrotechnice, Metody numeryczne w elektrotechnice


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:

W obrębie klasycznych metod numerycznych możemy wyróżnić m.in. takie zagadnienia jak:

(każdy wynik musi podlegać weryfikacji!)

Reprezentacja liczb w maszynie cyfrowej

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

0x01 graphic
, 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:

0x01 graphic
, może być zatem przechowywana w słowie o długości D + 1 > 5 bitów jako

0x01 graphic

Liczby rzeczywiste (zmiennopozycyjne):

0x01 graphic
, gdzie ½ < m < 1, m-mantysa, i-cecha

0x01 graphic

W maszynie cyfrowej mantysa zapisywana jest w t-bitach mt

0x01 graphic

w wyniku zaokrąglenia do t-bitów mantysy

0x01 graphic

Zmiennopozycyjna reprezentacja liczby rzeczywistej x oznaczana jest symbolem rd(x) i jest równa

0x01 graphic

0x01 graphic

!

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:

0x01 graphic
, gdzie 0x01 graphic

liczbę 2-t nazywamy dokładnością maszynową

przykład: Liczba 18,5 daje się przedstawić w postaci:

0x01 graphic

może być więc przechowywana w słowie o długości d + 1 > 11 bitów o t > 6 bitach mantysy:

0x01 graphic

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:

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:

0x01 graphic

0x01 graphic
0x01 graphic

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 0x01 graphic
, 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:

0x01 graphic
, dla i = 1, 2, 3, ...

(niewiadoma musi być po obu stronach równania)

0x01 graphic

0x01 graphic

0x01 graphic

0x01 graphic
, współczynnik przed funkcją (tutaj 1) musi być < 1 dlatego:

0x01 graphic

0x01 graphic

0x01 graphic
- iteracja i-ta

np. 0x01 graphic
takiego równania nie da się rozwiązać bo 7/2

ale można zrobić tak:

0x01 graphic

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ż 0x01 graphic
, to

0x01 graphic
, a ciąg po > 0

0x01 graphic
, przy i = 1, 2, 3, ... jest zbieżny do 1

Przykład:

rozwiązanie iteracyjne

0x01 graphic

0x01 graphic
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:

0x01 graphic

W przypadku dokładnych działań, gdy 0x01 graphic
i 0x01 graphic
, więc qi < 0 dla i = 1, 2, 3, ... jeśli pi>1/3. Oznacza to, że maleje wartość bezwzględna różnicy 0x01 graphic
odpowiednio dużego numeru iteracji i.

Dla obliczeń wykonywanych na maszynie cyfrowej, mamy:

0x01 graphic
, gdzie 0x01 graphic
, k = 1, 2, 3, ...

Dla liczby 0x01 graphic
uzyskamy 0x01 graphic

Jeżeli yi jest już dostatecznie dobrym przybliżeniem liczby 0x01 graphic
, 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

0x01 graphic

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:

0x01 graphic

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

0x01 graphic
dla j = 1, 2, ..., m

We wzorze tym, czynnikiem określającym wrażliwość yj na bezwzględną zmianę δxi jest 0x01 graphic

Analogiczny wzór można wyprowadzić dla przenoszenia się błędów względnych. Jeśli 0x01 graphic
dla j = 0,1, 2, ..., m i 0x01 graphic
dla i = 1, 2, ..., n to

0x01 graphic
, gdzie 0x01 graphic
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 0x01 graphic

0x01 graphic
= 0x01 graphic
, gdzie ε jest

największą z wartości ε1, ε2 i ε3, gdzie 0x01 graphic
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), 0x01 graphic
, 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)

0x01 graphic

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 0x01 graphic
0x01 graphic

wyrażenie to nazywamy wielomianem uogólnionym.

Wprowadzając macierz bazową

0x01 graphic
i macierz współczynników 0x01 graphic

mamy 0x01 graphic

warunek W(x0) = Y0, W(x1) = Y1, ..., W(xn) = Yn

można zapisać w postaci układu równań liniowych X · A = Y, gdzie

0x01 graphic
0x01 graphic

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

0x01 graphic

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ć:

0x01 graphic
dodatkowo musi być spełniony warunek:

0x01 graphic

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

0x01 graphic

Macierz X-1 dla bazy wielomianowej 0x01 graphic
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 0x01 graphic
, wówczas 0x01 graphic
, a macierz Φ przyjmuje postać

0x01 graphic

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

0x01 graphic

przyjmuje się funkcje bazowe w postaci

0x01 graphic

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:

0x01 graphic

współczynniki a0 ... an tego wielomianu wyznaczamy z równania:

X · A = Y, przy czym macierz X ma postać:

0x01 graphic

Macierz posiada tylko główną przekątną niezerową w związku z tym układ równań X · A = Y rozwiązuje się natychmiastowo

0x01 graphic

Można więc wielomian interpolacyjny Lagrange'a zapisać w postaci ułamka:

0x01 graphic

lub krócej

0x01 graphic
, 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

0x01 graphic

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

0x01 graphic

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 0x01 graphic

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

0x01 graphic

z którego wyznacza się wartości nieznanych współczynników a0, a1, a2, ..., an

0x01 graphic

Ostatecznie otrzymujemy I wzór interpolacyjny Newtona w postaci

0x01 graphic

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 0x01 graphic

0x01 graphic

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 0x01 graphic
.

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 0x01 graphic

Współczynniki a0, a1, a2, ..., an wyznaczamy jak poprzednio i dochodzimy do wzoru:

0x01 graphic

Przykład: Obliczyć T log p dla T = 1,65 z dokładnością do Δ2

0x01 graphic

0x01 graphic
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

0x01 graphic

Otrzymujemy zatem wielomian interpolacyjny w postaci

0x01 graphic

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

0x01 graphic
, gdzie i = 0, 1, ..., 2n

Czyli

0x01 graphic

Warunek interpolacji prowadzi do układu równań liniowych w postaci

0x01 graphic

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

0x01 graphic

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:

0x01 graphic

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:

0x01 graphic

Tworzymy macierz X:

0x01 graphic

Elementy macierzy X mnożymy przez 2/7, otrzymaną macierz transponujemy i obliczamy współczynniki wzoru interpolacyjnego

0x01 graphic

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:

Dokonując aproksymacji funkcji f(x) musimy rozwiązać dwa ważne problemy:

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)

0x01 graphic

Jako funkcje bazowe stosowane są:

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ą:

0x01 graphic

0x01 graphic

0x01 graphic

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ć

0x01 graphic

Przy czym współczynniki ai są tak określone, aby wyrażenie

0x01 graphic
0x01 graphic
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

0x01 graphic

0x01 graphic

gdzie Rj jest odchyleniem w punkcie xj.

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

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

otrzymujemy układ m+1 równań o niewiadomych ai zwany układem normalnym

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

0x01 graphic

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ć

0x01 graphic

gdzie

0x01 graphic
0x01 graphic
0x01 graphic

Aproksymacja wielomianowa

Jeżeli za funkcje bazowe przyjmiemy ciąg jednomianów (xi), i = 0, 1, ..., n to układ normalny przyjmuje postać

0x01 graphic
, k = 0, 1, ..., n

co po zmianie kolejności sumowania daje

0x01 graphic

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

0x01 graphic
0x01 graphic

przy czym 0x01 graphic

0x01 graphic
, gdzie: 0x01 graphic
, 0x01 graphic

0x01 graphic

0x01 graphic
0x01 graphic

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

0x01 graphic

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

0x01 graphic

przy czym

0x01 graphic
0x01 graphic

Analogicznie ciąg funkcyjny

0x01 graphic

nazywamy ortogonalnym na zbiorze punktów x0, x1, ..., xn jeśli

0x01 graphic
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

0x01 graphic

Załóżmy, że znamy n+1 równoodległych punktów xi (xi = x0+ih, i = 0, 1, ..., n). Za pomocą przekształcenia liniowego

0x01 graphic

przeprowadzimy te punkty w kolejne liczby całkowite od 0 do n poszukujemy ciągu wielomianów

0x01 graphic

(dolny indeks oznacza stopień i m ≤ n) spełniających warunek ortogonalności

0x01 graphic
dla j≠k

przy czym

0x01 graphic

gdzie

0x01 graphic

Często używa się unormowanego ciągu wielomianów spełniających warunek

0x01 graphic
, gdzie k = 0, 1, ..., m

0x01 graphic

Co po przekształceniu daje nam wzór na wielomiany Grama, zwane też wielomianami Czebyszewa stopni k = 0, 1, ..., m w postaci

0x01 graphic
, gdzie k = 0, 1, ..., m

Wzór aproksymacyjny oparty na wielomianach Grama ma postać:

0x01 graphic

gdzie,

0x01 graphic
oraz

0x01 graphic
0x01 graphic

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:

0x01 graphic

0x01 graphic

0x01 graphic
0x01 graphic

0x01 graphic
0x01 graphic

0x01 graphic
0x01 graphic

0x01 graphic
0x01 graphic

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:

0x01 graphic
0x01 graphic

Obliczamy ci i si oraz bi= ci / si i ze wzoru:

0x01 graphic

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 0x01 graphic

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:

0x01 graphic
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ć:

0x01 graphic

0x01 graphic

0x01 graphic

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:

0x01 graphic
, n<L

Współczynniki aj i bj wielomianu wyznaczamy tak, aby suma kwadratów różnic

0x01 graphic
była minimalna.

Korzystając z warunku ortogonalności otrzymujemy rozwiązanie układu normalnego

0x01 graphic

w postaci:

0x01 graphic

dla j=1, 2, …, n

0x01 graphic

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 0x01 graphic
, gdzie β = 0, 1, ..., n-1 można przybliżać wielomianem trygonometrycznym w postaci

0x01 graphic
, j - czynnik urojony

lub

0x01 graphic

gdzie

Θ = 1 0x01 graphic
dla n parzystych

Θ = 1 0x01 graphic
dla n nieparzystych

0x01 graphic

Aproksymacja za pomocą funkcji sklejonych stopnia 3 (funkcji giętych)

Każdą funkcję 0x01 graphic
można przedstawić w postaci kombinacyjnej liniowej

0x01 graphic
, a ≤ x ≤ b oraz gdzie 0x01 graphic
są funkcjami określonymi wzorem:

0x01 graphic

gdzie 0x01 graphic
, 0x01 graphic

Aproksymacja na przedziale <a,b>

Współczynniki ci dobieramy tak, aby odchylenie kwadratowe dane poniższym wzorem było jak najmniejsze:

0x01 graphic

w celu wyznaczenia minimum funkcji I=I(c-1 , c0 , … , cn+1) obliczamy pochodne:

0x01 graphic
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

0x01 graphic
j= -1, 0, 1, … ,n+1

Funkcje 0x01 graphic
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:

0x01 graphic
j = -1, 0, 1, … , n+1,

gdzie 0x01 graphic

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:

0x01 graphic

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:

0x01 graphic
, gdzie

0x01 graphic

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

0x01 graphic

Układ ten można zapisać także w postaci macierzowej A · X = B, gdzie

0x01 graphic

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.

0x01 graphic

Niewiadome można wtedy obliczyć 0x01 graphic
, aii ≠ 0, i = 0, 1, ..., n.

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

0x01 graphic

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

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

PRZYKŁAD:

0x01 graphic

zmienną x3 mamy daną wprost - x3 = 3

0x01 graphic

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ń

0x01 graphic

0x01 graphic
0x01 graphic

to można wykazać, że 0x01 graphic
, 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 0x01 graphic
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.

0x01 graphic

Za pomocą wzorów określających nowe współczynniki

0x01 graphic
i=2,3,...,n j=2,3,...,n+1

Jeśli 0x01 graphic
, to odejmiemy drugie równanie układu pomnożone przez 0x01 graphic
od i-tego równania układu i=3,4,...,n. W wyniku otrzymujemy układ:

0x01 graphic

A nowe współczynniki wyrażone są wzorami

0x01 graphic
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

0x01 graphic

0x01 graphic

Algorytm dochodzenia od układu trójkątnego do układu równań liniowych

0x01 graphic

Przykład

0x01 graphic
0x01 graphic

Obliczamy elementy macierzy C1

0x01 graphic

druki krok s=2

0x01 graphic
0x01 graphic

co daje nam macierz C1:

0x01 graphic

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

0x01 graphic

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.

0x01 graphic

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:

0x01 graphic

ponieważ

0x01 graphic

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.

0x01 graphic

0x01 graphic

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

0x01 graphic

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:

0x01 graphic

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ą

0x01 graphic

rozwiązaniem dokładnym jest X = [4,59, 7,49, 4,20, 3,44]

Ten sam układ w metodzie Gaussa - Siedla ma postać:

0x01 graphic

a kolejne przybliżenia:

0x01 graphic

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:

0x01 graphic

czyli mamy:

0x01 graphic

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:

0x01 graphic

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 0x01 graphic
, 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:

0x01 graphic

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:

0x01 graphic
, gdzie0x01 graphic
- oznacza moduł wektora λ

wówczas równanie charakterystyczne jest zbieżne gdy:

0x01 graphic

Można wykazać, że:

0x01 graphic

ostatecznie dla metod iteracyjnych można zapisać twierdzenie:

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

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 0x01 graphic
, ani też od wartości danego wektora Z, a wystarczy jedynie spełnienie jednego z następujących warunków:

0x01 graphic
norma wierszowa

0x01 graphic
norma kolumnowa

0x01 graphic
norma energetyczna

Kryterium „stopu” metod iteracyjnych

Obliczenia iteracyjne trwają do momentu kiedy zostanie spełniony warunek

0x01 graphic

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

0x01 graphic
, 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

0x01 graphic

który możemy zapisać w postaci wektorowej

f(x)=0

gdzie:

0x01 graphic
0x01 graphic

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 0x01 graphic
jest błędem pierwiastka przybliżonego 0x01 graphic
. Skoro istnieje ciągła pochodna funkcji f(x) możemy zapisać:

0x01 graphic

Zastępując błąd 0x01 graphic
przyrostem 0x01 graphic
i porównując prawą stronę powyższego wyrażenia do zera otrzymujemy równanie liniowe:

0x01 graphic

Wzór ten stanowi zapis rekurencyjny dla metody Newtona w postaci macierzowej.

Pochodna f'(x) jest macierzą Jakobiego

0x01 graphic

Przyjmując, że 0x01 graphic
jest macierzą nieosobliwą możemy równanie liniowe przekształcić do postaci:

0x01 graphic

stąd przyjmując dowolną wartość 0x01 graphic
otrzymujemy ciąg kolejnych przybliżeń pierwiastka równania 0x01 graphic
uzyskanych metodą Newtona.

Metoda iteracji

Dany jest układ równań nieliniowych z n niewiadomymi

0x01 graphic

który możemy zapisać w postaci wektorowej

x=g(x)

gdzie:

0x01 graphic
0x01 graphic

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 0x01 graphic
, w rezultacie otrzymujemy ciąg kolejnych przybliżeń 0x01 graphic
pierwiastka a równania.

Metoda siecznych

Rozpatrujemy układ równań nieliniowych w postaci

0x01 graphic

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:

0x01 graphic

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

0x01 graphic

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ń

  1. kwadratury z ustalonymi węzłami

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ę

0x01 graphic
, gdzie 0x01 graphic
są dowolnymi punktami pośrednimi

Jeżeli ciąg {SN} dla 0x01 graphic
jest zbieżny do tej samej granicy S przy każdym przedziale odcinka [a ; b] takim, że 0x01 graphic
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:

0x01 graphic

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

0x01 graphic

0x01 graphic
, gdzie E(f) - reszta kwadratury

0x01 graphic

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.

0x01 graphic

Całka oznaczona Riemanna

Interpretacją graficzną całki oznaczonej jest pole powierzchni pomiędzy osią OX a krzywą. Rysunek obok to ilustruje

0x01 graphic

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

0x01 graphic
(1) 0x01 graphic
(2)

Można wykazać, że suma ciągu S(f) jest zbieżna do całki oznaczonej I(f) dla 0x01 graphic

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

0x01 graphic

W celu przybliżonego obliczenia całki oznaczonej funkcji f(x) na skończonym przedziale 0x01 graphic
, 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ń

0x01 graphic
(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

0x01 graphic

gdzie:

0x01 graphic
dla wówczas=(0, 1, …, n)

n - oznacza również stopień wielomianu interpolacyjnego

0x01 graphic
(5), gdzie:

0x01 graphic
(6)

wówczas

0x01 graphic

0x01 graphic

ξ - dowolny punkt pośredni z przedziału (x,xi)

gdzie:

0x01 graphic
oraz 0x01 graphic
(9)

ostatecznie wzór kwadratury Newtona - Cotesa oraz reszta kwadratury

0x01 graphic
(10)

0x01 graphic
(11)

Twierdzenie 1

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

0x01 graphic

wówczas wzory kwadratur dla n parzystych mają postać

0x01 graphic

oraz dla nieparzystych mają postać:

0x01 graphic
(13)

Kwadratury z ustalonymi węzłami - wzór trapezów

Wyznaczmy wzór kwadratury dla n = 1

0x01 graphic

0x01 graphic

gdzie:

0x01 graphic

0x01 graphic

wówczas:

0x01 graphic

0x01 graphic
(16)

Ostatecznie otrzymujemy

0x01 graphic

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

0x01 graphic

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 0x01 graphic

wzór parabol (Simpsona)

0x01 graphic

wzór Bessela

0x01 graphic

Interpretacja graficzna

wzór parabol (Simpsona)

0x01 graphic

wzór Bessela

0x01 graphic

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

0x01 graphic

Jeżeli dla skończonego [a ; b] przedziału wybierzemy zbiór punktów węzłowych {x0, …, xn} takich, że :

0x01 graphic

gdzie:

0x01 graphic
dla i=(0, …, n)

oraz

0x01 graphic

W oparciu o Twierdzenie 1 można również wyznaczyć wzory kwadratur otwartych. Wówczas wzory kwadratur dla n parzystych mają postać

0x01 graphic

oraz dla n nieparzystych mają postać:

0x01 graphic

wzór prostokątów

0x01 graphic

0x01 graphic

0x01 graphic

0x01 graphic

Interpretacja geometryczna

wzór prostokątów

0x01 graphic

wzór trapezów

0x01 graphic

Kwadratury Newtona-Cotesa - porównanie

Przykładowe wyniki obliczeń całki

0x01 graphic

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

0x01 graphic

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

0x01 graphic

0x01 graphic
, wtedy 0x01 graphic
0x01 graphic

wówczas wyniki:

0x01 graphic

oraz oszacowanie błędu obliczeń:

0x01 graphic

rozwiązanie dokładne:

0x01 graphic

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

0x01 graphic
0x01 graphic

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:

0x01 graphic
, dla 0x01 graphic
, 0x01 graphic

mamy wtedy:

0x01 graphic

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:

0x01 graphic
(28)

gdzie błąd przybliżenia 0x01 graphic
(29)

oraz 0x01 graphic

0x01 graphic

złożony wzór trapezów:

0x01 graphic
(30)

gdzie błąd przybliżenia 0x01 graphic
(31)

oraz 0x01 graphic

0x01 graphic

złożony wzór parabol - dla m parzystych:

0x01 graphic
(32)

0x01 graphic
(33)

gdzie błąd przybliżenia 0x01 graphic
(34)

oraz 0x01 graphic

0x01 graphic

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 0x01 graphic
(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: 0x01 graphic
, przy 0x01 graphic

inaczej 0x01 graphic

0x01 graphic
, 0x01 graphic
, 0x01 graphic

Przykład:

Obliczyć poniższą całkę:

Rozwiązanie dokładne

0x01 graphic

Metoda Simpsona m=1, h=2

0x01 graphic

błąd = -3,17143

Metoda Simpsona m=1, h=1

0x01 graphic

błąd = -0,26570

Metoda Simpsona m=4, h=0,5

0x01 graphic

błąd = -0,01807

Przykładowe wyniki obliczeń całki:

0x01 graphic

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.

0x01 graphic
0x01 graphic

Kwadratury Gaussa dobierają takie węzły, aby osiągnąć optymalne przybliżenie całki: 0x01 graphic

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:

0x01 graphic

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ż:

0x01 graphic
(64)

Zgodnie z założeniami, powyższe przybliżenie będzie dokładne tylko gdy stopień funkcji f(x) mniejszy lub równy 0x01 graphic
wówczas

0x01 graphic
(65)

gdzie a0, a1, a2, a3 - dowolne stałe

a więc możemy zapisać:

0x01 graphic

Poszukujemy x1, x2 oraz c1, c2, więc można zapisać:

0x01 graphic
(67)

0x01 graphic
(68)

0x01 graphic
(69)

0x01 graphic
(70)

otrzymaliśmy układ liniowych równań algebraicznych, którego rozwiązanie jest:

0x01 graphic
0x01 graphic
0x01 graphic
0x01 graphic

ostatecznie otrzymaliśmy wzór na przybliżenie całki dla funkcji danej wzorem (64):

0x01 graphic

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:

0x01 graphic

Dla naszych dalszych rozważań, zajmiemy się tylko ciągiem wielomianów Legendre'a

0x01 graphic
, który posiada właściwości:

1. Dla każdego n, Pn(x) jest wielomianem stopnia n

2. 0x01 graphic
, gdzie P(x) jest wielomianem stopnia wyższego (niższego???) od n

Wzory kilku pierwszych wielomianów Legendre'a:

0x01 graphic
0x01 graphic
0x01 graphic
0x01 graphic

oraz wzór rekurencyjny do wyznaczania wielomianów krzywych ??? oraz ich pochodnych:

0x01 graphic
0x01 graphic

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:

0x01 graphic

wówczas jeżeli P(x) jest dowolnym wielomianem stopnia mniejszego od 2n wtedy:

0x01 graphic
(75)

reszta dla stopnia większego od (2n-1): 0x01 graphic

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 0x01 graphic
. W celu obliczenia przybliżenia całki w dowolnym przedziale 0x01 graphic
dla którego f(x) jest ciągła należy wykonać przekształcenie poprzez podstawienie.

Dokonujemy przekształcenia dowolnego przedziału 0x01 graphic
do przedziału 0x01 graphic
.

0x01 graphic
0x01 graphic

Wówczas możemy zapisać:

0x01 graphic

po zastosowaniu kwadratury Gaussa-Lagrange'a otrzymujemy

0x01 graphic

gdzie xi - pierwiastek wielomianu stopnia n

ci - współczynnik ze wzoru (74)

przykład:

Obliczyć całkę wykorzystując kwadraturę Gaussa:

0x01 graphic

dokończyć ??? najpierw przekształcenia przedziału całkowania:

0x01 graphic

po zastosowaniu kwadratury Gaussa-Lagrange'a, n=2 otrzymujemy:

c1=1 c2=1 0x01 graphic
0x01 graphic

0x01 graphic

Po zastosowaniu kwadratury Gaussa-Lagrange'a, n=3 otrzymujemy:

0x01 graphic
0x01 graphic
0x01 graphic
0x01 graphic
0x01 graphic
0x01 graphic

0x01 graphic

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 0x01 graphic
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:

0x01 graphic

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:

0x01 graphic
, gdzie f(x):

przypadek 1*

0x01 graphic

przypadek 2*

0x01 graphic

0x01 graphic

0x01 graphic

0x01 graphic

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:

0x01 graphic

Metoda adaptacyjna:

0x01 graphic

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

Całki wielokrotne:

rozpatrzmy całkę podwójną:

0x01 graphic

x - stała i stosunek kwadratury do obliczonej całki

0x01 graphic

otrzymujemy:

0x01 graphic

m - liczba przedziałów przedziału [c(x) ; d(x)]

wykorzystanie z kolei kwadratury Simpsona do obliczenia całki:

0x01 graphic

ostatecznie:

0x01 graphic

np.: 0x01 graphic

stos. złoż. kwadratury Simpsona:

dla n=4 i m=2

0x01 graphic

0x01 graphic

0x01 graphic

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

0x01 graphic
(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

0x01 graphic
(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:

0x01 graphic
(131)

gdzie P jest prawdopodobieństwem, przy czym błąd metody monte Carlo można określić wzorem

0x01 graphic
(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:

0x01 graphic
(135)

Każda ze zmiennych xi zawiera się w przedziale 0x01 graphic
(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:

0x01 graphic
dla k=(1, 2, ..., n), 0x01 graphic
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:

0x01 graphic
(138),

gdzie

0x01 graphic
(139)

jest objętością m-wymiarowej kostki. Przybliżona wartość In całki (135) określana wzorem:

0x01 graphic
(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ą:

Numeryczne obliczanie pochodnej

Aby obliczyć pochodną funkcji f w punkcie x0, mamy

0x01 graphic
(1)

W celu wyznaczenia przybliżonej wartości pochodnej f'(x) dla małych wartości h możemy użyć wyrażenia

0x01 graphic
(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]

0x01 graphic
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:

0x01 graphic
(3)

oraz po uproszczeniu

0x01 graphic
(4)

0x01 graphic

Jeżeli dokonamy interpolacji funkcji f wielomianem stopnia 2, wówczas

0x01 graphic
(5)

oraz po uproszczeniu

0x01 graphic
I(6)

Za pomocą wielomianów wyższych stopni możemy wyznaczyć wzory na przybliżenie pochodnej.

0x01 graphic

Przykład: Obliczyć pochodną funkcji f(x)=ln(x) dla x0=1,8 stosując metodę w przód i wstecz stopnia 1.

0x01 graphic
0x01 graphic
,

gdzie:

0x01 graphic

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 0x01 graphic

dodając oba równania stronami otrzymujemy:

wzór (9)

Rozwiązując równanie ze względu na f”(x0), otrzymujemy:

0x01 graphic
(11)

oraz

0x01 graphic
, gdzie 0x01 graphic

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

0x01 graphic
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 różnicowe wielokrokowe:

Metoda Eulera

Rozpatrzmy zadanie znalezienia funkcji y=y(t), które dla 0x01 graphic
spełnia równanie różniczkowe zwyczajne pierwszego rzędu:

0x01 graphic
(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.

0x01 graphic
, gdzie 0x01 graphic
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ć:

0x01 graphic
(11), gdzie 0x01 graphic

Oznaczamy h=ti+1-ti, wówczas otrzymujemy:

0x01 graphic
(12)

Użyjemy podstawienia y'(t)=f(t, y), wówczas otrzymujemy:

0x01 graphic
(13)

Zapisując ωi≈y(ti) oraz pomijając błąd przybliżenia otrzymujemy ωi+1i+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ć:

0x01 graphic
(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ć:

0x01 graphic
, gdzie 0x01 graphic
(16)

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

0x01 graphic
(17)

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

0x01 graphic
, gdzie 0x01 graphic

Globalny błąd dyskretyzacji g(x)

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

Dla wybranego punktu ti możemy zapisać:

gii-y(ti) (20), wówczas

0x01 graphic
, gdzie L- liczba Lipschnitz`a

Lokalny błąd dyskretyzacji

0x01 graphic

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+1i+hf(ti, ωi), gdzie f(t,y)=y-t2+1

ωi+1i+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

0x01 graphic

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

0x01 graphic

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

0x01 graphic

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:

0x01 graphic
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 = 0x01 graphic
lub inaczej wi+1 - wi = h0 (ti,wi,h) (26)

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

0x01 graphic

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, α11, 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)

0x01 graphic

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)

0x01 graphic

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

0x01 graphic
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:

  1. równania różniczkowe eliptyczne