Procedury numeryczne

Metody numeryczne –

procedury

na podstawie [Marciniak et. al. 1997] oraz [Bronsztejn et. al. 2004]

dr inż. Paweł Zalewski Akademia Morska w Szczecinie

Interpolacja wielomianowa:

Zadanie interpolacji polega na znalezieniu pewnej funkcji, która „przybliża” daną funkcję f. Dla funkcji f znane są przy tym wartości f(x_i) lub wartości jej pochodnych w pewnej skończonej liczbie różnych punktów x_i, które nazywamy węzłami interpolacji.

Wartość wielomianu interpolacyjnego Lagrange’a L_n w punkcie x dla danej funkcji f oblicza się z wzoru:

 

_

 

_







x

x

x

x

x

f

x

L

Interpolacja wielomianowa:

Wartość wielomianu interpolacyjnego w punkcie x dla danej funkcji f można obliczyć również algorytmem Neville’a na podstawie zależności rekurencyjnych:

 













i

n

j

i

x

x

P

x

x

P

x

x

P

n

i

x

f

P

,...,

1

,

0

,

,...,

1

,

0

,

,...,

1

,

0

,



















gdzie wartości x_i oraz f(x_i) (i = 0, 1, ... , n) są dane. Wartość elementu P_n,n jest poszukiwaną wartością wielomianu w punkcie x.

Aproksymacja wielomianowa:

Przykładem jest aproksymacja średniokwadratowa wielomianem, która dla danej funkcji f(x) polega na określeniu wielomianu p(x) stopnia m, który w przypadku danych dyskretnych minimalizuje sumę:

gdzie wartości x_i oraz f(x_i) (i = 0, 1, ... , n) są dane, m < n, a w(x) oznacza funkcję wagową. Współczynniki wielomianu aproksymującego funkcję

f(x):





_

 

_

 













_



f

x

a

x

a

a

a

H

,

,...,

   





 





p

x

w

x

x

f

x

a

x

a

a





...



można wyznaczyć poszukując minimum poniższego funkcjonału przy danych wartościach x_i oraz f(x_i) :

Aproksymacja wielomianowa:

 

x

a

 

x

x

k

m

f

a

H

,...,

1

,

0

,

0

















_







_

_

 

i

k

a

y

k

m

z

,...,

1

,

0

,

,







Powyższy układ równań liniowych można rozwiązać np. metodą eliminacji Gaussa lub metodami iteracyjnymi (przybliżonymi).

Wprowadzając oznaczenia:

 

 

 

i

k

x

i

k

m

z

m

k

x

x

f

k

y

,...,

1

,

0

,

,

,

,...,

1

,

0

,













Całkowanie numeryczne:

Wiele całek nie można przedstawić za pomocą funkcji elementarnych. Do ich obliczenia stosuje się metody przybliżone (przykładowo dla obliczenia

pola powierzchni pod krzywą). W celu wyznaczenia przybliżonej wartości całki:

 

_

   

_

 







x

f

A

dx

x

f

x

w

dx

x

F

funkcję podcałkową F(x) przedstawia się w postaci iloczynu w(x) f(x), gdzie funkcja w(x), zwana funkcją wagową zazwyczaj wynosi 1, a funkcję

f(x) przybliża się wielomianem. Wówczas:

 















F

x

dx

a

b

b

a

,

gdzie współczynniki A_k zależą od zastosowanej metody numerycznej obliczeń.

Całkowanie numeryczne:

W metodzie Simpsona dany przedział [a, b] dzieli się na równe

podprzedziały. W przypadku trzech podprzedziałów:

  

a

,

x



x

,

x

 

,

x

,

x

 

,

x

,

x

  



x

,

b

 



 



  



2

,

1

,

0

,

2

4

3



















i

x

x

h

I

x

f

h

x

f

x

f

h

dx

x

f

w każdym podprzedziale stosuje się wzór Simpsona:

Następnie każdy z przedziałów dzieli się znowu na trzy części i otrzymuje się drugie przybliżenie całki. Procedura ta jest kontynuowana do chwili, gdy różnica pomiędzy dwoma kolejnymi przybliżeniami całki jest dostatecznie mała lub, gdy długości podprzedziałów są mniejsze od:

N

n

a

b





,

Układy równań liniowych:

Metody dokładne rozwiązywania układu równań liniowych:

,

b

Ax



gdzie macierz kwadratowa A stopnia n i wektor b  Rn są dane,

teoretycznie pozwalają uzyskać rozwiązanie dokładne:

,

1

b

A

x



gdzie x  Rn_{, o ile macierz kwadratowa A jest nieosobliwa. Ponieważ}

wartości elementów macierzy A i wektora b nie są reprezentowane dokładnie (rozdzielczość typów liczbowych w operacjach zmiennoprzecinkowych w danym środowisku programowym), i podczas obliczeń występują błędy zaokrągleń, więc rozwiązanie obliczone metodą „dokładną” może różnić się od teoretycznego rozwiązania dokładnego.

Układy równań liniowych:

W przypadku układu równań liniowych z macierzą trójkątną górną:

n n nn n n n n

b

x

a

b

x

a

x

a

b

x

a

x

a

x

a

















...

,

...

,

...

jeżeli macierz układu nie jest osobliwa, tj. gdy a_ii ≠ 0 dla każdego

i = 1, 2,..., n to niewiadome x_i są obliczane na podstawie zależności:

1

,...,

1

,

,











_i

_n

_n

a

x

a

b

x

Układy równań liniowych:

W przypadku układu równań liniowych z macierzą trójkątną dolną:

n n nn n n

x

a

x

a

x

b

a

b

x

a

x

a

b

x

a















...

...

,

,

jeżeli macierz układu nie jest osobliwa, tj. gdy a_ii ≠ 0 dla każdego

i = 1, 2,..., n to niewiadome x_i są obliczane na podstawie zależności:

n

i

a

x

a

b

x

,

1

,

2

,...,









Następnie w macierzy A(k) _{zamieniamy wiersze o numerach i}

k oraz k,

a w wektorze b(k) składowe o takich samych numerach. Odpowiada to

lewostronnemu przemnożeniu macierzy A(k) _{i wektora b}(k) _{przez macierz}

permutacji P(k)_{. W otrzymanym} układzie równań:

Układy równań liniowych:

Rozwiązywanie układu równań liniowych metodą eliminacji Gaussa

z częściowym wyborem elementu podstawowego:

Macierz A jest rozkładana na iloczyn macierzy LU, gdzie L oznacza macierz trójkątną dolną z jedynkami na głównej przekątnej, a U macierz trójkątną górną. Krok k (k = 1, 2 , ... , n-1) operacji rozkładu macierzy A można opisać następująco:

Znajdujemy element  k taki, że: k k i

a

a

a



max

b

x

A



Układy równań liniowych:

Rozwiązywanie układu równań liniowych metodą eliminacji Gaussa

z częściowym wyborem elementu podstawowego:

eliminujemy niewiadomą x_k z wierszy o numerach k + 1, k + 2, .... , n korzystając z wzorów:            

,

,

...

,

2

,

1

,

,

...

,

2

,

1

,

,

n

k

k

j

n

k

k

i

b

l

b

b

a

l

a

a





















a

a

l



Układy równań liniowych:

Rozwiązywanie układu równań liniowych metodą eliminacji Gaussa

z częściowym wyborem elementu podstawowego:

po n-1 krokach otrzymamy układ równań postaci:



































1

...

...

...

...

...

...

0

...

1

0

...

0

1

0

...

0

0

1

l

l

l

l

l

l

L

,

Pb

L

Ux



w którym P=P(n-1)_P(n-2)_...P1_{, U=A}(n) _oraz:

Powyższy układ równań jest rozwiązywany za pomocą zależności:

Pb

L

b

n

n

i

x

u

b

x

ˆ

;

1

,...,

1

,

,

ˆ

_

_

_







Układy równań liniowych:

Wyznaczenie macierzy odwrotnej do macierzy rzeczywistej:

W celu wyznaczenia macierzy odwrotnej do macierzy rzeczywistej A:

,

...

...

...

...

...

...

...

...

...



































a

a

a

a

a

a

a

a

a

a

a

a

a

a

a

a

A

macierz A jest najpierw rozkładana na iloczyn LU metodą eliminacji Gaussa. Po wykonaniu tej operacji otrzymujemy:

gdzie P oznacza macierz permutacji.

P

L

U

A

LU

PA

_



Układy równań liniowych:

Wyznaczenie macierzy odwrotnej do macierzy rzeczywistej:

Elementy macierzy L-1 _{oraz elementy macierzy U}-1 są określone wzorami:

,

,

...

,

1

,

,

,

...

,

2

,

1

,

,

,

...

,

1

,

,

,

...

,

2

,

1

,

n

j

j

i

n

j

u

u

u

u

n

j

j

i

n

j

l

l

l

l









































.

gdy

,

1

,

gdy

,

0

j

i

j

i



Układy równań liniowych:

Metody iteracyjne rozwiązywania układu równań liniowych:

),

1

(

b

Ax



gdzie macierz kwadratowa A stopnia n i wektor b  Rn są dane, generują

ciąg kolejnych przybliżeń wektora x.

Rozpoczęcie procesu iteracyjnego wymaga podania początkowych przybliżeń wszystkich składowych tego wektora.

Jedną z metod iteracyjnych rozwiązywania układu równań liniowych z nieosobliwą macierzą kwadratową A jest metoda Jacobiego. W metodzie tej macierz A jest przekształcana na sumę trzech macierzy:

,

U

D

L

A







gdzie L oznacza macierz trójkątną dolną, D – macierz diagonalną, a U macierz trójkątną górną.

Układy równań liniowych:

Rozwiązywanie układu równań liniowych metodą iteracyjną (przybliżoną) Jacobiego:

Uwzględniając rozkład macierzy A, układ równań (1) można zapisać w postaci:









(

2

)

,

b

D

x

U

L

D

x

b

x

U

L

Dx

















z czego wynika następujący proces iteracyjny:









,

,

b

x

U

L

Dx

b

x

U

D

L















jeżeli promień spektralny macierzy –D-1_{(L+U) jest mniejszy od 1, to}

Układy równań liniowych:

Rozwiązywanie układu równań liniowych metodą iteracyjną (przybliżoną) Jacobiego:

Z zależności (2) wynika, że (k+1)-sze przybliżenie i-tej składowej rozwiązania jest określone wzorem:





,

0

lub

0

,

,

max









x

x

x

x

x

x



),

3

(

,

...

,

2

,

1

,

n

i

a

b

x

a

x











przy czym a_ii≠ 0. Proces iteracyjny kończy się, gdy:

gdzie:

,

max

x

x



Równania i układy równań nieliniowych:

Równanie skalarne z jedną niewiadomą może być zapisane w postaci:

 

x



0

,

f

gdzie f: X → Y, przy czym X, Y  R. Jeśli równanie to nie może być

rozwiązane analitycznie, to do znalezienia jego pierwiastków stosujemy zwykle metodę iteracyjną, za pomocą której można otrzymać przybliżoną wartość pierwiastków, w przypadku, gdy rozwiązanie równania istnieje. Stosując metodę regula falsi znajduje się wartość pierwiastka równania skalarnego z jedną niewiadomą leżącego w przedziale [a, b] takim, że:

   

a

f

b



0

Jeden krok metody regula falsi polega na zastąpieniu bieżącego przedziału [a, b] przedziałem mniejszym, którym jest jeden z przedziałów [a, z] lub [z, b], gdzie:

Równania i układy równań nieliniowych:

     

f

b

f

a

a

b

b

f

b

z









Wybór odpowiedniego przedziału zależy od spełnienia jednego z warunków f(a)f(z) ≤ 0 lub f(z)f(b) ≤ 0.

Wartość pierwiastka równania skalarnego z jedną niewiadomą leżącego

w przedziale [a, b] można również znaleźć za pomocą metody siecznych

(interpolacji liniowej), która polega na zastosowaniu zależności:

Równania i układy równań nieliniowych:



b

a



b



b

a



a



0

,

1



i



0

,

1



     

,

2

,

3

,

...

.











k

x

f

x

f

x

x

x

f

x

x

Za początkowe przybliżenia x₁ i x₂ w procesie iteracyjnym przyjąć można np. wartości (zależnie od rozdzielczości podziału przedziału [a, b]):

uporządkowane w taki sposób, że:

 

x

f

 

x

Różniczkowanie numeryczne:

Przybliżone wartości pochodnej funkcji f(x) mogą być wyznaczone z

wielomianu aproksymującego p(x). Ponieważ różnica f(x) - p(x) może być bardzo mała, a różnica f’(x) – p’(x) bardzo duża, podejście takie może spowodować powstanie zbyt dużego błędu. Gdy jest wymagana duża dokładność obliczania pochodnej zaleca się stosować metodę Romberga. W metodzie tej dla danego początkowego kroku h₀ wyznacza się wartości:

...

,

1

,

0

,

2





h

i

h



 



,

2

h

h

x

f

h

x

f

T









następnie korzystając z wzoru:

i

k

T

T

T

T

,

1

,

2

,

...

,

1

4









Różniczkowanie numeryczne:

 

x

i

k

f

T

T

T

T

T

T

T

i,

'

,

gdy

dla

danego

22 21 20 11 10 00



















Obliczenia kończą się, gdy:







T

T

gdzie ε oznacza zadaną z góry dokładność.

Jeżeli dla założonego i oraz k powyższy warunek nie będzie spełniony to dany krok h₀ dzieli się na połowę i cały proces obliczeń powtarza się. Gdy po mN iteracjach zadana dokładność ε nie jest osiągnięta to za rozwiązanie przyjmuje się ostatnio obliczone przybliżenie T_i,k.

Różniczkowanie numeryczne:

Rozwiązanie zagadnienia początkowego dla pojedynczego równania

różniczkowego to inaczej obliczenie przybliżonej wartości y(x₁)

rozwiązania skalarnego równania:

 

x

y

f

dx

dy

y

'





,

 

x

y

y



Jedną z metod numerycznych prowadzących do rozwiązania tego zagadnienia jest metoda Rungego-Kutty z automatycznym doborem kroku całkowania.

Różniczkowanie numeryczne:

Wartość y(x₁) oblicza się stosując zależność:



   

2

2



6

K

K

K

K

h

x

y

h

x

y













 





,

,

2

,

2

,

2

,

2

,

,

hK

y

h

x

f

K

K

h

y

h

x

f

K

K

h

y

h

x

f

K

y

x

f

K



















_

_















_

_





Różniczkowanie numeryczne:

Na podstawie przedstawionej zależności wyznacza się dwie wartości przybliżonego rozwiązania w punkcie x+h, gdzie h oznacza krok całkowania, a mianowicie: y(x+h,h) przy użyciu kroku h i y(x+h,h) stosując dwukrotnie krok h/2. Jeżeli:









,

(

4

)

2

,

,

,

max

2

,

,





























 













 





h

h

x

y

h

h

x

y

h

h

x

y

h

h

x

y

gdzie ε, η > 0 oznaczają wartości charakteryzujące dokładność rozwiązania, to końcową przybliżoną wartość rozwiązania w punkcie x+h oblicza się z wzoru:





,

2



,



,

h

h

x

y

h

h

x

y

h

h

x

y

h

x

y















 









 





Różniczkowanie numeryczne:

W przypadku, gdy nierówność (4) nie zostanie spełniona, to zmienia się krok całkowania h.

Ponieważ całkowanie od punktu x₀ do punktu x₁ wykonuje się w pewnej liczbie kroków, w praktyce błąd względny otrzymanego przybliżonego rozwiązania y(x₁) może być większy od ε.