Wyklad mn no 3 piątek

background image
background image

Rozważymy przypadek z wielokrotną prawą stroną.
Należy rozwiązać układ równań:

0

0;

;

1

x

2

x

-

5

3;

-

;

3

x

x

2

x

0

5;

;

2

x

4

x

4

-

;

2

-

;

0

x

x

2

x

4

4

3

4

3

2

3

2

3

2

1

Tworzymy macierz dołączoną:

0

0

1

2

1

0

0

5

3

3

1

2

1

0

0

5

2

0

4

1

0

4

2

0

0

1

2

4

Dzielimy pierwszy wiersz przez a

11

=4

background image

0

0

1

2

1

0

0

5

3

3

1

2

1

0

0

5

2

0

4

1

0

1

5

.

0

0

0

25

.

0

5

.

0

1

Drugi wiersz dzielimy przez a

22

=1 i eliminujemy wyraz a

32

=1

0

0

1

2

1

0

0

5

8

5

1

2

0

0

0

5

2

0

4

1

0

1

5

.

0

0

0

25

.

0

5

.

0

1

Trzeci wiersz dzielimy przez a

33

=2 i eliminujemy wyraz a

43

=-1

background image

5

.

2

4

5

.

1

5

.

2

0

0

0

5

.

2

4

5

.

2

5

.

0

1

0

0

0

5

2

0

4

1

0

1

5

.

0

0

0

25

.

0

5

.

0

1

0

0

1

2

1

0

0

5

8

5

1

2

0

0

0

5

2

0

4

1

0

1

5

.

0

0

0

25

.

0

5

.

0

1

i ostatecznie dzielimy czwarty wiersz przez a

44

=2.5

1

6

.

1

6

.

0

1

0

0

0

5

.

2

4

5

.

2

5

.

0

1

0

0

0

5

2

0

4

1

0

1

5

.

0

0

0

25

.

0

5

.

0

1

background image

Eliminujemy wyrazy znajdujące się nad główną przekątną
w czwartej kolumnie:

1

6

.

1

6

.

0

1

0

0

0

5

.

2

4

5

.

2

5

.

0

1

0

0

0

5

2

0

4

1

0

1

5

.

0

0

0

25

.

0

5

.

0

1

mnożąc wyrazy czwartego wiersza przez –a

34

i dodając do wiersza

trzeciego mamy:

1

6

.

1

6

.

0

1

0

0

0

2

2

.

3

2

.

2

0

1

0

0

0

5

2

0

4

1

0

1

5

.

0

0

0

25

.

0

5

.

0

1

background image

Eliminujemy kolumnę trzecią:

1

6

.

1

6

.

0

1

0

0

0

2

2

.

3

2

.

2

0

1

0

0

0

5

2

0

4

1

0

1

5

.

0

0

0

25

.

0

5

.

0

1

i mamy:

1

6

.

1

6

.

0

1

0

0

0

2

2

.

3

2

.

2

0

1

0

0

8

8

.

7

8

.

6

0

0

1

0

5

.

1

3

.

0

55

.

0

0

0

5

.

0

1

background image

Ostatnia eliminacja drugiej kolumny:

i mamy:

1

6

.

1

6

.

0

1

0

0

0

2

2

.

3

2

.

2

0

1

0

0

8

8

.

7

8

.

6

0

0

1

0

5

.

2

6

.

3

85

.

2

0

0

0

1

1

6

.

1

6

.

0

1

0

0

0

2

2

.

3

2

.

2

0

1

0

0

8

8

.

7

8

.

6

0

0

1

0

5

.

1

3

.

0

55

.

0

0

0

5

.

0

1

background image

0

0;

;

1

x

2

x

-

5

3;

-

;

3

x

x

2

x

0

5;

;

2

x

4

x

4

-

;

2

-

;

0

x

x

2

x

4

4

3

4

3

2

3

2

3

2

1

W kolumnach 5, 6, 7 otrzymujemy rozwiązania odpowiednio dla
przypadku I, II i III.

1

6

.

1

6

.

0

1

0

0

0

2

2

.

3

2

.

2

0

1

0

0

8

8

.

7

8

.

6

0

0

1

0

5

.

2

6

.

3

85

.

2

0

0

0

1

6

.

0

2

.

2

8

.

6

85

.

2

X

I

6

.

1

2

.

3

8

.

7

6

.

3

X

II

1

2

8

5

.

2

X

III

background image

30

10

0

5

x

x

x

x

0

2

.

0

1

4

1

10

0

2

.

0

2

0

5

1

4

1

2

0

4

3

2

1

Rozważmy układ równań

który jest układem już rozwiązanym po zmianie pierwszego
i ostatniego wiersza w układzie:

5

10

0

30

x

x

x

x

4

1

2

0

1

10

0

2

.

0

2

0

5

1

0

2

.

0

1

4

4

3

2

1

Nie może układ wierszy wpłynąć na istnienie rozwiązania!

background image

Metoda Gaussa z częściowym wyborem elementu

podstawowego

30

10

0

5

x

x

x

x

0

2

.

0

1

4

1

10

0

2

.

0

2

0

5

1

4

1

2

0

4

3

2

1

Największy element w kolumnie pierwszej znajduje się we wierszu

czwartym i dlatego zamieniamy te dwa wiersze miejscami:

5

10

0

30

x

x

x

x

4

1

2

0

1

10

0

2

.

0

2

0

5

1

0

2

.

0

1

4

4

3

2

1

background image

Przeprowadzamy eliminację:

157894737

.

8

578947368

.

3

021052632

.

1

0

0

421052632

.

8

021052632

.

1

989473684

.

9

0

0

578947368

.

1

421052632

.

0

010526316

.

0

1

0

5

.

7

0

05

.

0

25

.

0

1

0

2

.

0

1

4

/

1

5

4

1

2

0

5

.

8

1

99

.

9

05

.

0

0

5

.

7

2

05

.

0

75

.

4

0

5

.

7

0

05

.

0

25

.

0

1

background image

2971549

.

7

474582663

.

3

0

0

0

842992624

.

0

102212856

.

0

1

0

0

578947368

.

1

421052632

.

0

010526316

.

0

1

0

5

.

7

0

05

.

0

25

.

0

1

100152913

.

2

1

0

0

0

842992624

.

0

102212856

.

0

1

0

0

578947368

.

1

421052632

.

0

010526316

.

0

1

0

5

.

7

0

05

.

0

25

.

0

1

Jeszcze raz rzut oka na otrzymany układ równań:

1

.

2

x

843

.

0

x

102

.

0

x

5789

.

1

x

421

.

0

x

0105

.

0

x

5

.

7

x

0

x

05

.

0

x

25

.

0

x

4

4

3

4

3

2

4

3

2

1

background image

3897439

.

2

1

0

0

0

5987301

.

0

0

1

0

0

5788529

.

2

0

0

1

0

1147767

.

8

0

0

0

1

100152913

.

2

1

0

0

0

628329997

.

0

0

1

0

0

45660828

.

2

0

0

1

0

185835002

.

7

0

0

25

.

0

1

100152913

.

2

1

0

0

0

628329997

.

0

0

1

0

0

46322228

.

2

0

010526316

.

0

1

0

5

.

7

0

5

.

0

25

.

0

1

i redukcja górnej trójkątnej:

i ostatnia kolumna zawiera
rozwiązanie.

background image

Rozkład LU. Metoda Croute’a.
Rozkład na macierze trójkątne

Dana jest macierz A i przedstawiamy ją w postaci:

A=LU

gdzie macierz L jest macierzą dolną trójkątną:

55

54

53

52

51

44

43

42

41

33

32

31

22

21

11

l

l

l

l

l

0

l

l

l

l

0

0

l

l

l

0

0

0

l

l

0

0

0

0

l

L

background image

lub ogólnie:

 

i

j

dla

obliczane

i

j

dla

0

u

u

ij

ij

U

Macierz U górna trójkątna:

55

45

44

35

34

33

25

24

23

22

15

14

13

12

11

u

0

0

0

0

u

u

0

0

0

u

u

u

0

0

u

u

u

u

0

u

u

u

u

u

U

lub ogólnie:

 



j

i

dla

obliczane

j

i

dla

1

j

i

dla

0

l

l

ij

ij

L

background image

Jeżeli

A=LU

, to dla układu równań

AX=b

mamy:

b

LY

Y

UX

b

LUX

Rozwiązanie układu LY=b z dolną macierzą trójkątną jest łatwe:

 

1

i

1

j

j

ij

i

ii

i

11

1

1

y

l

b

l

1

y

l

b

y

i=2,3,...,N

background image

i rozwiązanie równania UX=Y z górną macierzą trójkątną
jest łatwe:

 

N

1

i

j

j

ij

i

ii

i

NN

N

N

x

u

y

u

1

x

u

y

x

i=N-1,N-2,...,1

Duża zaleta:
Znając rozkład LU
możemy go wykorzystać wielokrotnie dla
różnych prawych stron.

background image

Obliczanie wyrazów macierzy L i U

NN

kN

kk

N

2

k

2

22

N

1

k

1

12

11

1

N

,

N

Nj

1

N

1

k

,

k

1

k

21

u

0

0

0

0

0

.

.

.

.

.

.

u

.

u

.

0

0

.

.

.

.

.

.

u

.

u

.

u

0

u

.

u

.

u

u

1

l

.

l

.

l

.

.

.

.

.

.

0

.

1

l

.

l

.

.

.

.

.

.

0

0

0

0

1

l

0

0

0

0

0

1

w wyniku mnożenia obu macierzy mamy macierz B=[b

ij

]

Zaczynamy kolejno:
pierwszy wiersz macierzy L
razy k-ta kolumna macierzy U:

k

1

k

1

u

b

k-ty wiersz macierzy L razy pierwsza kolumna macierzy U:

1

k

11

1

k

b

u

l

background image

NN

kN

kk

N

2

k

2

22

N

1

k

1

12

11

1

N

,

N

Nj

1

N

1

k

,

k

1

k

21

u

0

0

0

0

0

.

.

.

.

.

.

u

.

u

.

0

0

.

.

.

.

.

.

u

.

u

.

u

0

u

.

u

.

u

u

1

l

.

l

.

l

.

.

.

.

.

.

0

.

1

l

.

l

.

.

.

.

.

.

0

0

0

0

1

l

0

0

0

0

0

1

k-ty wiersz macierzy L razy j-ta (jk) kolumna macierzy U:

kj

kj

1

k

i

1

i

ij

ki

b

u

u

l

j-ty wiersz (j>k) macierzy L razy k-ta kolumna macierzy U:

jk

k

i

1

i

ik

ji

b

u

l

background image

ponieważ musi zachodzić B=A, czyli b

ij

=a

ij

dla (i,j=1,2,...,N)

stąd otrzymujemy kolejno:

Pierwszy wiersz macierzy U:

k

1

k

1

a

u

pierwsza kolumna macierzy L:

11

1

k

1

k

u

a

l

k-ty wiersz macierzy U:

1

k

i

1

i

ij

ki

kj

kj

u

l

a

u

dla j=k,k+1,...,N

k-ta kolumna macierz L:

1

k

i

1

i

ik

ji

jk

kk

jk

u

l

a

u

1

l

dla j=k+1,k+2,...,N

background image

Przykład:

10

1

0

0

0

4

8

4

0

0

1

8

2

1

0

0

3

1

5

2

0

0

2

1

4

A

Zgodnie z:

k

1

k

1

a

u pierwszy wiersz macierzy U:

.

0

0

0

0

.

.

0

0

0

.

.

.

0

0

.

.

.

.

0

0

0

2

1

4

U

background image

Pierwsza kolumna macierzy L zgodnie z

10

1

0

0

0

4

8

4

0

0

1

8

2

1

0

0

3

1

5

2

0

0

2

1

4

A

11

1

k

1

k

u

a

l

gdzie u

11

=4

1

.

.

.

0

0

1

.

.

0

0

0

1

.

0

0

0

0

1

5

.

0

0

0

0

0

1

L

background image

10

1

0

0

0

4

8

4

0

0

1

8

2

1

0

0

3

1

5

2

0

0

2

1

4

A

drugi wiersz macierzy U zgodnie ze wzorem:

1

i

1

i

ij

i

2

j

2

j

2

u

l

a

u

j=2,3,4,5

.

0

0

0

0

.

.

0

0

0

.

.

.

0

0

0

3

2

5

.

5

0

0

0

2

1

4

U

1

.

.

.

0

0

1

.

.

0

0

0

1

.

0

0

0

0

1

5

.

0

0

0

0

0

1

L

background image

10

1

0

0

0

4

8

4

0

0

1

8

2

1

0

0

3

1

5

2

0

0

2

1

4

A

j=3,4,5

.

0

0

0

0

.

.

0

0

0

.

.

.

0

0

0

3

2

5

.

5

0

0

0

2

1

4

U

1

.

.

0

0

0

1

.

0

0

0

0

1

18182

.

0

0

0

0

0

1

5

.

0

0

0

0

0

1

L

1

i

1

i

2

i

ji

2

j

22

2

j

u

l

a

u

1

l

Druga kolumna macierzy L:

background image

10

1

0

0

0

4

8

4

0

0

1

8

2

1

0

0

3

1

5

2

0

0

2

1

4

A

trzeci wiersz macierzy U zgodnie ze wzorem:

2

i

1

i

ij

i

3

j

3

j

3

u

l

a

u

j=3,4,5

.

0

0

0

0

.

.

0

0

0

1

54545

.

8

36364

.

2

0

0

0

3

2

5

.

5

0

0

0

2

1

4

U

1

.

.

0

0

0

1

.

0

0

0

0

1

18182

.

0

0

0

0

0

1

5

.

0

0

0

0

0

1

L

background image

10

1

0

0

0

4

8

4

0

0

1

8

2

1

0

0

3

1

5

2

0

0

2

1

4

A

.

0

0

0

0

.

.

0

0

0

1

54545

.

8

36364

.

2

0

0

0

3

2

5

.

5

0

0

0

2

1

4

U

1

.

0

0

0

0

1

692308

.

1

0

0

0

0

1

18182

.

0

0

0

0

0

1

5

.

0

0

0

0

0

1

L

j=4,5

2

i

1

i

3

i

ji

3

j

33

3

j

u

l

a

u

1

l

trzecia kolumna macierzy L:

background image

10

1

0

0

0

4

8

4

0

0

1

8

2

1

0

0

3

1

5

2

0

0

2

1

4

A

czwarty wiersz macierzy U zgodnie ze wzorem:

3

i

1

i

ij

i

4

j

4

j

4

u

l

a

u

j=4,5

background image

1

.

0

0

0

0

1

692308

.

1

0

0

0

0

1

18182

.

0

0

0

0

0

1

5

.

0

0

0

0

0

1

L

.

0

0

0

0

30769

.

2

46158

.

6

0

0

0

1

54545

.

8

36364

.

2

0

0

0

3

2

5

.

5

0

0

0

2

1

4

U

background image

10

1

0

0

0

4

8

4

0

0

1

8

2

1

0

0

3

1

5

2

0

0

2

1

4

A

j=5

3

i

1

i

4

i

ji

4

j

44

4

j

u

l

a

u

1

l

czwarta kolumna macierzy L:

background image

1

15476

.

0

0

0

0

0

1

692308

.

1

0

0

0

0

1

18182

.

0

0

0

0

0

1

5

.

0

0

0

0

0

1

L

.

0

0

0

0

30769

.

2

46158

.

6

0

0

0

1

54545

.

8

36364

.

2

0

0

0

3

2

5

.

5

0

0

0

2

1

4

U

background image

Dla sprawdzenia czy nie popełniliśmy błędu obliczamy: B=LU

1

15476

.

0

0

0

0

0

1

692308

.

1

0

0

0

0

1

18182

.

0

0

0

0

0

1

5

.

0

0

0

0

0

1

L

35714

.

10

0

0

0

0

30769

.

2

46158

.

6

0

0

0

1

54545

.

8

36364

.

2

0

0

0

3

2

5

.

5

0

0

0

2

1

4

U

10

1

0

0

0

4

8

4

0

0

1

8

2

1

0

0

3

1

5

2

0

0

2

1

4

B

background image

10

1

0

0

0

4

8

4

0

0

1

8

2

1

0

0

3

1

5

2

0

0

2

1

4

A

10

1

0

0

0

4

8

4

0

0

1

8

2

1

0

0

3

1

5

2

0

0

2

1

4

B

Mając macierz A=LU możemy rozwiązać równanie LUX=b
dla dowolnego wektora prawej strony.

background image

Znając rozkład LU macierzy łatwo obliczyć wyznacznik
główny |A
| macierzy A=LU. Mamy:

U

L

LU

A

ale

1

L

a

N

1

i

ii

u

U

i ostatecznie:

N

i

1

i

ii

u

A

background image

10

1

0

0

0

4

8

4

0

0

1

8

2

1

0

0

3

1

5

2

0

0

2

1

4

A

35714

.

10

0

0

0

0

30769

.

2

46158

.

6

0

0

0

1

54545

.

8

36364

.

2

0

0

0

3

2

5

.

5

0

0

0

2

1

4

U

3480

A

background image

Obliczanie macierzy odwrotnej

Dana macierz:

10

1

0

0

0

4

8

4

0

0

1

8

2

1

0

0

3

1

5

2

0

0

2

1

4

A

Macierz odwrotna : AA

-1

=1

i A

-1

A=1

Oznaczając: X=A

-1

mamy N układów N równań liniowych:
AX=1

Metoda Gaussa - Jordana

background image

Dla określenia macierzy odwrotnej X mamy równanie:

1

0

0

0

0

0

1

0

0

0

0

0

1

0

0

0

0

0

1

0

0

0

0

0

1

x

x

x

x

x

x

x

x

x

x

x

x

x

x

x

x

x

x

x

x

x

x

x

x

x

10

1

0

0

0

4

8

4

0

0

1

8

2

1

0

0

3

1

5

2

0

0

2

1

4

55

54

53

52

51

45

44

43

42

41

35

34

33

32

31

25

24

23

22

21

15

14

13

12

11

background image

Zapisujemy w postaci tablicy dołączonej:

1

0

0

0

0

10

1

0

0

0

0

1

0

0

0

4

8

4

0

0

0

0

1

0

0

1

8

2

1

0

0

0

0

1

0

0

3

1

5

2

0

0

0

0

1

0

0

2

1

4

i procedura eliminacji Gaussa – Jordana:

1

0

0

0

0

10

1

0

0

0

0

1

0

0

0

4

8

4

0

0

0

0

1

0

0

1

8

2

1

0

0

0

0

1

5

.

0

0

3

2

5

.

5

0

0

0

0

0

25

.

0

0

0

5

.

0

25

.

0

1

background image

1

0

0

0

0

10

1

0

0

0

0

1

0

0

0

4

8

4

0

0

0

0

1

18182

.

0

090909

.

0

1

5454

.

8

3636

.

2

0

0

0

0

0

18182

.

0

090909

.

0

0

54545

.

0

36364

.

0

1

0

0

0

0

045455

.

0

22727

.

0

0

13636

.

0

40909

.

0

0

1

Ponieważ pierwsze dwie kolumny już nie ulegną zmianie
dlatego ze względu na oszczędność miejsca zostaną usunięte

1

0

0

0

0

10

1

0

0

0

0

1

0

0

0

4

8

4

0

0

0

0

1

0

0

1

8

2

1

0

0

0

0

1

5

.

0

0

3

2

5

.

5

0

0

0

0

0

25

.

0

0

0

5

.

0

25

.

0

1

background image

1

0

0

0

0

10

1

0

0

1

0

0

0

4

8

4

0

0

1

18182

.

0

090909

.

0

1

5454

.

8

3636

.

2

0

0

0

18182

.

0

090909

.

0

0

54545

.

0

36364

.

0

0

0

0

045455

.

0

22727

.

0

0

13636

.

0

40909

.

0

1

0

0

0

0

10

1

0

0

1

6923

.

1

3077

.

0

15385

.

0

3077

.

2

4616

.

6

0

0

0

42308

.

0

076925

.

0

038462

.

0

42308

.

0

6154

.

3

1

0

0

15385

.

0

15385

.

0

076923

.

0

15385

.

0

76925

.

0

0

0

0

17308

.

0

076924

.

0

21154

.

0

17308

.

0

6154

.

1

0

Pomijamy pierwszą kolumnę

background image

1

15476

.

0

2619

.

0

04762

.

0

02381

.

0

357

.

10

0

0

15476

.

0

2619

.

0

04762

.

0

02381

.

0

35714

.

0

1

0

55952

.

0

52379

.

0

09524

.

0

047621

.

0

7143

.

1

0

0

11905

.

0

47617

.

0

19048

.

0

095239

.

0

42858

.

0

0

0

25

.

0

24999

.

0

000001348

.

0

25

.

0

75

.

0

0

Pomijamy pierwszą kolumnę:

1

0

0

0

0

10

1

0

1

6923

.

1

3077

.

0

15385

.

0

3077

.

2

4616

.

6

0

0

42308

.

0

076925

.

0

038462

.

0

42308

.

0

6154

.

3

0

0

15385

.

0

15385

.

0

076923

.

0

15385

.

0

76925

.

0

0

0

17308

.

0

076924

.

0

21154

.

0

17308

.

0

6154

.

1

background image

096553

.

0

014943

.

0

025287

.

0

0045979

.

0

0022989

.

0

1

034483

.

0

14942

.

0

25287

.

0

045978

.

0

022989

.

0

0

16552

.

0

5339

.

0

48044

.

0

087358

.

0

04368

.

0

0

041381

.

0

11265

.

0

036779

.

0

18851

.

0

094254

.

0

0

072415

.

0

23879

.

0

23102

.

0

0034471

.

0

24828

.

0

0

i otrzymujemy macierz odwrotną:

1

15476

.

0

2619

.

0

04762

.

0

02381

.

0

357

.

10

0

15476

.

0

2619

.

0

04762

.

0

02381

.

0

35714

.

0

0

55952

.

0

52379

.

0

09524

.

0

047621

.

0

7143

.

1

0

11905

.

0

47617

.

0

19048

.

0

095239

.

0

42858

.

0

0

25

.

0

24999

.

0

000001348

.

0

25

.

0

75

.

0

background image

Sprawdzamy poprawność obliczonej macierzy odwrotnej
obliczając AA

-1

10

1

0

0

0

4

8

4

0

0

1

8

2

1

0

0

3

1

5

2

0

0

2

1

4

A

096553

.

0

014943

.

0

025287

.

0

0045979

.

0

0022989

.

0

034483

.

0

14942

.

0

25287

.

0

045978

.

0

022989

.

0

16552

.

0

5339

.

0

48044

.

0

087358

.

0

04368

.

0

041381

.

0

11265

.

0

036779

.

0

18851

.

0

094254

.

0

072415

.

0

23879

.

0

23102

.

0

0034471

.

0

24828

.

0

background image

1

AX

3

.

1

1

0

1

.

0

0

4

.

0

2

.

1

2

.

5

04

.

0

36

.

0

4

.

0

3

.

3

4

.

1

01

.

0

09

.

0

4

.

0

3

5

.

2

02

.

2

3

.

0

1

.

0

1

1

.

2

56

.

0

4

.

1

10

5

Macierz odwrotną można również obliczyć korzystając z
rozkładu LU

Niech A=LU

mamy rozwiązać N układów N równań algebraicznych:

LUX=1

oznaczając:

Y=UX

mamy:

LY=1

background image

Postępowanie jest proste:

Krok pierwszy – rozwiązujemy N - krotnie układ N równań
z dolną macierzą trójkątną L
wyznaczając Y:

LY=1

Krok drugi – rozwiązujemy N – krotnie układ N równań
z górną macierzą trójkątną U
wyznaczając macierz odwrotną
A

-1

=X:

UX=Y

background image

Dana macierz:

10

1

0

0

0

4

8

4

0

0

1

8

2

1

0

0

3

1

5

2

0

0

2

1

4

A

i

1

15476

.

0

0

0

0

0

1

692308

.

1

0

0

0

0

1

18182

.

0

0

0

0

0

1

5

.

0

0

0

0

0

1

L

35714

.

10

0

0

0

0

30769

.

2

46158

.

6

0

0

0

1

54545

.

8

36364

.

2

0

0

0

3

2

5

.

5

0

0

0

2

1

4

U

background image

LY=1

Równanie

jest

1

0

0

0

0

0

1

0

0

0

0

0

1

0

0

0

0

0

1

0

0

0

0

0

1

y

y

y

y

y

y

y

y

y

y

y

y

y

y

y

y

y

y

y

y

y

y

y

y

y

1

15476

.

0

0

0

0

0

1

692308

.

1

0

0

0

0

1

18182

.

0

0

0

0

0

1

5

.

0

0

0

0

0

1

55

54

53

52

51

45

44

43

42

41

35

34

33

32

31

25

24

23

22

21

15

14

13

12

11

background image

Macierz odwrotna do dolnej trójkątnej też jest macierzą dolną

trójkątną i w przypadku macierzy L główna przekątna

to 1 czyli

1

0

0

0

0

0

1

0

0

0

0

0

1

0

0

0

0

0

1

0

0

0

0

0

1

1

y

y

y

y

0

1

y

y

y

0

0

1

y

y

0

0

0

1

y

0

0

0

0

1

1

15476

.

0

0

0

0

0

1

692308

.

1

0

0

0

0

1

18182

.

0

0

0

0

0

1

5

.

0

0

0

0

0

1

54

53

52

51

43

42

41

32

31

21

background image

Pozostałe wyrazy macierzy Y wyznaczamy rozpoczynając
od pierwszej kolumny i kolejno następne:

1

15476

.

0

2619

.

0

047619

.

0

023809

.

0

0

1

6923

.

1

3077

.

0

15385

.

0

0

0

1

18182

.

0

09091

.

0

0

0

0

1

5

.

0

0

0

0

0

1

1

L

Y

Pozostaje do rozwiązania równanie:

UX=Y

background image

55

54

53

52

51

45

44

43

42

41

35

34

33

32

31

25

24

23

22

21

15

14

13

12

11

x

x

x

x

x

x

x

x

x

x

x

x

x

x

x

x

x

x

x

x

x

x

x

x

x

35714

.

10

0

0

0

0

30769

.

2

46158

.

6

0

0

0

1

54545

.

8

36364

.

2

0

0

0

3

2

5

.

5

0

0

0

2

1

4

1

15476

.

0

2619

.

0

047619

.

0

023809

.

0

0

1

6923

.

1

3077

.

0

15385

.

0

0

0

1

18182

.

0

09091

.

0

0

0

0

1

5

.

0

0

0

0

0

1

Startujemy kolejno od pierwszej kolumny kolejno do piątej,
a niewiadome w kolumnach wyznaczamy od ostatniej tj. x

Nk

background image

096552

.

0

014942

.

0

025287

.

0

0045977

.

0

0022988

.

0

034483

.

0

14942

.

0

25287

.

0

045978

.

022989

.

0

16552

.

0

53389

.

0

48045

.

0

087359

.

0

043679

.

0

04138

.

0

11264

.

0

03678

.

0

18851

.

0

094253

.

0

072415

.

0

23878

.

0

23103

.

0

003448

.

0

24828

.

0

1

X

A

Dla porównania macierz odwrotna obliczona

metodą Gaussa - Jordana

096553

.

0

014943

.

0

025287

.

0

0045979

.

0

0022989

.

0

034483

.

0

14942

.

0

25287

.

0

045978

.

0

022989

.

0

16552

.

0

5339

.

0

48044

.

0

087358

.

0

04368

.

0

041381

.

0

11265

.

0

036779

.

0

18851

.

0

094254

.

0

072415

.

0

23879

.

0

23102

.

0

0034471

.

0

24828

.

0

background image

Metody iteracyjne

Metoda iteracji prostej:

b

Ax

Macierz A przedstawiamy w postaci:

G

L

D

A

gdzie jeżeli

NN

Nk

2

N

1

N

kN

kk

2

k

1

k

N

2

k

2

22

21

N

1

k

1

12

11

a

.

a

.

a

a

.

.

.

.

.

.

a

.

a

.

a

a

.

.

.

.

.

.

a

.

a

.

a

a

a

.

a

.

a

a

A

to

background image

0

.

a

.

a

a

.

.

.

.

.

.

0

.

0

.

a

a

.

.

.

.

.

.

0

.

0

.

0

a

0

.

0

.

0

0

D

Nk

2

N

1

N

2

k

1

k

21

NN

kk

22

11

a

.

0

.

0

0

.

.

.

.

.

.

0

.

a

.

0

0

.

.

.

.

.

.

0

.

0

.

a

0

0

.

0

.

0

a

L

0

.

0

.

0

0

.

.

.

.

.

.

a

.

0

.

0

0

.

.

.

.

.

.

a

.

a

.

0

0

a

.

a

.

a

0

G

kN

N

2

k

2

N

1

k

1

12

background image

i równanie:

b

Ax

piszemy:

ponieważ

G

L

D

A

b

Gx

Lx

Dx

Ax

lub

x

G

D

b

Lx

i mamy:

x

G

D

L

b

L

x

1

1

Otrzymujemy algorytm:

n

1

1

1

n

x

G

D

L

b

L

x

Warunek zbieżności

1

G

D

L

1

background image

Przykład:

5

10

0

30

x

x

x

x

4

1

2

0

1

10

0

2

.

0

2

0

5

1

0

2

.

0

1

4

4

3

2

1

0

1

2

0

0

0

0

2

.

0

0

0

0

1

0

0

0

0

D

4

0

0

0

0

10

0

0

0

0

5

0

0

0

0

4

L

0

0

0

0

1

0

0

0

2

0

0

0

0

2

.

0

1

0

G

25

.

0

0

0

0

0

1

.

0

0

0

0

0

2

.

0

0

0

0

0

25

.

0

L

1

background image

podstawiając

n

1

1

1

n

x

G

D

L

b

L

x

25

.

1

1

0

5

.

7

b

L

1

0

1

2

0

1

0

0

2

.

0

2

0

0

1

0

2

.

0

1

0

25

.

0

0

0

0

0

1

.

0

0

0

0

0

2

.

0

0

0

0

0

25

.

0

G

D

L

1

0

25

.

0

5

.

0

0

1

.

0

0

0

02

.

0

4

.

0

0

0

2

.

0

0

05

.

0

25

.

0

0

G

D

L

1

background image

Niech zerowe przybliżenie

0

0

0

0

x

0

dla pierwszego mamy:

25

.

1

1

0

5

.

7

x

1

Podstawiamy do równania:

n

1

1

1

n

x

G

D

L

b

L

x

i znajdujemy

1

725

.

2

45

.

7

x

2

Badamy normę:

475

.

2

25

.

0

175

.

0

2

05

.

0

x

x

1

1

2

024

.

2

x

x

2

1

2

2

x

x

1

2

background image

Po 9-ciu krokach mamy:

x

9

8.112

2.565

0.602

2.383

x

10

8.111

2.576

0.599

2.382

=

X

d

8.1147767
2.5788529
0.5987301
2.3897439

Rozwiązanie dokładne:

Normy:

016

.

0

001

.

0

003

.

0

011

.

0

001

.

0

x

x

1

9

10

01149

.

0

x

x

2

9

10

011

.

0

x

x

9

10

background image

x

19

8.114772

2.578831

0.598735

2.389734

=

X

d

8.1147767
2.5788529
0.5987301
2.3897439

x

20

8.114771

2.578848

0.598731

2.389732

000024

.

0

x

x

1

19

20

00001761

.

0

x

x

2

19

20

000017

.

0

x

x

19

20

background image

Metoda Gaussa-Seidela

n

1

1

1

n

Gx

D

L

b

D

L

x

Przykład:

5

10

0

30

x

x

x

x

4

1

2

0

1

10

0

2

.

0

2

0

5

1

0

2

.

0

1

4

4

3

2

1

0

1

2

0

0

0

0

2

.

0

0

0

0

1

0

0

0

0

D

4

0

0

0

0

10

0

0

0

0

5

0

0

0

0

4

L

0

0

0

0

1

0

0

0

2

0

0

0

0

2

.

0

1

0

G

25

.

0

0

0

0

0

1

.

0

0

0

0

0

2

.

0

0

0

0

0

25

.

0

L

1

background image

A=(L+D)

-1

n

1

n

Gx

A

Ab

x

lub w postaci równań

1

n

,

3

1

n

,

2

1

n

,

4

n

,

4

1

n

,

1

1

n

,

3

n

,

4

1

n

,

1

1

n

,

2

n

,

3

n

,

2

1

n

,

1

x

25

.

0

x

5

.

0

25

.

1

x

x

1

.

0

x

02

.

0

1

x

x

4

.

0

x

2

.

0

x

x

05

.

0

x

25

.

0

5

.

7

x

A

0.25
0.05

5 10

3


0.026

0

0.2

0

0.1

0
0

0.1

0.025

0
0

0

0.25

background image







0

0

0

0

x

0

Zerowe przybliżenie

x

1

7.5
1.5
0.85

1.7875





x

2

7.8325
2.2815

0.6646

2.2246





7365

.

1

x

x

1

1

2

x

4

8.093477

2.55647

0.603687

2.377313





x

5

8.108933
2.572712

0.60009


2.386333





04432

.

0

x

x

1

4

5

background image

x

9

8.114744
2.578818

0.598738

2.389725





x

10

8.114768
2.578843

0.598732

2.389739





5

1

9

10

10

889

.

6

x

x

x

14

8.114777
2.578853

0.59873


2.389744





x

15

8.114777
2.578853

0.59873


2.389744





7

1

14

15

10

071

.

1

x

x

background image

D L

G

(

) x

15

30

4.35972

10

8

10

5

Sprawdzian:

5

10

0

30

x

x

x

x

4

1

2

0

1

10

0

2

.

0

2

0

5

1

0

2

.

0

1

4

4

3

2

1

background image

Interpolacja funkcji

Dane wartości funkcji y

n

w punktach x

n

, gdzie n=0,1,2, ....N-1.

x

y

x

0

y

0

x

n

y

n

x

N-1

y

N-1

background image

Interpolacja wielomianowa

Twierdzenie

Istnieje dokładnie jeden wielomian stopnia co najwyżej N (N>=0),

który w punktach x

0

, x

1

,...,x

N-1

przyjmuje wartości y

0

,y

1

,...,y

N-1

.

Wzór interpolacyjny Lagrange'a:

)

x

(

y

....

)

x

(

y

)

x

(

y

)

x

(

W

1

N

1

N

1

1

0

0

n

gdzie

1

-

0,1,...N

=

k

dla

)

x

(

k

jest wielomianem stopnia co najwyżej N.

background image

Z warunku interpolacyjnego:

1

-

N

0,1,....,

=

k

dla

y

)

x

(

W

k

k

N

powyższy układ N równań można najprościej rozwiązać
przyjmując dla wielomianów

k

(x) następujące warunki :

i

k

dla

1

i

k

dla

0

)

x

(

i

k

jako wielomian

k

(x) należy wybrać taki, który ma miejsca

zerowe we wszystkich punktach interpolacji

z wyjątkiem punktu x

k

, w którym funkcja ma wartość 1

,

1

N

1

k

1

k

1

0

x

,...,

x

,

x

,...,

x

,

x

Rozwiązaniem jest wielomian :

background image

Rozwiązaniem jest wielomian:

)

x

x

)...(

x

x

)(

x

x

)...(

x

x

)(

x

x

(

A

)

x

(

1

N

1

k

1

k

1

0

k

z warunku:

otrzymuje się:

1

)

x

(

k

k

)

x

x

)...(

x

x

)(

x

x

)...(

x

x

)(

x

x

(

1

A

1

N

k

k

1

k

k

1

k

1

k

0

k

Wielomian Lagrange'a przyjmuje postać:

)

x

x

)...(

x

x

)(

x

x

)...(

x

x

)(

x

x

(

)

x

x

)...(

x

x

)(

x

x

)...(

x

x

)(

x

x

(

y

)

x

(

W

1

N

k

1

k

k

1

k

k

1

k

0

k

1

N

1

k

1

k

1

0

k

N

background image

Ocena błędu interpolacji:

)

x

(

f

sup

M

)

x

x

(

)

x

(

gdzie

)

x

(

)!

1

N

(

M

)

x

(

W

)

x

(

f

)

1

N

(

]

b

,

a

[

x

1

N

N

k

0

k

k

1

N

1

N

1

n

N

background image

Przykład 1.

Zbudować wielomian interpolacyjny dla funkcji exp(x)
w przedziale [1,2] bazując na 5 węzłach interpolacyjnych.

Wybierzmy węzły równomiernie czyli

25

.

0

4

1

2

x

x

i

1.0

1.25

1.50

1.75

2.0

y

i

2.7182
8

3.4903
4

4.4816
9

5.7546 7.3890

6

mamy:

background image

Wielomian Lagrange’a jest:





























































75

.

1

2

5

.

1

2

25

.

1

2

1

2

75

.

1

x

5

.

1

x

25

.

1

x

1

x

38906

.

7

2

75

.

1

5

.

1

75

.

1

25

.

1

75

.

1

1

75

.

1

2

x

5

.

1

x

25

.

1

x

1

x

7546

.

5

2

5

.

1

75

.

1

5

.

1

25

.

1

5

.

1

1

5

.

1

2

x

75

.

1

x

25

.

1

x

1

x

48169

.

4

2

25

.

1

75

.

1

25

.

1

5

.

1

25

.

1

1

25

.

1

2

x

75

.

1

x

5

.

1

x

1

x

49034

.

3

2

1

75

.

1

1

5

.

1

1

25

.

1

1

2

x

75

.

1

x

5

.

1

x

25

.

1

x

71828

.

2

)

x

(

W

4

background image

lub































75

.

1

x

5

.

1

x

25

.

1

x

1

x

817

.

78

2

x

5

.

1

x

25

.

1

x

1

x

53

.

245

2

x

75

.

1

x

25

.

1

x

1

x

83

.

286

2

x

75

.

1

x

5

.

1

x

1

x

92

.

148

2

x

75

.

1

x

5

.

1

x

25

.

1

x

995

.

28

)

x

(

W

4

Wyniki obliczeń przedstawiono na wykresie:

background image

1

1.2

1.4

1.6

1.8

2

2

3.2

4.4

5.6

6.8

8

w x

( )

expx

( )

x

Dla lepszej oceny wykres błędu względnego:

100

)

x

exp(

)

x

exp(

)

x

(

w

)

x

(

eps

4

background image

1

1.2

1.4 1.6

1.8

2

0

0.002

0.004

0.006

0.008

0.01

eps x

( )

x

background image

Przykład 2.
W wyniku pomiarów zdjęto pierwotną krzywą magnesowania
B=F(H). Zbudować wielomian interpolacyjny Lagrange'a
dla zakresu 0<=H <=3000A/m.

H[A/m]

0

50

100 200 500 100

0

150

0

200

0

3000

B[T]

0

0.7

5

1.5 1.8 1.9

5

2.0 2.02 2.03 2.035

Kolejne wielomiany

k

(H) dla k=0,1,...8 są:

  

























3000

0

2000

0

1500

0

1000

0

3000

H

2000

H

1500

H

1000

H

500

0

200

0

100

0

50

0

500

H

200

H

100

H

50

H

H

0

lub po obliczeniu mianownika mamy:

background image

 













3000

H

2000

H

1500

H

1000

H

500

H

200

H

100

H

50

H

10

2222

.

2

H

22

0

 











3000

H

2000

H

1500

H

1000

H

500

H

200

H

100

H

H

10

4784

.

7

H

22

1

 











3000

H

2000

H

1500

H

1000

H

500

H

200

H

50

H

H

10

2019

.

7

H

22

2

 











3000

H

2000

H

1500

H

1000

H

500

H

100

H

50

H

H

10

1198

.

2

H

22

3

 











3000

H

2000

H

1500

H

1000

H

200

H

100

H

50

H

H

10

9753

.

1

H

23

4

 











3000

H

2000

H

1500

H

500

H

200

H

100

H

50

H

H

10

924

.

2

H

24

5

background image

 











3000

H

2000

H

1000

H

500

H

200

H

100

H

50

H

H

10

7366

.

6

H

25

6

 











3000

H

1500

H

1000

H

500

H

200

H

100

H

50

H

H

10

9965

.

9

H

26

7

 











2000

H

1500

H

1000

H

500

H

200

H

100

H

50

H

H

10

8554

.

1

H

26

8

i wielomian aproksymacyjny jest

 

 

 

 

 

 

 

 

 

 

H

035

.

2

H

03

.

2

H

02

.

2

H

2

H

95

.

1

H

8

.

1

H

5

.

1

H

75

.

0

H

0

H

B

8

7

6

5

4

3

2

1

0

lub

 

 

 

 

 

 

 

 

 

H

035

.

2

H

03

.

2

H

02

.

2

H

2

H

95

.

1

H

8

.

1

H

5

.

1

H

75

.

0

H

B

8

7

6

5

4

3

2

1

background image

0

600 1200 1800 2400 3000

4000

2800

1600

400

800

2000

B H

( )

H

Otrzymany wynik jest niemożliwy do przyjęcia!!!

background image

Interpolacja liniowa odcinkami:

H[A/m]

0

50

100 200 500 100

0

150

0

200

0

3000

B[T]

0

0.7

5

1.5 1.8 1.9

5

2.0 2.02 2.03 2.035

 

0

50

0

H

75

.

0

50

0

50

H

0

H

B

1

dla

50

H

0

lub po wykonaniu działań:

 

H

015

.

0

H

B

1

50

H

0

dla

i podobnie:

 

50

H

03

.

0

100

H

015

.

0

H

B

2

dla

100

H

50

 

100

H

018

.

0

200

H

015

.

0

H

B

3

dla

200

H

100

 

200

H

0065

.

0

500

H

006

.

0

H

B

4

500

H

200

dla

background image

 

500

H

004

.

0

1000

H

0039

.

0

H

B

5

dla

1000

H

500

 

1000

H

00404

.

0

1500

H

004

.

0

H

B

6

dla

1500

H

1000

 

1500

H

00406

.

0

2000

H

00404

.

0

H

B

7

dla

2000

H

1500

 

2000

H

002035

.

0

3000

H

00203

.

0

H

B

8

dla

3000

H

2000

background image

0

600 1200 1800 2400 3000

0

0.6

1.2

1.8

2.4

3

BaH

( )

H

B(H)

background image

0

100 200 300 400 500

0

0.4

0.8

1.2

1.6

2

BaH

( )

B H

( )

H

Porównanie Ba(H) – interpolacja liniowa
B(H) – wielomian 8-go stopnia

background image

Optymalny dobór węzłów interpolacji.

Dobrać węzły interpolacji tak aby kres górny wielomianu

)

x

(

sup

1

n

]

b

,

a

[

x

był jak najmniejszy.

Rozwiązanie otrzymuje się za pomocą wielomianów Czebyszewa.

Są to wielomiany zdefiniowane na przedziale x [-1,1]

i są zdefiniowane :

)

x

arccos

n

cos(

)

x

(

T

n

Przykładowe wykresy dla n=1,2,3,4:

background image

1

0.5

0

0.5

1

1

0.67

0.33

0

0.33

0.67

1

T 1 x

(

)

T 2 x

(

)

T 3 x

(

)

T 4 x

(

)

x

background image

Wielomiany spełniają następujące związki:

x

3

x

4

)

x

(

T

1

x

2

)

x

(

T

.

2,3,4,....

=

n

dla

)

x

(

T

)

x

(

xT

2

)

x

(

T

x

)

x

cos(arccos

)

x

(

T

1

)

x

(

T

3

3

2

2

2

n

1

n

n

1

0

Każdy z wielomianów ma n różnych pierwiastków określonych
zależnością:

1

-

n

0,1,2,...,

=

m

gdzie

)

n

2

1

m

2

cos(

x

m

w przedziale [-1,1].

background image

Dowodzi się, że jeżeli dla przedziału [-1,1] dobrać pierwiastki

zgodnie z zależnością
określającą pierwiastki wielomianu Czebyszewa to zachodzi:

)

n

2

1

m

2

cos(

x

m

n

1

n

]

1

,

1

[

x

1

n

n

1

n

2

1

)

x

(

sup

)

x

(

T

2

1

)

x

(

stąd wynika

background image

czyli ocena błędu w przedziale [-1,1] jest:

)!

1

n

(

2

M

)

x

(

W

)

x

(

f

n

1

n

n

Problem jest rozwiązany w przedziale [-1,1].

Aby go rozwiązać w przedziale [a,b] należy dokonać odwzorowania

przedziału [a,b] na przedział [-1,1].

Niech

]

1

,

1

[

z

a

]

b

,

a

[

x

otrzymujemy

a

b

a

x

2

1

z

-1

1

z

a

b

x

background image

i stąd mamy:

)

1

z

)(

a

b

(

5

.

0

a

x

a

b

a

x

2

1

z

Dla przedziału [a,b] należy dla optymalnej interpolacji
wybrać punkty według zależności:

1

-

n

0,1,2,...,

=

k

dla

1

)

n

2

1

k

2

cos(

)

a

b

(

5

.

0

a

x

k





Ocena błędu przyjmuje postać:

1

n

2

1

n

1

n

n

2

)

a

b

(

)!

1

n

(

M

)

x

(

W

)

x

(

f

background image

Przykład
Interpolacja funkcji exp(x) na przedziale [0,10] przy pięciu węzłach
i ich optymalnym wyborze.

Mamy a=0, b=10 i zgodnie ze wzorem:

1

-

n

0,1,2,...,

=

k

dla

1

)

n

2

1

k

2

cos(

)

a

b

(

5

.

0

a

x

k





mamy: n=4 oraz

0,1,2,3,4

=

k

dla

1

)

10

1

k

2

cos(

5

x

k





obliczając: x

0

=9.7553; x

1

=7.9389; x

2

=5.0; x

3

=2.0611;

x

4

=0.24472

i odpowiednio y=exp(x) mamy: y

0

=17245; y

1

=2804.3; y

2

=148.41;

y

3

=7.8546; y

4

=1.2773

background image

Kolejne wielomiany Lagrange’a są:

 







24472

.

0

x

0611

.

2

x

5

x

9389

.

7

x

0015821

.

0

x

0

 







24472

.

0

x

0611

.

2

x

5

x

7553

.

9

x

0041422

.

0

x

1

 







24472

.

0

x

0611

.

2

x

9389

.

7

x

7553

.

9

x

0051201

.

0

x

2

 







24472

.

0

x

5

x

9389

.

7

x

7553

.

9

x

0041429

.

0

x

3

 







0611

.

2

x

5

x

9389

.

7

x

7553

.

9

x

0015822

.

0

x

4

Wielomian interpolacyjny W(x) jest:

 

 

 

 

 

 

x

2773

.

1

x

8546

.

7

x

41

.

148

x

3

.

2804

x

17245

x

W

4

3

2

1

0

background image

0

2

4

6

8

10

110

4

2000

6000

1.410

4

2.210

4

310

4

w z

( )

expz

( )

z

background image

Błąd względny:

 

 

 

 

100

x

epx

x

epx

x

W

x

eps

0

2

4

6

8

10

0

60

120

180

240

300

eps z

( )

z

background image

Funkcje sklejane (spline)

Dane wartości funkcji y

n

w punktach x

n

, gdzie n=0,1,2, ....N-1.

x

y

a=x

0

y

0

x

k

y

k

b=x

N

y

N

przy czym: a=x

0

<x

1

<...<x

k

<...<x

N

=b

s

m

(x)

background image

Funkcję s

m

(x) określoną na przedziale [a,b] nazywamy

funkcją sklejaną stopnia m1, jeżeli:

1. s

m

(x) jest funkcją co najwyżej stopnia m na każdym

z podprzedziałów [x

k

,x

k+1

] gdzie k=0,1,...,N-1

2. funkcja s

m

(x) posiada ciągłą pochodną do rzędu m-1

w przedziale [a,b].

Najpowszechniej stosowane funkcje sklejane 3-go stopnia.

Z drugiego warunku wynika, że funkcja s

3

(x) musi mieć

ciągłą drugą pochodną w przedziale [a,b].

background image

Niech:

 

k

k

3

M

x

s

dla k=0,1,...,N

Druga pochodna będzie ciągła, jeżeli dla przedziału [x

k

,x

k+1

]

przyjmiemy:

 

1

k

k

k

k

1

k

k

1

k

k

3

x

,

x

x

dla

x

x

M

x

x

M

x

s



oraz k=0,1,...,N-1

gdzie

k

1

k

k

x

x

Całkując mamy:

 

1

k

k

k

k

2

k

1

k

k

2

1

k

k

3

x

,

x

x

dla

A

2

x

x

M

2

x

x

M

x

s

background image

i całkując pierwszą pochodną mamy:

 

1

k

k

k

k

k

k

3

k

1

k

k

3

1

k

k

3

x

,

x

x

dla

B

x

x

A

6

x

x

M

6

x

x

M

x

s

Z warunków interpolacji mamy:

 

N

0,1,...,

k

dla

y

x

s

k

k

3

Dla każdego z N przedziałów piszemy:

 

y

B

x

x

A

6

x

x

M

x

s

y

B

6

x

x

M

x

s

1

k

k

k

1

k

k

k

3

k

1

k

1

k

1

k

3

k

k

k

3

k

1

k

k

k

3

background image

ale

k

1

k

k

x

x

czyli

 

y

B

x

x

A

6

x

x

M

x

s

y

B

6

x

x

M

x

s

1

k

k

k

1

k

k

k

3

k

1

k

1

k

1

k

3

k

k

k

3

k

1

k

k

k

3

możemy zapisać:

 

y

B

A

6

M

x

s

y

B

6

M

x

s

1

k

k

k

k

2

k

1

k

1

k

3

k

k

2

k

k

k

3

stąd mamy:

6

M

y

B

2

k

k

k

k

i

6

M

M

y

y

A

k

k

1

k

k

k

1

k

k

background image

Pozostaje tylko zapewnić ciągłość pierwszej pochodnej

w węzłach interpolacji x

k

(k=1,2,...,N-1) czyli mamy N-1

równań:

k

3

k

3

x

s

x

s

stąd:

1

-

N

1,2,...,

k

dla

A

2

x

x

M

2

x

x

M

A

2

x

x

M

2

x

x

M

k

k

2

k

k

1

k

k

2

k

1

k

k

1

k

1

k

2

1

k

k

k

1

k

2

k

k

1

k

Biorąc pod uwagę, że

k

1

k

k

x

x

i

6

M

M

y

y

A

k

k

1

k

k

k

1

k

k

background image

mamy:

1

k

1

k

k

k

k

1

k

k

k

k

1

k

k

1

k

1

k

y

y

y

y

6

M

3

M

6

M

dla k=1,2,...,N-1

mamy N-1 równań, a niewiadomych M

k

: N+1.

Potrzebne dwa dodatkowe równania, które przyjmuje się
w jednej z trzech postaci:

1. Dana jest pochodna funkcji w punkcie a i b czyli:

 

 

N

3

0

3

p

b

s

i

p

a

s

background image

2. Dana jest druga pochodna na obu końcach przedziału, czyli:

 

 

N

3

0

3

d

b

s

i

d

a

s



3. Funkcja jest okresowa o okresie b-a i wtedy mamy:

 

 

 

 

 

 

b

s

a

s

b

s

a

s

b

s

a

s

3

3

3

3

3

3





W przypadku 1 mamy więc:

0

0

0

2

0

0

1

0

2

0

1

0

0

3

p

A

2

x

x

M

2

x

x

M

x

a

s

lub

0

0

0

1

0

1

0

p

y

y

6

M

M

2

background image

i podobnie dla drugiego końca przedziału mamy:

N

1

N

1

N

2

1

N

N

N

1

N

2

N

N

1

N

N

3

p

A

2

x

x

M

2

x

x

M

b

x

s

lub

1

N

1

N

N

N

1

N

N

1

N

y

y

p

6

M

2

M

i otrzymujemy układ N+1 równań:

0

0

0

1

0

1

0

p

y

y

6

M

M

2

1

k

1

k

k

k

k

1

k

k

k

k

1

k

k

1

k

1

k

y

y

y

y

6

M

3

M

6

M

dla k=1,2,...,N-1

1

N

1

N

N

N

1

N

N

1

N

y

y

p

6

M

2

M

background image

Dla drugiego typu warunków mamy:

 

 

N

3

0

3

d

b

s

i

d

a

s





czyli z definicji:

 

1

k

k

k

k

1

k

k

1

k

k

3

x

,

x

x

dla

x

x

M

x

x

M

x

s



wynika, że

0

0

d

M

i

N

N

d

M

i N-1 równań:

1

k

1

k

k

k

k

1

k

k

k

k

1

k

k

1

k

1

k

y

y

y

y

6

M

3

M

6

M

dla k=1,2,...,N-1

Zapiszemy otrzymane równania w postaci wygodniejszej
do obliczeń:

background image

N

N

N

1

N

N

k

1

k

k

k

k

1

k

k

1

2

1

1

1

0

1

0

1

0

0

0

f

M

c

M

b

......

..........

..........

f

M

c

M

b

M

a

..

..........

..........

..........

f

M

c

M

b

M

a

f

M

c

M

b

Takie układy rozwiązujemy techniką „wymiatania”.

Niech

1

1

1

0

K

M

L

M

Podstawiając do drugiego równania i wyznaczając M

1

mamy:

1

1

1

1

1

1

1

1

1

2

2

1

b

L

a

K

a

f

b

L

a

M

c

M

background image

lub oznaczając:

1

1

1

1

1

1

2

1

1

1

2

2

b

L

a

K

a

f

K

b

L

a

c

L

możemy zapisać w postaci:

2

2

2

1

K

M

L

M

Wyznaczmy związek rekurencyjny spełniany przez L

k

i K

k

.

Zakładając, że

k

k

k

1

k

K

M

L

M

i podstawiając do k-go równania mamy:

k

k

k

k

k

k

k

k

k

1

k

k

k

b

L

a

K

a

f

b

L

a

M

c

M

i oznaczając:

background image

k

k

k

k

k

k

1

k

k

k

k

k

1

k

b

L

a

K

a

f

K

b

L

a

c

L

mamy:

1

k

1

k

1

k

k

K

M

L

M

a więc współczynniki L

k

,K

k

potrafimy wyznaczyć z zależności

rekurencyjnej:

k

k

k

k

k

k

1

k

k

k

k

k

1

k

b

L

a

K

a

f

K

b

L

a

c

L

jeżeli znamy wartości startowe L

1

, K

1

.

background image

Ale porównując:

1

1

1

0

K

M

L

M

z pierwszym równaniem:

N

N

N

1

N

N

k

1

k

k

k

k

1

k

k

1

2

1

1

1

0

1

0

1

0

0

0

f

M

c

M

b

......

..........

..........

f

M

c

M

b

M

a

..

..........

..........

..........

f

M

c

M

b

M

a

f

M

c

M

b

mamy:

0

0

1

b

c

L

i

0

0

1

b

f

K

i na podstawie związku rekurencyjnego:

k

k

k

k

k

k

1

k

k

k

k

k

1

k

b

L

a

K

a

f

K

b

L

a

c

L

wyznaczmy kolejno dla k=1,2,...N-1

background image

Ostatnie równanie przyjmuje postać:

N

N

N

1

N

K

M

L

M

i podstawiając do ostatniego równania:

N

N

N

1

N

N

k

1

k

k

k

k

1

k

k

1

2

1

1

1

0

1

0

1

0

0

0

f

M

c

M

b

......

..........

..........

f

M

c

M

b

M

a

..

..........

..........

..........

f

M

c

M

b

M

a

f

M

c

M

b

mamy:

N

N

N

N

N

N

N

b

L

c

b

K

f

M

background image

Znając M

N

wyznaczamy ze wzoru:

1

k

1

k

1

k

k

K

M

L

M

kolejne M

k

przyjmując k=N-1,N-2,...,0

Przykład:
Dane są wartości funkcji:

x 0 0.3141

6

0.62832 0.9424

8

1.2566
4

1.570
8

y 0 0.3090

2

0.58779 0.8090

2

0.9510
6

1

i jej pierwsze pochodne:

0

y

1

y

5708

.

1

x

0

x

background image

Współczynniki M

0

, M

1

, M

2

, M

3

, M

4

, M

5

wyznaczamy z układu

równań:

0

0

0

1

0

1

0

p

y

y

6

M

M

2

1

k

1

k

k

k

k

1

k

k

k

k

1

k

k

1

k

1

k

y

y

y

y

6

M

3

M

6

M

dla k=1,2,3,4

4

4

5

5

4

5

4

y

y

p

6

M

2

M

gdzie ze względu na równomierny podział mamy:

31416

.

0

4

3

2

1

0

i podstawiając dane mamy:

background image

97520

.

2

M

2

M

65980

.

5

M

M

4

M

81417

.

4

M

M

4

M

49801

.

3

M

M

4

M

83898

.

1

M

M

4

M

31261

.

0

M

M

2

5

4

5

4

3

4

3

2

3

2

1

2

1

0

1

0

Z pierwszego równania mamy:

5

.

0

L

1

1563

.

0

K

1

i

i na mocy związku rekurencyjnego:

k

k

k

k

k

k

1

k

k

k

k

k

1

k

b

L

a

K

a

f

K

b

L

a

c

L

26796

.

0

L

26804

.

0

L

26923

.

0

L

28571

.

0

L

5

4

3

2

22889

.

1

K

07362

.

1

K

80875

.

0

K

48077

.

0

K

5

4

3

2

background image

Z ostatniego równania mamy:

00828

.

1

M

5

i ze związku rekurencyjnego:

1

k

1

k

1

k

k

K

M

L

M

wyznaczamy:

00004

.

0

M

31252

.

0

M

58888

.

0

M

81665

.

0

M

95871

.

0

M

0

1

2

3

4

a następnie stałe A

k

i B

k

ze wzorów:

6

M

y

B

2

k

k

k

k

6

M

M

y

y

A

k

k

1

k

k

k

1

k

k

i

k=0,1,2,3,4

background image

15838

.

0

A

45957

.

0

A

71612

.

0

A

93455

.

0

A

1

A

4

3

2

1

0

96683

.

0

B

82245

.

0

B

59748

.

0

B

31416

.

0

B

0

B

4

3

2

1

0

i ostatecznie mamy:

0,0.31416

x

dla

 

3

3

0

3

x

1658

.

0

x

31416

.

0

00002

.

0

x

x

s

62832

0.31416,0.

x

dla

 

3

3

1

3

31416

.

0

x

3141

.

0

x

62832

.

0

1658

.

0

31416

.

0

x

93455

.

0

31416

.

0

x

s

background image

94248

0.62832,0.

x

dla

 

3

3

2

3

62832

.

0

x

43325

.

0

x

94248

.

0

31241

.

0

62832

.

0

x

71612

.

0

59748

.

0

x

s

25664

0.94248,1.

x

dla

 

3

3

3

3

94248

.

0

x

50861

.

0

x

25664

.

1

43325

.

0

94248

.

0

x

45957

.

0

82245

.

0

x

s

5708

1.25664,1.

x

dla

 

3

3

4

3

25664

.

1

x

53491

.

0

x

5708

.

1

50861

.

0

25664

.

1

x

15838

.

0

96683

.

0

x

s

background image

0

0.31 0.63 0.94 1.26 1.57

0

0.33

0.67

1

s x

( )

sin x

( )

x

background image

0

0.31 0.63 0.94 1.26 1.57

0.01

0.004

0.002

0.008

0.014

0.02

eps x

( )

x

Błąd względny

 

 

 

 

x

sin

x

sin

x

s

x

eps

3

background image

B – funkcje sklejane (B – spline)

Podział przedziału [a,b] jest równomierny, czyli

N

a

b

Funkcja sklejana s

3

(x) jest przyjmowana w postaci:

 

 

1

N

k

1

k

k

k

3

x

B

a

x

s

gdzie funkcje B

k

(x) tworzą bazę przestrzeni s

3

(x) i mają

postać:

background image

x

y

a=x

0

y

0

x

k

=k

y

k

b=N

y

N

s

m

(x)

 

2

k

2

-

k

2

k

1

k

3

2

k

1

k

k

3

1

k

2

1

k

1

k

2

3

k

1

-

k

3

1

k

2

1

k

1

k

2

3

1

k

2

-

k

3

2

k

3

k

x

,

x

x

dla

0

x

,

x

x

dla

x

x

x

,

x

x

dla

x

x

3

x

x

3

x

x

3

x

,

x

x

dla

x

x

3

x

x

3

x

x

3

x

,

x

x

dla

x

x

1

x

B

background image

0

1

2

3

4

5

6

0

1

2

3

4

B 3 t

( )

t

x

B

k

(x)

x

k

x

k-1

x

k-2

x

k-3

x

k+1

x

k+2

x

k+3

background image

x

k-2

x

k-1

x

k

x

k+1

x

k+2

B

k

(x

)

0

1

4

1

0

0

3/

0

-

3/

0

0

6/

2

-

12/

2

6/

2

0

 

x

B

k

 

x

B

k



Na mocy tabeli w węzłach interpolacji mamy:

k=0

0

1

0

1

y

a

a

4

a

k=1

1

2

1

0

y

a

a

4

a

...............................................

k=i

i

1

i

i

1

i

y

a

a

4

a

.................................................

k=N

N

1

N

N

1

N

y

a

a

4

a

Mamy N+1 równań

i N+3 niewiadomych.

Dwa dodatkowe równania

z warunków brzegowych.

background image

Dla warunku pierwszego rodzaju:

 

 

N

3

0

3

p

b

s

i

p

a

s

mamy na mocy tablicy – pierwsze równanie:

0

1

1

p

a

3

a

3

i ostatnie równanie:

N

1

N

1

N

p

a

3

a

3

Dla warunku drugiego rodzaju:

 

 

N

3

0

3

d

b

s

i

d

a

s



mamy na mocy tablicy – pierwsze równanie:

0

1

2

0

2

1

2

d

a

6

a

12

a

6

background image

i ostatnie równanie:

N

1

N

2

N

2

1

N

2

d

a

6

a

12

a

6

Dla warunku trzeciego rodzaju:

 

 

 

 

 

 

b

s

a

s

b

s

a

s

b

s

a

s

3

3

3

3

3

3





równanie:

0

1

0

1

y

a

a

4

a

i ostatnie:

N

1

N

N

1

N

y

a

a

4

a

mają identyczne prawe strony, gdyż ze względu na okresowość
y

0

=y

N

i dlatego zamiast ostatniego równania piszemy:

1

N

N

1

N

1

0

1

a

a

4

a

a

a

4

a

i pozostałe dwa warunki dają równania:

1

N

1

N

1

1

a

3

a

3

a

3

a

3

background image

1

N

2

N

2

1

N

2

1

2

0

2

1

2

a

6

a

12

a

6

a

6

a

12

a

6

Rozwiązując powyższe trzy równania mamy:

1

N

1

N

0

1

N

1

a

a

a

a

a

a

a więc układ równań przyjmuje postać:

k=0

0

1

0

1

N

y

a

a

4

a

k=1

1

2

1

0

y

a

a

4

a

...............................................

k=i

i

1

i

i

1

i

y

a

a

4

a

.................................................

k=N-1

1

N

0

1

N

2

N

y

a

a

4

a

Jak rozwiązać otrzymany

układ równań metodą

„wymiatania” pokażemy

na przykładzie

background image

Dana jest elipsa o równaniu:

1

y

5

x

2

2

lub w postaci tablicy:

k

0

1

2

3

4

5

x

5

4.9240

4

4.6984

6

3.830

22

2.5

0

y

0

0.1736

5

0.3420

2

0.642

79

0.8660

3

1

k

6

7

8

9

10

x

-2.5

-

3.8302

2

-

4.6984

6

-

4.9240

4

-5

y

0.8660

3

0.6427

9

0.3420

2

0.1736

5

0

background image

k

11

12

13

14

15

x

-

4.9240

4

-

4.69846

-

3.8302

2

-2.5

0

y

-

0.1736

5

-

0.34202

-

0.6427

9

-

0.8660

3

1

k

16

17

18

19

20

x

2.5

3.8302

2

4.6984

6

4.9240

4

5

y

-

0.8660

3

-

0.6427

9

-

0.3420

2

-

0.1736

5

0

Interpolację krzywej zamkniętej możemy wykonać przyjmując
przedstawienie parametryczne:

 

 

21

k

1

k

k

k

t

B

a

t

x

i

 

 

21

k

1

k

k

k

t

B

b

t

y

i mamy układ równań:

background image

83022

.

3

a

a

4

a

69846

.

4

a

a

4

a

92404

.

4

a

a

4

a

5

a

a

4

a

92404

.

4

a

a

4

a

69846

.

4

a

a

4

a

83022

.

3

a

a

4

a

5

.

2

a

a

4

a

0

a

a

4

a

5

.

2

a

a

4

a

83022

.

3

a

a

4

a

69846

.

4

a

a

4

a

92404

.

4

a

a

4

a

5

a

a

4

a

14

13

12

13

12

11

12

11

10

11

10

9

10

9

8

9

8

7

8

7

6

7

6

5

6

5

4

5

4

3

4

3

2

3

2

1

2

1

0

1

0

19

92404

.

4

a

4

a

a

69846

.

4

a

a

4

a

83022

.

3

a

a

4

a

5

.

2

a

a

4

a

0

a

a

4

a

5

.

2

a

a

4

a

19

18

0

19

18

17

18

17

16

17

16

15

16

15

14

15

14

13

Rozwiązanie metodą „wymiatania”:

k

19

k

k

k

1

k

A

a

B

a

L

a

Podstawiając do k-go równania:

k

1

k

k

1

k

x

a

a

4

a

mamy:

k

k

k

k

19

k

k

1

k

k

L

4

A

x

L

4

a

B

L

4

a

a

background image

ale

1

k

19

1

k

1

k

1

k

k

A

a

B

a

L

a

i porównując z

k

k

k

k

19

k

k

1

k

k

L

4

A

x

L

4

a

B

L

4

a

a

mamy wzory rekurencyjne:

k

1

k

L

4

1

L

k

k

1

k

L

4

B

B

i

k

k

k

1

k

L

4

A

x

A

Wartości startowe L

1

, A

1

i B

1

wyznaczamy porównując wzór:

1

19

1

1

1

0

A

a

B

a

L

a

z pierwszym równaniem:

25

.

1

a

25

.

0

a

25

.

0

a

19

1

0

mamy: L

1

=-0.25, B

1

=-0.25 i A

1

=1.25

background image

i z dokładnością do 5 cyfr mamy:

L

2

=-0.26667

L

3

=-0.26786

L

4

=-0.26794

L

5

=L

6

=...=L

19

=-0.26795

B

2

=0.066667 B

15

=-2.44585e-9

B

3

=-0.017857 B

16

=6.55364e-10

B

4

=4.78469e-3 B

17

=-1.75604e-10

B

5

=-1.28205e-3 B

18

=4.7053e-11

B

6

=3.43525e-4 B

19

=-1.26078e-11

B

7

=-9.20471e-5

B

8

=2.4664e-5

B

9

=-6.60869e-6

B

10

=1.77079e-6

B

11

=-4.74482e-7

B

12

=1.27137e-7

B

13

=-3.40663e-8

B

14

=9.128037e-9

W równaniu

92404

.

4

a

4

a

a

19

18

0

B

19

poprawia 4

background image

A

2

=0.97974 A

14

=-0.76330

A

3

=0.99608 A

15

=-0.46535

A

4

=0.75939 A

16

=0.12469

A

5

=0.46640 A

17

=0.63646

A

6

=-0.12497 A

18

=0.85576

A

7

=-0.63639 A

19

=1.02965

A

8

=-0.85579

A

9

=-1.02964

A

10

=-1.04350

A

11

=-1.06014

A

12

=-1.03533

A

13

=-0.98153

Podstawiając do równania:

92404

.

4

a

4

a

a

19

18

0

mamy:

19

19

0

19

19

B

L

4

a

A

92404

.

4

a

czyli

0

19

19

19

a

D

C

a

gdzie

19

19

19

19

B

L

4

A

92404

.

4

C

19

19

19

B

L

4

1

D

background image

i zakładając:

0

1

k

1

k

1

k

a

D

C

a

i podstawiając do związku rekurencyjnego:

1

k

19

1

k

1

k

1

k

k

A

a

B

a

L

a

mamy:

0

k

k

k

a

D

C

a

gdzie

19

1

k

1

k

1

k

k

1

k

19

1

k

1

k

1

k

k

D

B

D

L

D

A

C

B

C

L

C

dla k=18,17,...,0

ze startowymi C

19

, D

19

z zależności:

19

19

19

19

B

L

4

A

92404

.

4

C

i

19

19

19

B

L

4

1

D

background image

C

19

=1.04350 D

19

=-0.26795 C

4

=0.46504 D

4

=3.70097e-4

C

18

=0.75004 D

18

=0.07180 C

3

=0.63978 D

3

=-1.38122e-3

C

17

=0.65479 D

17

=-0.01924 C

2

=0.80608 D

2

=5.15478e-3

C

16

=0.46101 D

16

=5.15478e-3 C

1

=0.83436 D

1

=-0.01924

C

15

=1.16148e-3 D

15

=-1.38122e-3 C

0

=0.78054 D

0

=0.07180

C

14

=-0.46566 D

14

=3.70097e-4

C

13

=-0.63853 D

13

=-9.91696e-5 i z równania:

C

12

=-0.81044 D

12

=2.65816e-5

C

11

=-0.81817 D

11

=-7.15657e-6

C

10

=-0.84091 D

10

=2.04473e-6 mamy:

C

9

=-0.81818 D

9

=-1.02237e-6

C

8

=-0.81042 D

8

=2.04473e-6

C

7

=-0.63861 D

7

=-7.15657e-6 czyli a

0

=a

20

=0.84091

C

6

=-0.46537 D

6

=2.65816e-5

C

5

=8.33928e-5 D

5

=-9.91696e-5

0

0

0

0

a

D

C

a

0

0

0

D

1

C

a

background image

i ze związku rekurencyjnego:

0

k

k

k

a

D

C

a

mamy:

dla k=1,2,...,19

a

1

=a

21

=0.81818 a

13

=-0.63861

a

2

=0.81042 a

14

=-0.46535

a

3

=0.63861 a

15

=0.00000

a

4

=0.46535 a

16

=0.46535

a

5

=-0.00000 a

17

=0.63861

a

6

=-0.46535 a

18

=0.81042

a

7

=-0.63861 a

-1

=a

19

=0.81818

a

8

=-0.81042

a

9

=-0.81818

a

10

=-0.84091

a

11

=-0.81818

a

12

=-0.81042

background image

0 2 4 6 8 10 12 14 16 18 20

5

4

3

2

1

0

1

2

3

4

5

v s

( )

xd s

( )

s

x(s)

 

 

21

k

1

k

k

k

s

B

a

s

x

 

 

10

s

cos

5

s

x

d

x

d

(s)

background image

Równanie krzywej dla współrzędnej x jest:

 

 

10

s

cos

5

s

x

d

dla parametru

20

,

0

s

Błąd między wielomianem interpolacyjnym a funkcją x(s)
definiujemy:

 

 

 

 

 

 

 

 



0

s

x

if

s

x

s

x

0

s

x

if

s

x

s

x

s

x

s

d

d

d

d

d

background image

0 2 4 6 8 10 12 14 16 18 20

0

0.2

0.4

0.6

0.8

 s

( )

s

background image

0 2 4 6 8 10 12 14 16 18 20

1

0.8

0.6

0.4

0.2

0

0.2

0.4

0.6

0.8

1

w s

( )

yd s

( )

s

 

 

21

k

1

k

k

k

s

B

b

s

y

 

 

10

s

sin

s

y

d

y(s)

y

d

(s)

background image

5

4

3

2

1

0

1

2

3

4

5

1

0.8

0.6

0.4

0.2

0

0.2

0.4

0.6

0.8

1

w s

( )

yd s

( )

v s

( ) xd s

( )

x(s),x

d

(s)

y(s)

y

d

(s)

elipsa otrzymana

w rezultacie

interpolacji

i

dokładna

background image

Aproksymacja

Dążenie do minimalizacji normy.

Przykłady stosowanych norm:

 

 

 

 

b

,

a

l

i

przestrzen

w

)

)]

x

(

f

[

(

f

b

,

a

L

i

przestrzen

w

dx

)

x

(

f

)

x

(

w

f

b

,

a

L

i

przestrzen

w

dx

)

x

(

f

f

b

,

a

C

i

przestrzen

w

)

x

(

f

sup

f

N
2

N

0

i

2

1

2

i

w

2,

2

1

2

b

a

w

,

2

2

2

1

2

b

a

2

0

]

b

,

a

[

x









background image

Zadanie aproksymacji polega na minimalizacji normy:

)

x

(

g

)

x

(

F

Niech

K

0

n

n

n

)

x

(

h

a

)

x

(

g

i dane są wartości funkcji :

i

i

F

)

x

(

F

w punktach i=0,2,...P.
Niech będzie zastosowana norma z wagą:





 

P

0

i

2

K

0

n

i

n

n

i

i

)

x

(

h

a

F

w

i szukamy minimum sumy ze względu na współczynniki a

n

.

background image

Przykład:

x

0

0.2

0.5

0.6

1

y

0.1

0.2

0.4

0.45

0.4

2

2

1

0

x

a

x

a

a

)

x

(

g

Przyjmujemy wielomian aproksymujący w postaci:

przyjmując funkcję wagową równą jedności otrzymujemy:

 

2

2

1

0

2

2

1

0

2

2

1

0

2

2

1

0

2

0

2

1

0

a

a

a

4

.

0

36

.

0

a

6

.

0

a

a

45

.

0

25

.

0

a

5

.

0

a

a

4

.

0

04

.

0

a

2

.

0

a

a

2

.

0

)

a

1

.

0

(

)

a

,

a

,

a

(

d

background image

Szukamy ekstremum funkcji d(a

0

,a

1

,a

2

) i przyrównując do zera

pierwsze pochodne względem a

0

, a

1

i a

2

otrzymujemy:

67

.

0

a

1937

.

1

a

349

.

1

a

65

.

1

a

d

91

.

0

a

2365

.

1

a

65

.

1

a

3

.

2

a

d

55

.

1

a

65

.

1

a

3

.

2

a

5

a

d

3

1

0

2

2

1

0

1

2

1

0

0

Rozwiązanie powyższych równań ma postać:

151

.

0

a

493

.

0

a

204

.

0

a

0

1

2

background image

Interpolacja z wagą

 

2

2

1

0

2

2

1

0

2

2

1

0

2

2

1

0

2

0

2

1

0

w

a

a

a

4

.

0

5

.

0

a

36

.

0

a

6

.

0

a

45

.

0

a

25

.

0

a

5

.

0

a

4

.

0

a

04

.

0

a

2

.

0

a

2

.

0

5

.

0

a

1

.

0

25

.

0

a

,

a

,

a

d

Po obliczeniu ekstrmum mamy: a

0

=0.319

a

1

=0.158

a

2

=-0.041

Wielomian aproksymujący jest:

 

2

w

x

041

.

0

x

158

.

0

319

.

0

x

W

background image

x

0

0.2

0.5

0.6

1

y

0.1

0.2

0.4

0.45

0.4

y

ap

0.151 0.241

4

0.346

5

0.373

3

0.44

1

x

0.5

0

0.2

0.1

0.2

0.4

y

y

y

ap

background image

Wielomiany trygonometryczne

n

1

k

k

k

0

n

)

kx

sin

b

kx

cos

a

(

2

a

)

x

(

F

aproksymacja funkcji okresowej na dyskretnym równoodległym
zbiorze punktów:

i=0,1,2,..., 2L-1

L

i

x

i

ciągły:

L

n

)

jx

sin

b

jx

cos

a

(

a

2

1

)

x

(

y

n

1

j

i

i

0

n

a współczynniki wyznacza się z równania:

background image

1

L

2

0

i

2

i

n

i

)

x

(

y

)

x

(

f

Różniczkując względem a

k

otrzymujemy następujące równania

dla wyznaczania współczynników:

0

kx

cos

jx

sin

b

jx

cos

a

2

a

)

x

(

f

i

1

L

2

0

i

n

1

j

i

j

i

j

0

i





0

kx

sin

jx

cos

0

=

k

=

j

dla

2L

0

k

=

j

dla

L

k

j

dla

0

kx

cos

jx

cos

1

L

2

0

l

l

l

1

L

2

0

i

i

i



Mamy tożsamości:

Na mocy powyższych tożsamości mamy:

background image

L

n

,...,

1

,

0

j

L

lj

cos

)

x

(

F

L

1

a

1

L

2

0

l

l

j

Podobnie wyznaczamy współczynniki b

j

z równania:

0

kx

sin

jx

sin

b

jx

cos

a

2

a

)

x

(

F

l

1

L

2

0

l

n

1

j

l

j

l

j

0

l





korzystając z tożsamości:

k

=

j

dla

L

k

j

dla

0

kx

sin

jx

sin

1

L

2

0

l

l

l

otrzymujemy:

L

n

,...,

1

,

0

j

L

jl

sin

)

x

(

F

L

1

b

1

L

2

0

i

l

j

background image

Szereg zespolony.

Dana jest funkcja określona przez podanie jej wartości f

n

w punktach:

N

n

2

x

n

gdzie n=0,1,2,...,N-1. Aproksymujemy funkcję wielomianem
trygonometrycznym postaci:

1

i

gdzie

)

ikx

exp(

c

)

x

(

w

1

N

0

k

k

Otrzymujemy N równań dla wyznaczenia współczynników c

k

:

1

N

0

k

n

k

n

n

)

ikx

exp(

c

f

)

x

(

w

background image

Rozwiązanie ostatniego układu równań czyli współczynniki c

k

są określane równaniami:

1

N

0

n

n

n

p

)

ipx

exp(

f

N

1

c

Idea szybkie transformaty Fouriera tzw. FFT

Fast Fourier Transform

Ponieważ

N

n

2

x

n

więc

N

pn

2

i

exp

ipx

exp

n

Oznaczmy:

N

i

2

epx

w

background image

Zauważmy, że

Re

Im

w=w

N+1

N

2

w

p

=w

p+N

Każda całkowita potęga w leży na okręgu jednostkowym
i co więcej jeżeli wykładnik p potęgi w

p

jest większy od N

to punkty się nakrywają. Na tym spostrzeżeniu bazuje FFT.

background image

Piszemy:

1

N

0

n

n

pn

p

f

w

N

1

c

Możemy zapisać w postaci macierzowej:

Oznaczając:

1

N

,...,

1

,

0

n

dla

N

f

g

n

n

 

 

 

n

pn

p

g

w

c

gdzie

 

2

1

N

1

N

0

1

N

1

0

0

0

0

pn

w

w

w

.

.

.

w

w

w

w

w

w

w

Ponieważ w

0

=1 więc nie będziemy pisać zerowej kolumny i wiersza.

background image

Dalej mamy związki:

N

i

2

epx

w

 

k

k

N

w

k

N

i

2

exp

k

N

i

2

exp

N

N

i

2

exp

k

N

N

i

2

exp

w

 

 





czyli

 

k

k

N

w

w

a więc wiersze i kolumny: 1 i N-1
2 i N-2
..............
k i N-k
..............
N/2-1 i N/2+1 dla N parzystych
(N-1)2 i (N+1)/2 dla N nieparzystych

są sprzężone

.

background image

W praktyce najczęściej stosowane N=2

M

.

Jeżeli liczba węzłów interpolacyjnych mniejsza od 2

M

,

to uzupełniamy zerami.

N=8

w

w

2

w

3

w

4

=-

1

(w

*

)

3

(w

*

)

2

w

*

w

2

w

4

w

6

w

8

=1 (w

*

)

6

(w

*

)

4

(w

*

)

2

w

3

w

6

w

9

w

12

=

-1

(w

*

)

9

(w

*

)

6

(w

*

)

3

w

4

w

8

w

12

w

16

=

1

(w

*

)

12

(w

*

)

8

(w

*

)

4

w

5

w

10

w

15

w

20

=

-1

(w

*

)

15

(w

*

)

10

(w

*

)

5

w

6

w

12

w

18

w

24

=

1

(w

*

)

18

(w

*

)

12

(w

*

)

6

w

7

w

14

w

21

w

28

=

-1

(w

*

)

21

(w

*

)

14

(w

*

)

7

8

i

2

epx

w

 

8

i

2

epx

w

*

background image

lub inaczej

a b c

-

1

c

*

b

*

a

*

b -

1

b

*

1 b -1 b

*

c b

*

a

-

1

a

*

b c

*

-1 1 -1 1 -1 1 -1

c

*

b a

*

-

1

a b

*

c

b

*

-

1

b 1 b

*

-1 b

a

*

b

*

c

*

-

1

c

b a

7

6

5

4

3

2

1

0

f

f

f

f

f

f

f

f

dla otrzymania tablicy mnożników wystarczy obliczyć połowę

pierwszego wiersza!!!

np. a=a

r

+ia

i

oraz a

*

=a

r

-ia

i

czyli np.

af

1

+a

*

f

7

=a

r

(f

1

+f

7

)+ia

i

(f

1

-f

7

)

i podobnie dla innych

operacji.

background image

Wykorzystanie przedstawionych uproszczeń pozwala w stosunku

do zwykłego algorytmu zawierającego N

2

działań zespolonych

zmniejszyć ich liczbę dla N=2

M

do 2NM

1

2

3

4

5

6

0

1200

2400

3600

4800

6000

2

2 m

m2

m 1

m

background image

Rozwiązywanie równań algebraicznych

f(x)=0

Metoda bisekcji

Przykład:

0

1

x

4

x

3

x

-1 0

-0.5

-0.25

-0.125

-

0.1875

f(x) -4 1

-

1.125

-

0.01562

5

0.4980

45

0.2430

8

x

-

0.21875

-

0.23437

5

-

0.242187

5

-

0.24609375

f(x

)

0.11453

2

0.04962

54

0.017044

5

0.00072103

7

background image

x

-

0.248046
9

-

0.2465820
3

-

0.2463379

-

0.246215
8

f(x) -

0.007449
1

-

0.0013209
8

-

0.0002999
3

0.000210

57

Zaleta metody: Jeżeli pierwiastek istnieje, to go znajdziemy.

Wada metody: Duża liczba obliczeń

Regula falsi.

Założenia:

a) funkcja ma w przedziale [a,b] tylko jeden pierwiastek
i zachodzi f(a)f(b)<0,
b) jest funkcja jest klasy C

2

[a,b], pierwsza i druga pochodna

nie zmieniają znaku na przedziale [a,b].

background image

Funkcja spełniająca powyższe założenia musi mieć w otoczeniu

miejsca zerowego jeden z następujących przebiegów:

f(a)

a

b

f(b)

x

y

f(a)

a

b

f(b)

x

y

f(a)

a

b

f(b)

x

y

f(a)

a

b

f(b)

x

y

background image

Przebieg obliczeń metodą regula falsi:

x

y

a

b

f(a)

f(b)

x

1

f(x

1

)

x

2

analitycznie: ustalamy koniec z
warunku f(x

1

)f(a)<0 lub f(x

1

)f(b)<0

Prowadzimy prostą:

   

   

a

f

b

f

a

f

x

f

a

b

a

x

background image

   

   

a

f

b

f

a

f

x

f

a

b

a

x

ale f(x

1

)=0 stąd

 

   

a

f

b

f

a

f

a

b

a

x

1

lub

     

a

f

a

f

b

f

a

b

a

x

1

Dla n-tej iteracji mamy b=x

n-1

i podstawiając mamy:

    

a

f

a

f

x

f

a

x

a

x

1

n

1

n

n

background image

Ocena błędu dla dostatecznie małego przedziału [x

n-1

,x

n

]

można przyjąć jako:

)

x

(

f

)

x

(

f

)

x

(

f

x

x

x

x

n

1

n

n

1

n

n

n

d

Metoda regula falsi jest zbieżna dowolnej funkcji ciągłej
na przedziale [a,b].
Poszukiwanie pierwiastka zostaje zakończone jeżeli:

1

n

n

x

x

Metoda jest wolno zbieżna.

Przykład:

0

1

x

4

x

3

background image

x

-1 0

-0.2

-

0.236641

2

-

0.24423354

f(x) -4 1 0.19

2

0.040134

27

0.00849729

2

     

 

   

2

.

0

x

4

4

1

1

0

1

x

a

f

a

f

b

f

a

b

a

x

1

1

1

Ponieważ f(-1)=-4, a f(x

1

)=0.192,

więc stałym punktem będzie x=-1

x

-

0.24583563

2

-

0.2461749

2

-

0.24624682

9

f(x) 0.00180035

9

0.0003816

03

0.00008089

1

x

-

0.24626207

2

f(x) 0.00001714

7

w metodzie bisekcji potrzebowaliśmy

14 kroków

ocena błędu:

)

x

(

f

)

x

(

f

)

x

(

f

x

x

x

x

n

1

n

n

1

n

n

n

d

background image

)

x

(

f

)

x

(

f

)

x

(

f

x

x

x

x

n

1

n

n

1

n

n

n

d

ocena błędu:

6

d

10

1

.

4

246262072

.

0

x

Metoda siecznych

Przepis:

)

x

(

f

)

x

(

f

)

x

(

f

x

x

x

x

n

1

n

n

1

n

n

n

1

n

Przykład:

0

1

x

4

x

3

x

-1

0

-0.2

-

0.24752475

2

-

0.24625643

9

f(x) -4

1 0.192

-

0.00526448

1

0.00004070

3

w regula falsi potrzeba 8 kroków

background image

x

-

0.2462661

7

f(x)

0.907E-8

w 6-tym kroku

Koniecznie trzeba obliczać f(x

n

) i jeżeli zaczyna narastać

należy zawęzić przedział i powtórzyć obliczenia.

Niebezpieczeństwo znalezienia fałszywego pierwiastka.

Metoda szybsza niż reguła falsi.

a

b

x

1

Pierwsza iteracja musi startować

z punktów spełniających warunek:

f(a)f(b)<0

background image

Metoda Newtona - Raphsona

Niech małe w mamy:

  

 

 

2

x

f

x

f

x

f

x

f

2

Pomijając małe drugiego rzędu

2

mamy, że f(x+)=0,

jeżeli

 

 

x

f

x

f

Graficznie:

x

y

x

n

n

 

n

n

x

f

tg

Równanie prostej stycznej
w punkcie x

n

jest:

  

  

n

n

n

x

f

x

x

x

f

y

x

n+1

background image

Prosta:

  

  

n

n

n

x

f

x

x

x

f

y

przechodzi przez zero, czyli y=0, w punkcie x

n+1

i mamy:

 

 

n

n

n

1

n

x

f

x

f

x

x

Przykład:

0

1

x

4

x

3

 

x

0

-0.25

-
0.24626865

7

-
0.24626617

2

f(x) 1 -

0.01562

5

-
0.00001039

-0.2E-10

W 3 krokach dokładność osiągana w metodzie siecznych

w 5 krokach

background image

W obliczeniach numerycznych pochodną najczęściej oblicza się
numerycznie:

Metoda Newtona – Raphsona jest zbieżna kwadratowo, tzn.

 

 

i

i

2
i

1

i

x

f

2

x

f



  

  

h

x

f

h

x

f

x

f

i

i

i

„Pechowe” przypadki:

x

f(x)

x

0

x

1

x

2

rozbieżna

Zmniejszyć przedział

[x

d

,x

0

]

x

d

background image

cykl

x

f(x)

x

1

=x

3

=...

x

2

=x

4

=...

x

d

Budując procedurę należy się zabezpieczyć przed taką możliwością.

Wystartować z punktu x

1

znajdującego się bliżej x

d

Pierwiastki wielokrotne:

 

 

0

x

f

i

0

x

f

d

d

Przy pierwiastkach wielokrotnych badać funkcję:

)

x

(

f

)

x

(

f

)

x

(

g

background image

Pierwiastki zespolone

Przykład

0

4

x

x

2

3

5

3

1

1

3

5

200

140

80

20

40

100

f x

( )

0

x

Szukamy zespolonych pierwiastków metodą Newtona - Raphsona

background image

n

2
n

2
n

3
n

n

1

n

x

2

x

3

4

x

x

x

x

Jako punkt startowy musimy wybrać liczbę zespoloną:

x

0

=i

gdzie

1

i

i

2309

.

1

8462

.

0

x

e

6056

.

3

e

1623

.

3

i

i

2

3

i

3

i

x

1

69

.

213

i

43

.

198

i

1

x

2

=-0.50605+1.21289i

x

3

=-0.49471+1.32934i

x

4

=-0.50119+1.32409i

x

5

=-0.500059+1.322855i

x

6

=-0.5+1.32275i

błąd=-0.00083198+0.00043738i

x

d

=-0.5+1.322288i

background image

Układy równań nieliniowych

Dany jest układ równań:

0

x

,...,

x

,

x

f

..........

..........

..........

0

x

,...,

x

,

x

f

0

x

,...,

x

,

x

f

n

2

1

n

n

2

1

2

n

2

1

1

Dla skrócenia zapisu wprowadzamy oznaczenia:

 

X

f

x

,...,

x

,

x

f

k

n

2

1

k

oraz

)

x

,...,

x

,

x

(

f

.

.

)

x

,...,

x

,

x

(

f

)

x

,...,

x

,

x

(

f

)

X

(

F

n

2

1

n

n

2

1

2

n

2

1

1

background image

i równanie zapisujemy krótko:

 

0

X

F

Metoda iteracji prostej

Równanie:

 

0

X

F

zapisujemy w postaci:

 

X

G

X

i procedura iteracji prostej ma postać:

1

n

n

X

G

X

Stosowana szczególnie w przypadkach jeżeli mamy dobre

przybliżenie początkowe. Sytuacja taka występuje np. w

przypadku małej zmiany parametrów równania.

background image

Przykład:

1

y

2

x

1

y

x

2

2

2

2

którego rozwiązaniem jest: x

1

=1; y

1

=0 oraz x

2

=-1; y

2

=0

Szukamy rozwiązania układu po małej zmianie parametrów:

1

y

01

.

2

x

95

.

0

y

x

2

2

2

2

mamy schemat iteracyjny:

01

.

2

x

1

y

y

95

.

0

x

2

1

n

n

2

1

n

n

Jako startowy punkt wybieramy: x

0

=1; y

0

=0 i mamy:

background image

n

0

1

2

3

4

x

n

1

0.9746

8

0.97468

0.9618

34

0.96183

4

y

n

0

0

0.15771

8

0.1577

18

0.19300

6

n

5

6

7

8

x

n

0.9553

79

0.95537

9

0.95215

09

0.952151

y

n

0.1930

06

0.20834

7

0.20834

7

0.215573

6

n

9

10

11

12

x

n

0.95054

09

0.950541 0.949738

9

0.949739

y

n

0.21557

34

0.2190799

4

0.219079

7

0.220803

6

background image

Z przedstawionych obliczeń widać, że metoda jest wolno

zbieżna i dlatego stosowana tylko w przypadkach, gdy

znamy bardzo dobrze zerowe przybliżenie. Zastosowanie

w równaniach różniczkowych.

Metoda Newtona - Raphsona

Rozwijamy funkcję f

k

(X) w szereg Taylora w otoczeniu

punktu X

i

:

background image

0

)

x

x

(

x

f

...

)

x

x

(

x

f

)

X

(

f

...

..........

..........

..........

..........

..........

..........

..........

..........

0

)

x

x

(

x

f

...

)

x

x

(

x

f

)

X

(

f

.........

..........

..........

..........

..........

..........

..........

..........

0

)

x

x

(

x

f

...

)

x

x

(

x

f

)

X

(

f

i

,

n

n

X

X

n

n

i

,

1

1

X

X

1

n

i

n

i

,

n

n

X

X

n

k

i

,

1

1

X

X

1

k

i

k

i

,

n

n

X

X

n

1

i

,

1

1

X

X

1

1

i

1

i

i

i

i

i

i

Dla uproszczenia zapisu wprowadzamy macierz Jacobiego
zdefiniowaną następująco:

background image

i

X

X

n

n

2

n

1

n

n

2

2

2

1

2

n

1

2

1

1

1

i

x

f

.

.

x

f

x

f

.

.

.

.

.

.

.

.

.

.

x

f

.

.

x

f

x

f

x

f

.

.

x

f

x

f

)

X

(

J

i w postaci macierzowej możemy krótko zapisać układ równań:

   

0

X

X

J

X

F

i

i

i

gdzie oznaczono:

i

1

i

i

X

X

X

i rozwiązując symbolicznie mamy:

   

i

i

1

i

1

i

X

F

X

J

X

X

background image

Przykład

B x

4

y

4

0.0221347267

1.7154013895

10

3



B x

3

y

3

0.4054281364

0.0757925914

y

3

1.6034305983

x

3

2.547524882

B x

2

y

2

2.0774058055
0.4332616428

y

2

1.2783983143

x

2

2.6115212513

B x

1

y

1

1.0181697533

0.2113862264

y

1

0.6296790316

x

1

0.8775588736

Nowa wartość startowa

x

n

y

n

x

n 1

y

n 1

J x

n 1

y

n 1

1

(

)

B x

n 1

y

n 1



n

1 10





y

0

1



x

0

1



B x y

(

)

f1 x y

(

)

f2 x y

(

)



J x y

(

)

pf1xx y

(

)

pf2xx y

(

)

pf1yx y

(

)

pf2yx y

(

)



pf2yx y

(

)

1

x y

sinx

( ) siny

( )



pf2xx y

(

)

cosx

( ) cos y

( )

1

x y



pf1yx y

(

)

expx y

(

) sin5 x

(

)

siny

( )



pf1xx y

(

)

expx y

(

) sin5 x

(

) 5 cos5 x

(

)

(

)



f2 x y

(

)

cos y

( ) sinx

( )

ln x y



f1 x y

(

)

expx y

(

) sin5 x

(

)

cos y

( )



background image

B x

10

y

10

7.6327832943

10

16

1.0061396161

10

16









y

10

1.5326556001

x

10

2.5104053429

B x

7

y

7

7.6327832943

10

16

1.0061396161

10

16









y

7

1.5326556001

x

7

2.5104053429

B x

6

y

6

1.429500962

10

11

2.5677376891

10

13









y

6

1.5326556001

x

6

2.5104053429

B x

5

y

5

6.4971110534

10

6

3.05616917610

6









y

5

1.5326533005

x

5

2.5104046859

B x

4

y

4

0.0221347267

1.7154013895

10

3

y

4

1.5304688912

x

4

2.5085770863

B x

3

y

3

0.4054281364

0.0757925914

y

3

1.6034305983

x

3

2.547524882

B x

2

y

2

2.0774058055
0.4332616428

y

2

1.2783983143

x

2

2.6115212513

B x

1

y

1

1.0181697533

0.2113862264

y

1

0.6296790316

x

1

0.8775588736

background image

B v

6

w

6

2.8284190988

10

4

4.0146896267

10

6









w

6

1.0538411892

v

6

0.0359317811

B v

5

w

5

3.082730668

10

3

6.399377386710

6









w

5

1.0541200235

v

5

0.0361162013

B v

4

w

4

0.0321949207

9.5125042042

10

4

w

4

1.0561052291

v

4

0.038130869

B v

3

w

3

0.2812551155

0.0206871382

w

3

1.097928869

v

3

0.0523869108

B v

2

w

2

1.4623634971

0.0241516982

w

2

0.8706086893

v

2

0.0653226558

B v

1

w

1

1.251877214

0.5666309304

w

1

0.9185509088

v

1

0.2566102428

v

n

w

n

v

n 1

w

n 1

PJ v

n 1

w

n 1

h

1

( )

B v

n 1

w

n 1



w

0

0.1



v

0

0.1



h 0.1



PJ x y

 h

(

)

f1 x h

 y

(

) f1 x y

(

)

h

f2 x h

 y

(

) f2 x y

(

)

h

f1 x y h

(

) f1 x y

(

)

h

f2 x y h

(

) f2 x y

(

)

h



Obliczenia przy numerycznie liczonej pochodnej

h=0.1

background image

B v

10

w

10

2.0727367989

10

8

1.9211454649

10

10









w

10

1.0538173605

v

10

0.0359127014

B v

7

w

7

2.6276351633

10

5

1.6937433446

10

7









w

7

1.053819747

v

7

0.0359144545

B v

6

w

6

2.8284190988

10

4

4.0146896267

10

6









w

6

1.0538411892

v

6

0.0359317811

B v

5

w

5

3.082730668

10

3

6.3993773867

10

6









w

5

1.0541200235

v

5

0.0361162013

B v

4

w

4

0.0321949207

9.5125042042

10

4

w

4

1.0561052291

v

4

0.038130869

B v

3

w

3

0.2812551155

0.0206871382

w

3

1.097928869

v

3

0.0523869108

B v

2

w

2

1.4623634971

0.0241516982

w

2

0.8706086893

v

2

0.0653226558

B v

1

w

1

1.251877214

0.5666309304

w

1

0.9185509088

v

1

0.2566102428

v

n

w

n

v

n 1

w

n 1

PJ v

n 1

w

n 1

h

1

(

)

B v

n 1

w

n 1




Document Outline


Wyszukiwarka

Podobne podstrony:
Wyklad mn no 8 piątek
Wyklad mn no 7 piątek
Wyklad mn no 4 piątek
Wyklad mn no 5 piątek
Wyklad mn no 6 piątek
Wyklad mn no 8 piątek
Wyklad mn no 1
Wyklad mn 2
Wyklad mn 9
Wyklad mn 16
Wyklad mn 9
Wyklad mn 3
Wyklad mn 6
Wyklad mn 12
Wyklad mn 10
Wyklad mn 6
Wyklad mn 15
Wyklad mn 8
Wyklad mn 5

więcej podobnych podstron