Metody numeryczne dla inżynierów [PDF]

Zitiervorschau

THE COLLEGE OF COMPUTER SCIENCE

E d w a rd K ą c k i, A n d rz e j M a ło le p s z y A lic ja R o m a n o w ic z

M eto d y num eryczne dla inżynierów

W yższa S z k o ła In fo r m a ty k i i v Łodzi

Spis treści 1. Informacje wstępne 1.1. Reprezentacje liczb

7 7

1.2. Rachunek zmiennopozycyjny

9

1.3. Algorytmy numeryczne

10

1.4. Numeryczna stabilność

algorytmów

2. Obliczanie wartości funkcji 2.1. Obliczanie wartości wielomianu 2.2. Obliczanie wartości wielomianu trygonometrycznego

14 16 16 19

3. Aproksymacja 3.1. Sformułowanie zagadnienia

23 23

3.2. Wielomiany ortogonalne

26

3.3. Wielomiany ortogonalne na zbiorze dyskretnym

29

3.4. Aproksymacja w przestrzeniach Hilberta

30

3.5. Aproksymacja na zbiorze dyskretnym

32

3.6. Przykłady

34

4. Interpolacja 4.1. Interpolacja wielomianowa 4.2. Wielomian interpolacyjny Lagrange’a

37 37 38

4.3. W zór interpolacyjny Lagrange’a dla węzłów równoodległych

40

4.4. Metoda Aitkena

41

4.5. Wzór interpolacyjny Newtona

44

4.6. Wzór interpolacyjny Hermite’a

46

4.7. Przykłady

48

5. Rozwiązywanie układów algebraicznych równań liniowych 5.1. Metody dokładne

49 49

5.1.1. Metoda eliminacji Gaussa

49

5.1.2. Oszacowanie błędu

56

5.1.3. Przykłady

59

5.2. Metody iteracyjne 5.2.1. Metoda iteracji prostej

63 64 3

5.2.2. Metoda Seidla 6

. Rozwiązywanie równań nieliniowych oraz ich układów 6.1. Metody wyznaczania pierwiastków w zadanym

69 70

przedziale

70

6.1.1. Metoda połowienia przedziału

70

6.1.2. Metoda cięciw zwana również regułą falsi

73

6.1.3. Metoda stycznych zwana również metodą Newtona 6.1.4. Metoda mieszana 6.2. Rozwiązywanie układów równań nieliniowych

76 77 80

6.2.1. Metoda Newtona

80

6.2.2. Metoda iteracji

84

7. Całkowanie numeryczne 7.1. Kwadratury liniowe

89 89

7.2. Kwadratury Newtona-Cotesa

91

7.3. Złożone kwadratury Newtona-Cotesa

94

7.4. Kwadratury Gaussa

98

7.5. Obliczanie całek niewłaściwych i całek z osobliwościami

8

104

7.6. Obliczanie numeryczne całek wielokrotnych

110

. Rozwiązywanie równań różniczkowych zwyczajnych 8.1. Metody różnicowe

112 112

8.1.1. Metody jednokrokowe

112

8.1.2. Liniowe metody wielokrokowe

114

8.1.3. Metoda współczynników nieoznaczonych

117

8

.1.4. Metody iteracyjne

8.1.5. Obliczanie wartości wstępnych 8.2. Metody Rungego-Kutty

119 122 123

8.3. Równania n-tego rzędu i układy równań różniczkowych zwyczajnych

128

8.3.1. Metoda Rungego-Kutty

128

8.3.2. Metody ekstrapolacyjne

133

8.3.3. Metody interpolacyjne

135

8.3.4. Metoda różnicowa zwykła i ulepszona

138

8.4. Metody wariacyjne 4

142

9. Rozwiązywanie równań różniczkowych cząstkowych 9.1. Metody różnicowe

151 151

9.1.1. Ogólne uwagi o dyskretyzacji i funkcjach siatkowych

152

9.1.2. Równanie Laplace’a

157

9.1.3. Równanie Poissona

163

9.1.4. Równanie Helmholtza

164

9.1.5. Równanie przewodnictwa

165

9.1.6. Równanie fali płaskiej

171

9.1.7. Równanie telegrafistów

173

9.1.8. Ulepszona metoda różnicowa

174

9.2. Stabilność schematów różnicowych

175

9.2.1. Uwagi wstępne

175

9.2.2. Stabilność schematów różnicowych dla równania typu parabolicznego

178

9.2.3. Stabilność schematu różnicowego dla równania typu hiperbolicznego

186

9.3. Metody Galerkina

190

9.4. Przykłady

197

9.5. Metoda elementu skończonego

205

10. Rozwiązywanie równań całkowych 10.1. Metoda kolejnych przybliżeń 10.2. Metoda sum skończonych 11. Metody Monte Carlo 11.1. Rozwiązywanie równań z jedną niewiadomą 11.2. Obliczanie całek oznaczonych

212 212 215 217 217 220

5

1. Inform acje w stępne 1.1 R eprezentacja liczb W komputerach liczby są reprezentowane przez skończoną liczbę cyfr ich rozwinięć pozycyjnych o pewnej podstawie rozwinięcia. Najczęściej stosowanymi podstawami rozwinięć liczb są podstawy 2, 8 lub 16. Arytmetyka tych rozwinięć jest podobna, dlatego też możemy ograniczyć się do arytmetyki dwójkowej inaczej zwanej binarną. W systemie dwójkowym mamy dwie cyfry 0 lub 1. W komputerach na reprezentację liczby przeznacza się słowo o skończonej długości np. d + 1 bitów (cyfr dwójkowych). Sposób rozwinięcia liczby zależy od jej typu. Liczby całkowite przedstawiane są w sposób stałopozycyjny. Dowolną liczbę całkowitą I w systemie dwójkowym możemy przedstawić w postaci

1=

e' 2‘ ’

( 1-1)

i= 0

gdzie zn - jest znakiem liczby /, e, - są jej cyframi dwójkowymi, przy czym jeśli I * 0, to en ^ Q . W komputerze na reprezentację liczby stało pozycyjnej jest przeznaczone słowo o długości d + 1 bitów i jeden bit jest przeznaczony na znak liczby, zatem liczba całkowita I może być przedstawiona w postaci (1.1) jeśli / e [ - 2 d + 1, 2d - 1]. Oznacza to, że n < d i mówimy wtedy, że liczba jest reprezentowana dokładnie. P r z y k ła d 1.1. Liczba I = 18 w systemie dwójkowym może być przedstawiona w postaci / = 1 ■24 + 0 • 2 3 + 0- 2 2 + 1 • 2 1 + 0 • 2°. Może więc być przechowywana w słowie o długości d + 1 > 5 bitów jako

000...010010 -------- V---------'

2«

ci cyfr rozwinięcia

Liczby rzeczywiste są przedstawione w sposób zmiennopozycyjny. Dowolną liczbę rzeczywistą x ^ 0 (zero jest na ogół zapamiętane jako słowo o wszystkich bitach zerowych) możemy przedstawić w postaci x = zn ■2Cm ,

( 1 .2 )

gdzie zn - jest znakiem liczby, c - jest liczbą całkowitą zwaną cechą, m - jest liczbą rzeczywistą z przedziału [ y , l)

zwaną mantysą. Liczba 7

rzeczywista jest więc przedstawiona za pomocą dwu grup bitów. Jeśli liczba x jest przechowywana w jednym słowie o długości d + 1 bitów, wtedy jeden bit przeznaczony jest na znak liczby, t bitów na reprezentację mantysy m, a pozostałe d - t bitów na reprezentację cechy c, przy czym jeden z bitów tej reprezentacji musi być przeznaczony na jej znak. Mantysa jako liczba rzeczywista z przedziału przedstawiona w rozwinięcia

systemie dwójkowym w

[y ,l) ,

może być

postaci nieskończonego

oo

(1.3)

m = Y ue- i 2 ' > gdzie e_, =

1

oraz e_; równe

0

lub

1

dla /> 1 .

Jednak na jej zapamiętanie jest przeznaczonych t bitów, może być zatem zapamiętane tylko jej t początkowych cyfr dwójkowych, tzn. mantysa zostaje prawidłowo zaokrąglona do t cyfr w następujący sposób: (1.4) 7

=1

W wyniku zaokrąglenia, a więc zastąpienia m (1.3) przez mt (1.4) mamy \m -

< -i- • 2 1.

(1.5)

Natomiast reprezentacja zmiennopozycyjna liczby rzeczywistej x danej wzorem ( 1 .2 ) oznaczana jest symbolem rd{x) i jest równa rd(x) = zn • 2 Cmt .

(1.6)

Ze wzorów (1.2), (1.4) i (1.6 ) otrzymujemy, że dla x # 0 rd(x) - x < x co oznacza, że rd(x) = x (l + £■),

gdzie | e |
0

Funkcja y 4 (x ) przyjmuje zatem postać

y 4 (x ) = 2 + ^ x - ^ - ( 2 x 2 - l ) - - j ( 4 x

3

- 3x)- y (

8

x4 -

8 * 2

+ l) ,

skąd po uporządkowaniu mamy /

\

32

4

16

3

8

7

(x)= y x - y x - y x dla x e [ - l ,

56

46

\ r

/

+ —x + — =7r(x+l)jl-x

y

1]

Błąd aproksymacji w danym przypadku jest określony wzorem (3.46). Wynosi on || / - y4 || , przy czym i / - > d i = v i i / i r + i b 4 r - 2 ( / > > ’4) =

i

£ [ f ’ T t\ fet, ( T , . T „ )

= J — M 8

15

=0,25.

35

P r zy k ła d

3.2.

Funkcję

mamy podaną w postaci tablicy

/(x )

wartości

i

Xi

fi

1

0,5

0,62

2

1,0

1,64

3

1,5

2,58

4

2,0

3,70

5

2,5

5,02

6

3,0

6,04

Aproksymujemy tę funkcję funkcją liniową. Przyjmujemy więc we wzorze (2.37) m = 2 oraz g, (x ) = 1, g 2 (x ) = x .Podstawiając te wartości do wzoru (3.37) otrzymamy

j( x ) = Ą + Ą x

Natomiast układ normalny (3.39) będzie miał postać: n A , 77

+

Y , Xi =

Z

/=1

/=1

A, ¿ x ,. + ^ ¿ x

1=1

n

i=1

,2

/ ' ’

= ¿ / x,..

!=1

A po uwzględnieniu wartości z podanej tabeli

6,00

/t, + 10,50/1, = 19,60

10,50/1, + 22,75/1, = 43,89

36

Rozwiązując ten układ równań otrzymujemy

/t,=

14 945 26,25

57 54 A t = ^ ^ = 2,192, 2 26,25

= -0 ,5 6 9 ,

oraz y (x) - - 0 , 5 6 9 + 2,192 x

Błąd aproksymacji w danym przypadku wynosi

i n

in

f - Z A f t ( X ,)

=

X (/;. + 0,569 - 2,192X )2 = 0,258

4. loterpolacja 4.1 interpolacja wielomianowa Dane są wartości y 0, y l } . . . , y „

funkcji f ( x )

w pewnych punktach

x 0 , x , , . . . ,x n , zwanych węzłami. Mogą być również pochodnych

funkcji

/ (x)

w

węzłach.

dane wartości

Zasadniczym

interpolacji jest określenie przybliżonych wartości funkcji

zadaniem / (x)

w

punktach nie będących węzłami, oraz oszacowaniem błędu wartości przybliżonych. Rozpatrzmy przypadek interpolacji za pomocą wielomianu. Wyznaczymy zatem wielomian w „ ( x ) stopnia n-tego, który w węzłach x0 , x , , . . . , x „

przyjmuje dane wartości y 0, y i , . . . , y „

funkcji / ( x ) .

Dla wyznaczenia wielomianu n

% ( x ) = X ak %k

(4 -1)

k =0

należy

rozwiązać

następujący

układ

n

+ 1

równań

liniowych

z

niewiadomymi współczynnikami ciQ,al ,... ,an :

37

X0 an + X0 1an~1 ■+ x"

Przy

założeniu,

+ x"

1 an_,

XX

+ X 'n~l a n -

że

x; ^ xy

h a0 - 3;0 >

+ •• • + n 0 = y l ,

(4.2)

1 + - ‘ - + «0 = 3 ;»-

dla

każdego

i& j ,

wyznacznik

charakterystyczny układu (4.2) jest różny od zera: r"

i

r ”“ 1 • ■ ■ x xa ,,_1 i

•••

1

xi 1 1

a zatem istnieje jedno rozwiązanie wielomian interpolacyjny (4.1).

,n

układu

równań

(4.2) i jeden

Na ogół dla punktów nie będących węzłami różnica f (x ) - w „ ( x ) = e (x ) jest różna od zera. Oszacowanie reszty

e (x )

(4.3) dla

x ^ x;

stanowi

oszacowanie błędu interpolacji. W ogólnym przypadku wielomian (4.1) może być przedstawiony jako kombinacja liniowa wielomianów g k { x ) co najwyżej n-tego stopnia określonego typu:

wA

x ) = Y j ak S k ix )

(4 -4 )

k= 0

4.2 W ielomian interpolacyjny Lagrange’a Wzór interpolacyjny Lagrange’a umożliwia prostą konstrukcję wielomianu (4.1) stopnia n, który w n + 1 węzłach x 0 ,x 1 , . . . , x n przyjmuje zadane wartości y 0 = / ( x 0) , y 1 = / ( x , ) , ... ,y n = f ( x a) . Wzór interpolacyjny Lagrange’a opiera się na wzorze (4.4) gdzie wielomiany g k ( x ) są stopnia n-tego i oznaczone przez !k ( x ) . W ielomiany lk ( x ) konstruuje się w ten sposób, aby wzór (4.4) postaci: 38

w,X x ) = Y J y k h ( x ) ■

(4-5)

k=O

Będzie tak wtedy, gdy h ( x ) = 5 kj»

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

gdzie f0

dla k ^ i ,

* “ {l

dla k = j ,

’

(4'6)

jest symbolem Kroneckera. Ponieważ chcemy, aby wielomian lk ( x ) był stopnia n-tego i spełniał warunek (3.6) musi być wyrażony wzorem:

fi (*-*/ ) h (*)= 4 ^ n

■

(4.7)

( xk - X i )

1=0

j*k

Wielomian (4.7) jest jedynym możliwym wielomianem stopnia n-tego spełniającym warunek Kroneckera i nie spełnia tych warunków żaden wielomian stopnia wyższego. Podstawiając (4.7) do (4.5) otrzymujemy, że wzór interpolacyjny Lagrange’a fu n k c ji) ma postać

0 ( * “ */) 1= 0

j * k ______________

-------- .

k=0

fl

(4.8)

( X ~ Xi )

1=0

J*k

nie wymaga zatem rozwiązania układu równań (4.2). Wzór interpolacyjny (4.8) ma sens przy założeniu, że xi ^ x k dla każdego i ^ k , a dla postaci:

x ^ xi (i =

0

, 1 ,.. . ,n )

można go również zapisać w

39

w.

ykp A x )

00 1

(4.9)

=

*=° ( x - x k ) p n ( x k )

gdzie

A 7 \ dPn z7» ( x a-) = \ dx s

/>» 0 0 = f i ( * “ * / ) ’ 7=0

Błąd (4.3) przybliżonej wartości funkcji / ( x )

(4.10)

obliczonej za pomocą

wzoru interpolacyjnego (4.8) można oszacować korzystając ze wzoru e (x ) =

P n (X ^

(4.11)

(n + 1)! gdzie Ę, leży w przedziale rozpiętym na węzłach x0,x l , ... ,xn i zależy od x. Wzór (4.11) wyraża także błąd w punktach będących węzłami, ponieważ P„ ( x,) = 0

dla i - 0 , 1 , . . . , « .

Jeżeli węzły uporządkujemy w ten sposób, że xi > xy dla każdego z > y , to

[ x 0,xn]

jest

przedziałem,

w

którym

leży

punkt

Ę, .

Oszacowanie błędu (4.11) określa relacja \Pn(X)\m hP„(X ) M — < e ( x ) \ < J------------ 1 ---J-----------1 (77 + 1 ) !

1

1

(4.12)

(77 + 1 ) !

gdzie 777

-

mm x e ( x 0 ,x„

przy założeniu, że funkcja / ( x )

M =

max x e { x 0 ,x„)

ma ciągłą pochodną /7-tego rzędu w

przedziale domkniętym [ x 0 ,x n ] . 4.3 W zór interpolacyjny Lagramge’a dla węzłów równoodległych Bardzo często węzły są rozłożone równomiernie w pewnym przedziale, tzn. odległości pomiędzy sąsiednimi węzłami są stałe i wynoszą 40

xi+, - xi = h > O

(4.13)

dla każdego i = 0 ,1 ,... ,n , wówczas wprowadzając nową zmienną i = iL z 5 L ,

(4.14)

h

możemy wzór Lagrange’a (4.8) sprowadzić do następującej postaci:

ń ( w ) ------------ .

(4-15)

n 7=0

j*k

gdzie t e [ 0 , 77] . Wielomian p n ( x ) będzie miał postać

# , ( 0 = 11 {t ~ j ) , 7=0

natomiast wzór (4.11) przyjmie postać

(77

-—

+ 1)!

- / ' " " U (4.15)

gdzie Ł, e[0,/7] .

4.4

¡VSetoda Aitkena

Wielomian interpolacyjny Lagrange’a n-tego stopnia można tworzyć metodą Aitkena Jest to metoda iteracyjna polegająca na interpolacji liniowej dwóch różnych wielomianów interpolacyjnych stopnia k < n zbudowanych (77

na

dowolnych

k +1

podzbiorach

węzłów

danego

+ l)-elem entow ego zbioru węzłów, w rezultacie czego otrzymujemy

wielomian interpolacyjny stopnia o jeden wyższego, tzn. stopnia k + 1 . Przez

w, o

( x)

oznaczmy

wielomian

interpolacyjny

n-tego

stopnia zbudowany na n + 1 różnych węzłach

xig ,xh , . . . ,jc,. , gdzie

dla

pierwszego

k A / x;. * x(. .

Wielomian

Lagrange’a

stopnia 41

zbudowany na węzłach xk ,xt

i przyjmujący w tych węzłach zadane

wartości y k ,y, zgodnie z zależnością (4.8) ma postać X - X, wk i( x ) = y k

gdzie

xk - X,

X - X,

+ yr

k , l e {0,1,... ,n }

i

x, - xk

xk - X,

x — x,

y,

x - xk

yk

,

(4.16)

k * l . Wielomian interpolacyjny stopnia

drugiego tworzymy przez interpolację liniową dwóch wielomianów stopnia pierwszego wkl ( x ) i wkm( x ) , gdzie k ^ m , k & l i / ^ m , w następujący sposób:

1 wkim( x ) =

x.„ - x.

x-x,

Wkl( x )

x - x „,

wh,Xx )

(4.18)

Wielomian stopnia r + 1 tworzymy przez interpolację liniową dwóch wielomianów stopnia r zbudowanych na zbiorach wielomiany w. .

(. i>5( x )

i wi o (

r + 1 węzłów. Są to

p { x ) , przy czym wymienione

zbiory węzłów różnią się tylko jednym elementem, a mianowicie xs ^ xp . Sposób tworzenia wielomianu ( r + 1 )-ego stopnia określa zależność X - X s. w, ( X P ~ Xs )

JtP

W: i W-/0 ,l,■

,

(X ) P

Wyznaczając wielomian interpolacyjny n-tego stopnia

(x)

.

wQ1

(4.19)

„(x)

zbudowany na węzłach x0 , x , ,.. . ,x n o wartościach w węzłach równych y 0,y n ... ,y n , obliczamy kolejno wielomiany wQk( x ) d i a k = 1,2, ... , n, w0lk ( x ) dla k = 2, 3, ... , n, wm2k ( x ) dla k = 1, 2, ... , n, itd. aż do uzyskania wielomianu stopnia n. Metoda Aitkena ze względu na swój iteracyjny charakter jest szczególnie przydatna przy komputerowej interpolacji wielomianami. Na rysunku 4.1 jest przedstawiony schemat blokowy algorytmu interpolacji wielomianami Lagrange’a n-tego stopnia skonstruowany na podstawie iteracyjnej metody Aitkena. Poniżej przedstawiamy program do tego schematu napisany w języku C. 42

S TA R T

Rys.4.1 43

/* Interpolacja wielomianami Lagrange'a */ /* n-tego stopnia, na podstawie */ /* iteracyjnej metody Aitkena */ void interpolacja(void) { /* x[], y[] wartości funkcji w węzłach */ for (i=l;in-i-l;k++) { l=k+i; y[k]=(x[k]*y[k+l]-x[l]*y[k] )/ (x [k] x [l] ) ;

4.5 Wzór interpolacyjny Newtona Dane

są wartości

y 0,y x, ... ,y„

funkcji

/(x)

wn+1

węzłach

x0,x,, ... ,x n . Wzór interpolacyjny Newtona n-tego stopnia zbudowany na n + 1 węzłach ma postać w „ ( x ) = ® ( x 0) + ( x - x 0) M = y 0, O ( x 0, x , ) =

yi ~ y 0

Xi - x0 ^ , O ( x 1,x2) - O ( x 0, x 1) O ( x 0,X;,x2) = --------------------------------- , x2 - x0

, ® ( X , , x2, ... , xn) - O ( X 0 , X j , ... , X H_ j ) O ( x 0,xl 5 ... , x j = .

44

(4.21)

Zastosowanie wzoru interpolacyjnego (4.20) wymaga wcześniejszego policzenia ilorazów różnicowych (4.21), które obliczamy kolejno według następującego schematu:

®(*„) ® ( x 0,x ,) 3>(x,)

® ( x 0, x , , x 2) ® ( x , , x 2)

® ( x 2)

0 ( x , ,x2 ,x3) ® ( x 2,x3)

0 ( x 3)

:

:

O ( x 0, x , , ... , x „)

® ( x b_2,x „_1,x „)

(4.22) Ilorazy różnicowe (4.21) można rozszerzać o następne węzły i ilorazy różnicowe, zwiększając w rezultacie stopień wielomianu interpolacyjnego aż do zagwarantowania odpowiedniej dokładności (4.3) określonej dla wzoru interpolacyjnego Newtona przez następującą zależność: e ( x ) = f ( x ) - w „ ( x ) = p „ ( x ) & ( x 0, X i , ... , x „ , x ) ,

(4.23)

gdzie p n ( x ) jest wielomianem określonym wzorem (4.10). Jeżeli węzły

x 0, x , , .. . ,x„

są rozłożone równomiernie, tzn. że

odległość pomiędzy sąsiednimi węzłami są stałe i wynoszą /? > 0 , to po wprowadzeniu do wzoru (4.20) nowej zmiennej t danej wzorem (4.14) otrzymamy w rezultacie następującą postać wzoru interpolacyjnego Newtona: v-i A* w, » ;,(0 = X k-0 k ■

-

7

( 4- 24)

gdzie: 45

k

a0( t ) = 1,

ak ( t ) = ] ^ [ ( t - m )

d la k > 0 ,

(4.25)

m=O

A/c>»0 jest

/c-tą różnicą progresywną funkcji

y = / (x )

określoną

wzorem (4.26)

Wielomian interpolacyjny wyrażony wzorem (4.24) jest różnicową postacią wzoru Newtona (4.20).

4.6 W ielomian interpolacyjny Hermite’a Rozpatrzmy przypadek interpolacji za pomocą wielomianu stopnia

m = n + r + 1 , który w

wartości

y 0, y j, - - - , y „

funkcji

węzłach

x 0 , x , ,.. . ,x n

f ( x ) , oraz

w

wm( x ) przyjmuje

pewnych węzłach

x0,x 1, ... ,xr gdzie r < n , przyjmuje wartości y'o ,y 'i

y ’r pochodnej

funkcji f ' ( x ) . Wzór interpolacyjny Hermite’a opiera się na wzorze n w ,,t ( x

) = Z y/t hk ( x )

+ Z y* hk ( x ) »

k =0

gdzie hk ( x ) i hk ( x ) są wielomianami stopnia Aby

wzór

(4.27)

wartości y 0 ,y } , . . . ,y„

(4 -27)

k=0

w

węzłach

funkcji

/(x),

przyjmował wartości y ’o,y'i,...,y'r być spełnione warunki

m= n + r + 1.

x0,xl s . . . , x B oraz węzłach

przyjmował x0,x l , . . . , x r

pochodnej funkcji f ( x ) , powinny

[ 1 - (x - x k) ( / ’t J x k) + l ’r (x k) \ lkn( x ) lkr(x ), dla k = 0,1,... ,r,

hk ( x ) = ( x — xk ) l kl. ( x ) lkn ( x ),

dla k = 0,1,... ,r. (4-29)

46

Błąd tej metody wyraża się wzorem

£ (-t) = 7 ^ T ^ / ' " ł " I>(i ) , (n + r + 2)!

(4.30)

gdzie £ leży w przedziale rozpiętym na węzłach x0 ,x ,, ... ,x n i zależy od x. Wielomiany oparte na n+1 węzłach są określone wzorami (4.7) i (4.10), a wielomiany p n(x) i /te(x) są określone tymi samymi wzorami opartymi na r+ 1 węzłach.

hk ( x j ) = 8kj,

dla k , j = 0,1,... ,/?,

K ( * / ) = 0,

dla k = 0, 1,... ,r; j = 0, 1,... , 77 , (4.28)

dla k = 0, 1,... , 77; j = 0, 1,... , 7%

K ( X j) = 0, f

dla k , j = 0, 1,... ,7'.

K ( xj ) = 5kj>

gdzie J /(/ jest symbolem Kroneckera (4.6). Będzie tak wtedy, gdy W zór (4.27) jest nazywany zmodyfikowanym wzorem interpolacyjnym Herm ite’a. Jeśli położymy r=n, otrzymamy

,U ) = Ż .y ,A ( * )

w „,

.

(4-31)

.

/r=0

7=0

gdzie K ( x ) = [ 1 - 2(x - x k ) /; (x /f) ] l 2 k (x)

(4.32)

hk ( x ) = ( x - x k ) l j ( x ) , dla /c = 0 ,1 ,... , 77. Błąd (4.30) przyjmie w tym przypadku postać

=

(4.33)

(277 + 2 ) !

Wielomian określony wzorem interpolacyjnym Herm ite’a.

(4.31)

nazywany

jest

wzorem

47

4.7 Przykłady

P r z y k ła d 4.1. Wyznaczymy wielomian interpolacyjny Lagrange’a

funkcji / ( x ) = ex w przedziale [0,2,0,5] korzystając z następujących danych / ( 0 , 2 ) = 1,2214 , / ( 0 , 4 ) = 1,4918 , / ( 0 , 5 ) = 1,6487. Korzystając ze wzoru (4.8) mamy w2 ( jc)

= 1,2214

-0 ,2 -(-0 ,3 )

0,2-(-0,1) = 0,724 x 2 + 0,9175* + 1,009.

Z wyprowadzonego wzoru interpolacyjnego obliczymy np. w2 (0,3), a więc przybliżoną wartość e0,3 , oraz oszacujemy błąd. Mamy zatem w (0,3) = 1,3493 , przy czym błąd (4.10) wynosi

Wielkość

eą

w

przedziale

[0,2,0,5]

spełnia

nierówność

1,2214 < & < 1,6487 . Stąd mamy oszacowanie błędu e(0,3) 0,0004 < e ( 0 , 3 ) < 0,0006 .

48

5 Rozwiązywanie układów algebraicznych równań liniow ych 5.1 Metody dokładne 5.1.1 Metoda eliminacji Gaussa Dany jest układ algebraiczny n równań liniowych z n niewiadomymi, o stałych współczynnikach: au x, + an x2+ ••• + aht x„ = 6,,

ci~,\ Xj + ci11 X-,+ •■•+ ciln Xn= b-,, (5.1) alń x, + a n2 x2+ — + am x„ = b n, który w postaci macierzowej wyraża się wzorem A X = B ,

(5.2)

gdzie

q22

, . Qn\

an2

••

B =

am_

.. o

#21

'x ,'

b\

■ n ... a2n

, 1

A =

Ć7j2

l

a\ i

x=

x2

_x«_

Zakładamy, że A jest macierzą nieosobliwą. Oznacza to, że wyznacznik charakterystyczny układu równań (5.1) jest różny od zera: det A

0

(5.3)

a zatem układ równań (5.1) ma dokładnie jedno rozwiązanie. Istnieje wiele modyfikacji metody eliminacji Gaussa. Jedną z nich obecnie się zajmiemy. Polega ona na odpowiednio usystematyzowanym rugowaniu kolejnych niewiadomych w ten sposób, ażeby z danego układu równań z n niewiadomymi uzyskać n równań, z których każde będzie tylko z jedną kolejną niewiadomą x, ,x2,.. . ,x„ . Z warunku (5.3) wynika, że co najmniej jeden ze współczynników przy niewiadomej x, w układzie równań (5.1) jest różny od zera, 49

możemy zatem przyjąć, że au * 0 . Rugujemy niewiadomą xx z kolejnych równań począwszy od drugiego przez dodawanie stronami równania pierwszego pomnożonego stronami przez czynnik (5.4)

Yk1

dla

a.

k = 2 ,3 ,... ,n . Operację tę

możemy zapisać

macierzowo w

następujący sposób: Wj A X = W XB .

(5.5)

gdzie 1

0

• 0

1 0

■• 0

0

1

• • 0

Yn1 0

0

•■ 1

0

Y 21 Wx = Y3i

(5.6)

W rezultacie operacji (5.5) otrzymujemy równanie Ax X = B x , które po

1

1 J -i

a \3

0

c i 22

a {X)

a 23 w

• a 2/7 m

0

¿70 ) 32

a 33 w

■ a 3n w

x3

o 1

2

l

l

1

1

rozpisaniu ma postać

a n2 {X)

a unm3

• ann w

. xn .

=

¿ 3( , )

(5.7)

_y?

Obecnie eliminujemy niewiadomą x2 ze wszystkich kolejnych równań układu (5.7) począwszy od równania trzeciego, czego dokonujemy dodając stronami równanie drugie pomnożone stronami przez czynnik

Yki =

50

a k2 w (i) a 22

(5.8)

do każdego /c-tego równania dla k = 3 , 4 , . . . ,n . Operację tę możemy zapisać w postaci zależności W2 A, X = W.; 5, .

(5.9)

gdzie

1

0

0

0

o

0

1

0 0

o

0

/u

1

0

r 42 0

W2 =

o

0

(5.10)

o

1

o r„2 o o ••• i W rezultacie operacji (5.9) otrzymujemy równanie A2 X = B2 , które po rozpisaniu ma postać a\\

£7)2

ai3

«14

'•

0

a2 w2

a 2w3

■••

0

0

0

0

(2) aW (2) «43

aw4 u2 (2) «34

_0

0

(2) a»3

a(2) 44 (2) «„4

«1n

am 2n (2) • •• «3,1 (2) ■’ • « 1

•• • an(2) n

-

X, '

6,

x2

ę

x3

¿14 = -Jr Ą 4 = -

i

■Tworzymy macierz V4 według wzoru (5.18):

wykonujemy

operację

1

0

0

0

1

0

5 42 19 84

0

0

1

29 294

0

0

0

1

V4 A3 X = V4 B3 ,

otrzymując

zależność

A (3> X = B {3) o następującej postaci:

61

'2

1 -3 1_

(i

X

1 2

2

2

0

y

57 4

0

0

J_ 7

0

z

3 14

0

0

0

u

-2 1

0

-4 2

Dla dalszych przekształceń wyznaczamy współczynniki Sk3 według zależności (5.20) dla

£ = 1 ,2 , otrzymując

Su = 21

i

t>23 = - -y .

Piszemy obecnie macierz V3 według wzoru (5.22):

i w ykonujem y

A (2> X = B {2)

0

0

1

-

f

0

0

1 0

0

0

0 1

V2 A (3) X = V3 B (3)

otrzym ując

1 0

0"

X

r

0

2 2

0

0

y

21 2

0

0

2 7

0

z

22 17

0

0

0

u

-2 1

-4 2

Obecnie obliczamy współczynnik

yn

ze wzoru (5.20), otrzymując

i tworzymy macierz V2 według wzoru (5.22):

"l j

0 o"

0

1 0

0

0

0

1_

a następnie wykonujemy operację

0

V2 A (2) X = V2 B {2)

rezultacie zależność A (l} X = B m o następującej postaci:

62

zależność

następującej postaci:

"2

8U = j

21 0"

V = 3 0

operację o

"1

otrzymując w

o

I o

0

CN 1 0

0

0

y

4' 21

X

2

0

0

y

0

z

3 14

0

0

0

- 42

u

-2 1

skąd mamy wartości niewiadomych x = 2 , y = - 3 , z = j

i u= \ .

5.2 Wletody Steracyjne Jeżeli mamy układ równań (5.2), to macierz współczynników A może być nieduża ale pełna, tzn. posiada mało elementów zerowych. Dla takich układów równań stosujemy metody dokładne rozwiązywania układów równań. Często jednak jest tak, że macierz A jest rzadka ale bardzo duża, tzn. posiada mało elementów niezerowych, które zwykle leżą w pobliżu głównej przekątnej macierzy A. Takie macierze występują powszechnie przy rozwiązywaniu równań różniczkowych cząstkowych. Mamy wtedy trudności z umieszczeniem elementów macierzy A w pamięci komputera. W tym przypadku bardziej użyteczne są metody iteracyjne rozwiązywania układów algebraicznych równań liniowych. Dla takich macierzy metody iteracyjne mogą być lepsze także pod względem ilości wykonywanych operacji arytmetycznych, a więc i czasu obliczeń od metod dokładnych. Dla rozwiązywania układów algebraicznych równań liniowych rozpatrujemy metody iteracyjne liniowe zmiennych x, , x2, ,x„ nie tylko dlatego że rozpatrywane zagadnienie jest liniowe, ale także dlatego, że analiza nieliniowych zagadnień iteracyjnych jest trudniejsza, oraz metody takie są mało efektywne, gdyż wym agają wielokrotnego mnożenia macierzy przez wektory w każdym kroku iteracyjnym. Metody iteracyjne polegają na przekształceniu układu równań (5.2) do postaci X = CX + D , i tworzeniu

ciągu

kolejnych

przybliżeń

(5.30) X ] , X 2, X 3, ...

wektora

X

spełniającego równanie (5.30) tworzonych według wzoru X n+l = C X n + D .

(5.31)

Jeżeli proces iteracyjny (5.31) jest zbieżny, wtedy mamy lim X n = X , w-»co

(5.32)

gdzie X spełnia równanie (5.30). 63

5.2.1 Metoda iteracji prostej Dany układ równań liniowych (5.1) z n niewiadomymi x 15x2,

,x„ o

stałych współczynnikach i wyznaczniku charakterystycznym różnym od zera \A\ ^ 0. Równania porządkujemy w ten sposób, ażeby na głównej przekątnej macierzy współczynników wszystkie elementy były różne od zera, jest to możliwe, wobec założenia | ^ | ^ 0 . Tak uporządkowany układ (5.1) przekształcamy, otrzymując

dx,

x, = 0 x l + cux2 + c13x3 +

■• ■ + cXnxn +

x9 = c2jXj + 0 x-> + c23x3 +

• • ■ + c~,^nxn + d-,,

x3 = c31Xj + c32x2+ 0 x3 + ••• + c3nxn + d3, Xn = CnlX1+ C,m X2 + C„3^3 +

■'• + Cn,n-lX»-l + 0 Xn + d n»

co w zapisie macierzowym ma postać (5.30), gdzie '0

c=

C|2

C13 ■’

Cl„

^"23

C 2n C2 n

C21

0

C31

c32

0

Cn2

C„3

_C«1

...

0

,

D =

d] d2 d3

Xi X, ,

x =

x3

A,_

W metodzie iteracji prostej ciąg kolejnych przybliżeń X l , X 2, X 3, ... wektora X spełniającego równanie (5.30) tworzonych według wzoru rekurencyjnego (5.31), przy założeniu zbieżności procesu iteracyjnego, tzn. spełnieniu warunku (5.32). Warunkiem koniecznym i wystarczającym zbieżności procesu iteracyjnego określonego wzorem (5.31) przy dowolnym wektorze początkowym X 0 i dowolnym danym wektorze D jest, aby wszystkie wartości własne macierzy C były co do modułu mniejsze od jedności. W praktyce jest na ogół wygodniej korzystać z następującego twierdzenia o warunku wystarczającym zbieżności procesu iteracyjnego (5.31). Jeżeli norma macierzy C je s t mniejsza od jedności, to proces iteracyjny określony wzorem (5.31) je s t zbieżny. Z przytoczonych twierdzeń wynika, że zbieżność procesu iteracyjnego nie zależy ani od wartości początkowej wektora X 0 , ani też od wartości danego wektora D, a wystarczy jedynie spełnienie jednego z następujących warunków: 64

n

I < 1’

(5.33)

max X |cy |< 1 ’

(5.34)

max ZI i

cj

7=1

lub ii J

/=i

lub

(5.35)

gdzie ^

są elementami macierzy C.

Oczywiście proces iteracyjny (5.31) jest zbieżny, gdy dla każdego i , j = 1 ,2 ,... ,n , 1 Ca < —

(5.36)

W tym przypadku każdy z warunków (5.33), spełniony. Oszacowanie błędu n-tego przybliżenia

(5.34) i (5.35) jest

uzyskanego na drodze

iteracji (5.31) dla wartości wektora początkowego X 0 przy założeniu ||C|| < 1 określa relacja

(5.37) gdzie X jest dokładnym rozwiązaniem równania (5.30). W praktycznych oszacowaniach błędu e znacznie wygodniej jest korzystać z relacji uzależniającej oszacowanie od różnicy dwóch kolejnych przybliżeń wektora X. Z twierdzenia \x - x Ą < ią \x ~ x l

(5.38)

otrzymujemy następujące oszacowanie błędu: (5.39)

65

które dla ¡C||
0 ograniczającej z góry błąd e = U - X H\. /*

R o z w ią z y w a n ie u k ła d u rów nań lin io w y c h */ / * z n n ie w ia d o m y m i m e to d ą i t e r a c j i p ro s te j * / v o id i t e r a c j a p r o s t a ( v o id ) {

/* X [i+1]=C*X[i]+D */ in t i; f l o a t e_m ax; f o r ( i = l ; i< n ; i+ + ) y [i]= z [ i ] ; do { f o r ( i = l ; i < n ; i+ + ) x [i ] = y[i ]; f o r ( i = l ; i< n ; i+ + ) { s= 0 ; for (j=l;j = e p s ilo n ) ; }

66

START

Rys.5.2 P r zy k ła d 5.2. niewiadom ym i

Dany

je st

układ

1

1

10

4

x= — y—

równań

liniowych

z

trzem a

7

z+ — , 4

1 2 V=-n ZH 19 , 2

5

5

z = —1 x + -1 y - 11. 4 5 Układ ten w zapisie m acierzowym ma postać

U = A U + B, gdzie 67

Sprawdzamy, czy jest spełniony warunek wystarczający zbieżności procesu iteracyjnego określonego wzorem

U„+i = A U n + B , przyjmując normę macierzy postaci

Ml = maxŻ h ] i=i

i normę wektora

IMI = S W /=!

W danym przypadku

= E k l = ° . 75/=i

Rozważany proces iteracyjny jest zbieżny, ponieważ |^4| < 1. Wartość początkową wektora U przyjmujemy

u =

tzn. że x0 = 4, y0 = 2 i z0 = 2. Wyznaczymy piąte przybliżenie wektora U i oszacujemy błąd postać

e = |{7 - 1/5||. W zór iteracyjny po rozpisaniu ma

Wyznaczamy pierwsze przybliżenie U 1 1 7 x, = — - 2 - - - 2 + - = 1,45, '1 0 4 4 1 2 19 y, = — • 4 + — • 2 + — = 6,6, 2 5 5 1 1 z, = — ■4 + —■2 - 1 = 0,4. W analogiczny sposób obliczamy kolejne przybliżenia wektora U, a mianowicie U2, U3, UĄ i U5: x2 = 2,31,

y 2 = 4,865,

z, = 0,6825,

x 3 = 2,04788,

_y3 = 5,228,

z3 = 0,5145,

x4 = 2,14418,

y 4 = 5,02974,

z4 = 0,55757,

x5 = 2,11358,

_y5 = 5,09512,

z5 = 0,54200.

Obecnie oszacujemy błąd e, korzystając ze wzoru (5.39). W tym celu wyznaczymy normę różnicy U5 - UĄ\ ||C/5- 6 / 4|| = |x 5-

x 4| +

|^5- ^ 4| + [z5-

z 4|

= 0,11155,

a następnie obliczamy

Ml \U5 - U 4\\ = 0,33465

i-

i zgodnie ze wzorem wyznaczonej wartości U5\

(5.39)

otrzymujemy

oszacowanie

błędu

e = \\U - U 5\\< 0,335. 5.2.2

Metoda Seidla

Dany jest układ równań liniowych (5.1) z n niewiadomymi o stałych współczynnikach i wyznaczniku charakterystycznym różnym od zera \A\ + 0. Układ ten przedstawiamy w postaci (5.30). W metodzie Seidla po obliczeniu dowolnej składowej x f wektora X korzystamy z niej przy obliczeniu składowej xf+1 jeszcze w tej samej k-tej iteracji według zależności: 69

* f , = x > r ci ‘ > + i ‘:i * r , , + 4 7=1 7=/+l

(5.41)

Warunkiem wystarczającym zbieżności metody iteracyjnej Seidla jest, aby ||C|| < 1, gdzie C = [cc] (patrz układ równań (5.30). W praktyce często przy ocenie zbieżności metody Seidla korzystamy z twierdzenia: Jeżeli macierz C je s t dodatnio określona, to metoda Seidla (5.41) je s t zbieżna dla dowolnego wektora początkowego X 0. Metoda Seidla jest szybciej zbieżna od metody iteracji prostej, przy założeniu zbieżności obu metod. Przy oszacowaniu błędu e k-tego przybliżenia uzyskanego metodą Seidla korzystamy z tych samych wzorów, co przy oszacowaniu błędu w metodzie iteracji prostej.

6

Rozwiązywanie y kład ów

równań

nieliniow ych

oraz

Ich

6.1 Metody wyznaczania pierwiastka równania w zadanym przedziale Dane jest równanie / ( * ) = 0,

(6.1)

gdzie f (x) jest funkcją ciągłą w przedziale [a ,b ] i f ( a ) f ( b ) < 0. 6.1.1 Metoda połowienia przedziału Wyznaczymy pierwiastek równania (6.1) leżący w przedziale [a ,b \ przez

dzielenie

przedziału

na

połowy

i

odrzucenie

części

nie

zawierającej pierwiastka. Jeżeli f ( \ ( a + b j) = 0, to ć, - \ j a + b) jest pierwiastkiem

równania.

Jeżeli

f { \ ( c i + b j ) ^ 0,

to

wybieramy

tę

połówkę przedziału [ a , \ ( a + b)~\ lub [\( c i + 6 ),6 ], na krańcach której funkcja f { x ) ma znak przeciwny i w ten sposób otrzymujemy nowy przedział

[a l ,b l \

zawierający pierwiastek równania (6.1). Obecnie

połowimy przedział [a x,b x~\ postępując tak samo jak z przedziałem [a ,b ] i w wyniku wielokrotnego powtórzenia tej procedury otrzymujemy ciąg przedziałów [a x,b x], [a 2,b2], ... , [a k ,bk ], ... o własności

70

/ ( Ą )f(bk) < O

dla k = 1, 2,...,

(6 .2 )

przy czym długości lk tych przedziałów tworzą ciąg geometryczny b- a 2k

Wobec ciągłości funkcji

/(x)

(6.3)

każdy z wymienionych przedziałów

\ak ,bk ] zawiera pierwiastek ¿(równania (6.1), a zatem mamy \\m a k = lim Z^,

k —>co

k —>co

Przyjmując przybliżone rozwiązanie xn , jako środek przedziału [a n,bn], mamy jednocześnie oszacowanie błędu e przybliżonego pierwiastka: (6.4) Na rysunku 6.1 jest przedstawiony schemat blokowy algorytmu wyznaczania przybliżonego pierwiastka równania (6.1) metodą połowienia przedziału. Poniżej przedstawiamy program w języku C dla algorytmu przedstawionego na rysunku 6.1.

/* Wyznaczanie przybliżonego pierwiastka */ /* równania nieliniowego metodą połowienia */ /* przedziału */ float połowienie(void) { /* a, b granice przedziału */ /* f(x)=0 rozwiązywane równanie */ float x; do { FA=f(a); x=(a+b)/2; FC=f(x); i f (FA*FC .

Rys.6.4 Metoda cięciw dla przypadku / ' ( * ) > 0 i f " { x ) < 0 dla x e[a,b] 75

Z założeń dotyczących funkcji

/ (x)

dla równoważnych dwóch

przypadków wynika, że istnieją liczby dodatnie m-\, M-\ ograniczające moduł pochodnej f ' { x ) w przedziale [a, 6 ], tzn. że 0 < n \ < | / ' ( x )| ^ M < +00 •

(6-9)

Błąd przybliżonego pierwiastka xk możemy oszacować również na podstawie następującej zależności:

i

l

M, - ni

i

i

m,

6.1.3

Metoda stycznych zwana również metodą Newtona

Przybliżoną wartość xk+i pierwiastka równania (6.1), leżącego w przedziale

[a, 6 ],

wyznaczymy,

dzieląc

przedział

[a, 6]

punktem

przecięcia osi Ox styczną do linii y = f { x ) w punkcie xk, przy czym o funkcji / (x) zakładamy, że ma drugą pochodną ciągłą w przedziale [a ,b ]

oraz że znak pierwszej i drugiej pochodnej są stałe w tym

przedziale. Ciąg kolejnych przybliżeń

pierwiastka równania (6.1) jest

określony przez wzór rekurencyjny /(* * )

Jeżeli dla każdego x e [ a , 6 ] mamy / ' ( x ) / " ( x ) < 0, to x 0 = a i ciąg {x^} jest ciągiem rosnącym ograniczonym z góry x0 < x, < ... < xk < x/c+1 < ... < ^ < b . Jeżeli dla każdego x e [a ,i)] mamy / ' ( x ) / " ( x ) > 0, to x0 = b i {x/;} jest ciągiem malejącym ograniczonym z dołu: a < %< ... < xk+l < x k < ... < x, < x0. Z przyjętych założeń wynika, że w przedziale (a,£>) istnieje dokładnie jeden pierwiastek ¿(równania (6.1), przy czym x/c -» ¿(.

76

Rys.6.5 Oszacowanie błędu przybliżonego pierwiastka xk określa wzór (6.10) lub zależność

e =

~ X k - \ ) 2>

\lv ~ X k \ -

(6 -1 2 )

gdzie 777, =

m in

|/ '( x ) |,

M2=

x s < a ,b >

max

|/" (x )|.

x e < a ,b >

Na rysunku 6.5 jest przedstawiona graficzna interpretacja metody stycznych przy założeniu, że dla każdego x e [a ,b ] f '( x ) > 0 i f "(x ) < 0. 6.1.4

Metoda mieszana

Przy założeniach, że / ( a ) f (b) < 0 i funkcja f(x) ma ciągłą pochodną w przedziale [a ,b ] wyznaczamy przybliżoną wartość xk pierwiastka równania (6.1) ze wzoru Xk + xk — .

gdzie

(6-13)

x k jest k-tym przybliżeniem pierwiastka obliczonym metodą

cięciw, a x k - metodą stycznych. Na rysunku 6.6 jest przedstawiona graficzna interpretacja metody mieszanej przy założeniach, że dla każdego x e [a ,6 ] f ' ( x ) > 0 i f " ( x ) > 0. 77

Jeżeli dla każdego x e [a ,b ] mamy f ( a ) f ( b ) < 0, to f ( X k - 1) x k = x k- i - — =— -— — - ( x ^ - a ) , J (x k- \ ) - f (a )

(6.14)

gdzie x 0 = b , oraz = = x k = x k- i

gdzie x 0 = a , przy czym ciągiem rosnącym:

/(* * _ ,) s------ , / ' ( * * - 1)

(6.15)

{x^} jest ciągiem malejącym, a {x k} jest

a < x i < ... < x k-\ < x k Jeżeli dla każdego x e [ a , 6 ]

x k = X t -1 -

x k+\ < x k < x k-\ < ... < X] < b . f ' { x ) f " { x ) > 0, to

/(x*-i)

-

— — — (x*_, —Z»), f { x k- \ ) - f { b )

(6.17)

gdzie x 0 = a , oraz x k określone jest wzorem (6.15), przy czym x 0 = b i { x i} jest ciągiem rosnącym, a {x *} jest ciągiem malejącym

a

< x\ < . . .