Numerická matematika Banka řešených příkladů

Numerická matematika Banka řešených příkladů Radek Kučera, Pavel Ludvík, Zuzana Morávková Katedra matematiky a deskriptivní geometrie Vysoká škola báňská – Technická Univerzita Ostrava

K D

H

M G

ISBN 978-80-248-3894-6

OBSAH

1

Řešení nelineárních rovnic

5

2

Soustavy lineárních rovnic: přímé metody

19

3

Soustavy lineárních rovnic: iterační metody

30

4

Interpolace a aproximace funkcí

41

5

Numerické integrování a derivování

59

6

Obyčejné diferenciální rovnice: počáteční úlohy

68

Literatura

79

3

Předmluva Studijní materiály tvořící tato skripta jsou určeny převážně pro studenty kombinované i prezenční formy fakulty strojní Vysoké školy báňské – Technické univerzity Ostrava navštěvující předmět Numerická matematika. Naše skripta obsahují řešené příklady a jsou přirozený doplněk skript Radka Kučery Numerické metody, která jsou zaměřená na teoretický výklad základních partií numerické matematiky. Rádi bychom upozornili na webovou stránku http://mdg.vsb.cz/portal/nm, kde je umístěn nejen tento text, ale také řada dalších souvisejících studijních materiálů. Tento studijní text vznikl za finanční podpory projektu IRP-FRVŠ 158/2015 Inovace předmětu Numerická matematika na Fakultě strojní Vysoké školy báňské - Technické univerzitě Ostrava a Katedry matematiky a deskriptivní geometrie VŠB-TUO. Příjemně strávený čas s numerickou matematikou přeje kolektiv autorů.

4

KAPITOLA

1 ŘEŠENÍ NELINEÁRNÍCH ROVNIC

Příklad 1.1: Metodou půlení intervalu určete všechny kořeny rovnice x3 = ln(10 − x ) s přesností ε = 10−3 . f(x)=x3−ln(10−x) 70

60

50

40

30

20

10

0

−10

−20 −2

−1

0

1

2

3

4

Z grafu vidíme, že kořen (průsečík s osou x) leží v intervalu h1, 2i. Vytabelujeme si na tomto 5

KAPITOLA 1. ŘEŠENÍ NELINEÁRNÍCH ROVNIC intervalu hodnotu funkce f ( x ) = x3 − ln(10 − x ) a zjistíme, že znaménko funkčních hodnot se mění mezi 1.2 a 1.3. Funkce je spojitá na daném intervalu. Rovnice má jeden kořen na intervalu h1.2, 1.3i. x

1

1.1

1.2

1.3

1.4

1.5

1.6

1.7

1.8

f (x)

-1.1972

-0.8551

-0.4468

0.0337

0.5922

1.2349

1.9678

2.7967

3.7279

Počítáme kořen na intervalu h1.2, 1.3i, tedy a0 = 1.2, b0 = 1.3. Spočítáme první aproximaci x1 : b0 + a0 1.3 + 1.2 = = 1.25. 2 2

x1 =

Spočítáme funkční hodnoty funkce f ( x ) = x3 − ln(10 − x ) v bodech a0 , x1 , b0 : f ( a0 ) = −0.446,

f ( x1 ) = −0.2159,

f (b0 ) = 0.0337

a určíme interval h a1 , b1 i: f ( a0 ) f ( x1 ) 6< 0 ⇒ a1 = x0 = 1.25, Určíme chybu aproximace b Spočítáme další aproximaci

0 − a0

2

b1 = b0 = 1.3.

= 0.05 6< ε = 10−3 . Ve výpočtu pokračujeme dál.

x2 =

a1 + b1 1.25 + 1.3 = = 1.275 2 2

a určíme interval h a2 , b2 i: f ( a1 ) = −0.2159,

f ( x2 ) = −0.0935,

f ( a1 ) f ( x2 ) 6< 0 ⇒ a2 = x2 = 1.275 Určíme chybu aproximace b Spočítáme další aproximaci

1 − a1

2

x3 =

f (b1 ) = 0.0337, b2 = b1 = 1.3.

= 0.025 6< ε = 10−3 . Ve výpočtu pokračujeme dál. a2 + b2 1.275 + 1.3 = = 1.2875 2 2

a určíme interval h a3 , b3 i: f ( a2 ) = −0.0935,

f ( x3 ) = −0.0305,

f ( a2 ) · f ( x3 ) 6< 0 ⇒ a3 = x3 = 1.2875, Určíme chybu aproximace b Spočítáme další aproximaci

2 − a2

2

x4 =

f (b2 ) = 0.0337, b3 = b2 = 1.3.

= 0.0125 6< ε = 10−3 . Ve výpočtu pokračujeme dál. a3 + b3 1.2875 + 1.3 = = 1.2938 2 2 6

KAPITOLA 1. ŘEŠENÍ NELINEÁRNÍCH ROVNIC a určíme interval h a4 , b4 i: f ( a3 ) = −0.0305,

f ( x3 ) = 0.0014,

f (b3 ) = 0.0337,

f ( a3 ) f ( x3 ) < 0 ⇒ a4 = a3 = 1.2875, 3

b4 = x3 = 1.2938

3

Určíme chybu aproximace b −2 a = 0.0062 6< ε = 10−3 . Ve výpočtu pokračujeme dál. Spočítáme další aproximaci x4 : a4 + b4 1.2875 + 1.2938 = = 1.2906 2 2

x4 = Určíme interval h a5 , b5 i:

f ( a4 ) = −0.0305

f ( x4 ) = −0.0146

f ( a4 ) · f ( x4 ) 6< 0 ⇒ a5 = x4 = 1.2906, 4

f (b4 ) = 0.0014 b5 = b4 = 1.2938

4

Určíme chybu aproximace b −2 a = 0.0031 6< ε = 10−3 . Ve výpočtu pokračujeme dál. Spočítáme další aproximaci x5 : a5 + b5 1.2906 + 1.2938 = = 1.2922 2 2

x5 = Určíme interval h a6 , b6 i:

f ( a5 ) = −0.0146

f ( x5 ) = −0.0066

f ( a5 ) · f ( x5 ) > 0 ⇒ a6 = x5 = 1.2906, 4

f (b5 ) = 0.0014 b6 = b5 = 1.2938

4

Určíme chybu aproximace b −2 a = 0.0016 6< ε = 10−3 . Ve výpočtu pokračujeme dál. Spočítáme další aproximaci x4 : x6 = Určíme chybu aproximace b Vše si zapíšeme do tabulky:

1.2922 + 1.2938 a6 + b6 = = 1.293 2 2

6 − a6

2

= 0.0008 < ε = 10−3 . Je dosaženo zadané přesnosti.

i

ai

f ( ai )

xi

0 1 2 3 4 5 6

1.2 1.25 1.275 1.2875 1.2875 1.2906 1.2922

− − − − − − −

1.25 1.275 1.2875 1.2938 1.2906 1.2922 1.2930

f ( xi )

bi

f ( bi )

|bi − ai |/2

− 1.3 − 1.3 − 1.3 + 1.3 − 1.2938 − 1.2938 − 1.2938

+ + + + + + +

0.05 0.025 0.0125 0.0062 0.0031 0.0016 0.0008

Kořen rovnice je x = 1.293 ± 0.001.

7

KAPITOLA 1. ŘEŠENÍ NELINEÁRNÍCH ROVNIC Příklad 1.2: Metodou půlení intervalu určete všechny kořeny rovnice: x + sin( x ) − 2 = 0 s přesností ε = 10−2 . f(x)=x+sin(x)−2 2

1

0

−1

−2

−3

−4

−5

−6 −3

−2

−1

0

1

2

3

Z grafu vidíme, že kořen (průsečík s osou x) leží v intervalu h1, 2i. Vytabelujeme si na tomto intervalu hodnotu funkce f ( x ) = x + sin( x ) − 2 a zjistíme, že znaménko funkčních hodnot se mění mezi 1.1 a 1.2. Funkce je spojitá na daném intervalu. Rovnice má jeden kořen na intervalu h1.1, 1.2i. x

1

1.1

1.2

1.3

1.4

1.5

1.6

1.7

f (x)

-0.1585

-0.0088

0.1320

0.2636

0.3854

0.4975

0.5996

0.6917

Počítáme kořen na intervalu h1.1, 1.2i a tedy a0 = 1.1, b0 = 1.2,. Spočítáme počáteční aproximaci x0 : x0 =

b0 + a0 1.2 + 1.1 = = 1.15 2 2

Spočítáme funkční hodnoty funkce f ( x ) = x + sin( x ) − 2 v bodech a0 , x0 , b0 . f ( a0 ) = −0.0088

f ( x0 ) = 0.0628

f (b0 ) = 0.1320

f ( a0 ) · f ( x0 ) < 0 ⇒ a1 = a0 = 1.1,

b1 = x0 = 1.15

A určíme interval h a1 , b1 i:

8

KAPITOLA 1. ŘEŠENÍ NELINEÁRNÍCH ROVNIC Spočítáme chybu aproximace b0 −2 a0 = 0.05 6< ε = 10−2 . Ve výpočtu pokračujeme dál. Spočítáme další aproximaci x1 : x1 =

1.1 + 1.15 a1 + b1 = = 1.125 2 2

Určíme interval h a2 , b2 i: f ( a1 ) = −0.0088

f ( x1 ) = 0.0273

f (b1 ) = 0.0628

f ( a1 ) · f ( x1 ) < 0 ⇒ a2 = a1 = 1.1

b2 = x1 = 1.125

Spočítáme chybu aproximace b1 −2 a1 = 0.025 6< ε = 10−2 . Spočítáme další aproximaci x2 : x2 =

1.1 + 1.125 a2 + b2 = = 1.1125 2 2


f ( x2 ) = 0.0093

f ( a2 ) · f ( x2 ) < 0 ⇒ a3 = a2 = 1.1,

f (b2 ) = 0.0273 b3 = x2 = 1.1125

Spočítáme chybu aproximace b2 −2 a2 = 0.0125 6< ε = 10−2 . Spočítáme další aproximaci x3 : x3 =

a3 + b3 1.1 + 1.1125 = = 1.1063 2 2


f ( x3 ) = 0.0003

f ( a3 ) · f ( x3 ) < 0 ⇒ a4 = a3 = 1.1, Spočítáme chybu aproximace Vše si zapíšeme do tabulky: i

ai

0 1 2 3

1.1 1.1 1.1 1.1

f ( ai )

b3 − a3 2

f (b3 ) = 0.0093 b4 = x3 = 1.1063

= 0.0062 < ε = 10−2 . Je dosaženo zadané přesnosti. xi

− 1.15 − 1.125 − 1.1125 − 1.1063

bi

f ( bi )

|bi − ai |/2

+ 1.2 + 1.15 + 1.125 + 1.1125

+ + + +

0.05 0.025 0.0125 0.0062

f ( xi )

Kořen rovnice je: x = 1.11 ± 0.01.

9

KAPITOLA 1. ŘEŠENÍ NELINEÁRNÍCH ROVNIC Příklad 1.3: Newtonovou metodou určete všechny kořeny rovnice: x3 = ln(10 − x ) s přesností ε = 10−6 . f(x)=x3−ln(10−x) 70

60

50

40

30

20

10

0

−10

−20 −2

−1

0

1

2

3

4

Z grafu vidíme, že kořen (průsečík s osou x) leží v intervalu h1, 2i. Vytabelujeme si na tomto intervalu hodnotu funkce f ( x ) = x3 − ln(10 − x ) a zjistíme, že znaménko funkčních hodnot se mění mezi 1.2 a 1.3. Funkce je spojitá na daném intervalu. Rovnice má jeden kořen na intervalu h1.2, 1.3i. x

1

1.1

1.2

1.3

1.4

1.5

1.6

1.7

1.8

f (x)

-1.1972

-0.8551

-0.4468

0.0337

0.5922

1.2349

1.9678

2.7967

3.7279

Spocítáme první a druhou derivaci: f ( x ) = x3 − ln(10 − x ),

f 0 ( x ) = 3 · x2 +

1 , 10 − x

f 00 ( x ) = 6 · x +

Ověříme předpoklady metody: f ( a) −0.4468 f 0 ( a) = 4.4336 = 0.1008 6< 0.1 = b − a

1 , (10 − x )2

a vidíme, že podmínka není splněna. Musíme zmenšit interval, na kterém hledáme kořen. 10

KAPITOLA 1. ŘEŠENÍ NELINEÁRNÍCH ROVNIC 3

f(x)=x −ln(10−x) 0.6

0.4

0.2

0

−0.2

−0.4

−0.6

−0.8

−1 1.1

1.15

1.2

1.25

1.3

1.35

1.4

1.45

Z grafu vidíme, ze kořen leží v intervalu h1.25, 1.3i. A opět začneme ověřovat předpoklady metody, tentokrát na intervalu h1.25, 1.3i: f ( a) −0.2159 f 0 ( a) = 4.8018 = 0.0450 < 0.05 = b − a f (b) 0.0337 f 0 (b) = 5.1849 = 0.0065 < 0.05 = b − a xi

f 0 ( xi )

f 00 ( xi )

1.25 1.26 1.27 1.28 1.29 1.3

4.8018 4.8772 4.9532 5.0299 5.1071 5.1849

7.5131 7.5731 7.6331 7.6932 7.7532 7.8132

Z tabelace první a druhé derivace na intervalu h1.25, 1.3i můžeme usoudit, že: f 0 ( x ) > 0 na h1.25, 1.3i f 00 ( x ) > 0 na h1.25, 1.3i . Jsou tedy splněny předpoklady Newtonovy metody. Zvolíme počáteční aproximaci x0 = b = 1.3. Spočítáme další aproximaci x1 : x1 = x0 −

f ( x0 ) 0.03367697433946 = 1.3 − = 1.29350485098864 f 0 ( x0 ) 5.18494252873563 11

KAPITOLA 1. ŘEŠENÍ NELINEÁRNÍCH ROVNIC a chybu aproximace | x1 − x0 | = 0.00649514901136 6< ε = 10−6 . Ve výpočtu pokračujeme dál. Spočítáme další aproximaci x2 : x2 = x1 −

f ( x1 ) 0.00016453367995 = 1.29350485098864 − = 1.29347280513989 0 f ( x1 ) 5.13432117883453

a chybu aproximace | x2 − x1 | = 0.00003204584875 6< ε = 10−6 . Ve výpočtu pokračujeme dál. Spočítáme další aproximaci x3 : x3 = x2 −

f ( x2 ) 0.00000000399178 = 1.29347280513989 − = 1.29347280436238 f 0 ( x2 ) 5.13407205040059

a chybu aproximace | x1 − x0 | = 0.00000000077751 = 7.7751 · 10−10 < ε = 10−6 . Je dosaženo zadané přesnosti. Vše si zapíšeme do tabulky: i

xi

| x i − x i −1 |

0 1.3 — 1 1.29350485098864 0.00649514901136 2 1.29347280513989 0.00003204584875 3 1.29347280436238 0.00000000077751

Kořen rovnice je:

x = 1.293473 ± 10−6

Příklad 1.4: Newtonovou metodou určete všechny kořeny rovnice: x − 4 cos2 ( x ) = 0 s přesností ε = 10−8 . Nejprve separujeme kořeny rovnice. K tomu si nakreslíme graf funkce na intervalu h0, 5i. Tento interval jsme zvolili proto, protože jsou na něm videt všechny kořeny rovnice.

12

KAPITOLA 1. ŘEŠENÍ NELINEÁRNÍCH ROVNIC f(x)=x−4cos2(x) 5

4

3

2

1

0

−1

−2

−3

−4

0

0.5

1

1.5

2

2.5

3

3.5

4

4.5

5

Z grafu vidíme, že naše rovnice má tři kořeny x1 ∈ h1, 2i , x2 ∈ h2, 3i , x3 ∈ h3, 4i. Všechny tři kořeny aproximujeme Newtonovou metodou tak, jak jsme to už učinili v příkladu 1.3. Nejdříve se tedy budeme zabývat prvním kořenem x1 . Nejdříve ho separujeme na menší interval o délce jední desetiny. Vytabelujeme si na intervalu h1, 2i hodnoty funkce f ( x ) = x − 4 cos2 ( x ) a zjistíme, že náš kořen x1 leží v intervalu h1, 1.1i. Spocítáme první a druhou derivaci: f ( x ) = x − 4 cos2 ( x ), f 0 ( x ) = 1 + 8 sin( x ) cos( x ) = 1 + 4 sin(2x ), f 00 ( x ) = 8 cos(2x )

Ověříme předpoklady metody na intervalu h1, 1.1i: f ( a) 1 − 4 cos2 (1) f 0 ( a) = 1 + 4 sin(2 · 1) = 0.04 < 0.1 = b − a f (b) 1.1 − 4 cos2 (1.1) f 0 (b) = 1 + 4 sin(2 · 1.1) = 0.07 < 0.1 = b − a

13

KAPITOLA 1. ŘEŠENÍ NELINEÁRNÍCH ROVNIC xi

f 0 ( xi )

f 00 ( xi )

1.00 1.01 1.02 1.03 1.04 1.05 1.06 1.07 1.08 1.09 1.10

4.6371 4.6031 4.5677 4.5308 4.4925 4.4528 4.4117 4.36937 4.3255 4.2804 4.2339

-3.3291 -3.4739 -3.6174 -3.7593 -3.8998 -4.0387 -4.1760 -4.3116 -4.4455 -4.5777 -4.7080

Z tabelace první a druhé derivace na intervalu h1, 1.1i můžeme usoudit, že: f 0 ( x ) > 0 na h1, 1.1i f 00 ( x ) < 0 na h1, 1.1i . Jsou tedy splněny předpoklady Newtonovy metody, že první ani druhá derivace nemění na tomto intervalu znaménko. Zvolíme počáteční aproximaci x0 = a = 1. Spočítáme další aproximaci x1 : x1 = x0 −

f ( x0 ) 1 − 4 cos2 (1) = 1 − = 1.0361655092 f 0 ( x0 ) 1 + 4 sin(2 · 1)

a chybu aproximace | x1 − x0 | = 0.0361655092 6< ε = 10−8 . Ve výpočtu pokračujeme dál. Spočítáme další aproximaci x2 : f ( x1 ) 1.0361655092 − 4 cos2 (1.0361655092) = 1.0361655092 − f 0 ( x1 ) 1 + 4 sin(2 · 1.0361655092) = 1.0366737657

x2 = x1 −


x3 = x2 −


x4 = x3 −

a chybu aproximace | x4 − x3 | = 0.00000000000001 < ε = 10−8 . Je dosaženo žádané přesnosti. Vše si zapíšeme do tabulky: 14

KAPITOLA 1. ŘEŠENÍ NELINEÁRNÍCH ROVNIC

Kořen rovnice je:

i

xi

| x i − x i −1 |

0 1 2 3 4

1 1.03616550917501 1.03667376568303 1.03667387601395 1.03667387601396

— 0.03616550917501 0.00050825650802 0.00000011033092 0.00000000000001

x1 = 1.03667388 ± 10−8

Podobně to uděláme se zbývajícími kořeny. Tabelací funkce f ( x ) = x − 4 cos2 ( x ) na intervalech h2, 3i , h3, 4i zjistíme, že kořen x2 ∈ h2.4, 2.5i a kořen x3 ∈ h3.5, 3.6i. Stejný postup jako výše provedeme nejdříve pro kořen x2 . Ověříme předpoklady metody na intervalu h2.4, 2.5i: f ( a) 2.4 − 4 cos2 (2.4) = f 0 ( a) 1 + 4 sin(2 · 2.4) = 0.08 < 0.1 = b − a f (b) 1 − 4 cos2 (1.1) f 0 (b) = 1 + 4 sin(2 · 1.1) = 0.02 < 0.1 = b − a xi

f 0 ( xi )

f 00 ( xi )

2.40 2.41 2.42 2.43 2.44 2.45 2.46 2.47 2.48 2.49 2.50

-2.9846 -2.9768 -2.9674 -2.9565 -2.9439 -2.9298 -2.9141 -2.8968 -2.8780 -2.8576 -2.8356

0.6999 0.8592 1.0181 1.1766 1.3346 1.4920 1.6489 1.8052 1.9607 2.1154 2.2692

Z tabelace první a druhé derivace na intervalu h1, 1.1i můžeme usoudit, že: f 0 ( x ) < 0 na h2.4, 2.5i f 00 ( x ) > 0 na h2.4, 2.5i . Jsou tedy splněny předpoklady Newtonovy metody, že první ani druhá derivace nemění na tomto intervalu znaménko. Zvolíme počáteční aproximaci x0 = a = 2.4. Spočítáme další aproximaci x1 : x1 = x0 −

f ( x0 ) 2.4 − 4 cos2 (2.4) = 2.4 − = 2.47538619175203 f 0 ( x0 ) 1 + 4 sin(2 · 2.4)

a chybu aproximace | x1 − x0 | = 0.0753861918 6< ε = 10−8 . Ve výpočtu pokračujeme dál. 15

KAPITOLA 1. ŘEŠENÍ NELINEÁRNÍCH ROVNIC Spočítáme další aproximaci x2 : 2.4753861918 − 4 cos2 (2.4753861918) f ( x1 ) = 2.4753861918 − f 0 ( x1 ) 1 + 4 sin(2 · 2.4753861918) = 2.47646766323749

x2 = x1 −

a chybu aproximace | x2 − x1 | = 0.0010814715 6< ε = 10−8 . Ve výpočtu pokračujeme dál. Spočítáme další aproximaci x3 : x3

2.4764676632 − 4 cos2 (2.4764676632) f ( x2 ) = 2.4764676632 − = x2 − 0 f ( x2 ) 1 + 4 sin(2 · 2.4764676632) = 2.4764680473

a chybu aproximace | x3 − x2 | = 0.0000003840 6< ε = 10−8 . Ve výpočtu pokračujeme dál. Spočítáme další aproximaci x4 : 2.4764680473 − 4 cos2 (2.4764680473) f ( x3 ) = 2.4764680473 − f 0 ( x3 ) 1 + 4 sin(2 · 2.4764680473) = 2.4764680473

x4 = x3 −

a chybu aproximace | x4 − x3 | = 0.00000000000005 < ε = 10−8 . Je dosaženo žádané přesnosti. Vše si zapíšeme do tabulky:

Kořen rovnice je:

i

xi

| x i − x i −1 |

0 1 2 3 4

2.4 2.47538619175203 2.47646766323749 2.47646804730806 2.47646804730811

— 0.07538619175203 0.00108147148546 0.00000038407057 0.00000000000005

x2 = 2.47646805 ± 10−8

Nyní ještě provedeme stejný postup pro kořen x3 . Ověříme předpoklady metody na intervalu h3.5, 3.6i: f ( a) 3.5 − 4 cos2 (3.5) f 0 ( a) = 1 + 4 sin(2 · 3.5) = 0.002 < 0.1 = b − a f (b) 3.6 − 4 cos2 (3.6) = f 0 (b) 1 + 4 sin(2 · 3.6) = 0.09 < 0.1 = b − a

16

KAPITOLA 1. ŘEŠENÍ NELINEÁRNÍCH ROVNIC xi

f 0 ( xi )

f 00 ( xi )

3.50 3.51 3.52 3.53 3.54 3.55 3.56 3.57 3.58 3.59 3.60

3.6279 3.6877 3.7464 3.8040 3.8605 3.9158 3.9700 4.0230 4.0748 4.1253 4.1746

6.0312 5.9249 5.8162 5.7052 5.5919 5.4763 5.3586 5.2387 5.1168 4.9928 4.8668

Z tabelace první a druhé derivace na intervalu h1, 1.1i můžeme usoudit, že: f 0 ( x ) > 0 na h3.5, 3.6i f 00 ( x ) > 0 na h3.5, 3.6i . Jsou tedy splněny předpoklady Newtonovy metody, že první ani druhá derivace nemění na tomto intervalu znaménko. Zvolíme počáteční aproximaci x0 = a = 3.6. Spočítáme další aproximaci x1 : x1 = x0 −

f ( x0 ) 3.6 − 4 cos2 (3.6) = 2.4 − = 3.5081850213 f 0 ( x0 ) 1 + 4 sin(2 · 3.6)


x2 = x1 −


x3 = x2 −


x4 = x3 −

a chybu aproximace | x4 − x3 | = 0.0000000007 < ε = 10−8 . Je dosaženo žádané přesnosti. Vše si zapíšeme do tabulky: 17

KAPITOLA 1. ŘEŠENÍ NELINEÁRNÍCH ROVNIC i

xi

| x i − x i −1 |

0 1 2 3 4

3.6 3.50818502125718 3.50217696363258 3.50214739193511 3.50214739121355

— 0.09181497874282 0.00600805762461 0.00002957169746 0.00000000072156

Kořen rovnice je:

x3 = 3.50214739 ± 10−8

Kořeny rovnice jsou: x1 = 1.03667388 ± 10−8

x2 = 2.47646805 ± 10−8

18

x3 = 3.50214739 ± 10−8

KAPITOLA

2 SOUSTAVY LINEÁRNÍCH ROVNIC: PŘÍMÉ METODY

Příklad 2.1: Je dána matice 

 −1 −3 2 A =  −6 −19 10  . 3 9 −5

Vypočtěte LU-rozklad A = LU pomocí Gaussovy eliminační metody bez výběru hlavního prvku. Využijte nalezeného LU-rozkladu k určení řešení soustavy lineárních rovnic Ax = b, kde b = (−9, −59, 28)> , inverzní matice A−1 a determinantu det A. Gaussovou eliminační metodou dostáváme:       −6 3 −1 −3 2 −1 −3 2 −1 −3 2 −6 −19 10  ← ∼  0 −1 −2 −1 −2 −+ 0 ∼  0 3 9 −5 ←−−−− + 0 0 1 0 0 1 ← −+

Matice U je výsledkem Gaussovy eliminační metody, matici L sestavíme z multiplikátorů (s opačnými znaménky), které jsme použili během eliminace:     1 0 0 −1 −3 2 L =  6 1 0  , U =  0 −1 −2  . −3 0 1 0 0 1 Při řešení Ax = b nejprve spočteme y z Ly = b a potom z Ux = y vypočítáme x:    y1 = −9  −9 y =  −5  , 6y1 +y2 = −59  −3y1 + y3 = 28 1 19

KAPITOLA 2. SOUSTAVY LINEÁRNÍCH ROVNIC: PŘÍMÉ METODY    − x1 −3x2 +2x3 = −9  2  x= − x2 −2x3 = −5 3 .  x3 = 1 1

Inverzní matici vypočítáme tak, že řešíme soustavy Aai = ei , i = 1, 2, 3, kde ei jsou sloupce jednotkové matice I = (e1 , e2 , e3 ). Získané vektory ai jsou sloupci inverzní matice, tj. A−1 = (a1 , a2 , a3 ). Všechny tři soustavy budeme řešit současně, přičemž jejich řešení opět rozložíme do dvou kroků:    y1 = 1 0 0  1 0 0 Y =  −6 1 0  , 6y1 +y2 = 0 1 0  −3y1 + y3 = 0 0 1 3 0 1    − x1 −3x2 +2x3 = 1 0 0  5 3 8 A−1 =  0 − 1 − 2  . − x2 −2x3 = −6 1 0  3 0 1 x3 = 3 0 1

Determinant det A vypočítáme jako součin diagonálních prvků matic L a U: det A = det L det U = (1 · 1 · 1) · ((−1) · (−1) · 1) = 1.

Příklad 2.2: Je dána matice



 −2 4 1 A =  0 −3 2  . −6 18 0

Vypočtěte LU-rozklad A = LU pomocí Gaussovy eliminační metody bez výběru hlavního prvku. Využijte nalezeného LU-rozkladu k určení řešení soustavy lineárních rovnic Ax = b, kde b = (19, 11, 42)> , inverzní matice A−1 a determinantu det A. Gaussovou eliminační metodou dostáváme:       −2 4 1 −2 4 1 −2 4 1 0 −3  0 −3 2 ← ∼  0 −3 2  −3 2 −+ 2 ∼  0 0 6 −3 ← 0 0 1 −6 18 0 ←−−− + −+

Matice U je výsledkem Gaussovy eliminační metody, matici L sestavíme z multiplikátorů (s opačnými znaménky), které jsme použili během eliminace: 

   1 0 0 −2 4 1 L= 0 1 0  , U =  0 −3 2  . 3 −2 1 0 0 1

Při řešení Ax = b nejprve spočteme y z Ly = b a potom z Ux = y vypočítáme x:    y1 = 19  19 y =  11  , y2 = 11  3y1 −2y2 +y3 = 42 7 20

KAPITOLA 2. SOUSTAVY LINEÁRNÍCH ROVNIC: PŘÍMÉ METODY    −2x1 +4x2 + x3 = 19  −4 x =  1 . −3x2 +2x3 = 11  x3 = 7 7

Inverzní matici vypočítáme tak, že řešíme soustavy Aai = ei , i = 1, 2, 3, kde ei jsou sloupce jednotkové matice I = (e1 , e2 , e3 ). Získané vektory ai jsou sloupci inverzní matice, tj. A−1 = (a1 , a2 , a3 ). Všechny tři soustavy budeme řešit současně, přičemž jejich řešení opět rozložíme do dvou kroků:    y1 = 1 0 0  1 0 0 Y =  0 1 0 , y2 = 0 1 0  3y1 −2y2 +y3 = 0 0 1 −3 2 1    −6 3 11 −2x1 +4x2 + x3 = 1 0 0  6   A−1 =  −2 1 32  . −3x2 +2x3 = 0 1 0  x3 = −3 2 1 −3 2 1

Determinant det A vypočítáme jako součin diagonálních prvků matic L a U: det A = det L det U = (1 · 1 · 1) · ((−2) · (−3) · 1) = 6.


 6 0 3 A =  −18 2 −5  . −4 −1 −6

Vypočtěte LU-rozklad A = LU pomocí Gaussovy eliminační metody bez výběru hlavního prvku. Využijte nalezeného LU-rozkladu k určení řešení soustavy lineárních rovnic Ax = b, kde b = (0, 8, −8)> , inverzní matice A−1 a determinantu det A. Gaussovou eliminační metodou dostáváme: 

     4 2 6 0 3 6 0 3 6 0 3 3 6=3 1 −18 2 −5 ← ∼ 0 2 ∼ 0 2 4  4  −+ 2 0 −1 −4 ← −4 −1 −6 ←−−− + 0 0 −2 −+

Matice U je výsledkem Gaussovy eliminační metody, matici L sestavíme z multiplikátorů (s opačnými znaménky), které jsme použili během eliminace:     1 0 0 6 0 3 1 0 , U =  0 2 L =  −3 4 . 1 2 0 0 −2 −3 −2 1

Nezapomínejte, že při výpočtu matice L Gaussovou eliminační metodou nesmíte měnit pořadí řádků, ani násobit řádky konstantou. Není proto možné vyhnout se práci se zlomky. V dalších úlohách se seznámite s LU-rozkladem s permutační maticí, který můžete použít v případech, kdy je prohození řádků nutné. 21

KAPITOLA 2. SOUSTAVY LINEÁRNÍCH ROVNIC: PŘÍMÉ METODY Při řešení Ax = b nejprve spočteme y z Ly = b a potom z Ux = y vypočítáme x:    y1 = 0  0 −3y1 +y2 = −8 y =  −8  ,  2 1 4 − 3 y1 − 2 y2 + y3 = 8    6x1 +3x3 = 0  1 x =  0 . +2x2 +4x3 = −8  −2x3 = 4 −2

Inverzní matici vypočítáme tak, že řešíme soustavy Aai = ei , i = 1, 2, 3, kde ei jsou sloupce jednotkové matice I = (e1 , e2 , e3 ). Získané vektory ai jsou sloupci inverzní matice, tj. A−1 = (a1 , a2 , a3 ). Všechny tři soustavy budeme řešit současně, přičemž jejich řešení opět rozložíme do dvou kroků:    y1 = 1 0 0  1 0 0 −3y1 +y2 = 0 1 0 Y =  3 1 0 ,  1 2 13 1 − 3 y1 − 2 y2 + y3 = 0 0 1 6 2 1  17   1 1 6x1 +3x3 = 1 0 0  24 8 4   A−1 =  11 +2x2 +4x3 = −6 1 0 1 1 . 3  1 1 −2x3 = 3 0 1 − 13 12 − 4 − 2 Determinant det A vypočítáme jako součin diagonálních prvků matic L a U: det A = det L det U = (1 · 1 · 1) · (6 · 2 · (−2)) = −24.


 −2 −4 2 0  −2 −6 6 −2  . A=  0 1 3 6  −1 0 0 2 Vypočtěte LU-rozklad A = LU pomocí Gaussovy eliminační metody bez výběru hlavního prvku. Využijte nalezeného LU-rozkladu k určení řešení soustavy lineárních rovnic Ax = b, kde b = (−4, −22, 29, 13)> , inverzní matice A−1 a determinantu det A. Práce s čtvercovou maticí se čtyřmi řádky je o něco náročnější než v případě matic s menšími rozměry, postup je nicméně totožný. Gaussovou eliminační metodou dostáváme:     −1 0 − 12 −2 −4 2 0 −2 −4 2 0 1  −2 −6 6 −2 ←  0 −2 4 −2 1 2     −+ ∼   0   1 −7 −4 ←−−−− + 0 1 −7 −4 ← −+ −1 0 0 2 0 2 −1 2 ←−−−−−− + ←−−− +     −2 −4 2 0 −2 −4 2 0  0 −2 4 −2  0 −2 4 −2   ∼ ∼ 3  0   0 −5 0 5 5 0 5 5  0 0 3 0 0 0 0 −3 ← −+ 22

KAPITOLA 2. SOUSTAVY LINEÁRNÍCH ROVNIC: PŘÍMÉ METODY Matice U je výsledkem Gaussovy eliminační metody, matici L sestavíme z multiplikátorů (s opačnými znaménky), které jsme použili během eliminace:     1 0 0 0 −2 −4 2 0  1   1 0 0   , U =  0 −2 4 −2  . L= 1  0 − 1 0   0 0 5 5  1 2

2

−1

3 5

0

1

0 0 −3

Opět mějte stále na paměti, že tento algoritmus hledání LU-rozkladu zapovídá prohazování řádků a násobení řádků konstantou. Při řešení Ax = b nejprve spočteme y z Ly = b a potom z Ux = y vypočítáme x:    y1 = −4  −4   −18  y1 + y2 = −22   y= 1  20  , = 29  − 2 y2 + y3   1 −15 −y2 + 35 y3 +y4 = 13 2 y1    −2x1 −4x2 +2x3 = −4  −3   2  −2x2 +4x3 −2x4 = −18   x=  −1  . +5x3 +5x4 = 20    −3x4 = −15 5

Inverzní matici vypočítáme tak, že řešíme soustavy Aai = ei , i = 1, 2, 3, 4, kde ei jsou sloupce jednotkové matice I = (e1 , e2 , e3 , e4 ). Získané vektory ai jsou sloupci inverzní matice, tj. A−1 = (a1 , a2 , a3 , a4 ). Všechny čtyři soustavy budeme řešit současně, přičemž jejich řešení opět rozložíme do dvou kroků:    1 0 0 0 y1 = 1 0 0 0      −1 1 0 0  y1 + y2 = 0 1 0 0  , Y=  −1 1  = 0 0 1 0  − 12 y2 +y3 1 0    2 2  3 1 6 3 7 y − y + y + y = 0 0 0 1 2 4 2 1 5 3 − 5 10 − 5 1    7 −8 − 17 − 27 0 0  −2x1 −4x2 +2x3 = 1 0  2    13 5   −2x2 +4x3 −2x4 = −1 1 0 0  −1  − 83 13 3 3 3  A =  . 1 1 1 2 2 1 +5x3 +5x4 = − 2 2 1 0  − − −   3 3 3 3     7 6 3 3 −3x4 = − 5 10 − 5 1 0 −2 − 2 − 12

Determinant det A vypočítáme jako součin diagonálních prvků matic L a U:

det A = det L det U = (1 · 1 · 1 · 1) · ((−2) · (−2) · 5 · (−3)) = −60. Příklad 2.5: Je dána matice 

 4 −3 0 3  −8 5 0 −4  . A=  12 −7 −1 8  8 −3 −1 1 Vypočtěte LU-rozklad A = LU pomocí Gaussovy eliminační metody bez výběru hlavního prvku. Využijte nalezeného LU-rozkladu k určení řešení soustavy lineárních rovnic

23

KAPITOLA 2. SOUSTAVY LINEÁRNÍCH ROVNIC: PŘÍMÉ METODY Ax = b, kde b = (−16, 31, −52, −33)> , inverzní matice A−1 a determinantu det A. Gaussovou eliminační metodou dostáváme:   4 −3 0 3 2 −3 −2  −8 5  0 −4 ← −+   12 −7 −1 8  ←−−− + 8 −3 −1 1 ←−−−−−− +    4 −3 0 3 4 0 −1 0 2 0  ∼ ∼ 0 0 −1 3 0 −1 0 0 −1 1 ← 0 −+



 4 −3 0 3 0 −1 0 2  2 3  ∼ 0 2 −1 −1 ← −+ 0 3 −1 −5 ←−−− +  −3 0 3 −1 0 2   0 −1 3  0 0 −2

Matice U je výsledkem Gaussovy eliminační metody, matici L sestavíme z multiplikátorů (s opačnými znaménky), které jsme použili během eliminace:     1 0 0 0 4 −3 0 3  −2  1 0 0  0 2   , U =  0 −1 . L=  3 −2 1 0   0 0 −1 3  2 −3 1 1 0 0 0 −2 Při řešení Ax = b nejprve spočteme y z Ly = b a potom z Ux = y vypočítáme x:    y1 = −16  − 16    −1  −2y1 +y2 = 31  y=  −6  , 3y1 −2y2 +y3 = −52    2y1 −3y2 +y3 +y4 = −33 2    4x1 −3x2 +3x4 = −16  −4   −1  − x2 +2x4 = −1   x=  3 . − x3 +3x4 = −6    −2x4 = 2 −1

Inverzní matici vypočítáme tak, že řešíme soustavy Aai = ei , i = 1, 2, 3, 4, kde ei jsou sloupce jednotkové matice I = (e1 , e2 , e3 , e4 ). Získané vektory ai jsou sloupci inverzní matice, tj. A−1 = (a1 , a2 , a3 , a4 ). Všechny čtyři soustavy budeme řešit současně, přičemž jejich řešení opět rozložíme do dvou kroků:    y1 = 1 0 0 0  1 0 0 0   2 1 −2y1 +y2 = 0 1 0 0  0 0  , Y=   3y1 −2y2 +y3 = 0 0 1 0  1 2 1 0   2y1 −3y2 +y3 +y4 = 0 0 0 1 3 1 −1 1   19 9 3 3  − − − 8 8 8 8 4x1 −3x2 +3x4 = 1 0 0 0        − 5 − 2 1 − 1 − x2 +2x4 = 2 1 0 0  A−1 =   − 11 − 7 1 − 3  . − x3 +3x4 = 1 2 1 0    2 2 2 2   −2x4 = 3 1 −1 1 3 1 1 − 2 − 2 2 − 12

Determinant det A vypočítáme jako součin diagonálních prvků matic L a U:

det A = det L det U = (1 · 1 · 1 · 1) · (4 · (−1) · (−1) · (−2)) = −8.

24

KAPITOLA 2. SOUSTAVY LINEÁRNÍCH ROVNIC: PŘÍMÉ METODY Příklad 2.6: Je dána matice 

−4 3  −12 10 A=  −12 11 8 −3

 1 −3 4 −11  . 6 −16  0 0

Vypočtěte LU-rozklad A = LU pomocí Gaussovy eliminační metody bez výběru hlavního prvku. Využijte nalezeného LU-rozkladu k určení řešení soustavy lineárních rovnic Ax = b, kde b = (5, 11, 1, −19)> , inverzní matice A−1 a determinantu det A. Gaussovou eliminační metodou dostáváme:    −3 −3 2 −4 3 1 −3 −4 3 1 −12 10 4 −11 ←  0 1 1   −+ ∼ −12 11 6 −16 ←−−−− +  0 2 3 8 −3 0 0 0 3 2 ←−−−−−−− +     −4 3 1 −3 −4 3 1 −3  0 1 1 −2  0 1 1 −2    ∼ ∼  0 0 1 −3  0 0 1 −3 1 0 0 −1 0 0 0 0 −3 ← −+

 −3 −2 −3 −2   −7 ← −+ −6 ←−−−− +

Matice U je výsledkem Gaussovy eliminační metody, matici L sestavíme z multiplikátorů (s opačnými znaménky), které jsme použili během eliminace:     1 0 0 0 −4 3 1 −3  3 1   0 0   , U =  0 1 1 −2  . L=  3 2  0 0 1 −3  1 0 

−2 3 −1 1

0 0 0 −3

Při řešení Ax = b nejprve spočteme y z Ly = b a potom z Ux = y vypočítáme x:    y1 = 5  5    −4  3y1 +y2 = 11  y=  −6  , 3y1 +2y2 +y3 = 1    −2y1 +3y2 −y3 +y4 = −19 −3    −4x1 +3x2 + x3 −3x4 = 5  −2    1  x2 + x3 −2x4 = −4  x=  −3  . x3 −3x4 = −6    −3x4 = −3 1

Inverzní matici vypočítáme tak, že řešíme soustavy Aai = ei , i = 1, 2, 3, 4, kde ei jsou sloupce jednotkové matice I = (e1 , e2 , e3 , e4 ). Získané vektory ai jsou sloupci inverzní matice, tj. A−1 = (a1 , a2 , a3 , a4 ). Všechny čtyři soustavy budeme řešit současně, přičemž jejich řešení opět rozložíme do dvou kroků:    y1 = 1 0 0 0  1 0 0 0   −3 3y1 +y2 = 0 1 0 0  1 0 0   Y=  3 −2 1 0  , 3y1 +2y2 +y3 = 0 0 1 0    −2y1 +3y2 −y3 +y4 = 0 0 0 1 14 −5 1 1 25

KAPITOLA 2. SOUSTAVY LINEÁRNÍCH ROVNIC: PŘÍMÉ METODY

−4x1 +3x2 + x3 −3x4 x2 + x3 −2x4 x3 −3x4 −3x4

= 1 0 1 = −3 = 3 −2 = 14 −5

0 0 1 1

   

0 0 0 1

  

A−1



1  − 21 12 − 14 4 1   −4 4 −2 3 3 3 3 . =  −11 3 0 −1  5 1 1 − 14 3 3 −3 −3

Determinant det A vypočítáme jako součin diagonálních prvků matic L a U: det A = det L det U = (1 · 1 · 1 · 1) · ((−4) · 1 · 1 · (−3)) = 12.

Příklad 2.7: Je dána matice



 1 −1 7 A =  −2 2 8 . −1 5 3

Vypočtěte LU-rozklad s permutační maticí, tj. rozklad PA = LU, pomocí Gaussovy eliminační metody s výběrem hlavního prvku. Při výpočtu LU-rozkladu s permutační maticí budeme postupovat takto: ˜ = A, P˜ = I a L˜ = I (= jednotková matice); • vytvoříme pomocné matice U ˜ provádíme dopředný chod Gaussovy eliminační metody s výběrem hlav• v matici U ního prvku; ˜ • v matici P˜ přehazujeme řádky stejně jako v matici U; • do matice L˜ zapíšeme v každé fázi multiplikátory (s opačnými znaménky) a při pře˜ přehodíme v L˜ řádky i sloupce; hození řádků v U ˜ L = L˜ a U = U. ˜ • nakonec dostáváme P = P, ˜ P˜ a L. ˜ Výpočet tedy zahájíme tím, že definujeme pomocné matice U, y       1 −1 7 ← 1 0 0 ← 1 − − ˜ ˜ ˜      U = −2 2 8 ← −, P= 0 1 0 ← −, L= 0 −1 5 3 0 0 1 0

y 

 0 0 ← −  1 0 ← −. 0 1

Po výběru hlavního prvku v první fázi (tj. v absolutní hodnoty největšího prvku v prvním sloupci): 

     1 − 12 0 1 0 1 0 0 −2 2 8 2 ˜ =  1 −1 7 ← U , P˜ =  1 0 0  , L˜ =  0 1 0  . −+ 0 0 1 0 0 1 −1 5 3 ←−−− + Po eliminaci v první fázi s multiplikátory m21 = 





1 2

a m31 = − 12 :





1 −2 2 8 0 1 0 ˜ =  0 0 11  ← U − , P˜ = 1 0 0 ← − , L˜ = − 12 1 0 4 −1 ← 0 0 1 ← − − 2

26

y  y   0 0 1 0 ← −. 0 1 ← −

KAPITOLA 2. SOUSTAVY LINEÁRNÍCH ROVNIC: PŘÍMÉ METODY Po výběru hlavního prvku ve druhé fázi (bez prohození druhého a třetího řádku bychom v eliminaci vůbec nemohli pokračovat!): 

     1 0 0 −2 2 8 0 1 0 ˜ =  0 4 −1 ˜ = 0 0 1 , L˜ =  1 1 0 . U 0 , P 2 − 12 0 1 0 0 11 1 0 0 ← −+ Eliminace ve druhé fázi je multiplikátorem m32 = 0, čili veskrze symbolická:      1 0 0 −2 2 8 0 1 0 ˜ =  0 4 −1 , P˜ = 0 0 1 , L˜ =  1 1 0 . U 2 − 12 0 1 0 0 11 1 0 0 

Výsledkem jsou matice ˜ U = U. ˜ ˜ L = L, P = P, Můžete si vyzkoušet, že skutečně platí PA = LU. Poznámka. Mohli byste považovat za nevýhodu, že získaný tvar ve skutečnosti není rozkladem matice A. Ten byste dostali jeho vynásobením obou stran rovnosti zleva maticí P−1 . Výpočet inverzní matice k permutační matici je naštěstí velice snadný, platí totiž P−1 = P> . Matici A proto můžeme rozložit na součin takto: A = P−1 LU = P> LU, neboli 

       1 0 0 1 −1 7 0 0 1 −2 2 8  −2 2 8 = 1 0 0 ·  1 1 0 ·  0 4 −1 . 2 − 12 0 1 −1 5 3 0 1 0 0 0 11 Příklad 2.8: Je dána matice  2 −1 −8 −2  −1 3 −9 6   . A=  −3 −1 6 3  1 −6 3 3 

Vypočtěte LU-rozklad s permutační maticí, tj. rozklad PA = LU, pomocí Gaussovy eliminační metody s výběrem hlavního prvku. Při výpočtu LU-rozkladu s permutační maticí budeme postupovat takto: ˜ = A, P˜ = I a L˜ = I (= jednotková matice); • vytvoříme pomocné matice U ˜ provádíme dopředný chod Gaussovy eliminační metody s výběrem hlav• v matici U ního prvku; ˜ • v matici P˜ přehazujeme řádky stejně jako v matici U; • do matice L˜ zapíšeme v každé fázi multiplikátory (s opačnými znaménky) a při pře˜ přehodíme v L˜ řádky i sloupce; hození řádků v U ˜ L = L˜ a U = U. ˜ • nakonec dostáváme P = P,

27

KAPITOLA 2. SOUSTAVY LINEÁRNÍCH ROVNIC: PŘÍMÉ METODY ˜ P˜ a L. ˜ Výpočet tedy zahájíme tím, že definujeme pomocné matice U, 





1 0 2 −1 −8 −2 ← − 0 1   ˜ =  −1 3 −9 6  , P˜ =  U 0 0  −3 −1 6  3 ← − 0 0 1 −6 3 3

y   0 0 ← 1 − 0 0 0  , L˜ =   0 1 0 ← − 0 1 0 

 y

0 1 0 0

 0 ← − 0  . 0 ← − 1

0 0 1 0

Po výběru hlavního prvku v první fázi (tj. v absolutní hodnoty největšího prvku v prvním sloupci):   − 13 23 13 0 −3 −1 6 3   −1 3 −9 6  ← +  − ˜ = ˜ = 0 , P U 1  2 −1 −8 −2 ←−−−− + 0 1 −6 3 3 ←−−−−−− + 

0 1 0 0

1 0 0 0

Po eliminaci v první fázi s multiplikátory m21 = − 13 , m31 = 

−3 −1 6 10  0 −11 3 ˜ = U  0 − 5 −4  3 0 − 19 5 3





3 0 0  5 ← −  ˜ = 0 1 , P  1 0 0 0 0 4 ← −



 0 0  , L˜ = 0 1

2 3

a m41 = 13 :



1 0  0 0



1 0  1 0 −  ← ˜ =  32 , L − 0 3 − 13 1 ← −

1 0 0 0

y  0 1 0 0

0 1 0 0

0 0 1 0

Po výběru hlavního prvku ve druhé fázi (toho docílíme prohozením druhého řádku):      −3 −1 6 3 1 0 0 0 1 0 1 5 10  0 − 19   − 19 19 5 4 0 0 0 1 3  − 3 1   ˜ ˜ ˜ U= , P= , L= 2  5 −+ 1 0 0 0  0 − 3 −4 0 ← − 3 0 10 1 0 1 0 0 0 −11 5 ←−−−− + 0 3 3 5 Při eliminaci ve druhé fázi jsme použili multiplikátory m32 = − 19 a m32 =

  −3 −1 6 3 0   0 − 19 5 4   3 ˜ = ˜ = 0 U , P  0 101 20  1 0 − 19 − 19  ← −  135 0 0 0 − 159 ← − 19 19 

0 0 0 1

 0 1  , L˜ = 0 ← − 0 ← −

1 0 0 0



1 − 1  3  2 − 3 1 3

10 19 :

0 1 5 19 − 10 19

0 0 1 0

 0 0  . 0 1

 y  0 0 −  ← .  0 1 ← − a čtvrtého  0 0 0 0   . 1 0 0 1 y  y   0 0 0 0  .  1 0 ← − 0 1 ← −

Po výběru hlavního prvku ve třetím kroku (prohozením třetího a čtvrtého řádku): 

−3 −1 6  0 − 19 5 3 ˜ = U  0 − 159  0 19 0

0

 3 4   135  19 

− 101 − 20 19 19

− 101 159

← −+



0 0 , P˜ =  0 1

0 0 1 0

1 0 0 0

Nakonec provedeme eliminaci koeficientem m43 = − 101 159 : 28

  1 0 0 1  −3 1 1  , L˜ =   10 0  13 − 19 5 0 − 23 19

 0 0 0 0   . 1 0 0 1

KAPITOLA 2. SOUSTAVY LINEÁRNÍCH ROVNIC: PŘÍMÉ METODY



−3 −1 6 19  0 − 5 3 ˜ = U  0 0 − 159  19 0

0

0

 3 4   ˜ = 135  , P 19  − 259 53



0 0  0 1

0 0 1 0

1 0 0 0

 0 1  , L˜ = 0 0

Výsledkem jsou matice ˜ U = U. ˜ ˜ L = L, P = P, Můžete si vyzkoušet, že skutečně platí PA = LU.

29



1 0 − 1 1  3  1 − 10  3 19 2 5 − 3 19

0 0 1 101 159

 0 0  . 0  1

KAPITOLA

3 SOUSTAVY LINEÁRNÍCH ROVNIC: ITERAČNÍ METODY

Příklad 3.1: Vyřešte soustavu lineárních rovnic

− x1 + x2 +3x3 = 0, − x1 +2x2 − x3 = 9, 3x1 + x2 + x3 = 10, pomocí Jacobiho iterační metody s přesností e = 10−2 . Než aplikujeme rekurentní vzorce pro Jacobiho metodu, musíme zajistit, aby posloupnost vektorů, kterou metoda vytváří, konvergovala ke správnému výsledku. Z teorie víme, že konvergence je zaručena, pokud je matice soustavy ostře diagonálně dominantní (v ostatních případech konvergovat může, ale také nemusí). Pro tento účel stačí přehodit první a třetí rovnici a následně přičíst první rovnici ke druhé rovnici. Dostaneme pak: 3x1 + x2 + x3 = 10, 2x1 +3x2 = 19, − x1 + x2 +3x3 = 0. Matice soustavy



 3 1 1  2 3 0 −1 1 3

je ostře diagonálně dominantní, neboť absolutní hodnoty prvků na diagonále jsou (ostře) větší než součty absolutních hodnot ostatních prvků na příslušných řádcích, tj. 3 > 1 + 1, 3 > 2 + 0 a 3 > 1 + 1. 30

KAPITOLA 3. SOUSTAVY LINEÁRNÍCH ROVNIC: ITERAČNÍ METODY Upravenou soustavu převedeme na iterační tvar x1 = x2 = x3 =

1 3 (10 − x2 − x3 ), 1 3 (19 − 2x1 ), 1 3 ( x1 − x2 )

a připsáním iteračních indexů dostaneme rekurentní vzorce pro Jacobiho metodu ( k +1)

= = =

x1

( k +1)

x2

( k +1)

x3

(k) (k) 1 3 (10 − x2 − x3 ), (k) 1 3 (19 − 2x1 ), (k) 1 (k) 3 ( x1 − x2 ).

Poznamenejme, že rekurentní vzorce můžeme zapsat také v maticovém tvaru, který je obzláště výhodný při řešení úlohy pomocí počítačového softwaru.  ( k +1)     (k)   10  x1 x1 0 − 31 − 13 3  ( k +1)   2   19  (k)   · = + 0 0 −  x2   x2   3  . 3 1 1 ( k +1) (k) −3 0 0 x3 x3 3 (0)

(0)

(0)

Protože výpočet konverguje pro libovolnou počáteční aproximaci x(0) = ( x1 , x2 , x3 )> , zvolíme pro jednoduchost x(0) = ( 1, 1, 1 )> . Dosadíme do pravé strany rekurentních vzorců a vypočteme x(1) = ( 2.6667, 5.6667, 0 )> . Jestliže takto pokražujeme dále, dostáváme x(2) = ( 1.4444, 4.5556, −1 )> , x(3) = ( 2.1481, 5.3704, −1.0370 )> , atd. Současně přitom zjišťujeme, zda platí

kx(k) − x(k−1) k R ≤ 10−3 , takže postupně počítáme

kx(1) − x(0) k R = max{|2.6667 − 1|, |5.6667 − 1|, |0 − 1|} = 4.6667, kx(2) − x(1) k R = max{|1.4444 − 2.6667|, |4.5556 − 5.6667|, | − 1 − 0|} = 1.2222, atd. Výpočet zapíšeme do tabulky: k 0 1 2 3 4 5 6 7 8 9 10

(k)

x1

1 2.6667 1.4444 2.1481 1.8889 2.0576 1.9767 2.0146 1.9931 2.0034 1.9982

(k)

x2

1 5.6667 4.5556 5.3704 4.9012 5.0741 4.9616 5.0155 4.9902 5.0046 4.9978

(k)

x3

1 0 −1.0000 −1.0370 −1.0741 −1.0041 −1.0055 −0.9950 −1.0003 −0.9990 −1.0004 31

k x( k ) − x( k −1) k R — 4.6667 1.2222 0.8148 0.4691 0.1728 0.1125 0.0540 0.0253 0.0143 0.0068

KAPITOLA 3. SOUSTAVY LINEÁRNÍCH ROVNIC: ITERAČNÍ METODY Ukončili jsem po jedenácté iteraci, protože kx(10) − x(9) k R = 0.0068 ≤ 10−2 . Požadovaná přesnost je ε = 10−2 , proto jednotlivé souřadnice výsledku zaokrouhlíme na 2 desetinná místa. Výsledek zapíšeme jako x1 = 2.00 ± 10−2 , x2 = 5.00 ± 10−2 , x3 = −1.00 ± 10−2 .

Příklad 3.2: Vyřešte soustavu lineárních rovnic x1 + x2 +5x3 = −6, 5x1 + x2 + x3 = −18, 7x1 +7x2 +7x3 = −14, pomocí Gaussovy-Seidelovy iterační metody s přesností ε = 10−2 . Z teorie víme, že konvergence Gaussovy-Seidelovy metody je zaručena, pokud je matice soustavy ostře diagonálně dominantní (v ostatních případech konvergovat může, ale také nemusí). Pro dosažení takové matice můžeme například nejprve umístit druhou rovnici na první pozici, třetí rovnici na druhou pozici a první rovnici na třetí pozici: 5x1 + x2 + x3 = −18, 7x1 +7x2 +7x3 = −14, x1 + x2 +5x3 = −6, Pak již jen stačí odečíst první a třetí rovnici od druhé rovnice: 5x1 + x2 + x3 = −18, x1 +5x2 + x3 = 10, x1 + x2 +5x3 = −6. Matice soustavy



 5 1 1 2 5 1 1 1 5

je ostře diagonálně dominantní, neboť absolutní hodnoty prvků na diagonále jsou (ostře) větší než součty absolutních hodnot ostatních prvků na příslušných řádcích, tj. 5 > 1 + 1, 5 > 2 + 1 a 5 > 1 + 1. Upravenou soustavu převedeme na iterační tvar x1 = x2 = x3 =

1 5 (−18 − x2 − x3 ), 1 5 (10 − x1 − x3 ), 1 5 (−6 − x1 − x2 )

a připsáním iteračních indexů dostaneme rekurentní vzorce pro Gaussovu-Seidelovu metodu ( k +1) (k) (k) x1 = 51 (−18 − x2 − x3 ), ( k +1)

x2

( k +1)

x3

( k +1)

(k)

= 15 (10 − x1 − x3 ), ( k +1) ( k +1) 1 = − 5 (−6 − x1 − x2 ). 32

KAPITOLA 3. SOUSTAVY LINEÁRNÍCH ROVNIC: ITERAČNÍ METODY (0)

(0)

(0)

Protože výpočet konverguje pro libovolnou počáteční aproximaci x(0) = ( x1 , x2 , x3 )> , zvolíme pro jednoduchost x(0) = ( 1, 1, 1 )> . Dosadíme do pravé strany rekurentních vzorců a vypočteme x(1) = ( −4, 2.6, −0.92 )> . Jestliže takto pokražujeme dále, dostáváme x(2) = ( −3.936, 2.9712, −1.007 )> , x(3) = ( −3.9928, 3, −1.0014 )> , atd. Současně přitom zjišťujeme, zda platí


kx(1) − x(0) k R = max{| − 4 − 1|, |2.6 − 1|, | − 0.92 − 1|} = 5, kx(2) − x(1) k R = max{| − 3.936 + 4|, |2.9712 − 2.6|, | − 1.007 + 0.92|} = 0.3712, kx(3) − x(2) k R = 0.0568, atd. Výpočet zapíšeme do tabulky: k 0 1 2 3 4

(k)

(k)

x1

x2

1 −4.0000 −3.9360 −3.9928 −3.9997

1 2.6000 2.9712 3.0000 3.0002

(k)

x3

1 −0.9200 −1.0070 −1.0014 −1.0001

k x( k ) − x( k −1) k R — 5.0000 0.3712 0.0568 0.0069

Ukončili jsem po jedenácté iteraci, protože kx(4) − x(3) k R = 0.0069 ≤ 10−2 . Požadovaná přesnost je ε = 10−2 , proto jednotlivé souřadnice výsledku zaokrouhlíme na 2 desetinná místa. Výsledek zapíšeme jako x1 = −4.00 ± 10−2 , x2 = 3.00 ± 10−2 , x3 = −1.00 ± 10−2 .

Příklad 3.3: Vyřešte soustavu lineárních rovnic 3x1 + x2 +2x3 = 12, 5x1 − x2 −4x3 = −11, −2x1 +4x2 −2x3 = 2, jak pomocí Jacobiho, tak pomocí Gaussovy-Seidelovy iterační metody s přesností ε = 10−2 . Všímejte si odlišností obou metod.

33

KAPITOLA 3. SOUSTAVY LINEÁRNÍCH ROVNIC: ITERAČNÍ METODY

Aby metody konvergovaly, je třeba převést soustavu rovnic do ekvivalentního tvaru, jehož matice soustavy je ostře diagonálně dominantní. Zdůrazněme, že tato postačující podmínka konvergence je pro obě metody společná! Všimněte si, že v matici soustavy nenajdeme na žádném řádku ostře dominantní prvek ve smyslu, že by byl v absolutní hodnotě (ostře) větší než součet absolutních hodnot prvků zbylých. Vhodné ekvivalentní úpravy se nám proto budou hledat obtížněji, protože však je matice soustavy regulární (má nenulový determinant - ověřte si!), takové úpravy nutně existují. Se soustavou budeme pracovat v maticovém tvaru, protože je to přehlednější. Po krátkém experimentování (ovšem s vědomím toho, co chceme získat) dostáváme následující úpravy: 

3  5 −2  3 ∼ 1 5

  12 1 2 3 1 1 ∼ 5 −1 −1 −4 −11 2 4 −2 1 5 ← −+   1 2 12 8 0 ← −+ ∼ 1 5 5 0 14  −1 −4 −11 5 −1 1

 12 2 −4 −11 ← − 14 0 ← −  −2 1 0 14  . −4 −11

V tuto chvíli jsme získali ostrou diagonální dominanci na prvním a druhém řádku. Na třetím řádku bychom měli rádi dominantní třetí prvek (tj. −4), přičemž nám vadí první prvek (tj. 5). Nabízí se odečíst od třetího řádku první (na první pozici pak bude místo 5 menší −3), jenže tím dominance nedosáhneme. Ukázalo se, že odečtení celého prvního řádku bylo příliš „hrubou“ operací. Zkusme tedy být jemnější a od třetího řádku odečíst 85 prvního řádku. Po provedení závěrečné kosmetické úpravy dostaneme: 

8 0 1 5 5 −1  8 0  ∼ 1 5 0 8



−2 1 0 14  −4 −11 ← −+  −2 1 0 14 . 22 93

− 58



 8 0 −2 1  14  0 ∼ 1 5  93 0 −1 − 11 − | · −8 4 8

Zapsáno opět v rovnicovém tvaru tedy:

8x1 −2x3 = 1, x1 +5x2 = 14, 8x2 +22x3 = 93. Matice soustavy



 8 0 −2 1 5 0  . 0 8 22

je samozřejmě ostře diagonálně dominantní, neboť absolutní hodnoty prvků na diagonále jsou (ostře) větší než součty absolutních hodnot ostatních prvků na příslušných řádcích, tj. 8 > 0 + 2, 5 > 1 + 0 a 22 > 0 + 8. 34

KAPITOLA 3. SOUSTAVY LINEÁRNÍCH ROVNIC: ITERAČNÍ METODY Upravenou soustavu převedeme do iteračního tvaru x1 = x2 = x3 =

1 8 (1 + 2x3 ), 1 5 (14 − x1 ), 1 (93 − 8x2 ). − 22

V tuto chvíli je třeba napsat rekurentní vzorce a ty jsou u Jacobiho i Gaussovy-Seidelovy metody odlišné. Nejprve rekurentní vzorce pro Jacobiho metodu: ( k +1)

(k)

= 18 (1 + 2x3 ), (k) = 51 (14 − x1 ), (k) 1 = − 22 (93 − 8x2 ).

x1

( k +1)

x2

( k +1)

x3

A dále rekurentní vzorce pro Gaussovu-Seidelovu metodu: ( k +1)

x1

( k +1)

x2

( k +1)

x3

(k)

= 18 (1 + 2x3 ), ( k +1) ), = 51 (14 − x1 ( k +1) 1 = − 22 (93 − 8x2 ). (0)

(0)

(0)

Výpočet konverguje pro libovolnou počáteční aproximaci x(0) = ( x1 , x2 , x3 )> , u obou metod tedy zvolíme pro například x(0) = ( 1, 1, 1 )> . Dosadíme do pravé strany rekurentních vzorců a vypočteme první prvky. V případě Jacobiho metody dostaneme 85 > x(1) = ( 38 , 13 5 , 22 ) . Po prvním kroku je hodnota chyby 85 kx(1) − x(0) k R = max{| 38 − 1|, | 13 5 − 1|, | 22 − 1|} =

63 22

> 10−2 .

V případě Gaussovy-Seidelovy metody to bude x(1) = (

3 8,

109 40 ,

178 55

)>

a chyba bude mít hodnotu 123 −2 55 > 10 . 10−2 . Výpočty zapíšeme

178 kx(1) − x(0) k R = max{| 83 − 1|, | 109 40 − 1|, | 55 − 1|} =

Výpočet opakujeme, dokud není chyba menší nebo rovna tabulek. Jacobiho metoda: Gaussova-Seidelova metoda: k 0 1 2 3 4

(k)

x1

1 0.3750 0.9341 0.9443 0.9444

(k)

x2

1 2.7250 2.6132 2.6111 2.6111

(k)

x3

1 3.2364 3.2770 3.2778 3.2778 35

k x( k ) − x( k −1) k R — 2.2364 0.5591 0.0102 0.0002

do

KAPITOLA 3. SOUSTAVY LINEÁRNÍCH ROVNIC: ITERAČNÍ METODY k 0 1 2 3 4 5 6

(k)

x1

1 0.3750 1.0909 0.9455 0.9341 0.9471 0.9445

(k)

x2

1 2.6000 2.7250 2.5818 2.6109 2.6132 2.6106

(k)

x3

k x( k ) − x( k −1) k R

1 3.8636 3.2818 3.2364 3.2884 3.2779 3.2770

— 2.8636 0.7159 0.1455 0.0521 0.0130 0.0026

Jacobiho metodu jsme ukončili po šesté iteraci a Gaussovu-Seidelovu po čtvrté iteraci. Požadovaná přesnost je ε = 10−2 , proto jednotlivé souřadnice výsledku zaokrouhlíme na 2 desetinná místa. Výsledky obou metod jsou totožné (všimněte si však, že na dalších desetinných místech se liší) a zapíšeme je jako x1 = 0.94 ± 10−2 , x2 = 2.61 ± 10−2 , x3 = 3.28 ± 10−2 .

Příklad 3.4: Vyřešte soustavu lineárních rovnic

−6x1 +2x2 + x3 + x4 2x1 +5x2 + x4 x1 +6x2 +5x3 2x1 + x2 − x3 +5x4

= −17, = −9, = −21, = 12.

pomocí Jacobiho iterační metody s přesností ε = 10−2 . Než aplikujeme rekurentní vzorce pro Jacobiho metodu, musíme zajistit, aby posloupnost vektorů, kterou metoda vytváří, konvergovala ke správnému výsledku. Z teorie víme, že konvergence je zaručena, pokud je matice soustavy ostře diagonálně dominantní (v ostatních případech konvergovat může, ale také nemusí). Pro tento účel stačí od třetí rovnice odečíst rovnici druhou a dostaneme:

−6x1 +2x2 + x3 + x4 2x1 +5x2 + x4 − x1 + x2 +5x3 − x4 2x1 + x2 − x3 +5x4 Matice soustavy



−6  2   −1 2

= −17, = −9, = −12, = 12.

 2 1 1 5 0 1   1 5 −1 1 −1 5

je ostře diagonálně dominantní, neboť absolutní hodnoty prvků na diagonále jsou (ostře) větší než součty absolutních hodnot ostatních prvků na příslušných řádcích, tj. 6 > 2 + 1 + 1, 5 > 2 + 0 + 1, 5 > 1 + 1 + 1 a 5 > 2 + 1 + 1.

36

KAPITOLA 3. SOUSTAVY LINEÁRNÍCH ROVNIC: ITERAČNÍ METODY Upravenou soustavu převedeme na iterační tvar x1 x2 x3 x4

= − 16 (−17 − 2x2 − x3 − x4 ), = 15 (−9 − 2x1 − x4 ), = 51 (−12 + x1 − x2 + x4 ), = 51 (12 − 2x1 − x2 + x3 )

a připsáním iteračních indexů dostaneme rekurentní vzorce pro Jacobiho metodu ( k +1)

x1

( k +1)

x2

( k +1)

x3

( k +1)

x4

= = = =

(k) (k) (k) 1 6 (17 + 2x2 + x3 + x4 ), (k) (k) 1 5 (−9 − 2x1 − x4 ), (k) (k) (k) 1 5 (−12 + x1 − x2 + x4 ), (k) (k) (k) 1 5 (12 − 2x1 − x2 + x3 ).

Poznamenejme, že rekurentní vzorce můžeme zapsat také v maticovém tvaru, který je obzláště výhodný při řešení úlohy pomocí počítačového softwaru.         1 1 1 ( k +1) (k) 17 0 x1 x1 3 6 6 6  ( k +1)   2      x  − 5 0 0 − 15   x (k)   − 95  2 2          (k+1)  =  1 − 1 0 1  ·  (k)  + − 12  . x  x3     5   5  5 4 3 ( k +1) (k) 12 2 1 1 −5 −5 5 0 x4 x4 5 (0)

(0)

(0)

(0)

Protože výpočet konverguje pro libovolnou počáteční aproximaci x(0) = ( x1 , x2 , x3 , x4 )> , zvolíme pro jednoduchost x(0) = ( 0, 0, 0, 0 )> . Dosadíme do pravé strany rekurentních vzorců a vypočteme x(1) = (

17 6,

− 59 , − 12 5,

12 5

)> .

Jestliže takto pokražujeme dále, dostáváme x(2) = ( 2.2333, −3.4133, −0.9933, 1.1467 )> , x(3) = ( 1.7211, −2.9227, −1.0413, 1.9907 )> , atd. Současně přitom zjišťujeme, zda platí

kx(k) − x(k−1) k R ≤ 10−2 , takže postupně počítáme 9 12 12 kx(1) − x(0) k R = max{| 17 6 − 2.2333|, | − 5 + 3.4133|, | − 5 + 0.9933|, | 5 − 1.1467|}

= 2.8333, kx(2) − x(1) k R = 1.6133, kx(3) − x(2) k R = 0.8440, atd. Výpočet zapíšeme do tabulky: 37

KAPITOLA 3. SOUSTAVY LINEÁRNÍCH ROVNIC: ITERAČNÍ METODY k 0 1 2 3 4 5 6 7 8

(k)

(k)

x1

0 2.8333 2.2333 1.7211 2.0173 2.0403 1.9842 1.9963 2.0036

(k)

x2

x3

0 −1.8000 −3.4133 −2.9227 −2.8866 −3.0245 −3.0073 −2.9914 −3.0003

0 −2.4000 −0.9933 −1.0413 −1.0731 −1.0017 −0.9959 −1.0040 −1.0007

(k)

x4

k x( k ) − x( k −1) k R

0 2.4000 1.1467 1.9907 2.0878 1.9558 1.9885 2.0086 1.9989

— 2.8333 1.6133 0.8440 0.2962 0.1379 0.0561 0.0201 0.0097

Ukončili jsem po osmé iteraci, protože kx(8) − x(7) k R = 0.0097 ≤ 10−2 . Požadovaná přesnost je ε = 10−2 , proto jednotlivé souřadnice výsledku zaokrouhlíme na 2 desetinná místa. Výsledek zapíšeme jako x1 = 2.00 ± 10−2 , x2 = −3.00 ± 10−2 , x3 = −1.00 ± 10−2 , x4 = 2.00 ± 10−2 . Příklad 3.5: Vyřešte soustavu lineárních rovnic

−6x1 +2x2 + x3 +2x4 x1 −5x2 − x3 +2x4 2x1 − x2 −5x3 + x4 5x1 −6x3 +7x4

= 16, = 3, = 21, = 54.

pomocí Gaussovy-Seidelovy iterační metody s přesností ε = 10−2 . Z teorie víme, že konvergence Gaussovy-Seidelovy metody je zaručena, pokud je matice soustavy ostře diagonálně dominantní (v ostatních případech konvergovat může, ale také nemusí). Pro dosažení takové matice soustavy stačí od čtvrté rovnice odečíst rovnici třetí a dostaneme: −6x1 +2x2 + x3 +2x4 = 16, x1 −5x2 − x3 +2x4 = 3, 2x1 − x2 −5x3 + x4 = 21, 3x1 + x2 − x3 +6x4 = 33. Matice soustavy



 −6 2 1 2  1 −5 −1 2    2 −1 −5 1 3 1 −1 6 je ostře diagonálně dominantní, neboť absolutní hodnoty prvků na diagonále jsou (ostře) větší než součty absolutních hodnot ostatních prvků na příslušných řádcích, tj. 6 > 2 + 1 + 2, 5 > 1 + 1 + 2, 5 > 2 + 1 + 1 a 6 > 3 + 1 + 1. Upravenou soustavu převedeme na iterační tvar x1 x2 x3 x4

= − 61 (16 − 2x2 − x3 − 2x4 ), = − 15 (3 − x1 + x3 − 2x4 ), = − 15 (21 − 2x1 + x2 − x4 ), = 16 (33 − 3x1 − x2 + x3 ) 38

KAPITOLA 3. SOUSTAVY LINEÁRNÍCH ROVNIC: ITERAČNÍ METODY a připsáním iteračních indexů dostaneme rekurentní vzorce pro Gaussovu-Seidelovu metodu ( k +1) (k) (k) (k) x1 = − 16 (16 − 2x2 − x3 − 2x4 ), ( k +1)

( k +1)

= − 15 (3 − x1

( k +1)

= − 51 (21 − 2x1

( k +1)

=

x2

x3 x4

(k)

( k +1)

( k +1) 1 6 (33 − 3x1

(k)

+ x3 − 2x4 ), ( k +1)

+ x2

( k +1)

− x2

(k)

− x4 ), ( k +1)

+ x3

) (0)

(0)

(0)

(0)

Protože výpočet konverguje pro libovolnou počáteční aproximaci x(0) = ( x1 , x2 , x3 , x4 )> , zvolíme pro jednoduchost x(0) = ( 1, 1, 1, 1 )> . Dosadíme do pravé strany rekurentních vzorců a vypočteme x(1) = ( −1.8333, −0.7667, −4.58, 5.7811 )> . Jestliže takto pokražujeme dále, dostáváme x(2) = ( −1.7585, 2.2767, −4.2025, 5.2994 )> , x(3) = ( −0.8417, 2.1919, −3.9152, 4.9030 )> , atd. Současně přitom zjišťujeme, zda platí


kx(1) − x(0) k R = max{| − 1.8333 − 1|, | − 0.7667 − 1|, | − 4.58 − 1|, |5.7811 − 1|} = 5.58, kx(2) − x(1) k R = 3.0434, kx(3) − x(2) k R = 0.9168, atd. Výpočet zapíšeme do tabulky: k 0 1 2 3 4 5 6 7

(k)

x1

(k)

x2

0 0 −1.8333 −0.7667 −1.7585 2.2767 −0.8417 2.1919 −0.9542 1.9534 −1.0187 1.9891 −1.0017 2.0052 −0.9982 2.0003

(k)

x3

0 −4.5800 −4.2025 −3.9152 −3.9918 −4.0081 −3.9998 −3.9993

(k)

x4

0 5.7811 5.2994 4.9030 4.9863 5.0098 5.0000 4.9992

k x( k ) − x( k −1) k R — 5.5800 3.0434 0.9168 0.2385 0.0645 0.0171 0.0049

Ukončili jsem po sedmé iteraci, protože kx(7) − x(6) k R = 0.0049 ≤ 10−2 . Požadovaná přesnost je ε = 10−2 , proto jednotlivé souřadnice výsledku zaokrouhlíme na 2 desetinná místa.

39

KAPITOLA 3. SOUSTAVY LINEÁRNÍCH ROVNIC: ITERAČNÍ METODY Výsledek zapíšeme jako x1 = −1.00 ± 10−2 , x2 = 2.00 ± 10−2 , x3 = −4.00 ± 10−2 , x4 = 5.00 ± 10−2 .

40

KAPITOLA

4 INTERPOLACE A APROXIMACE FUNKCÍ

Příklad 4.1: Pro uzly xi a funkční hodnoty yi dané následující tabulkou i=1

xi yi

i=2

i=3

i=4

i=5

−3 −1

0

2

3

−1

0

−3

1

0

sestavte interpolační polynom v Lagrangeově tvaru. Interpolační polynom v Lagrangeově tvaru je určen předpisem p ( x ) = y0 ϕ0 ( x ) + y1 ϕ1 ( x ) + y2 ϕ2 ( x ) + y3 ϕ3 ( x ) + y4 ϕ4 ( x ), kde ϕ0 ( x ), ϕ1 ( x ), ϕ2 ( x ), ϕ3 ( x ), ϕ4 ( x ) jsou polynomy Lagrangeovy báze dané úlohy: ϕ0 ( x ) = ϕ0 ( x ) = ϕ0 ( x ) = ϕ0 ( x ) = ϕ0 ( x ) =

( x + 1)( x − 0)( x − 2)( x − 3) 1 = x ( x + 1)( x − 2)( x − 3), (−3 + 1)(−3 − 0)(−3 − 2)(−3 − 3) 180 ( x + 3)( x − 0)( x − 2)( x − 3) 1 = − x ( x + 3)( x − 2)( x − 3), (−1 + 3)(−1 − 0)(−1 − 2)(−1 − 3) 24 ( x + 3)( x + 1)( x − 2)( x − 3) 1 = ( x + 3)( x + 1)( x − 2)( x − 3), (0 + 3)(0 + 1)(0 − 2)(0 − 3) 18 ( x + 3)( x + 1)( x − 0)( x − 3) 1 = − x ( x + 3)( x + 1)( x − 3), (2 + 3)(2 + 1)(2 − 0)(2 − 3) 30 ( x + 3)( x + 1)( x − 0)( x − 2) 1 = x ( x + 3)( x + 1)( x − 2), (3 + 3)(3 + 1)(3 − 0)(3 − 2) 72

Dohromady pak dostáváme 41

KAPITOLA 4. INTERPOLACE A APROXIMACE FUNKCÍ

p( x ) =

1 1 x ( x + 1)( x − 2)( x − 3) − ( x + 3)( x + 1)( x − 2)( x − 3)− 180 18 1 − x ( x + 3)( x + 1)( x − 2). 24

Příklad 4.2: Pro uzly xi a funkční hodnoty yi dané následující tabulkou i=1

xi yi

i=2

i=3

i=4

−3 −1

0

2

3

−1

0

−3

1

0

sestavte interpolační polynom v Newtonově tvaru. Interpolační polynom v Newtonově tvaru je určen předpisem p( x ) = y0 + f [ x1 , x0 ]( x − x0 ) + f [ x2 , x1 , x0 ]( x − x0 )( x − x1 ) +

+ f [ x3 , x2 , x1 , x0 ]( x − x0 )( x − x1 )( x − x2 ) + + f [ x4 , x3 , x2 , x1 , x0 ]( x − x0 )( x − x1 )( x − x2 )( x − x3 ), kde f [ x1 , x0 ] , f [ x2 , x1 , x0 ], f [ x3 , x2 , x1 , x0 ] a f [ x4 , x3 , x2 , x1 , x0 ] jsou poměrné diference 1.,2.,3. a 4. řádu. Poměrné diference prvního řádu. Spočítáme všechny poměrného diference prvního řádu podle vzorce: f [ x i +1 , x i ] =

f i +1 − f i x i +1 − x i

i = 0, 1, . . . , n − 1

pro

V tomto příkladě je počet uzlů n = 5, budou tedy čtyři poměrné diference prvního řádu: pro

i=0

pro

i=1

pro

i=2

pro

i=3

f1 − f0 x1 − x0 f2 − f1 f [ x2 , x1 ] = x2 − x1 f3 − f2 f [ x3 , x2 ] = x3 − x2 f − f3 f [ x4 , x3 ] = 4 x4 − x3 f [ x1 , x0 ] =

0−1 −1 = −1 − (−3) 2 −1 − 0 = = −1 0 − (−1) 0 − (−1) 1 = = 2−0 2 −3 − 0 = = −3 3−2

=

Poměrné diference druhého řádu. Spočítáme všechny poměrného diference druhého řádu podle vzorce: f [ x i +2 , x i +1 , x i ] =

f [ x i +2 , x i +1 ] − f [ x i +1 , x i ] x i +2 − x i

42

pro

i = 0, 1, . . . , n − 2

KAPITOLA 4. INTERPOLACE A APROXIMACE FUNKCÍ Poměrné diference druhého řádu budou tři: i=0

−1 − 12 f [ x2 , x1 ] − f [ x1 , x0 ] 1 f [ x2 , x1 , x0 ] = = =− x2 − x0 0 − (−3) 6

pro

i=1

1 f [ x3 , x2 ] − f [ x2 , x1 ] 1 2 − (−1) f [ x3 , x2 , x1 ] = = = x3 − x1 2 − (−1) 2

pro

i=2

f [ x4 , x3 , x2 ] =

pro

−3 − (− 21 ) 7 f [ x4 , x3 ] − f [ x3 , x2 ] = =− x4 − x2 3 − (−1) 6

Poměrné diference třetího řádu. Spočítáme všechny poměrného diference třetího řádu podle vzorce: f [ x i +3 , x i +2 , x i +1 ] − f [ x i +2 , x i +1 , x i ] x i +3 − x i

f [ x i +3 , x i +2 , x i +1 , x i ] =

pro

i = 0, 1, . . . , n − 3

Poměrné diference třetího řádu budou dvě: pro

i=0

f [ x3 , x2 , x1 , x0 ] =

1 − (− 61 ) 2 f [ x3 , x2 , x1 ] − f [ x2 , x1 , x0 ] = 2 = x3 − x0 2 − (−3) 15

pro

i=1

f [ x4 , x3 , x2 , x1 ] =

− 76 − 12 5 f [ x4 , x3 , x2 ] − f [ x3 , x2 , x1 ] =− = x4 − x1 3 − (−1) 12

Poměrné diference čtvrtého řádu. Spočítáme poměrnou diference čtvrtého řádu podle vzorce:

f [ x i +4 , x i +3 , x i +2 , x i +1 , x i ] =

f [ x i +4 , x i +3 , x i +2 , x i +1 ] − f [ x i +3 , x i +2 , x i +1 , x i ] x i +4 − x i

pro

i = 0, 1, . . . , n − 4

Poměrná diference čtvrého řádu bude jedna: pro

i=0

5 2 − 12 − 15 f [ x4 , x3 , x2 , x1 ] − f [ x3 , x2 , x1 , x0 ] 11 f [ x4 , x3 , x2 , x1 , x0 ] = = =− x4 − x0 3 − (−3) 120

Poměrných diference lze zapsat do tabulky: i

xi

yi

1.řád

2.řád

3.řád

4.řád

0

−3

1

− 12

− 16

11 − 120

1

−1

0

−1

2

0

−1

1 2

1 2 − 76

2 15 5 − 12

3

2

0

−3

4

3

−3 43

KAPITOLA 4. INTERPOLACE A APROXIMACE FUNKCÍ Interpolační polynom v Newtonově tvaru je určen předpisem: p( x ) = f 0 + f [ x1 , x0 ]( x − x0 ) + f [ x2 , x1 , x0 ]( x − x0 )( x − x1 ) +

+ f [ x3 , x2 , x1 , x0 ]( x − x0 )( x − x1 )( x − x2 ) + + f [ x4 , x3 , x2 , x1 , x0 ]( x − x0 )( x − x1 )( x − x2 )( x − x3 ), Pomocí hodnot vypočtených v prvním řádku tabulky sestavíme interpolační polynom: 1 2 11 1 x ( x − 2)( x + 1)( x + 3). p( x ) = 1 − ( x + 3) − ( x + 1)( x + 3) + x ( x + 1)( x + 3) − 2 6 15 120 Příklad 4.3: Pro uzly xi a funkční hodnoty yi dané následující tabulkou i=1

xi yi

i=2

i=3

−3 −1

1

−1

2

1

sestavte interpolační polynom: a) v základním tvaru, b) v Lagrangeově tvaru, c) v Newtonově tvaru. a) Interpolační polynom hledáme ve tvaru p( x ) = a0 + a1 x + a2 x2 . Dosazením do interpolačních požadavků p( xi ) = yi , i = 0, 1, 2 dostaneme p(−3) = 1 =⇒ a0 − 3a1 + 9a2 = 1, p(−1) = −1 =⇒ a0 − a1 + a2 = −1, p (1) = 2 =⇒ a0 + a1 + a2 = 2. Tyto rovnosti představují soustavu lineárních rovnic, kterou je možno zapsat také ve tvaru      1 −3 9 a0 1  1 −1 1   a1  =  −1  . 1 1 1 a2 2

Vyřešením dostaneme a0 = − 18 , a1 = 32 , a2 = 58 , takže hledaným interpolační polynom má tvar 1 3 5 p( x ) = − + x + x2 . 8 2 8 b) Interpolační polynom v Lagrangeově tvaru je určen předpisem p ( x ) = y0 ϕ0 ( x ) + y1 ϕ1 ( x ) + y2 ϕ2 ( x ), kde ϕ0 ( x ), ϕ1 ( x ), ϕ2 ( x ) jsou polynomy Lagrangeovy báze dané úlohy:

( x + 1)( x − 1) (−3 + 1)(−3 − 1) ( x + 3)( x − 1) ϕ1 ( x ) = (−1 + 3)(−1 − 1) ( x + 3)( x + 1) ϕ2 ( x ) = = (1 + 3)(1 + 1)

ϕ0 ( x ) =

44

1 ( x + 1)( x − 1), 8 1 = − ( x + 3)( x − 1), 4 1 ( x + 3)( x + 1). 8

=

KAPITOLA 4. INTERPOLACE A APROXIMACE FUNKCÍ Dohromady pak dostáváme p( x ) =

1 1 1 ( x + 1)( x − 1) + ( x + 3)( x − 1) + ( x + 3)( x + 1). 8 4 4

c) Interpolační polynom v Newtonově tvaru je určen předpisem p( x ) = y0 + f [ x1 , x0 ]( x − x0 ) + f [ x2 , x1 , x0 ]( x − x0 )( x − x1 ), kde f [ x1 , x0 ] a f [ x2 , x1 , x0 ] jsou poměrné diference 1. a 2. řádu. Výpočet poměrných diferencí je proveden v následující tabulce: i

xi

0 1 2

yi

1.řád

2.řád

−3 1 −1 −1 1 2

−1

5 8

3 2

Pomocí hodnot vypočtených v prvním řádku tabulky sestavíme interpolační polynom: 5 p( x ) = 1 − ( x + 3) + ( x + 3)( x + 1). 8

p( x ) = −

1 3 5 + x + x2 8 2 8

1 1 1 ( x + 1)( x − 1) + ( x + 3)( x − 1) + ( x + 3)( x + 1) 8 4 4 5 p( x ) = 1 − ( x + 3) + ( x + 3)( x + 1) 8 p( x ) =

Příklad 4.4: Aproximujte následující data i=1

i=2

i=3

i=4

i=5

xi

1

2.5

3

5

7

yi

0

1

6.5

6

15

přímkou ϕ( x ) = c1 x + c2 metodou nejmenších čtverců. Pro nalezení nejlepší aproximace ve smyslu nejmenších čtverců je třeba nejprve sestavit normální rovnice. Pokud hledáme aproximaci ve tvaru ϕ( x ) = c1 ϕ1 ( x ) + c2 ϕ2 ( x ), kde c1 , c2 ∈ R jsou vhodné koeficienty, a v bodech xi máme předepsané hodnoty yi , přičemž i = 1, . . . , n, pak mají normální rovnice tvar n

n

c1 ∑ ( ϕ1 ( xi )) + c2 ∑ ϕ1 ( xi ) · ϕ2 ( xi ) = i =1 n

2

i =1

∑ y i · ϕ1 ( x i )

n

i =1 n

i =1

i =1

c1 ∑ ϕ2 ( xi ) · ϕ1 ( xi ) + c2 ∑ ( ϕ2 ( xi ))2 = i =1

n

45

∑ y i · ϕ2 ( x i )

KAPITOLA 4. INTERPOLACE A APROXIMACE FUNKCÍ V našem případě je ϕ1 ( x ) = x a ϕ2 ( x ) = 1 (jedná se o konstantní funkci, která má v každém bodě hodnotu 1). Nyní již soustavu normálních rovnic pro neznámé koeficienty c1 , c2 snadno zapíšeme. 5

5

i =1 5

i =1 5

c1 ∑ xi2 + c2 ∑ xi = c1 ∑ x i + c2 ∑ 1 = i =1

i =1

5

∑ yi · xi

i =1 5

∑ yi

i =1

Vypočítáme jednotlivé součty: 5

∑ xi2

= 12 + 2.52 + 32 + 52 + 72 = 90.25

∑ xi

= 1 + 2.5 + 3 + 5 + 7 = 18.5

i =1 5

i =1 5

∑ yi · xi

= 0 · 1 + 1 · 2.5 + 6.5 · 3 + 6 · 5 + 15 · 7 = 157

i =1

5

∑ yi

= 0 + 1 + 6.5 + 6 + 15 = 28.5

i =1

A dosadíme do soustavy lineárních rovnic: 90.25c1 + 18.5c2 = 157 18.5c1 + 5c2 = 28.5 Soustavu vyřešíme a dostaneme koeficienty c1 = 2.3647, c2 = −3.0493. Nalezená funkce má rovnici ϕ( x ) = 2.3647x − 3.0493.

46

KAPITOLA 4. INTERPOLACE A APROXIMACE FUNKCÍ 15 zadana data nalezena funkce

10

5

0 1

2

3

4

5

6

7


i=2

i=3

i=4

i=5

i=6

i=7

i=8

i=9

i=10

xi

0.2

0.8

4.6

6.2

6.6

11.5

11.8

13.1

14.6

14.9

yi

29

11

-42

-36

-35

-59

-55

-46

-47

-46

funkcí ϕ( x ) = c1 sin x + c2 ln x metodou nejmenších čtverců. Pro nalezení nejlepší aproximace ve smyslu nejmenších čtverců je třeba nejprve sestavit normální rovnice. Pokud hledáme aproximaci ve tvaru ϕ( x ) = c1 ϕ1 ( x ) + c2 ϕ2 ( x ), kde c1 , c2 ∈ R jsou vhodné koeficienty, a v bodech xi máme předepsané hodnoty yi , přičemž i = 1, . . . , n, pak mají normální rovnice tvar n

n

i =1 n

i =1

c1 ∑ ( ϕ1 ( xi ))2 + c2 ∑ ϕ1 ( xi ) · ϕ2 ( xi ) = n

i =1 n

i =1

i =1

c1 ∑ ϕ2 ( xi ) · ϕ1 ( xi ) + c2 ∑ ( ϕ2 ( xi ))2 = i =1

n

∑ y i · ϕ1 ( x i )

∑ y i · ϕ2 ( x i )

V našem případě je ϕ1 ( x ) = sin x a ϕ2 ( x ) = ln x. Nyní již soustavu normálních rovnic pro neznámé koeficienty c1 , c2 snadno zapíšeme. 10

10

c1 ∑ sin xi + c2 ∑ sin xi ln xi = i =1 10

2

i =1

∑ yi · sin xi

10

i =1 10

i =1

i =1

c1 ∑ ln xi sin xi + c2 ∑ ln2 xi = i =1

10

47

∑ yi · ln xi

KAPITOLA 4. INTERPOLACE A APROXIMACE FUNKCÍ Vypočítáme jednotlivé součty: 10

∑ sin2 xi

= sin2 (0.2) + sin2 (0.8) + sin2 (4.6) + sin2 (6.2) + sin2 (6.6) + sin2 (11.5) +

i =1

+ sin2 (11.8) + sin2 (13.1) + sin2 (14.6) + sin2 (14.9) = 4.4748 10

∑ sin xi ln xi

= sin(0.2) · ln(0.2) + sin(0.8) · ln(0.8) + sin(4.6) · ln(4.6) +

i =1

+ sin(6.2) · ln(6.2) + sin(6.6) · ln(6.6) + sin(11.5) · ln(11.5) + + sin(11.8) · ln(11.8) + sin(13.1) · ln(13.1) + sin(14.6) · ln(14.6) + + sin(14.9) · ln(14.9) = 0.2505 10

∑ ln2 xi

= ln2 (0.2) + ln2 (0.8) + ln2 (4.6) + ln2 (6.2) + ln2 (6.6) + ln2 (11.5) +

i =1

+ ln2 (11.8) + ln2 (13.1) + ln2 (14.6) + ln2 (14.9) = 45.0191 10

∑ yi sin xi

= 29 · sin(0.2) + 11 · sin(0.8) − 42 · sin(4.6) − 36 · sin(6.2) −

i =1

−35 · sin(6.6) − 59 · sin(11.5) − 55 · sin(11.8) − 46 · sin(13.1) − −47 · sin(14.6) − 46 · sin(14.9) = 38.564 10

∑ yi ln xi

= 29 · ln(0.2) + 11 · ln(0.8) − 42 · ln(4.6) − 36 · ln(6.2) − 35 · ln(6.6) −

i =1

−59 · ln(11.5) − 55 · ln(11.8) − 46 · ln(13.1) − 47 · ln(14.6) −46 · ln(14.9) = −893.4086

A dosadíme do soustavy lineárních rovnic: 4.4748c1 + 0.2505c2 = 38.564 0.2505c1 + 45.0191c2 = −893.4086 Soustavu vyřešíme a dostaneme koeficienty c1 = 9.7321, c2 = −19.8993. Nalezená funkce má rovnici ϕ( x ) = 9.7321 sin x − 19.8993 ln x.

48


zadana data nalezena funkce

30

20

10

0

−10

−20

−30

−40

−50

2

4

6

8

10

12

14


i=2

i=3

i=4

i=5

xi

-3

-1

3

4

8

yi

2

1

5.5

4.4

0.5

pomocí metody nejmenších čtverců nejprve přímkou ϕ( x ) = c1 x + c2 a posléze funkcí ψ( x ) = d1 sin x + d2 cos x. Určete, která ze získaných funkcí aproximuje zadaná data lépe. Nejprve je třeba sestavit normální rovnice. V obou případech jde o aproximaci tvaru ϕ( x ) = c1 ϕ1 ( x ) + c2 ϕ2 ( x ), kde c1 , c2 ∈ R jsou vhodné koeficienty. Pokud máme v bodech xi předepsané hodnoty yi , přičemž i = 1, . . . , n, pak mají normální rovnice tvar n

n

i =1 n

i =1

c1 ∑ ( ϕ1 ( xi ))2 + c2 ∑ ϕ1 ( xi ) · ϕ2 ( xi ) = n

i =1 n

i =1

i =1

c1 ∑ ϕ2 ( xi ) · ϕ1 ( xi ) + c2 ∑ ( ϕ2 ( xi ))2 = i =1

n

∑ y i · ϕ1 ( x i )

∑ y i · ϕ2 ( x i )

Po zkušenostech získaných v předchozích úlohách jsme s normálními rovnicemi obeznámeni natolik, že se můžeme podívat, jak vypadá přepis do maticového tvaru: ! ! ! 2 c1 ∑in=1 ( ϕ1 ( xi )) ∑in=1 ϕ1 ( xi ) · ϕ2 ( xi ) ∑in=1 yi · ϕ1 ( xi ) · = 2 c2 ∑in=1 yi · ϕ2 ( xi ) ∑in=1 ϕ2 ( xi ) · ϕ1 ( xi ) ∑in=1 ( ϕ2 ( xi )) 49

KAPITOLA 4. INTERPOLACE A APROXIMACE FUNKCÍ V prvním případě je ϕ1 ( x ) = x, ϕ2 ( x ) = 1 a ve druhém ϕ1 ( x ) = sin x, ϕ2 ( x ) = cos x. Nyní již obě soustavy normálních rovnic snadno zapíšeme. Aproximace funkcí tvaru ϕ( x ) = c1 x + c2 : ! ! ∑5i=1 xi2 ∑5i=1 xi ∑in=1 yi xi c1 = · c2 ∑in=1 yi ∑5i=1 xi ∑5i=1 1 Po dosazení hodnot ze zadání dostaneme soustavu rovnic: 31.1 99 11 c1 = · 13.4 11 5 c2 Soustavu vyřešíme a dostaneme koeficienty (zaokrouhleno na 2 desetinná místa) c1 = 0.02, c2 = 2.63. Nalezená funkce má rovnici ϕ( x ) = 0.02 x + 2.63. Aproximace funkcí tvaru ψ( x ) = d1 sin x + d2 cos x: ! ∑5i=1 sin2 xi ∑5i=1 sin xi cos xi · ∑5i=1 cos xi sin xi ∑5i=1 cos2 xi

d1 d2

!

=

∑in=1 yi sin xi ∑in=1 yi cos xi

!

Po dosazení hodnot ze zadání dostaneme soustavu rovnic (zaokrouhleno na 2 desetinná místa): 2.24 −0.23 d1 −10.87 · = −0.23 2.76 d2 −7.63 Soustavu vyřešíme a dostaneme koeficienty (zaokrouhleno na 2 desetinná místa) c1 = −1.55, c2 = −3.7. Nalezená funkce má rovnici ψ( x ) = −1.55 sin x − 3.7 cos x. Jak určit, která z aproximací je lepší? Někdy se to zdá být zjevné z obrázku, ale rádi bychom uměli „dobrost“ aproximace měřit exaktně. Principem metody nejmenších čtverců je nalézt takovou funkci určitého tvaru, jejíž součet druhých mocnin (tj. čtverců) rozdílů funkčních hodnot a zadaných hodnot v předepsaných bodech je co nejmenší. Pro měření „dobrosti“ aproximace proto využijeme právě těchto součtů. Čím menší součet dostaneme, tím lepší je aproximace. V případě, že dostaneme soušet nulový, je aproximace ideální a předepsaných bodech nabývá předepsaných hodnot. Pro aproximaci ϕ( x ) dat xi , yi , kde i = 1, . . . , n, definujeme chybu jako ∑in=1 ( ϕ( xi ) − yi )2 . Chyba aproximace ϕ( x ) = 0.02x + 2.63 je proto 5

∑ (0.02xi + 2.63 − yi )2 = 18.91

i =1

Chyba aproximace ψ( x ) = −1.55 sin x − 3.7 cos x je proto 5

∑ (−1.55 sin xi − 3.7 cos xi − yi )2 = 13.53

i =1

50

KAPITOLA 4. INTERPOLACE A APROXIMACE FUNKCÍ Protože 13.53 < 18.91, je lepší aproximací dat funkce ψ( x ). 8 zadana data 1. nalezena funkce 2. nalezena funkce 6

4

2

0

−2

−4

−6 −3

−2

−1

0

1

2

3

4

5

6

7

8


i=2

i=3

i=4

i=5

i=6

i=7

xi

0.2

2

2.7

4.5

4.7

7.5

9.3

yi

0

3.4

3.8

5.8

6.4

6.6

7.1

funkcí ϕ( x ) = c1 x2 + c2 x + c3 metodou nejmenších čtverců. Pro nalezení nejlepší aproximace ve smyslu nejmenších čtverců je třeba nejprve sestavit normální rovnice. Pokud hledáme aproximaci ve tvaru ϕ( x ) = c1 ϕ1 ( x ) + c2 ϕ2 ( x ) + c3 ϕ3 ( x ), kde c1 , c2 , c3 ∈ R jsou vhodné koeficienty, a v bodech xi máme předepsané hodnoty yi , přičemž i = 1, . . . , n, pak mají normální rovnice tvar n

n

n

i =1

n

i =1 n

i =1 n

i =1

i =1

c1 ∑ ( ϕ1 ( xi )) + c2 ∑ ϕ1 ( xi ) · ϕ2 ( xi ) + c3 ∑ ϕ1 ( xi ) · ϕ3 ( xi ) = 2

i =1 n

c1 ∑ ϕ2 ( xi ) · ϕ1 ( xi ) + c2 ∑ ( ϕ2 ( xi ))2 + c3 ∑ ϕ2 ( xi ) · ϕ3 ( xi ) =

n

∑ y i · ϕ1 ( x i )

∑ y i · ϕ2 ( x i )

n

n

n

i =1 n

i =1

i =1

i =1

i =1

i =1

c1 ∑ ϕ3 ( xi ) · ϕ1 ( xi ) + c2 ∑ ( ϕ2 ( xi ))2 + c3 ∑ ( ϕ3 ( xi ))2 =

∑ y i · ϕ3 ( x i )

V našem případě je ϕ1 ( x ) = x2 , ϕ2 ( x ) = x a ϕ3 ( x ) = 1 (tj. funkce, která má v každém bodě hodnotu 1). Nyní již soustavu normálních rovnic pro neznámé koeficienty c1 , c2 , c3 snadno zapíšeme. 51


7

7

7

i =1 7

i =1 7

i =1 7

i =1 7

i =1 7

i =1 7

i =1 7

i =1 7

c1 ∑ xi4 + c2 ∑ xi3 + c3 ∑ xi2 = c1 ∑ xi3 + c2 ∑ xi2 + c3 ∑ xi = c1 ∑ xi2 + c2 ∑ xi + c3 ∑ 1 = i =1

i =1

i =1

7

∑ yi · xi2

∑ yi · xi

∑ yi

i =1


∑ xi4 = 0.24 + 24 + 2.74 + 4.54 + 4.74 + 7.54 + 9.34 = 11611.7589

i =1 7

∑ xi3 = 0.23 + 23 + 2.73 + 4.53 + 4.73 + 7.53 + 9.33 = 1448.71

i =1 7

∑ xi2 = 0.22 + 22 + 2.72 + 4.52 + 4.72 + 7.52 + 9.32 = 196.41

i =1 7

∑ xi = 0.2 + 2 + 2.7 + 4.5 + 4.7 + 7.5 + 9.3 = 30.9

i =1 7

∑ yi · xi2 = 0 · 0.22 + 3.4 · 22 + 3.8 · 2.72 + 5.5 · 4.52 + 6.4 · 4.72 + 6.6 · 7.54 +

i =1

+ 7.1 · 9.34 = 1285.457 7

∑ yi · xi = 0 · 0.2 + 3.4 · 2 + 3.8 · 2.7 + 5.5 · 4.5 + 6.4 · 4.7 + 6.6 · 7.5+

i =1

+ 7.1 · 9.3 = 188.77 7

∑ yi = 0 + 3.4 + 3.8 + 5.8 + 6.4 + 6.6 + 7.1 = 33.1

i =1

A dosadíme do soustavy lineárních rovnic: 11611.7589c1 + 1448.71c2 + 196.41c3 = 1285.457 1448.71c1 + 196.41c2 + 30.9c3 = 188.77 196.41c1 + 30.9c2 + 7c3 = 33.1 Soustavu vyřešíme a dostaneme koeficienty c1 = −0.124, c2 = 1.9131, c3 = −0.2375. Nalezená funkce má tvar: ϕ( x ) = −0.124x2 + 1.9131x − 0.2375.

52

KAPITOLA 4. INTERPOLACE A APROXIMACE FUNKCÍ 7


6

5

4

3

2

1

0

1

2

3

4

5

6

7

8

9


i=2

i=3

i=4

i=5

i=6

i=7

i=8

i=9

xi

-2.1

-1.6

-0.3

0.2

0.9

1.4

1.9

2.5

3.1

yi

4.4

4

3.2

3.3

3.6

3.8

4.7

6

7.9

funkcí ϕ( x ) = c1 e x + c2 x2 + c3 metodou nejmenších čtverců. Pro nalezení nejlepší aproximace ve smyslu nejmenších čtverců je třeba nejprve sestavit normální rovnice. Pokud hledáme aproximaci ve tvaru ϕ( x ) = c1 ϕ1 ( x ) + c2 ϕ2 ( x ) + c3 ϕ3 ( x ), kde c1 , c2 , c3 ∈ R jsou vhodné koeficienty, a v bodech xi máme předepsané hodnoty yi , přičemž i = 1, . . . , n, pak mají normální rovnice tvar n

n

i =1 n

i =1

n

c1 ∑ ( ϕ1 ( xi ))2 + c2 ∑ ϕ1 ( xi ) · ϕ2 ( xi ) + c3 ∑ ϕ1 ( xi ) · ϕ3 ( xi ) =

∑ y i · ϕ1 ( x i )

n

i =1 n

i =1 n

i =1

i =1

n

i =1 n

i =1

i =1

c1 ∑ ϕ2 ( xi ) · ϕ1 ( xi ) + c2 ∑ ( ϕ2 ( xi ))2 + c3 ∑ ϕ2 ( xi ) · ϕ3 ( xi ) = i =1

n

n

n

i =1

i =1

c1 ∑ ϕ3 ( xi ) · ϕ1 ( xi ) + c2 ∑ ( ϕ2 ( xi )) + c3 ∑ ( ϕ3 ( xi ))2 = 2

∑ y i · ϕ2 ( x i )

∑ y i · ϕ3 ( x i )

V našem případě je ϕ1 ( x ) = e x , ϕ2 ( x ) = x2 a ϕ3 ( x ) = 1 (tj. funkce, která má v každém bodě hodnotu 1). Nyní již soustavu normálních rovnic pro neznámé koeficienty c1 , c2 , c3 snadno zapíšeme.

53


9

c1 ∑ e

2xi

i =1 9

9

+ c2 ∑

xi2 e xi

i =1

9

+ c3 ∑ e

xi

9

=

9

i =1 9

i =1 9

i =1 9

i =1 9

i =1 9

c1 ∑ e xi xi2 + c2 ∑ xi4 + c3 ∑ xi2 = i =1

9

c1 ∑ e xi + c2 ∑ xi2 + c3 ∑ 1 = i =1

∑ y i · e xi

i =1

i =1

∑ yi · xi2

∑ yi

i =1


∑ e2xi

i =1 9

∑ xi2 exi

= e−4.2 + e−3.2 + e−0.6 + e0.4 + e1.8 + e2.8 + e3.8 + e5 + e6.2 = 710.4541 = −2.1 · e−2.1 − 1.6 · e−1.6 − 0.3 · e−0.3 + 0.2 · e0.2 + 0.9 · e0.9 +

i =1

+1.4 · e1.4 + 1.9 · e1.9 + 2.5 · e2.5 + 3.1 · e3.1 = 324.7119 9

∑ e xi

= e−2.1 + e−1.6 + e−0.3 + e0.2 + e0.9 + e1.4 + e1.9 + e2.5 + e3.1 = 49.8677

∑ xi4

= (−2.1)4 + (−1.6)4 + (−0.3)4 + 0.24 + 0.94 + 1.44 + 1.94 + 2.54 +

i =1 9

i =1

+3.14 = 174.9558 9

∑ xi2

= (−2.1)2 + (−1.6)2 + (−0.3)2 + 0.22 + 0.92 + 1.42 + 1.92 +

i =1

+2.52 + 3.12 = 29.34 9

∑ y i · e xi

= 4.4 · e−2.1 + 4 · e−1.6 + 3.2 · e−0.3 + 3.3 · e0.2 + 3.6 · e0.9 + 3.8 · e1.4 +

i =1

+4.7 · e1.9 + 6 · e2.5 + 7.9 · e3.1 = 311.8945 9

∑ yi · xi2

= 4.4 · (−2.1)2 + 4 · (−1.6)2 + 3.2 · (−0.3)2 + 3.3 · 0.22 + 3.6 · 0.92 +

i =1

+3.8 · 1.42 + 4.7 · 1.92 + 6 · 2.52 + 7.9 · 3.12 = 170.814 9

∑ yi

= 4.4 + 4 + 3.2 + 3.3 + 3.6 + 3.8 + 4.7 + 6 + 7.9 = 40.9

i =1

A dosadíme do soustavy lineárních rovnic: 710.4541c1 + 324.7119c2 + 49.8677c3 = 311.8945 324.7119c1 + 174.9558c2 + 29.34c3 = 170.814 49.8677c1 + 29.34c2 + 9c3 = 40.9 Soustavu vyřešíme a dostaneme koeficienty c1 = 0.0842 c2 = 0.3004 c3 = 3.0985. 54

KAPITOLA 4. INTERPOLACE A APROXIMACE FUNKCÍ Nalezená funkce má tvar ϕ( x ) = 0.0842e x + 0.3004x2 + 3.0985. zadana data nalezena funkce

7.5

7

6.5

6

5.5

5

4.5

4

3.5 −2

−1.5

−1

−0.5

0

0.5

1

1.5

2

2.5

3


i=2

i=3

i=4

i=5

i=6

xi

-5.2

-4.1

0.7

1.5

1.7

3.6

yi

-12.3

1.2

-5.6

3.7

3.2

7.8

funkcí ϕ( x ) = c1 x3 + c2 x2 + c3 x + c4 metodou nejmenších čtverců. Pro nalezení nejlepší aproximace ve smyslu nejmenších čtverců je třeba nejprve sestavit normální rovnice. Ty bychom mohli v případě aproximace tvaru ϕ( x ) = c1 ϕ1 ( x ) + c2 ϕ2 ( x ) + c3 ϕ3 ( x ) + c4 ϕ4 ( x ), kde c1 , c2 , c3 , c4 ∈ R jsou vhodné koeficienty, sestavit analogickým způsobem, jako v předchozích úlohách. Pokud je koeficientů větší počet, je výhodné uchýlit se k „triku“, který zde popíšeme. Soustavu normálních rovnic totiž můžeme elegantně zapsat v maticovém tvaru. Předpokládejme, že v bodech xi máme předepsané hodnoty yi , přičemž i = 1, . . . , n. Definujme následující objekty: y = ( y1 , y2 , . . . , y n ) > c = ( c1 , c2 , c3 , c4 ) >   3 3 x1 x2 . . . xn3  x2 x2 . . . x2   n X= 1 2   x1 x2 . . . x n  1

1

55

...

1

KAPITOLA 4. INTERPOLACE A APROXIMACE FUNKCÍ Pak můžeme soustavu normálních rovnic zapsat jako (není to samozřejmé, ale pravdivost tohoto tvrzení snadno ověříte roznásobením matic) XX > c = Xy Naším cílem je nalézt vektor c, čehož dosáhneme vynásobením rovnice zleva inverzní maticí k XX > , tj.

( XX > )−1 XX > c = ( XX > )−1 Xy c = ( XX > )−1 Xy Zbývá nám tedy v zásadě vytvořit matici X a dosadit do předchozího vzorce. Matice X vypadá následovně (její hodnoty jsou zaokrouhleny na 2 desetinná místa) 

−140.61 −68.92  27.04 16.81 X=  −5.20 −4.10 1.00 1.00

0.34 0.49 0.70 1.00

3.38 2.25 1.50 1.00

 4.91 46.66 2.89 12.96  1.70 3.60  1.00 1.00

Dosazením do vzorce získáme vektor c = (0.22, 0.30, −1.87, 1.54)> (opět zaokrouhlujeme na 2 desetinná místa). Hledaná funkce je proto ϕ( x ) = 0.22x3 + 0.30x2 − 1.87x + 1.54. Na závěr ještě poznamenejme, že uvedený postup lze s drobnou obměnou (matice X mění rozměry) použít pro aproximaci polynomy libovolných řádů. zadana data nalezena funkce

8 6 4 2 0 −2 −4 −6 −8 −10 −12 −5

−4

−3

−2

−1

56

0

1

2

3

KAPITOLA 4. INTERPOLACE A APROXIMACE FUNKCÍ Příklad 4.10: Aproximujte následující data i=1

i=2

i=3

i=4

i=5

i=6

i=7

i=8

xi

0.3

1.2

2.1

3.3

4.1

5.2

6.3

6.8

yi

21

7

3

6

13

40

55

61

funkcí ϕ( x ) = aebx metodou nejmenších čtverců. Zadaná funkce ϕ( x ) není ve tvaru aϕ1 ( x ) + bϕ2 ( x ), jak jsme byli zvyklí z předchozích úloh. Rozdíl spočívá v tom, že funkce není v proměnných a, b lineární. Naším cílem stále zůstává nalézt hodnoty a, b ∈ R takové, aby s co možná nejmenší odchylkou platilo aebxi = yi , pro i = 1, . . . , 8. Dohodli jsme se, že minimální odchylkou budeme rozumět takovou, při níž je minimální součet čtverců, tj. veličina 2 8 bxi ∑ ae − yi . i =1

Pokud bychom se pokusili obvyklým způsobem odvodit normální rovnice, nedostali bychom soustavu lineárních rovnic (zkuste si samostatně!), ale soustavu velmi obtížně řešitelnou. Pomůžeme si proto určitým trikem. Rovnosti, kterých bychom rádi dosáhli, zlogaritmujeme: ln( aebxi ) = ln yi , ln ebxi + ln a = ln yi , bxi + ln a = ln yi . Tím jsme dostali ekvivalentní soubor podmínek. Pokud zavedeme nové proměnné c1 = b, c2 = ln a a zi = ln yi pro i = 1, . . . , 8, vypadají podmínky takto c1 xi + c2 = zi , kde i = 1, . . . , 8. Tato nová úloha je ovšem již lineární! Jde v ní o to aproximovat data [ xi , zi ], i = 1, . . . , 8, funkcí tvaru ψ( x ) = c1 x + c2 , přičemž c1 c2 ∈ R. Toho můžeme docílit metodou nejmenších čtverců, jak jsme ji procvičovali v předchozích úlohách. Nyní již snadno zapíšeme příslušné normální rovnice: 8

8

i =1 8

i =1 8

c1 ∑ xi2 + c2 ∑ xi = c1 ∑ x i + c2 ∑ 1 = i =1

i =1

8

∑ zi · xi ,

i =1 8

∑ zi .

i =1

Po vyčíslení jednotlivých součtů dostaneme (zaokrouhlujeme na 2 desetinná místa) 146.61c1 + 29.30c2 = 94.37, 29.30c1 + 8.00c2 = 22.25. 57

KAPITOLA 4. INTERPOLACE A APROXIMACE FUNKCÍ Soustavu vyřešíme a dostaneme koeficienty c1 = 0.33, c2 = 1.58. Řešením pomocné úlohy je tedy funkce χ( x ) = 0.33x + 1.58. Pro nalezení řešení původní úlohy se vrátíme zpět k původním neznámým a, b, které zpětným dosazením snadno spočteme c1 = b ⇒ b = c1 = 0.33, c2 = ln a ⇒ a = ec2 = e1.58 = 4.87. Výsledná funkce má tvar ϕ( x ) = 4.87e0.33 x .

Poznámka: Výsledný postup je trikem i v tom smyslu, že nám obecně nedává nejlepší aproximaci původní úlohy ve smyslu nejmenších čtverců. Tu spočíst by bylo velmi obtížné. Nalezená funkce je však rozumným a dostupným odhadem takové aproximace. 60


50

40

30

20

10

1

2

3

4

58

5

6

KAPITOLA

5 NUMERICKÉ INTEGROVÁNÍ A DERIVOVÁNÍ

Příklad 5.1: Vypočtěte integrál

Z 3

−1

2

e x dx

složenou lichoběžníkovou formulí pro n = 8. Spočítáme si h, tj. velikost dílku při dělení intervalu h a, bi na 8 částí 3 − (−1) b−a = = 0.5 . n 8 Do tabulky si zapíšeme body dělení x0 , . . . , x8 a hodnotu integrované funkce v těchto bodech. h=

i

xi

0 1 2 3 4 5 6 7 8

-1 -0.5 0 0.5 1 1.5 2 2.5 3

y i = e xi

2

2.7183 1.284 1 1.284 2.7183 9.4877 54.5982 518.0128 8103.0839

59

KAPITOLA 5. NUMERICKÉ INTEGROVÁNÍ A DERIVOVÁNÍ Nakonec spočítáme přibližnou hodnotu integrálu: ! ! 7 y 0 + y n n −1 y0 + y8 + ∑ yi = 0.5 + ∑ yi = I0.5 = h 2 2 i =1 i =1 2.7183 + 8103.0839 = 0.5 + 1.284 + 1 + 1.284 + 2.7183 + 9.4877+ 2 8105.8022 + 588.385 = 2320.64305. +54.5982 + 518.0128) = 0.5 2 Výsledek je I0.5 = 2320.64305.


Z 3

−1

2

e x dx

složenou lichoběžníkovou formulí se zadanou přesností ε = 10−3 . Počítáme integrál složenou lichoběžníkovou formulí pro n = 2, 4, 8, 16, . . . (tj. násobky čísla 2), dokud je chyba větší než zadaná přesnost, tj. | Ih − I2h | > ε.

Výsledek:

n

h

Ih

| Ih − I2h |

2 4 8 16 32 64 128 256 512 1024 2048 4096 8192 16384

2 1

8111.2388 4111.2176 2320.6431 1688.9283 1508.6402 1461.7928 1449.9621 1446.9969 1446.2551 1446.0696 1446.0232 1446.0116 1446.0087 1446.0080

— 4000.0212 1790.5745 631.7148 180.2881 46.8474 11.8307 2.9652 0.7418 0.1855 0.0464 0.0116 0.0029 0.0007

1 2 1 4 1 8 1 16 1 32 1 64 1 128 1 256 1 512 1 1024 1 2048 1 4096

Z 3

−1

2

ε = 10−3 = 0.001

> 0.001 > 0.001 > 0.001 > 0.001 > 0.001 > 0.001 > 0.001 > 0.001 > 0.001 > 0.001 > 0.001 > 0.001 ≤ 0.001

e x dx = 1446.008 ± 0.001.

60

KAPITOLA 5. NUMERICKÉ INTEGROVÁNÍ A DERIVOVÁNÍ Příklad 5.3: Vypočtěte integrál Z 2 0

arctg2 ( x ) + x + 1 dx

složenou Simpsonovou formulí pro n = 16. Spočítáme si h velikost dílků dělení intervalu h a, bi, h=

b−a 2−0 = = 0.125 . n 16

Do tabulky si zapíšeme body dělení x0 , . . . , x16 a hodnotu integrované funkce v těchto bodech.

i

xi

yi = arctg2 ( xi ) + xi + 1

i

xi

yi = arctg2 ( xi ) + xi + 1

0 1 2 3 4 5 6 7 8

0 0.1250 0.2500 0.3750 0.5000 0.6250 0.7500 0.8750 1

1 1.1405 1.3100 1.5037 1.7150 1.9370 2.1641 2.3917 2.6169

9 10 11 12 13 14 15 16

1.1250 1.2500 1.3750 1.5000 1.6250 1.7500 1.8750 2

2.8376 3.0529 3.2624 3.4659 3.6636 3.8560 4.0432 4.2258

Nakonec spočítáme přibližnou hodnotu integrálu: I0.125 =

0.125 h ( y0 + y n + 4 · S L + 2 · SS ) = (y0 + y16 + 4 · (y1 + y3 + y5 + 3 3 +y7 + y9 + y11 + y13 + y15 ) + 2(y2 + y4 + y6 + y8 + y10 + y12 + y14 ) =

= 0.0417 · (5.2258 + 4 · 20.7797 + 2 · 18.1808) = 5.2002.

Výsledek: I0.125 = 5.2002.


Z 5 0

arctg2 ( x ) − 2x dx

složenou Simpsonovou formulí se zadanou přesností ε = 10−8 .

61

KAPITOLA 5. NUMERICKÉ INTEGROVÁNÍ A DERIVOVÁNÍ

Počítáme integrál složenou Simpsonovou formulí pro n = 2, 4, 8, 16, . . ., dokud je chyba větší než zadaná přesnost, tj. | Ih − I2h | > ε.

Výsledek:

n

h

Ih

| Ih − I2h |

2 4 8 16 32 64 128

2.5 1.25 0.625 0.3125 0.15625 0.078125 0.0390625

-18.7055080634 -18.8342068359 -18.8981149073 -18.9021335763 -18.9021508072 -18.9021508687 -18.9021508726

— 0.1286987725 0.0639080714 0.0040186690 0.0000172309 0.0000000615 0.0000000039

Z 5 0

ε = 10−3 = 0.001

> 10−8 > 10−8 > 10−8 > 10−8 > 10−8 ≤ 10−8

arctg2 ( x ) − 2x dx = −18.90215087 ± 10−8 .

Příklad 5.5: Pomocí dvojného přepočtu s extrapolací vypočtěte přibližnou hodnotu integrálu Z 2 ln(1 + x ) dx 1 1 + cos( x ) s přesností ε = 10−6 . Použijte složenou lichoběžníkovou formuli. Označme přibližnou hodnotu integrálu spočtenou pomocí složené lichoběžníkové formule při délce kroku h jako I (h). Metoda dvojného přepočtu s extrapolací nám pro případ složeného lichoběžníkového pravidla dává následující vzorce pro odhad integrálu a odhad chyby při délce kroku h: I (h) − I (2h) I (h) − I (2h) a E(h) = . 3 3 Pro n = 2, 4, 8, 16, . . . postupně počítáme pro jednotlivé kroky hodnoty I, pomocí nich dále určíme hodnoty E. Ve výpočtu pokračujeme, dokud je odhady chyby větší než zadaná přesnost, tj. | E(h)| > ε. Nakonec spočteme hodnotu I1 , kterou prohlásíme za odhad integrálu. Poznamenejme, že výhodou metody dvojného přepočtu s extrapolací oproti standardnímu postupu je její rychlost. I1 (h) = I (h) +

62

KAPITOLA 5. NUMERICKÉ INTEGROVÁNÍ A DERIVOVÁNÍ n

h

I (h)

| E(h)|

2 4 8 16 32 64 128 256 512

0.5 0.25 0.125 0.0625 0.03125 0.015625 0.0078125 0.00390625 0.001953125

1.01079488 0.96728553 0.95599330 0.95314108 0.95242615 0.95224729 0.95220257 0.95219139 0.95218860

— 0.01450312 0.00376408 0.00095074 0.00023831 0.00005962 0.00001491 0.00000373 0.00000093

ε = 10−6 = 0.000001

> 0.000001 > 0.000001 > 0.000001 > 0.000001 > 0.000001 > 0.000001 > 0.000001 ≤ 0.000001

Nyní spočteme závěrečnou zpřesněnou aproximaci: I1 (h) = 0.95218860 +

Výsledek:

0.95218860 − 0.95219139 = 0.95218767. 3

Z 2 ln(1 + x ) 1

1 + cos x

dx = 0.952188 ± 10−6 .

Příklad 5.6: Pomocí dvojného přepočtu s extrapolací vypočtěte přibližnou hodnotu integrálu Z 0.8 3 x +x dx −1 cos2 x s přesností ε = 10−7 . Použijte složenou Simpsonovu formuli.

Označme přibližnou hodnotu integrálu spočtenou pomocí složené Simpsonovy formule při délce kroku h jako I (h). Metoda dvojného přepočtu s extrapolací nám pro případ složeného lichoběžníkového pravidla dává následující extrapolační vzorce pro odhad integrálu a odhad chyby při délce kroku h: I (h) − I (2h) I (h) − I (2h) a E(h) = . 15 15 Pro n = 2, 4, 8, 16, . . . postupně počítáme pro jednotlivé kroky hodnoty I, pomocí nich dále určíme hodnoty E. Ve výpočtu pokračujeme, dokud je odhady chyby větší než zadaná přesnost, tj. | E(h)| > ε. Nakonec spočteme hodnotu I1 , kterou prohlásíme za odhad integrálu. Poznamenejme, že výhodou metody dvojného přepočtu s extrapolací oproti standardnímu postupu je její rychlost. Je užitečné uvědomit si, že algoritmus nemusíme začít u n = 2, ale libovolné vyšší mocniny 2. To má za následek nižší počet kroků vedoucí k dosažení potřebné přesnosti. Začněme tedy například s n = 8. I1 (h) = I (h) +

63

KAPITOLA 5. NUMERICKÉ INTEGROVÁNÍ A DERIVOVÁNÍ n

h

I (h)

| E(h)|

8 16 32 64 128 256

0.2250 0.1125 0.0563 0.0281 0.0140 0.0070

-0.89202663 -0.88094879 -0.87998697 -0.87991974 -0.87991541 -0.87991513

— 0.00073852 0.00006412 0.00000448 0.00000029 0.00000002

ε = 10−7 = 0.0000001

> 10−7 > 10−7 > 10−7 > 10−7 ≤ 10−7

Nyní spočteme závěrečnou zpřesněnou aproximaci: I1 (h) = −0.87991513 +

Výsledek:

−0.87991513 − (−0.87991541) = −0.87991511. 15

Z 0.8 3 x +x −1

cos2 x

dx = −0.8799151 ± 10−7 .

Ra x3 + x Poznámka: Protože je integrand f ( x ) = cos 2 x lichou funkcí, jsou integrály typu − a f ( x ) dx nulové pro libovolné a ∈ R. Před samotným numerickým výpočtem jsme mohli využít této vlastnosti (společně s aditivitou integrálu vůči integračním mezím): Z 0.8 3 x +x −1

cos2 x

dx =

Z −0.8 3 x +x −1

cos2 x

dx +

Z 0.8 3 x +x −0.8

cos2 x

dx =

Z 0.8 3 x +x −1

cos2 x

dx + 0 =

Z 0.8 3 x +x −1

cos2 x

dx

a integrovat funkci na menším intervalu. Požadované přesnosti bychom pak dosáhli o několik kroků dříve. Příklad 5.7: Pro funkci f ( x ) = x3 − 5x vypočtěte přibližně její první derivaci na intervalu h1, 2i pomocí vzorce f 0 (x) ≈

f ( x + h) − f ( x − h) 2h

s krokem h = 0.2. Pro přibližný výpočet derivace s krokem h = 0.25 nejdříve určíme uzly xi = ih, i = 0, . . . , 8 a v nich vypočítáme funkční hodnoty f i = f ( xi ), viz prní tři sloupce tabulky: Vzorec pro přibližný výpočet derivace můžeme zapsat ve tvaru f 0 ( xi ) ≈

f i +1 − f i −1 = f [ x i +1 , x i −1 ] , 2h

64

KAPITOLA 5. NUMERICKÉ INTEGROVÁNÍ A DERIVOVÁNÍ který nám umožní spočítat přibližné hodnoty derivace ve všech vnitřních uzlech. 1.43 − 5 · 1.4 − (13 − 5 · 1) = −0.64, 2 · 0.2 1.63 − 5 · 1.6 − (1.23 − 5 · 1.2) = 0.92, 2 · 0.2 1.83 − 5 · 1.8 − (1.43 − 5 · 1.4) = 2.72, 2 · 0.2 23 − 5 · 2 − (1.63 − 5 · 1.6) = 4.76. 2 · 0.2

0

f (1.2) ≈ f 0 (1.4) ≈ f 0 (1.6) ≈ f 0 (1.8) ≈

Přibližné hodnoty derivace zapíšeme do tabulky. xi

1

1.2

1.4

1.6

1.8

2

f 0 ( xi )

–

-0.64

0.92

2.72

4.76

–

Příklad 5.8: Pro funkci

arctg x 1 + x4 vypočtěte přibližně její první derivaci na intervalu h0, 2i pomocí vzorce f ( x ) = sin x2 +

f 0 (x) ≈

f ( x + h) − f ( x − h) 2h

s krokem h = 0.25 a h = 0.1. Výsledky porovnejte s hodnotami přesně vypočítané derivace. Přesná derivace je určena vzorcem 0

2

f ( x ) = 2x cos x +

1+ x 4 1+ x 2

− 4x3 arctg x (1 + x 4 )2

.

Pro přibližný výpočet derivace s krokem h = 0.25 nejdříve určíme uzly xi = ih, i = 0, . . . , 8 a v nich vypočítáme funkční hodnoty f i = f ( xi ), viz prní tři sloupce tabulky:

i=0 i=1 i=2 i=3 i=4 i=5 i=6 i=7 i=8

xi

fi

f [ x i +1 , x i −1 ]

f 0 ( xi )

ei

0.00 0.25 0.50 0.75 1.00 1.25 1.50 1.75 2.00

0.0000 0.3065 0.6838 1.0221 1.2342 1.2603 0.9402 0.1803 −0.6917

1.3676 1.4313 1.1008 0.4764 −0.5880 −2.1600 −3.2637 -

1.0000 1.4213 1.5165 1.1284 0.5452 −0.4570 −2.1948 −3.6746 −2.7254

0.0538 0.0852 0.0276 0.0688 0.1310 0.0347 0.4109 -

Vzorec pro přibližný výpočet derivace můžeme zapsat ve tvaru f 0 ( xi ) ≈

f i +1 − f i −1 = f [ x i +1 , x i −1 ] , 2h 65

KAPITOLA 5. NUMERICKÉ INTEGROVÁNÍ A DERIVOVÁNÍ který nám umožní spočítat přibližné hodnoty derivace ve všech vnitřních uzlech, např. f 0 (0.75) ≈

1.2342 − 0.6838 = 1.1008, 2 · 0.25

viz čtvrtý sloupec tabulky. Pátý sloupec obsahuje přesné hodnoty derivace, které jsou v posledním sloupci porovnány s hodnotami přibližnými ei = | f [ xi+1 , xi−1 ] − f 0 ( xi )|. Pro h = 0.1 je celý postup analogický, dostaneme však přesnější hodnoty, jak ukazují následující obrázky. h=0.25 2 1 0 derivace aproximace derivace

−1 −2 −3 −4

0

0.2

0.4

0.6

0.8

1

1.2

1.4

1.6

1.8

2

1.2

1.4

1.6

1.8

2

h=0.1 2 1 0 derivace aproximace derivace

−1 −2 −3 −4

0

0.2

0.4

0.6

0.8

1

Příklad 5.9: Pro funkci f ( x ) = sin x2 + arctg x vypočtěte přibližně její druhou derivaci na intervalu h0, 2i pomocí vzorce f 00 ( x ) ≈

f ( x + h) − 2 f ( x ) + f ( x − h) h2

s krokem h = 0.25 a h = 0.1. Výsledky porovnejte s přesnými hodnotami. Přesná druhá derivace je určena vzorcem f 00 ( x ) = 2 cos x2 − 4x2 sin x2 −

2x . (1 + x 2 )2

Pro přibližný výpočet s krokem h = 0.25 nejdříve určíme uzly xi = ih, i = 0, . . . , 8 a v nich vypočítáme funkční hodnoty f i = f ( xi ), viz prní tři sloupce tabulky: 66

KAPITOLA 5. NUMERICKÉ INTEGROVÁNÍ A DERIVOVÁNÍ

i=0 i=1 i=2 i=3 i=4 i=5 i=6 i=7 i=8

xi

fi

2 f [ x i +1 , x i , x i −1 ]

f 00 ( xi )

ei

0.00 0.25 0.50 0.75 1.00 1.25 1.50 1.75 2.00

0.0000 0.3074 0.7111 1.1768 1.6269 1.8960 1.7609 1.1307 0.3503

1.5388 0.9942 −0.2510 −2.8946 −6.4689 −7.9208 −2.4017 -

2.0000 1.5376 1.0504 −0.1225 −2.7853 −6.6139 −8.5430 −3.1737 10.6416

0.0012 0.0562 0.1285 0.1093 0.1450 0.6222 0.7720 -

Vzorec pro přibližný výpočet druhé derivace můžeme zapsat ve tvaru f 00 ( xi ) ≈

f i +1 − 2 f i + f i −1 = 2 f [ x i +1 , x i , x i −1 ] , h2

který nám umožní spočítat přibližné hodnoty derivace ve všech vnitřních uzlech, např. f 00 (0.75) ≈

1.626869 − 2 · 1.176804 + 0.711052 = −0.2510, 0.252

viz čtvrtý sloupec tabulky. Pátý sloupec obsahuje přesné hodnoty druhé derivace, které jsou v posledním sloupci porovnány s hodnotami přibližnými ei = |2 f [ xi+1 , xi , xi−1 ] − f 00 ( xi )|. Pro h = 0.1 je celý postup analogický, dostaneme však přesnější hodnoty, jak ukazují následující obrázky. h=0.25 15 druha derivace aproximace druhe derivace

10 5 0 −5 −10

0

0.2

0.4

0.6

0.8

1

1.2

1.4

1.6

1.8

2

1.4

1.6

1.8

2

h=0.1 15 druha derivace aproximace druhe derivace

10 5 0 −5 −10

0

0.2

0.4

0.6

0.8

1

67

1.2

KAPITOLA

6 OBYČEJNÉ DIFERENCIÁLNÍ ROVNICE: POČÁTEČNÍ ÚLOHY

Příklad 6.1: Řešte diferenciální rovnici s počáteční podmínkou y 0 = ( x + y )2 ,

y(0) = 0,

Eulerovou metodou na intervalu h0, 1.5i s předepsaným krokem h = 0.25. Nejprve si spočítáme počet dílků dělení intervalu h a, bi n=

b−a 1.5 − 0 = = 6. h 0.25

Z počáteční podmínky víme, že x0 = 0, , y0 = 0. Spočítáme x1 , y1 Eulerovou metodou x1 = x0 + h = 0 + 0.25 = 0.25, y1 = y0 + h( x0 + y0 )2 = 0 + 0.25(0 + 0)2 = 0. Dále spočítáme x2 , y2 , x3 , y3 x2 = x1 + h = 0.25 + 0.25 = 0.5, y2 = y1 + h( x1 + y1 )2 = 0 + 0.25(0.25 + 0)2 = 0.0156, x3 = x2 + h = 0.5 + 0.25 = 0.75, y3 = y2 + h( x2 + y2 )2 = 0.0156 + 0.25(0.5 + 0.0156)2 = 0.0821. Další přibližné hodnoty řešení x4 , y4 , . . . , x6 , y6 počítáme obdobně. 68

KAPITOLA 6. OBYČEJNÉ DIFERENCIÁLNÍ ROVNICE: POČÁTEČNÍ ÚLOHY Nalezli jsme přibližné hodnoty řešení: xi

0

0.25

0.5

0.75

1

1.25

1.5

yi

0

0

0.0156

0.0821

0.2552

0.6491

1.5507

1.6

1.4

1.2

1

0.8

0.6

0.4

0.2

0 0

0.5

1

1.5

Příklad 6.2: Řešte diferenciální rovnici s počáteční podmínkou y0 =

1 − 0.1 y, x2 + 1

y(−2) = −1,

Eulerovou metodou na intervalu h−2, 3i s předepsaným krokem h = 0.5. Nejprve si spočítáme počet dílků dělení intervalu h a, bi n=

b−a 3 − (−2) = = 10 . h 0.5

Z počáteční podmínky víme, že x0 = −2, , y0 = −1. Spočítáme x1 , y1 Eulerovou metodou: x1 = x0 + h = −2 + 0.5 = −1.5, ! 1 1 y1 = y0 + h − 0.1 y0 = −1 + 0.5 − 0.1 (−1) = −0.85. (−2)2 + 1 x02 + 1

69

KAPITOLA 6. OBYČEJNÉ DIFERENCIÁLNÍ ROVNICE: POČÁTEČNÍ ÚLOHY Dále spočítáme x2 , y2 , x3 , y3 x2 = x1 + h = −1.5 + 0.5 = −1, ! 1 1 y2 = y1 + h − 0.1 (−0.85) = − 0.1 y1 = −0.85 + 0.5 (−1.5)2 + 1 x12 + 1 x3 y3

= −0.6537, = x2 + h = −1 + 0.5 = −0.5, ! 1 1 = y1 + h − 0.1 y2 = −0.6537 + 0.5 − 0.1 (−0.6537) = (−1)2 + 1 x22 + 1 = −0.371.

Další přibližné hodnoty řešení x4 , y4 , . . . , x10 , y10 počítáme obdobně. Nalezli jsme přibližné hodnoty řešení: xi

-2

-1.5

-1

-0.5

0

yi

-1

-0.85

-0.6537

-0.371

0.0476 0.5

1

1.5

2

2.5

3

0.5452

0.9179

1.122

1.2198

1.2588

1.2648

1

1.5

2

2.5

0.8

0.6

0.4

0.2

0

−0.2

−0.4

−0.6

−0.8

−1 −1

−0.5

0

0.5

70

3

KAPITOLA 6. OBYČEJNÉ DIFERENCIÁLNÍ ROVNICE: POČÁTEČNÍ ÚLOHY

Příklad 6.3: Řešte diferenciální rovnici s počáteční podmínkou y 0 = ( x + y )2 ,

y(0) = 0,

metodou Runge–Kutta 4. řádu na intervalu h0, 1.5i s předepsaným krokem h = 0.25. Nejprve si spočítáme počet dílků dělení intervalu h a, bi n=

b−a 1.5 − 0 = = 6. h 0.25

Z počáteční podmínky víme, že x0 = 0, , y0 = 0. Spočítáme x1 , y1 Runge–Kutta 4. řádu k1 = h · f ( x0 , y0 ) = h · ( x0 + y0 )2 = 0.25 · (0 + 0)2 = 0, h k1 h k1 2 k 2 = h · f x0 + , y0 + = h · x0 + + y0 + = 2 2 2 2 0 2 0.25 +0+ = 0.0039, = 0.25 · 0 + 2 2 h k2 h k2 2 = k 3 = h · f x0 + , y0 + = h · x0 + + y0 + 2 2 2 2 0.25 0.0039 2 = 0.25 · 0 + +0+ = 0.0040, 2 2 k4 = h · f ( x0 + h, y0 + k3 ) = h · ( x0 + h + y0 + k3 )2 =

x1 y1

= 0.25 · (0 + 0.25 + 0 + 0.0040)2 = 0.0161, = x0 + h = 0 + 0.25 = 0.25 1 = y0 + · ( k 1 + 2 · k 2 + 2 · k 3 + k 4 ) = 6 1 = 0 + · (0 + 2 · 0.0039 + 2 · 0.0040 + 0.0161) = 0.0053. 6

Dále spočítáme x2 , y2 k1 = h · f ( x1 , y1 ) = h · ( x1 + y1 )2 = 0.25 · (0.25 + 0.0053)2 = 0.0163, h k1 h k1 2 k 2 = h · f x1 + , y1 + = h · x1 + + y1 + = 2 2 2 2 0.25 0.0163 2 = 0.25 · 0.25 + + 00053 + = 0.0377, 2 2 h k2 h k2 2 k 3 = h · f x1 + , y1 + = h · x1 + + y1 + = 2 2 2 2 0.25 0.0377 2 = 0.25 · 0.25 + + 0.0053 + = 0.0398, 2 2

71

KAPITOLA 6. OBYČEJNÉ DIFERENCIÁLNÍ ROVNICE: POČÁTEČNÍ ÚLOHY k4 = h · f ( x1 + h, y1 + k3 ) = h · ( x1 + h + y1 + k3 )2 =

= 0.25 · (0.25 + 0.25 + 0.0053 + 0.0398)2 = 0.0743, = x1 + h = 0.25 + 0.25 = 0.5, 1 = y1 + · ( k 1 + 2 · k 2 + 2 · k 3 + k 4 ) = 6 1 = 0.0053 + · (0.0163 + 2 · 0.0377 + 2 · 0.0398 + 0.0743) = 0.0463. 6

x2 y2


0

0.25

0.5

0.75

1

1.25

1.5

yi

0

0.0053

0.0463

0.1816

0.5572

1.7534

10.6766

12

10

8

6

4

2

0 0

0.5

1

1.5

Příklad 6.4: Řešte diferenciální rovnici s počáteční podmínkou y0 = sin( x ) · y − x3 + y,

y(1) = 0,

Rungeovou–Kuttovou metodou na intervalu h1, 2i s předepsaným krokem h = 0.2. Nejprve si spočítáme počet dílků dělení intervalu h a, bi n=

b−a 2−1 = = 5. h 0.2

Z počáteční podmínky víme, že x0 = 1, , y0 = 0. 72

KAPITOLA 6. OBYČEJNÉ DIFERENCIÁLNÍ ROVNICE: POČÁTEČNÍ ÚLOHY Spočítáme x1 , y1 Runge–Kutta 4. řádu k1 = h · f ( x0 , y0 ) = h · sin ( x0 ) · y0 − x03 + y0 = 0.2 · sin (1) · 0 − 13 + 0 = k2

k3

k4

= −0.2, h k1 = h · f x0 + , y0 + = 2 2 ! h k1 h 3 k1 = h · sin x0 + · y0 + − x0 + = + y0 + 2 2 2 2 ! 0.2 −0.2 0.2 3 −0.2 = 0.2 · sin 1 + · 0+ − 1+ = −0.304, +0+ 2 2 2 2 h k2 = h · f x0 + , y0 + = 2 2 ! h k2 h 3 k2 = = h · sin x0 + · y0 + − x0 + + y0 + 2 2 2 2 ! 0.2 −0.304 0.2 3 −0.304 = 0.2 · sin 1 + · 0+ − 1+ +0+ = 2 2 2 2 = −0.3237, = h · f ( x0 + h, y0 + k3 ) = 3 = h · sin ( x0 + h) · (y0 + k3 ) − ( x0 + h) + y0 + k3 = = 0.2 · sin (1 + 0.2) · (0 − 0.3237) − (1 + 0.2)3 + 0 − 0.3237 = −0.4707,

x1 = x0 + h = 1 + 0.2 = 1.2, 1 y1 = y0 + · ( k 1 + 2 · k 2 + 2 · k 3 + k 4 ) = 6 1 = 0 + · (−0.2 + 2 · (−0.304) + 2 · (−0.3237) − 0.4707) = −0.321. 6 Dále spočítáme x2 , y2

x13

k1 = h · f ( x1 , y1 ) = h · sin ( x1 ) · y1 − + y1 = = 0.2 · sin (1.2) · (−0.321) − 13 − 0.321 = k2

= −0.4696, h k1 = = h · f x1 + , y1 + 2 2 ! h k1 h 3 k1 = = h · sin x1 + · y1 + − x1 + + y1 + 2 2 2 2 0.2 −0.4696 0.2 3 = 0.2 · sin 1.2 + · −0.321 + − 1.2 + − 0.321+ 2 2 2 −0.4696 + = −0.6577, 2 73

KAPITOLA 6. OBYČEJNÉ DIFERENCIÁLNÍ ROVNICE: POČÁTEČNÍ ÚLOHY k3 =

= =

k4 =

=

h k2 h · f x1 + , y1 + = 2 2 ! h k2 h 3 k2 h · sin x1 + · y1 + − x1 + + y1 + = 2 2 2 2 −0.6577 0.2 3 0.2 · −0.321 + − 1+ 0.2 · sin 1.2 + + 0+ 2 2 2 −0.6577 + = −0.6946, 2 h · f ( x1 + h, y1 + k3 ) = h · sin ( x1 + h) · (y1 + k3 ) − ( x1 + h)3 + y1 + k3 = 3 0.2 · sin (1 + 0.2) · (0 − 0.3237) − (1 + 0.2) + 0 − 0.3237 = −0.9521,

x2 = x1 + h = 1.2 + 0.2 = 1.4, 1 y2 = y1 + · ( k 1 + 2 · k 2 + 2 · k 3 + k 4 ) = 6 1 = −0.321 + · (−0.4696 + 2 · (−0.6577) + 2 · (−0.6946) − 0.9521) = −1.0087. 6 Další přibližné hodnoty řešení x3 , y3 , . . . , x5 , y5 počítáme obdobně. Nalezli jsme přibližné hodnoty řešení: xi

1

1.2

1.4

1.6

1.8

2

yi

0

-0.3210

-1.0087

-2.3257

-4.6586

-8.5351

0

−1

−2

−3

−4

−5

−6

−7

−8

−9 1

1.1

1.2

1.3

1.4

1.5

74

1.6

1.7

1.8

1.9

2

KAPITOLA 6. OBYČEJNÉ DIFERENCIÁLNÍ ROVNICE: POČÁTEČNÍ ÚLOHY

Příklad 6.5: Řešte diferenciální rovnici s počáteční podmínkou y0 =

y2 ln( x ) − y , x

y(1) = 1,

Eulerovou metodou i Rungeovou–Kuttovou metodou na intervalu h1, 2i s předepsaným krokem h = 0.2. Řešení porovnejte. Nejprve spočítáme diferenciální rovnice Eulerovou metodou. Spočítáme počet dílků dělení intervalu h a, bi 2−1 b−a = = 5. n= h 0.2 Z počáteční podmínky víme, že x0 = 1, , y0 = 1. Spočítáme x1 , y1 Eulerovou metodou x1 = x0 + h = 1 + 0.2 = 1.2, y2 ln x0 − y0 12 ln 1 − 1 y1 = y0 + h 0 = 1 + 0.2 = 0.8. x0 1 Dále spočítáme x2 , y2 , x3 , y3 x2 = x1 + h = 1.2 + 0.2 = 1.4, y2 ln x1 − y1 0.82 ln 1.2 − 0.8 = 0.8 + 0.2 y2 = y1 + h 1 x1 1.2 = 0.6861, x3 = x2 + h = 1.4 + 0.2 = 1.6, y2 ln x2 − y2 0.68612 ln 1.4 − 0.6861 y3 = y2 + h 2 = 0.6861 + 0.2 x2 1.4 = 0.6107,

Další přibližné hodnoty řešení x4 , y4 , x5 , y5 počítáme obdobně. Nalezli jsme přibližné hodnoty řešení: xi

1

1.2

1.4

1.6

1.8

2

yi

1

0.8

0.6861

0.6107

0.5563

0.5147

Nyní vyřešíme diferenciální rovnice Rungeovou–Kuttovou metodou. Nejprve si spočítáme počet dílků dělení intervalu h a, bi: n=

b−a 2−1 = = 5. h 0.2

Z počáteční podmínky víme, že x0 = 1, , y0 = 1.

75

KAPITOLA 6. OBYČEJNÉ DIFERENCIÁLNÍ ROVNICE: POČÁTEČNÍ ÚLOHY Spočítáme x1 , y1 Runge–Kutta 4. řádu: y20 ln( x0 ) − y0 12 ln(1) − 1 = 0.2 · = −0.2 x0 1 2 h k1 ln x + y + 0 0 2 2 − y0 + h k h · f x0 + , y0 + 1 = h · 2 2 x0 + 2h 2 0.2 1 − 0.2 ln 1 + 0.2 2 2 − 1− 2 0.2 · = −0.1496 1 + 0.2 2 2 k2 h y + ln x + 0 0 2 2 − y0 + h k2 h · f x0 + , y0 + = h· 2 2 x0 + 2h 2 0.1496 0.2 0.1496 1− 2 ln 1 + 2 − 1 − 2 0.2 · = −0.1534 1 + 0.2 2

k 1 = h · f ( x0 , y0 ) = h · k2 =

=

k3 =

=

k1 2

k2 2

=

=

(y0 + k3 )2 ln ( x0 + h) − (y0 + k3 ) = x0 + h (1 − 0.1534)2 ln (1 + 0.2) − (1 − 0.1534) 0.2 · = −0.1193 1 + 0.2 x0 + h = 1 + 0.2 = 1.2 1 y0 + · ( k 1 + 2 · k 2 + 2 · k 3 + k 4 ) = 6 1 1 + · (−0.2 + 2 · (−0.1496) + 2 · (−0.1534) − 0.1193) = 0.8458 6

k4 = h · f ( x0 + h, y0 + k3 ) = h ·

= x1 = y1 =

=

Dále spočítáme x2 , y2 : y21 ln( x1 ) − y1 0.84582 ln(1.2) − 0.8458 = −0.1192 = 0.2 · x1 1.2 2 k1 k1 h y + ln x + − y + 1 1 1 2 2 2 h k h · f x1 + , y1 + 1 = h · = h 2 2 x1 + 2 2 0.1192 0.2 0.1192 0.8458 − 2 ln 1.2 + 2 − 0.8458 − 2 0.2 · = −0.0960 1.2 + 0.2 2 2 k2 k2 h y + ln x + − y + 1 1 1 2 2 2 h k2 h · f x1 + , y1 + = h· = 2 2 x1 + 2h 2 0.0960 0.8458 − 0.0960 ln 1.2 + 0.2 2 2 − 0.8458 − 2 0.2 · = −0.0970 1.2 + 0.2 2

k 1 = h · f ( x1 , y1 ) = h ·

k2 =

=

k3 =

=

76

(y1 + k3 )2 ln ( x1 + h) − (y1 + k3 ) h · f ( x1 + h, y1 + k3 ) = h · = x1 + h (0.8458 − 0.0970)2 ln (1.2 + 0.2) − (0.8458 − 0.0970) 0.2 · = −0.08 1.2 + 0.2 x1 + h = 1.2 + 0.2 = 1.4 1 y1 + · ( k 1 + 2 · k 2 + 2 · k 3 + k 4 ) = 6 1 0.8458 + · (−0.1192 + 2 · (−0.0960) + 2 · (−0.0970) − 0.08) = 0.7482 6

k4 =

= x2 = y2 =

=


1

1.2

1.4

1.6

1.8

2

yi

1

0.8458

0.7482

0.6803

0.6298

0.5906

Nalezli jsme přibližné hodnoty řešení: xi

1

1.2

1.4

1.6

1.8

2

Euler yi

1

0.8

0.6861

0.6107

0.5563

0.5147

RK yi

1

0.8458

0.7482

0.6803

0.6298

0.5906

1 Euler RK

0.9

0.8

0.7

0.6

0.5

0.4 1

1.1

1.2

1.3

1.4

1.5

77

1.6

1.7

1.8

1.9

2

LITERATURA

Základní doporučená literatura [1] GILAT, Amos a Vish SUBRAMANIAM. Numerical methods for engineers and scientists: an introduction with applications using MATLAB. 2nd ed. Hoboken, N.J.: Wiley, 2011, xvi, 495 p. ISBN 9780470565155. [2] KUČERA, Radek. Numerické metody. 1. vyd. Ostrava: VŠB - Technická univerzita, 2006, 132 s. ISBN 80-248-1198-7. Dostupné z: http://mdg.vsb.cz/portal/nm/nm.pdf. [3] WOODFORD, Chris a Chris PHILLIPS. Numerical methods with worked examples: Matlab edition. 2nd ed. New York: Springer, 2012, x, 256 p. ISBN 9789400713666.

Doplňková literatura [4] FIEDLER, Miroslav. Speciální matice a jejich použití v numerické matematice. 1. vyd. Praha: Státní nakladatelství technické literatury, 1981, 266 s. Teoretická knižnice inženýra. [5] HAMMING, R. Numerical methods for scientists and engineers. 2nd ed. New York: Dover, 1973, ix, 721 p. ISBN 0486652416. [6] MÍKA, Stanislav. Numerické metody algebry. 2., nezm. vyd. Praha: SNTL, 1985, 169 s. Matematika pro vysoké školy technické. [7] PŘIKRYL, Petr. Numerické metody matematické analýzy: vysokoškolská příručka pro vysoké školy technické. 2. nezm. vyd. Praha: Státní nakladatelství technické literatury, 1988, 187 s. Matematika pro vysoké školy technické. [8] RALSTON, Anthony. Základy numerické matematiky: příručka pro univerzity ČSR. 2. čes. vyd. Praha: Academia, 1978, 635, [1] s. 78

[9] SÜLI, Endre a David MAYERS. An introduction to numerical analysis. New York: Cambridge University Press, 2003, x, 433 p. ISBN 0521007941. ˇ [10] TEODORESCU, P, Nicolae-Doru STANESCU a Nicolae PANDREA. Numerical analysis with applications in mechanics and engineering. Hoboken, NJ: Wiley, 2013, xi, 633 pages. ISBN 9781118077504. [11] WALLISCH, Pascal. MATLAB for neuroscientists: an introduction to scientific computing in MATLAB. Burlington: Elsevier, 2009, xiv, 384 s., [8] s. obr. příl. ISBN 978-0-12-374551-4. [12] VITÁSEK, Emil. Numerické metody. Praha: SNTL, 1987.

79

Název: Numerická matematika: Banka řešených příkladů Katedra: Katedra matematiky a deskriptivní geometrie Autoři: Radek Kučera, Pavel Ludvík, Zuzana Morávková Místo, rok, vydání: Ostrava, 2016, 1. vydání Počet stran: 79 Vydala: Vysoká škola báňská-Technická univerzita Ostrava Neprodejné ISBN 978-80-248-3894-6

Numerická matematika Banka řešených příkladů

Recommend Documents