Aplikace pro numerické ˇreˇsen´ı matematických úloh

(1)

Aplikace pro numerick´ e ˇ reˇ sen´ı matematick´ ych ´ uloh

Diplomov´ a pr´ ace

Studijn´ı program: N2612 – Elektrotechnika a informatika Studijn´ı obor: 1802T007 – Informaˇcn´ı technologie Autor pr´ace: Bc. Zdenˇek Kybl

Vedouc´ı pr´ace: RNDr. Dana ˇCern´a, Ph.D.

(2)

Applications for the numerical solution of mathematical problems

Diploma thesis

Study programme: N2612 – Electrotechnology and informatics Study branch: 1802T007 – Information technology

Author: Bc. Zdenˇek Kybl

Supervisor: RNDr. Dana ˇCern´a, Ph.D.

(3)

(4)

(5)

(6)

Abstrakt

Diplomová práce si klade za c´ıl vytvoˇren´ı reˇserˇse vybraných nu- merických metod a zhotoven´ı aplikace, jeˇz slouˇz´ı zejména jako di- daktická pom˚ucka student˚um pˇri studiu problematiky numerické matematiky. Teoretická ˇcást je rozdˇelena do ˇsesti kapitol, pˇriˇcemˇz v kaˇzdé kapitole jsou charakterizovány hlavn´ı principy numerických metod jednoho odvˇetv´ı numerické matematiky. Práce postupnˇe seznamuje s algoritmy zabývaj´ıc´ımi se aproximac´ı a interpolac´ı funkce, numerickou integrac´ı a derivac´ı, ˇreˇsen´ım nelineárn´ıch rovnic, metodami pro ˇreˇsen´ı soustav lineárn´ıch rovnic a s algoritmy slouˇz´ıc´ımi pro výpoˇcet vlastn´ıch ˇc´ısel a vlastn´ıch vektor˚u reálných symetrických matic. V praktické ˇcásti je nejprve pˇredstavena aplikace implementovaná v jazyce C# z pohledu jej´ıho návrhu, kdy jsou bl´ıˇze pˇredstaveny vˇsechny vrstvy aplikace. Pozdˇeji je ilustrována interakce jednotlivých vrstev aplikace a dˇen´ı v pozad´ı aplikace pˇri jej´ım uˇz´ıván´ı uˇzivatelem. Pro zajiˇstˇen´ı dostateˇcné didaktické úrovnˇe vyuˇz´ıvá aplikace nástroj˚u, pomoc´ı nichˇz docház´ı k zobrazen´ı nejen z´ıskaného ˇreˇsen´ı, ale i postupu, který vedl k jeho dosaˇzen´ı. V pˇr´ıpadˇe aproximace a interpolace funkce, ˇreˇsen´ı nelineárn´ıch rovnic a nume- rické integrace je didaktická úroveˇn umocnˇena grafickou interpre- taci úlohy a jej´ıho ˇreˇsen´ı. Aplikace dále obsahuje sadu cviˇcných

´

uloh a podporuje exporty do dalˇs´ıch formát˚u. Pro distribuci aplikace byly zhotoveny webové stránky.

Kl´ıˇ cov´ a slova

Numerická matematika, aproximace funkc´ı, numerická integrace, numerická derivace, nelineárn´ı rovnice, soustavy lineárn´ıch rovnic, vlastn´ı ˇc´ısla a vektory symetrických matic, didaktická aplikace, postup výpoˇctu, cviˇcné úlohy, export

(7)

Abstract

The aim of the thesis is a review of numerical methods and ma- nufacturing applications, which are mainly used as a didactic aid for students studying the problems of numerical mathematics. The theoretical part is divided into six chapters, where each chapter outlines the main principles of numerical methods, one branch of numerical mathematics. Work gradually introduces algorithms dealing with approximations and interpolation functions, numerical integration and differentiation, solution of nonlinear equations, methods for solving systems of linear equations and algorithms ser- ving for calculating eigenvalues and eigenvectors of real symmetric matrices. The practical part is firstly introduced by application im- plemented in language C# in terms of its design, which introduces each application layer in more details. Then the interaction of the layers of the application during its use is illustrated. To ensure suffi- cient levels of educational uses the application uses tools for viewing not only the results, but also the process leading to their achie- vement. In the case of approximation and interpolation functions, solving nonlinear equations and numerical integration a didactic le- vel is enhanced by graphical interpretation of the examples and its solutions. The application also includes a set of training tasks and supports exports to other formats. The websites were made for dis- tribution of the application.

Keywords

Numerical methods, approximation of function, numerical integration, numerical differentiation, nonlinear equations, linear equations, eigenvalues and eigenvectors of symmetric matrices, didactic applications, calculation procedure, practice tasks, export

(8)

Podˇ ekov´ an´ı

Rád bych podˇekoval vedouc´ı mé práce RNDr. Danˇe ˇCerné, Ph.D.

za odborné veden´ı, cenné rady, trpˇelivost a ochotu, kterou mi v pr˚ubˇehu zpracován´ı diplomové práce vˇenovala.

(9)

Obsah

Seznam zkratek . . . 10

Seznam obr´azk˚u. . . 11

Seznam tabulek . . . 11

Uvod´ 13 1 Aproximace funkce 14 1.1 Aproximace interpolaˇcn´ım polynomem . . . 14

1.2 Interpolace spline funkcemi. . . 17

1.3 Metoda nejmenˇs´ıch ˇctverc˚u. . . 21

2 Numerick´a derivace 23 2.1 Derivace pomoc´ı interpolace . . . 23

2.2 Richardsonova extrapolace . . . 25

3 Numerick´a integrace funkc´ı 27 3.1 Newton-Cotesovy vzorce . . . 27

3.2 Metoda poloviˇcn´ıho kroku . . . 34

3.3 Rombergova kvadratura . . . 35

3.4 Gaussova kvadratura . . . 36

4 Neline´arn´ı rovnice 42 4.1 Metoda p˚ulen´ı intervalu . . . 42

4.2 Metoda regula falsi . . . 43

4.3 Metoda seˇcen . . . 43

4.4 Newtonova metoda . . . 44

4.5 Steffensenova metoda . . . 44

4.6 Halleyova metoda . . . 44

4.7 Sturmova posloupnost . . . 45

5 Line´arn´ı rovnice 46 5.1 Z´akladn´ı pojmy . . . 46

5.2 Pˇr´ım´e metody . . . 48

5.3 Iteraˇcn´ı metody . . . 51

5.4 Soustavy s obd´eln´ıkovou matic´ı . . . 54

(10)

6 Vlastn´ı ˇc´ısla a vlastn´ı vektory matic 56

6.1 Z´akladn´ı pojmy . . . 56

6.2 C´ˇasteˇcn´y probl´em vlastn´ıch ˇc´ısel . . . 59

6.3 Upln´´ y probl´em vlastn´ıch ˇc´ısel . . . 63

7 Implementace aplikace 72 7.1 Implementaˇcn´ı vrstva aplikace . . . 72

7.2 Kontroln´ı a v´ypoˇcetn´ı vrstva . . . 80

7.3 Prezentaˇcn´ı vrstva aplikace . . . 84

7.4 Ilustrace interakce vrstev aplikace . . . 90

7.5 Distribuce a instalace aplikace . . . 96

Z´avˇer 98

Literatura 100

A Uk´azky aplikace 105

B Obsah DVD 116

(11)

Seznam zkratek

N Mnoˇzina pˇrirozených ˇc´ısel R Mnoˇzina reálných ˇc´ısel C Mnoˇzina komplexn´ıch ˇc´ısel f (x) Hodnota funkce f v bodˇe x

f⁰(x) Hodnota prvn´ı derivace funkce f v bodˇe x

b

R

a

f (x)dx Integr´al funkce f na intervalu ha, bi lim_x→0f Limita funkce f pro x jdouc´ı k nule δ_ij Kroneckerovo delta

∀ Pro vˇsechna

∃ Existuje

∈ Je prvkem

Cⁿ⁺¹(I) Prostor (n + 1) spojit´ych derivac´ı na intervalu I Π_n Prostor polynom˚u stupnˇe nejv´yˇse n

Π˜_n Prostor normovaných polynom˚u stupnˇe nejvýˇse n A^T Transponovaná matice k matici A

A⁻¹ Inverzn´ı matice k matici A

A⁺ Pseudoinverzn´ı matice k matici A h(A) hodnost matice A

I Jednotková matice det(A) Determinant matice A λ Vlastn´ı ˇc´ıslo matice ha, bi Uzavˇrený interval (a, b) Otevˇrený interval

ρ(A) Spektráln´ı polomˇer matice A σi i-té singulárn´ı ˇc´ıslo

kxk Euklidovsk´a norma vektoru x

|a| Absolutn´ı hodnota a

AJAX Asynchronous JavaScript and XML BMP Bit Mapped Picture

GIF Graphic Interchange Format GUI Graphic User Interface

HTML HyperText Markup Language JPG Joint Photographic Experts Group LIFO Last In First Out

MathML Mathematical Markup Language MSIE Microsoft Internet Explorer PDF Portable Document Format PNG Portable Network Graphics RPN Reverzn´ı polsk´a notace TIFF Tag Image File Format XML Extensible Markup Language

XSD XML Schema

(12)

Seznam obr´ azk˚ u

6.1 Odvozen´ı Householderovy matice . . . 65

6.2 Odvozen´ı matice rovinn´e rotace . . . 68

7.1 Zjednoduˇsen´y diagram n´avrhu aplikace . . . 79

7.2 Uk´azka pˇrevodu v´yrazu do MathML . . . 87

7.3 Zjednoduˇsený postup výpoˇctu úlohy obdéln´ıkovým pravidlem s následným zobrazen´ı pomoc´ı Awesomia . . . 91

7.4 Diagram znázorˇnuj´ıc´ı naˇcten´ı cviˇcné úlohy . . . 92

7.5 Uk´azka postupu generov´an´ı HTML pro Gaussovu eliminaci a jeho zobrazen´ı pomoc´ı Awesomia . . . 95

7.6 Uk´azka postupu pˇri exportu do PDF . . . 96

(13)

Seznam tabulek

2.1 Richardsonova extrapolace . . . 26 3.1 Uzly a v´ahy Gaussovy kvadratury . . . 41 7.1 Priorita oper´ator˚u v metodˇe RPN. . . 82

(14)

Uvod ´

Skoro kaˇzdý ˇclovˇek modern´ıho svˇeta dennˇe vyuˇz´ıvá technických pom˚ucek ˇci vymoˇzenost´ı, aniˇz by si uvˇedomoval, ˇze za jejich objeven´ım stoj´ı velice ˇcasto r˚uzná odvˇetv´ı matematiky. Pˇri vývoji tˇechto modern´ıch výdobytk˚u ˇcasto naráˇz´ıme na velice sloˇzité matematické modely a metody. V pˇredpoˇc´ıtaˇcové éˇre bylo naprosto nepˇredstavitelné provádˇet velké mnoˇzstv´ı poˇcetn´ıch operac´ı, a proto se odborn´ıci aplikované matematiky snaˇzili nalézt ˇreˇsen´ı analytickým zp˚usobem, pomoc´ı nˇehoˇz doˇslo k redukci poˇctu provádˇených operac´ı. Ovˇsem z´ıskán´ı tohoto pˇresného ˇreˇsen´ı je mnohdy nemoˇzné, ˇci velice sloˇzité a pracné. Naˇstˇest´ı se ukazuje, ˇze pro reálné vyuˇzit´ı nám ve vˇetˇsinˇe pˇr´ıpad˚u postaˇc´ı pouze ˇreˇsen´ı pˇribliˇzné. A v tuto chv´ıli pˇricház´ı na ˇradu numerická matematika.

Numerická matematika zaznamenala prudkého rozmachu aˇz s rozvojem poˇc´ıtaˇc˚u, kdy doˇslo k výraznému zlevnˇen´ı poˇcetn´ı operace. Proto mohly pˇrej´ıt do popˇred´ı výpoˇcty ˇcasto cyklického charakteru, v nichˇz m˚uˇzeme uvaˇzovat dˇr´ıve nemyslitelné poˇcty operac´ı.

I pˇres nesporné vyuˇzit´ı numerické matematiky a matematiky obecnˇe pˇri ˇreˇsen´ı reálných problém˚u (jmenujme napˇr´ıklad odvˇetv´ı stroj´ırenstv´ı, teorie obvod˚u ˇci sta- vitelstv´ı) nen´ı o studium této problematiku pˇr´ıliˇs velký zájem. Z tohoto d˚uvodu si pˇredkládaná diplomová práce klade nejprve za úkol seznámit ˇctenáˇre s vybranými numerickými metodami z oblasti aproximace funkce, numerického výpoˇctu urˇcitého integrálu a derivace, s postupy pro ˇreˇsen´ı nelineárn´ıch rovnic ˇci metodami pro ˇreˇsen´ı soustav lineárn´ıch rovnic a v neposledn´ı ˇradˇe také s algoritmy pro výpoˇcet vlastn´ıch ˇc´ısel a vlastn´ıch vektor˚u reálných symetrických matic.

Na základˇe této teoretické stati si diplomová práce klade za c´ıl vytvoˇren´ı di- daktické aplikace, která na rozd´ıl od jiˇz existuj´ıc´ıch matematických program˚u, bude kromˇe samotného ˇreˇsen´ı uˇzivateli interpretovat i postup výpoˇctu, pomoc´ı nˇehoˇz bylo ˇreˇsen´ı dosaˇzeno. Implementovaná aplikace by mˇela dále obsahovat soubor cviˇcných

´

uloh vˇcetnˇe podpory pro vykreslován´ı graf˚u a export˚u dat. Pro snazˇs´ı distribuci aplikace pˇr´ıpadnému ˇctenáˇri bude v rámci ˇreˇsen´ı diplomové práce zhotovena webová stránka, která bude obsahovat zhotovenou aplikaci.

Pˇredkládaná diplomová práce by tedy mˇela slouˇzit zejména jako studijn´ı opora a souˇcasnˇe jako výuková pom˚ucka pro zájemce zabývaj´ıc´ı se problematikou nume- rické matematiky.

(15)

1 Aproximace funkce

V prvn´ı kapitole se budeme zabývat metodami pro interpolaci a aproximaci funkce. Princip aproximace funkce spoˇc´ıvá v nahrazen´ı jisté funkce f jinou funkc´ı φ, která je v jistém smyslu p˚uvodn´ı funkci podobná.

Aproximuj´ıc´ı funkce φ by mˇela být co moˇzná nejjednoduˇsˇs´ı, a zároveˇn by mˇela být zadaným bod˚um f (x_i) pro i = 0, . . . ., n co nejbl´ıˇze.

Vyuˇzit´ı aproximac´ı funkce je pomˇernˇe r˚uznorodé. Pokud napˇr´ıklad chceme na poˇc´ıtaˇci vypoˇc´ıtat funkˇcn´ı hodnotu jisté funkce, tak výpoˇcet této hodnoty se ˇcasto dˇeje právˇe pomoc´ı aproximac´ı funkc´ı, kdy je vstupn´ı funkce f nahrazena jistým polynomem P , a to zejména z d˚uvodu snadné a hlavnˇe rychlé práce s mnohoˇcleny.

Polynomy jsou totiˇz pomˇernˇe snadno vyˇc´ısliteln´e, jejich derivace i integr´aly jsme schopni snadno a rychle vypoˇc´ıtat.

Dalˇs´ı oblast´ı, které se budeme vˇenovat pozdˇeji a kde m˚uˇzeme vidˇet vyuˇzit´ı aproximac´ı funkce, je výpoˇcet integrálu, kdy opˇet nahrazujeme vstupn´ı funkci f jistým polynomem P .

Nyn´ı vyvstává otázka, jak nalézt funkci, která bude aproximovat p˚uvodn´ı funkci.

Existuje nˇekolik typ˚u metod pro jej´ı urˇcen´ı. Prvn´ım typem, se kter´ym se budeme bl´ıˇze seznamovat, je aproximace interpolaˇcn´ım polynomem.

1.1 Aproximace interpolaˇ cn´ım polynomem

O interpolaci mluv´ıme tehdy, je-li ´ukolem stanovit hodnotu funkce f v bodech leˇz´ıc´ıch mezi dvˇema tabulkov´ymi body.

Snaˇz´ıme se nalézt funkci, která v bodech x₀, x₁, x₂, . . . , x_nnabývá hodnoty f (x₀), f (x₁), f (x₂), . . . , f (x_n). Body x₀, x₁, x₂, . . . , x_n nazýváme uzlové body.

1.1.1 Lagrange˚ uv interpolaˇ cn´ı polynom

Lagrangeova interpolace je aproximace polynomem L_n, pro který plat´ı, ˇze splˇnuje základn´ı úlohu interpolace.

Definice 1.1.1. Uvaˇzujme n + 1 bod˚u, které jsou navzájem r˚uzné. Hledáme polynom L_n stupnˇe nejvýˇse n takový, ˇze plat´ı f (x_i) = L_n(x_i), i = 0, . . . , n. Tento problém budeme oznaˇcovat jako základn´ı úlohu interpolace. Polynom Ln se nazývá Lagrange˚uv interpolaˇcn´ı polynom.

(16)

Vˇeta 1.1.1. Mˇejme dány body [x_i, f (x_i)], i = 0, . . . , n. Pak existuje právˇe jeden interpolaˇcn´ı polynom L_n stupnˇe nejvýˇse n takový, ˇze L_n(x_i) = f (x_i), i = 0, . . . , n.

D˚ukaz. Pˇredpokládejme, ˇze existuj´ı dva polynomy L_na R_nstupnˇe nejvýˇse n takové, ˇze L_n(x_i) = R_n(x_i) = f (x_i) pro i = 0, . . . , n. Ukáˇzeme, ˇze jsou tyto dva polynomy shodné.

Poloˇzme Q_n = L_n− R_n. Potom Q_n(x_i) = L_n(x_i) − R_n(x_i) = 0. Polynom Q_n je polynom stupnˇe n, který má n + 1 nulových bod˚u. Podle základn´ı vˇety algebry je tento polynom identicky roven nule, a tedy L_n≡ R_n.

Hledáme polynom stupnˇe nejvýˇse n , který splˇnuje podm´ınku interpolace. Máme tedy n + 1 bod˚u, kterými mus´ı graf hledaného polynomu procházet. Tuto úlohu m˚uˇzeme ˇreˇsit jako soustavu lineárn´ıch rovnic, kdy bychom zjistili hodnoty koeficient˚u hledaného polynomu (pomoc´ı tzv. Vandermondovy matice). Nicménˇe tento zp˚usob ˇreˇsen´ı se pˇr´ıliˇs nevyuˇz´ıvá a lze se mu vyhnout pomoc´ı Lagrangeova inter- polaˇcn´ıho polynomu.

Vˇeta 1.1.2. Lagrange˚uv interpolaˇcn´ı polynom lze vyj´adˇrit ve tvaru L_n(x) =

n

X

i=0

f (x_i)g_i(x), (1.1)

kde g_i(x) = Q

j=0,j6=i x−xj

xi−x_j.

Vˇeta 1.1.3. Necht’ f ∈ Cⁿ⁺¹(I), kde I je nejmenˇs´ı interval obsahuj´ıc´ı x₀, x₁, x₂, . . . , x_n−1, x_n, x^∗ a x₀, x₁, x₂, . . . , x_n jsou navz´ajem r˚uzn´e uzly.

Pak ∀ x^∗ ∈ I ∃ ξ ∈ I, pro kter´e plat´ı:

f (x^∗) − L_n(x^∗) = f⁽ⁿ⁺¹⁾(ξ)ω_n+1(x^∗)

(n + 1)! , (1.2)

kde ω_n+1(x) = (x − x₀)(x − x₁) . . . (x − x_n).

D˚ukaz. Pˇredpokl´adejme, ˇze xi = x^∗. Potom

f (x_i) − L_n(x_i) = fⁿ⁺¹(ξ)ω_n+1(x_i)

(n + 1)!. (1.3)

Z vlastnosti interpolace plyne

f (x_i) − L_n(x_i) = 0. (1.4)

Pokud x_i 6= x^∗, potom definujme funkci

F (x) = f (x) − L_n(x) − tω_n+1(x), (1.5) kde x ∈ I, t ∈ R.

(17)

Funkce F (x) má n + 1 nulových bod˚u (uzlové body x_i). Hledáme takové t, aby byla splnˇena rovnost

F (x^∗) = 0. (1.6)

To splˇnuje

t = f (x^∗) − Ln(x^∗)

ω_n+1(x^∗) . (1.7)

Funkce F má tedy n + 2 nulových bod˚u. Z Rolleovy vˇety plyne, ˇze F⁰ má alespoˇn n + 1 nulových bod˚u. F⁰⁰ má alespoˇn n nulových bod˚u. F⁽ⁿ⁺¹⁾ má alespoˇn jeden nulový bod ξ,

F⁽ⁿ⁺¹⁾(ξ) = 0. (1.8)

Protoˇze L⁽ⁿ⁺¹⁾n ≡ 0, dostaneme

F⁽ⁿ⁺¹⁾(ξ) = f⁽ⁿ⁺¹⁾(ξ) − 0 − t(n + 1)!. (1.9) Dosad´ıme za promˇennou t

0 = F⁽ⁿ⁺¹⁾(ξ) = f⁽ⁿ⁺¹⁾(ξ) −f (x^∗) − L_n(x^∗)

ω_n+1(x^∗) (n + 1)! (1.10) a vyj´adˇr´ıme v´yslednou chybu

f (x^∗) − Ln(x^∗) = f⁽ⁿ⁺¹⁾(ξ)ω_n+1(x^∗)

(n + 1)! . (1.11)

V´ıce informac´ı nalezneme napˇr´ıklad v [1] nebo [2]

1.1.2 Newton˚ uv interpolaˇ cn´ı polynom

Tento polynom je pro n + 1 uzlových bod˚u interpolaˇcn´ım polynomem stupnˇe nejvýˇse n. Bude se tedy jednat o nový zp˚usob zápisu Lagrangeova interpolaˇcn´ıho polynomu. Tento polynom budeme hledat ve tvaru

N_n(x) = a₀+ a₁(x − x₀) + a₂(x − x₀)(x − x₁) + . . . (1.12) + a_n(x − x₀)(x − x₁) . . . (x − x_n−1).

Newton˚uv interpolaˇcn´ı polynom tedy mus´ı vyhovovat podm´ınce interpolace

N_n(x_i) = f (x_i), i = 0, 1, . . . , n. (1.13) Koeficienty a_i, i = 0, . . . , n, vyj´adˇr´ıme pomoc´ı pomˇern´ych diferenc´ı.

Definice 1.1.2. Pomˇerná diference nultého ˇrádu je definována:

f [x_i] = f (x_i), i = 0, . . . , n. (1.14)

(18)

Pomˇerná diference prvn´ıho ˇrádu je definována:

f [x_i, x_i+1] = f [x_i+1] − f [x_i]

x_i+1− x_i , i = 0, . . . , n − 1. (1.15) Pomˇerná diference k-tého ˇrádu je definována rekurentnˇe:

f [xi, xi+1, . . . , xi+k] = f [x_i+1, . . . , x_i+k] − f [x_i, x_i+1, . . . , x_i+k−1]

x_i+k− x_i . (1.16)

Newton˚uv interpolaˇcn´ı polynom m˚uˇzeme zapsat pomoc´ı pomˇern´ych diferenc´ı:

N_n(x) = f [x₀] + f [x₀, x₁](x − x₀) + f [x₀, x₁, x₂](x − x₀)(x − x₁) + . . .

+ . . . + f [x₀, x₁, . . . , x_n](x − x₀)(x − x₁) . . . (x − x_n−1). (1.17) Newtonova interpolace má oproti Lagrangeovˇe interpolaci podstatnou výhodu, která spoˇc´ıvá v tom, ˇze je výpoˇcetnˇe ménˇe nároˇcné pˇridat dalˇs´ı bod s jeho funkˇcn´ı hodnotou, protoˇze nˇekteré výpoˇcty z˚ustanou beze zmˇeny (napˇr´ıklad pˇredchoz´ı koeficienty a_k se nezmˇen´ı). Bliˇzˇs´ı informace jsou k nalezen´ı napˇr´ıklad v [1].

1.2 Interpolace spline funkcemi

Pro intervaly vˇetˇs´ı délky je pouˇzit´ı interpolaˇcn´ıho polynomu n´ızkého stupnˇe mnohdy nepˇresné a pouˇzit´ı interpolaˇcn´ıho polynomu vyˇsˇs´ıch stupˇn˚u nevhodné vzhle- dem k vlastnosti, kdy interpolaˇcn´ı polynom na kraj´ıch interval˚u nepˇr´ıjemnˇe osciluje.

Vhodnˇejˇs´ı je vyuˇzit´ı interpolaˇcn´ıch spline funkc´ı, kter´e jsou po ˇc´astech polynomy.

Princip spline interpolace spoˇc´ıvá v rozdˇelen´ı intervalu na podintervaly. Na kaˇzdém z tˇechto podinterval˚u poté budeme konstruovat obecnˇe jiný polynom.

Mezi nejˇcastˇeji vyuˇz´ıvané spliny patˇr´ı lineárn´ı a kubický [3].

1.2.1 Line´ arn´ı interpolaˇ cn´ı spline

Definice 1.2.1. (Lineárn´ı interpolaˇcn´ı spline) Lineárn´ım splinem nazýváme funkci φ(x), která je spojitá na intervalu hx₀, x_ni a na kaˇzdém podintervalu hx_i, x_i+1i, i = 0, . . . , n − 1 je polynomem prvn´ıho stupnˇe.

Lineárn´ı interpolaˇcn´ı spline, který procház´ı uzlovými body, tj. φ(x_i) = f (x_i), na podintervalu hx_i, x_i+1i, i = 0, . . . , n−1 m˚uˇzeme zkonstruovat následuj´ıc´ım zp˚usobem

φ_i(x) = f (x_i) + f (x_i+1) − f (x_i)

h_i (x − x_i) , (1.18)

kde h_i = x_i+1− x_i.

Grafem tohoto spline je lomen´a ˇc´ara. V´ıce informac´ı lze z´ıskat napˇr´ıklad v [4].

(19)

1.2.2 Kvadratick´ y interpolaˇ cn´ı spline

Definice 1.2.2. (Kvadratický interpolaˇcn´ı spline) Kvadratickým splinem nazýváme funkci φ(x), která je spojitá na intervalu hx_0,x_ni, pˇriˇcemˇz na kaˇzdém podintervalu hx_i, x_i+1i, i = 0, . . . , n − 1, je polynomem druhého stupnˇe. Tyto polynomy na sebe v uzlových bodech hladce navazuj´ı – maj´ı spojitou prvn´ı derivaci.

Kvadratický interpolaˇcn´ı spline procházej´ıc´ı uzlovými body, tj. φ(x_i) = f (x_i), na podintervalu hx_i, x_i+1i, i = 0, . . . , n − 1, m˚uˇzeme zkonstruovat následuj´ıc´ım zp˚usobem

φ_i(x) = d_i+1− d_i

2(x_i+1− x_i)(x − x_i)²+ d_i(x − x_i) + f (x_i), (1.19) kde d_i+1 = 2^{f (x}_xⁱ⁺¹^{)−f (x}ⁱ⁾

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

Hodnoty di vycházej´ı z podm´ınky prvn´ı spojité derivace ve vnitˇrn´ıch bodech kvadratického splinu. Jednotlivé spline na intervalech tedy ve výsledku tvoˇr´ı hladkou kˇrivku.

Jelikoˇz se budeme zabývat pˇrirozeným kvadratickým splinem, klademe d1 = 0.

V´ıce informac´ı týkaj´ıc´ı se konstrukce kvadratického interpolaˇcn´ıho splinu m˚uˇzeme nalézt v [5].

1.2.3 Kubick´ y interpolaˇ cn´ı spline

Interpolace pomoc´ı kubických splin˚u je nejˇcastˇeji vyuˇz´ıvanou spline interpolac´ı, protoˇze podle [3] se ukazuje, ˇze právˇe tato po ˇcástech polynomiáln´ı interpolace do- sahuje nejlepˇs´ıch výsledk˚u.

Definice 1.2.3. Kubickým splinem nazveme funkci φ(x), která má na intervalu hx_0,x_ni dvˇe spojité derivace a na kaˇzdém podintervalu hx_i, x_i+1i, i = 0, . . . , n − 1, je polynomem tˇret´ıho stupnˇe.

Konstrukce pˇrirozen´eho kubick´eho spline

Mˇejme v uzlových bodech x_i dány funkˇcn´ı hodnoty f (x_i), i = 0, . . . , n, funkce f . Naˇs´ım úkolem je sestrojit kubický interpolaˇcn´ı polynom, který splˇnuje základn´ı

´

ulohu interpolace.

Z definice kubického polynomu φ(x) = ax³ + bx² + cx + d plyne, ˇze kubický polynom je urˇcen ˇctyˇrmi koeficienty. Máme-li n + 1 bod˚u, budeme m´ıt n interval˚u.

Z toho vyplývá, ˇze φ(x) je na intervalu hx₀, x_ni urˇcena 4n parametry. Podm´ınky spojitosti φ(x), φ⁰(x) a φ⁰⁰(x) ve vnitˇrn´ıch bodech intervalu hx₀, x_ni dávaj´ı dalˇs´ıch 3n − 3 podm´ınek. Interpolaˇcn´ı podm´ınky nám daj´ı n + 1 podm´ınek. Celkem tedy máme 4n − 2 podm´ınek pro 4n neznámých. Protoˇze konstruujeme pˇrirozený kubický spline, zbylé dvˇe podm´ınky dopln´ıme tak, ˇze poloˇz´ıme druhé derivace v krajn´ıch bodech intervalu hx₀, x_ni rovny nule.

Pokud budeme vycházet z pˇredpokladu, ˇze φ(x) je kubický polynom, pak φ⁰(x) je kvadratický polynom a φ⁰⁰(x) je lineárn´ı polynom, který procház´ı body [xi, Mi]

(20)

a [x_i+1, M_i+1], kde M_i = f⁰⁰(x_i) a M_i+1 = f⁰⁰(x_i+1) jsou tzv. momenty splinu. Jelikoˇz je pˇr´ımka jednoznaˇcnˇe urˇcena dvˇema body, m˚uˇzeme poloˇzit φ⁰⁰(x) = L₁(x) s nulovou chybou, a zároveˇn m˚uˇzeme odvodit vztahy pro výpoˇcet kubické interpolaˇcn´ı spline funkce:

φ⁰⁰_i(x) = L₁(x) = f⁰⁰(x_i)x − x_i+1

x_i− x_i+1 + f⁰⁰(x_i+1) x − x_i

x_i+1− x_i (1.20)

= M_ix − x_i+1

−h_i + M_i+1x − x_i

h_i = −M_ix − x_i+1

h_i + M_i+1x − x_i

h_i , (1.21)

φ⁰_i(x) = −M_i 2hi

(x − x_i+1)² +M_i+1 2hi

(x − x_i)²+ A_i, (1.22)

φ_i(x) = −M_i

6h_i (x − x_i+1)³+M_i+1

6h_i (x − x_i)³+ A_i(x − x_i) + B_i, (1.23) kde h_i = x_i+1− x_i a A_i a B_i jsou integraˇcn´ı konstanty.

Nyn´ı urˇc´ıme integraˇcn´ı konstanty, pˇriˇcemˇz vyuˇzijeme skuteˇcnosti, ˇze funkce φ_i(x) mus´ı na intervalu hx_i, x_i+1i splˇnovat podm´ınky interpolace:

φi(xi) = f (xi) a φi(xi+1) = f (xi+1).

Dostaneme

φ_i(x_i) = −M_i

6h_i (x_i − x_i+1)³+M_i+1

6h_i (x_i− x_i)³+ A_i(x_i− x_i) + B_i (1.24)

= −M_i

6h_i (−h_i)³ + B_i (1.25)

a m˚uˇzeme vyj´adˇrit

Bi = f (xi) −Mi

6 hi2

. (1.26)

Analogicky dostaneme φ_i(x_i+1) = − M_i

6h_i(x_i+1− x_i+1)³+ M_i+1

6h_i (x_i+1− x_i)³+ A_i(x_i+1− x_i) + B_i (1.27)

= M_i+1 6hi

h³_i + A_ih_i+ B_i (1.28)

= M_i+1

6 h²_i + Aihi+ f (xi) −M_i 6 hi2

. (1.29)

Z tˇechto vztah˚u vyj´adˇr´ıme

A_i = f (x_i+1) − f (x_i) hi

+M_i− M_i+1

6 h_i. (1.30)

Nyn´ı urˇc´ıme momenty splinu M_i. Protoˇze konstruujeme pˇrirozen´y kubick´y spline, tak M₀ = M_n = 0. Momenty splinu urˇc´ıme pomoc´ı dalˇs´ı podm´ınky, kterou mus´ı

(21)

funkce φ_i(x) splˇnovat. Poˇzadujeme, aby derivace funkce φ_i(x) byla zleva i zprava spojit´a, tj. φ⁰_i−1(x_i) = φ⁰_i(x_i).

Jelikoˇz v´ıme, ˇze

φ⁰_i(x) = −M_i

2h_i(x − x_i+1)² +M_i+1

2h_i (x − x_i)²+ A_i. (1.31) Potom

φ⁰_i−1(xi) = −M_i−1

2h_i−1(xi− xi)²+ M_i

2h_i−1(xi− xi−1)²+ Ai−1 (1.32)

= M_i 2hi−1

(h_i−1)²+ A_i−1 (1.33)

= M_i

2h_i−1(h_i−1)²+f (x_i) − f (x_i−1)

h_i−1 +M_i−1− M_i

6 h_i−1 (1.34)

= M_i−1+ 2M_i

6 h_i−1+f (x_i) − f (x_i−1)

h_i−1 . (1.35)

Analogicky dostaneme φ⁰_i(x_i) = −M_i

2hi

(x_i− x_i+1)²+ M_i+1 2hi

(x_i− x_i)²+ A_i (1.36)

= −M_i

2h_i(−hi)²+ Ai (1.37)

= −M_i

2 h_i+ f (x_i+1) − f (x_i)

h_i + M_i− M_i+1

6 h_i (1.38)

= −Mi+1+ 2Mi

6 h_i+f (xi+1) − f (xi)

h_i . (1.39)

Nyn´ı m˚uˇzeme napsat n´asleduj´ıc´ı rovnost:

M_i−1+ 2M_i

6 h_i−1+f (x_i) − f (x_i−1)

h_i−1 = −M_i+1+ 2M_i

6 h_i+f (x_i+1) − f (x_i)

h_i , (1.40) M_i−1+ 2M_i

6 h_i−1+M_i+1+ 2M_i

6 h_i = f (x_i+1) − f (x_i)

h_i −f (x_i) − f (x_i−1)

h_i−1 . (1.41) Pˇri pˇredpokladu ekvidistantn´ıho dˇelen´ı intervalu plat´ı: h_i−1= h_i.

M_i−1+ 2M_i

6 h_i+M_i+1+ 2M_i

6 h_i = f (x_i+1) − f (x_i)

h_i −f (x_i) − f (x_i−1)

h_i ,(1.42) M_i−1+ 4M_i+ M_i+1

6 h_i = f (x_i+1) − 2f (x_i) + f (x_i−1)

h_i . (1.43)

Pˇri praktickém výpoˇctu interpolace pomoc´ı kubické spline funkce postupujeme v n´ıˇze naznaˇcených kroc´ıch.

Nejdˇr´ıve vypoˇcteme pomoc´ı rovnic (1.41) momenty splinu Mi, i = 1, . . . , n − 1.

V praxi se pˇri v´ypoˇctu moment˚u splinu pˇri pˇredpokladu konstrukce pˇrirozen´eho

(22)

kubického splinu vyuˇz´ıvá Gaussova eliminace aplikovaná na tˇr´ıdiagonáln´ı matici pro ˇreˇsen´ı n − 1 rovnic o n − 1 neznámých.

Soustava m´a tvar







h0+h1

3 h1

. .. ...6 . ..

hi−1

6

hi−1+hi

3

hi+1

. .. . ..6 . ..

hn−2

6

hn−2+hn−1

3











 M₁

... Mi−1

M_i M_i+1

... M_n−1







=





 g₁

... gi−1

g_i g_i+1

... g_n−1







, (1.44)

kde gi = ^{f (x}ⁱ⁺¹_h^{)−f (x}ⁱ⁾

i − ^{f (x}ⁱ^{)−f (x}_h ⁱ⁻¹⁾

i−1 .

Po v´ypoˇctu moment˚u splinu mus´ıme urˇcit integraˇcn´ı konstanty A_i, i = 1, . . . , n, a B_i, i = 1, . . . , n, podle vztah˚u (1.26) a (1.30).

V tuto chv´ıli jiˇz známe vˇsechny neznámé a m˚uˇzeme vypoˇc´ıtat kubické spliny na jednotlivých intervalech dosazen´ım do vztahu (1.23). Bliˇzˇs´ı informace lze nalézt ve [2] a [6].

1.3 Metoda nejmenˇ s´ıch ˇ ctverc˚ u

Metoda nejmenˇs´ıch ˇctverc˚u jiˇz nen´ı interpolaˇcn´ı metoda. Jej´ı princip m˚uˇzeme popsat tak, ˇze zadanými body [x_i, y_i] pro i = 1, . . . , n prokládáme funkci φ tak, aby souˇcet druhých mocnin rozd´ılu hodnot y_i a funkˇcn´ıch hodnot φ(x_i) byl minimáln´ı.

Odvozen´ı metody nejmenˇs´ıch ˇctverc˚u pro line´arn´ı regresi

Uvaˇzujme mnoˇzinu n bod˚u [x_i, y_i] pro i = 1, . . . , n. Vzdálenost bod˚u od pˇr´ımky y = ax + b m˚uˇzeme vyjádˇrit jednoduchým zp˚usobem s_i = |ax_i + b − y_i|, kde i = 1, . . . , n.

Smyslem t´eto metody je minimalizovat funkci S(a, b) = Pn

i=1(ax_i + b − y_i)². Protoˇze je známo, ˇze lokáln´ı extrém diferencovatelné funkce m˚uˇze nastat pouze ve stacionárn´ım bodˇe, vyuˇzijeme pro minimalizaci funkce S parciáln´ı derivace

∂S

∂a = 2

n

X

i=1

(ax_i+ b − y_i)x_i = 2 a

n

X

i=1

x²_i + b

n

X

i=1

x_i−

n

X

i=1

x_iy_i

!

, (1.45)

∂S

∂b = 2

n

X

i=1

(axi+ b − yi) = 2 a

n

X

i=1

xi+ b

n

X

i=1

1 −

n

X

i=1

yi

!

. (1.46)

Proto nyn´ı poloˇz´ıme tyto parci´aln´ı derivace rovny nule:

2 a

n

X

i=1

x²_i + b

n

X

i=1

x_i −

n

X

i=1

x_iy_i

!

= 0, (1.47)

2 a

n

X

i=1

x_i+ b

n

X

i=1

1 −

n

X

i=1

y_i

!

= 0.

(23)

Tyto rovnice uprav´ıme jednoduchými úpravami na následuj´ıc´ı tvar,

a

n

X

i=1

x²_i + b

n

X

i=1

x_i =

n

X

i=1

x_iy_i, (1.48)

a

n

X

i=1

xi+ bn =

n

X

i=1

yi.

ˇc´ımˇz z´ıskáme soustavu rovnic s neznámými a, b.

Pokud tedy chceme naj´ıt aproximaˇcn´ı polynom prvn´ıho ˇr´adu φ(x) = ax + b, mus´ıme zjistit hodnotu koeficient˚u a, b. Tyto koeficienty zjist´ıme vyˇreˇsen´ım soustavy (1.48).

Protoˇze vˇsechny hlavn´ı subdeterminanty matice jsou kladné, je kvadratická forma pozitivnˇe definitn´ı. Podle Sylvestrova kritéria je tedy nalezený stacionárn´ı bodem minima.

Pro aproximaci obecn´ym polynomem φ(x) = Pm

i=0b_ixⁱ ˇreˇs´ıme n´asleduj´ıc´ı soustavu rovnic:

nb₀+

n

X

i=1

x_ib₁+

n

X

i=1

x²_ib₂ + . . . +

n

X

i=1

x^m_i b_m =

n

X

i=1

y_i, (1.49)

n

X

i=1

x_ib₀+

n

X

i=1

x²_ib₁+ . . . +

n

X

i=1

x^m+1_i b_m =

n

X

i=1

x_iy_i, ...

n

X

i=1

x^m_i b₀+

n

X

i=1

x^m+1_i b₁+ . . . +

n

X

i=1

x^2m_i b_m =

n

X

i=1

x^m_i y_i.

Dalˇs´ı informace lze nal´ezt v [7].

(24)

2 Numerick´ a derivace

Ve druhé kapitole této práce se budeme zabývat metodami pro urˇcen´ı numerické derivace funkce f v bodˇe.

Definice 2.0.1. Funkce f m´a v bodˇe x₀ derivaci, je-li definov´ana v okol´ı bodu x₀ a existuje limita

h→0lim

f (x0+ h) − f (x0)

h . (2.1)

Tuto limitu naz´yv´ame derivac´ı funkce f v bodˇe x₀ a znaˇc´ıme ji f⁰(x₀).

Metody pro výpoˇcet numerické derivace v bodˇe vycházej´ı pˇr´ımo z definice nebo napˇr´ıklad z myˇslenky nahrazen´ı funkce f v okol´ı bodu x₀ interpolaˇcn´ım polynomem (funkci lze nahradit napˇr´ıklad i aproximac´ı z´ıskanou metodou nejmenˇs´ıch ˇctverc˚u, nebo ˇCebyˇsevovými polynomy) [8].

2.1 Derivace pomoc´ı interpolace

V následuj´ıc´ıch odstavc´ıch ukáˇzeme, jak m˚uˇzeme odvodit vztahy pro výpoˇcet derivace funkce v bodˇe pomoc´ı interpolace.

Funkci f m˚uˇzeme aproximovat Lagrangeov´ym interpolaˇcn´ım polynomem f (x) = L_n(x) + 1

(n + 1)!f⁽ⁿ⁺¹⁾(ξ)ω_n+1(x), (2.2) f (x) =

n

X

i=0

f (x_i)l_i(x) + 1

(n + 1)!f⁽ⁿ⁺¹⁾(ξ)ω_n+1(x). (2.3) Pro chybu Lagrangeovy interpolace plat´ı

f (x) − L_n(x) = 1

(n + 1)!f⁽ⁿ⁺¹⁾(ξ)ω(x), (2.4) kde ωn+1(x) =

n

Q

i=0

(x − xi).

Zderivujeme-li ω_n+1(x) podle promˇenn´e x, dostaneme ω_n+1⁰ (x) =

n

X

i=0 n

Y

j=0,j6=i

(x − x_j). (2.5)

(25)

Nyn´ı provedeme derivaci funkce f f⁰(x) =

n

X

i=0

f (x_i)l_i⁰(x) + 1

(n + 1)!f⁽ⁿ⁺¹⁾(ξ)

n

Y

j=0,j6=i

(x − x_j)

!

(2.6) a pro derivaci funkce f v bodˇe x₀ tedy plat´ı

f⁰(x₀) =

n

X

i=0

f (x_i)l_i⁰(x₀) + 1

(n + 1)!f⁽ⁿ⁺¹⁾(ξ)

n

Y

j=1

(x₀− x_j)

!

. (2.7)

Podle vztahu (2.7) m˚uˇzeme tedy konstruovat vztahy pro v´ypoˇcet derivace funkce f (x) v bodˇe x0. Nyn´ı odvod´ıme vztah pro derivaci v bodˇe pomoc´ı interpolaˇcn´ıho polynomu prvn´ıho stupnˇe.

Funkci f nahrad´ıme v okol´ı bodu x₀ interpolaˇcn´ım polynomem prvn´ıho stupnˇe.

Tento polynom pot´e zderivujeme podle promˇenn´e x.

Sestroj´ıme interpolaˇcn´ı polynom prvn´ıho stupnˇe pro body [x_i, f (x_i)] a [x_i + h, f (x_i + h)],

L₁(x) = f (x_i)x − (x_i+ h)

x_i− (x_i+ h) + f (x_i+ h) x − x_i

(x_i+ h) − x_i (2.8)

= f (x_i)x − (x_i+ h)

−h + f (x_i+ h)x − x_i

h . (2.9)

Polynom L₁(x) zderivujeme podle promˇenn´e x. Dostaneme L⁰₁(x) = 1

h(f (x_i+ h) − f (x_i)) (2.10) f⁰(x₀) = f (x₀+ h)–f (x₀)

h − 1

2f⁽²⁾(ξ)h. (2.11)

Pokud stejným zp˚usobem vyuˇzijeme pro urˇcen´ı prvn´ı derivace funkce f (x) v bodˇe x₀ interpolaˇcn´ı polynom druhého ˇrádu dostáváme

L⁰₂(x) = 1

2h(f (xi+ h) − f (xi− h)), (2.12) f⁰(x₀) = f (x₀+ h)–f (x₀− h)

2h − 1

6f⁽³⁾(ξ)h². (2.13) Jestliˇze L2 zderivujeme jeˇstˇe jednou dostáváme pˇredpis pro aproximaci druhé derivace funkce f v bodˇe x₀

f⁰⁰(x₀) = f (x₀+ h) − 2f (x₀) + f (x₀− h)

h² . (2.14)

Pˇri praktickém výpoˇctu numerické derivace funkce f (x) v bodˇe x₀, mus´ıme dbát na vliv zaokrouhlovac´ıch chyb, které mohou m´ıt podstatnou mˇerou vliv na výsledek.

Jmenovatelé uvedených vzorc˚u obsahuj´ı krok h, který mus´ı být jakýmsi kompromi- sem mezi dostateˇcnˇe pˇresnou aproximac´ı, která vyˇzaduje dostateˇcnˇe malý krok h, a zaokrouhlovac´ı chybou, která se naopak pˇri malém h zvyˇsuje. Bliˇzˇs´ı informace v [3]

nebo [9].

(26)

2.2 Richardsonova extrapolace

Richardsonova extrapolace je technika, která je v praxi hojnˇe vyuˇz´ıvaná. Setkáme se s n´ı napˇr´ıklad v numerické integraci v Rombergovˇe kvadratuˇre nebo napˇr´ıklad v Bulirsch-Stoerovu algoritmu, jenˇz se zabývá ˇreˇsen´ım obyˇcejných diferenciáln´ıch rovnic. Jej´ı podstata vycház´ı z pˇredpokladu, ˇze ze dvou pˇribliˇzných výsledk˚u m˚uˇzeme pomoc´ı lineárn´ı kombinace vypoˇc´ıtat tˇret´ı, který bude pˇresnˇejˇs´ı, pˇriˇcemˇz tato nová hodnota se nacház´ı mimo interval, ohraniˇcený pˇredcházej´ıc´ımi dvˇema hodnotami (proto hovoˇr´ıme o extrapolaci).

Odvozen´ı Richardsonovy extrapolace

Mˇejme funkci R a krok h. Pˇredpokládejme, ˇze funkci R je moˇzné vyjádˇrit moc- ninnou ˇradou

R (h) = a₀+ a₁h + a₂h²+ a₃h³+ a₄h⁴+ . . . . (2.15) Jestliˇze h < 1, pak je R(h) dobrou aproximac´ı ˇclenu a₀. Lepˇs´ı aproximaci ovˇsem z´ısk´ame, kdyˇz urˇc´ıme hodnotu funkce R s krokem ^h₂

R h 2

= a₀+ a₁h

2 + a₂ h 2

2

+ a₃ h 2

3

+ . . . .

| {z }

O(h²)

. (2.16)

Pomoc´ı vhodné lineárn´ı kombinace m˚uˇzeme vyjádˇrit ˇclen a0 s chybou O(h²). Ve vyjádˇren´ı (2.15) a (2.16) zanedbáme dalˇs´ı ˇcleny rozvoje

R (h) − a0 = a1h + O h² , (2.17)

R h 2

− a₀ = a₁h

2 + O h² . (2.18)

Rovnici (2.17) odeˇcteme od dvojn´asobku rovnice (2.18) R₂(h) = 2R h

2

− R (h) = a₀+ a⁽²⁾₂ (h)²+ a⁽²⁾₃ (h)³+ . . . .

| {z }

O(h²)

. (2.19)

Pomoc´ı vhodné lineárn´ı kombinace vztah˚u, z nichˇz oba aproximuj´ı hodnotu funkce R s chybou O(h), jsme vyjádˇrili vztah pro a₀ s chybou O(h²).

Podle [13] plat´ı, ˇze |a⁽²⁾_i | < |a_i|. R₂(h) je lepˇs´ı aproximace neˇz R(h) a pro mal´e h je i lepˇs´ı aproximac´ı neˇz R(^h₂).

Obecn´y vztah pro Richardsonovu extrapolaci m˚uˇzeme zapsat pomoc´ı n´asleduj´ıc´ı iteraˇcn´ı formule:

R_j+1(h) = R_j h 2

+ Rj h

2 − Rj(h)

2^j− 1 , kde j = 1, 2, . . . . (2.20) Výpoˇcet m˚uˇzeme uspoˇrádat do následuj´ıc´ı tabulky:

(27)

R(h)

R(^h₂) R₂(h)

R(^h₄) R₂(^h₂) R₃(h)

R(^h₈) R₂(^h₄) R₃(^h₂) R₄(h)

R(₁₆^h) R₂(^h₈) R₃(^h₄) R₄(^h₂) R₅(h) ... ... ... ... ... . ..

↑ ↑ ↑ ↑ ↑

O(h) O(h²) O(h³) O(h⁴) O(h⁵)

Tabulka 2.1: Richardsonova extrapolace

Pˇri výpoˇctu postupujeme po ˇrádc´ıch, pˇriˇcemˇz prvek v prvn´ım sloupci vypoˇcteme pomoc´ı základn´ı metody. Ostatn´ı prvky v ˇrádku vypoˇcteme pomoc´ı vztahu (2.20).

Výpoˇcet ukonˇc´ıme ve chv´ıli, kdy je rozd´ıl dvou diagonáln´ıch prvk˚u tabulky menˇs´ı neˇz poˇzadovaná pˇresnost. V´ıce napˇr´ıklad v [10], [11], [12] nebo [13].

(28)

3 Numerick´ a integrace funkc´ı

Ve dalˇs´ı kapitole se budeme zab´yvat numerickou integrac´ı.

Je-li funkce f (x) na intervalu ha, bi spojitá a známe jej´ı primitivn´ı funkci F (x), m˚uˇzeme hodnotu urˇcitého integrálu Rb

a f (x)dx urˇcit pomoc´ı Newton - Leibnizova vztahu

Z b a

f (x)dx = F (b) − F (a). (3.1)

V praxi ovˇsem ˇcasto nem˚uˇzeme primitivn´ı funkci F (x) naj´ıt analyticky nebo je analytický výpoˇcet integrálu pˇr´ıliˇs sloˇzitý a pracný. V tuto chv´ıli vyuˇzijeme metody numerické kvadratury, kdy funkci f (x) obvykle nahrazujeme jednoduˇsˇs´ı aproximuj´ıc´ı funkc´ı φ(x)(nejˇcastˇeji polynomem), a hodnotu integrálu f (x) poloˇz´ıme pˇribliˇznˇe hodnotˇe integrálu φ(x) [14].

Definice 3.0.1. Vyj´adˇreme integr´al ve tvaru

I(f ) = Z b

a

f (x)dx ≈ I_h(f ) = Z b

a

φ(x)dx =

n

X

i=0

α_if (x_i) (3.2)

kde Ih(f ) = Pn

i=0αif (xi) se nazývá kvadraturn´ı formule, αi jsou koeficienty kvadraturn´ı formule, které nezávis´ı na funkci f (x), x_i jsou uzly kvadraturn´ı formule, pˇriˇcemˇz x_i ∈ ha, bi. Chybu kvadraturn´ı formule R_n(f ) vyjádˇr´ıme R_n(f ) = I(f ) − Ih(f ).

U kvadraturn´ıch formul´ı nás zaj´ımá, jakou pˇresnost daná kvadraturn´ı formule pˇri aproximaci urˇcitého integrálu má. ˇRekneme, ˇze kvadraturn´ı formule je ˇrádu n, kdyˇz integruje polynomy stupnˇe n s nulovou chybou a polynomy stupnˇe n + 1 integruje s nenulovou chybou [2]. To nás vede k následuj´ıc´ı definici.

Definice 3.0.2. ˇR´ad kvadraturn´ı formule I_h(f ) = Pn

i=0α_if (x_i) je maxim´aln´ı ˇc´ıslo m ∈ N ∪ {0} takov´e, ˇze Rb

ap(x)dx =Pn

i=0αip(xi), kde p(x) je polynom stupnˇe m.

3.1 Newton-Cotesovy vzorce

Velkou skupinou v oblasti numerick´e integrace funkc´ı pˇredstavuj´ı Newton-Cotesovy vzorce. Tyto kvadraturn´ı vzorce pˇredpokl´adaj´ı ekvidistantn´ı dˇelen´ı intervalu integrace.

(29)

Newton-Cotesovy vzorce lze rozdˇelit do dvou podskupin. Prvn´ı skupinu oznaˇcujeme jako otevˇrené a druhou jako uzavˇrené. Rozd´ıl mezi tˇemito skupinami spoˇc´ıvá v pˇr´ıstupu ke krajn´ım bod˚um intervalu integrace. Uzavˇrené vzorce tyto body povaˇzuj´ı za uzly kvadratury, otevˇrené nikoli. Uzly kvadratury otevˇrených vzorc˚u jsou urˇceny symet- ricky podle stˇredu intervalu.

Newton-Cotesovy kvadraturn´ı vzorce aproximuj´ı hodnotu integr´alu pomoc´ı nahrazen´ı funkce f (x) Lagrangeov´ym interpolaˇcn´ım polynomem L_n(x). V´ıce napˇr´ıklad v [2], [15].

Odvozen´ı Newton-Cotesov´ych vzorc˚u

Nyn´ı odvod´ıme obecný tvar Newton - Cotesových vzorc˚u. Hodnotu integrálu funkce f na intervalu ha, bi aproximujeme pomoc´ı Lagrangeova interpolaˇcn´ıho polynomu L_n(x).

Z b a

f (x)dx ≈ Z b

a

L_n(x)dx =

= Z b

a n

X

i=0

f (x_i)

n

Y

j=0,j6=i

x − x_j x_i− x_j

| {z }

li(x)

dx =

n

X

i=0

f (x_i) Z b

a

l_i(x)dx

| {z }

αi

, (3.3)

kde a = x₀ < x₁ < x₂ < . . . < x_n = b. Urˇc´ıme koeficienty α_i. Poloˇzme h = ^b−a_n a x₀ = a + 0h, x₁ = a + h, . . . , x_n= a + nh. (3.4) Nyn´ı provedeme substituci

x = a + th, dx = hdt . (3.5)

Pouˇzit´ım substituce doˇslo ke zmˇenˇe integraˇcn´ıch mez´ı t = x − a

h . (3.6)

Pokud za promˇennou x dosad´ıme doln´ı mez intervalu integrace dostáváme jako novou doln´ı mez nulovou hodnotu. Dosad´ıme-li do vztahu (3.6) za promˇennou x výraz b = a + nh, dostáváme pro novou horn´ı mez integrálu hodnotu n:

α_i = h Z n

0 n

Y

j=0,j6=i

t − j i − j

dt . (3.7)

Z Newton - Cotesových vzorc˚u lze pomˇernˇe jednoduˇse odvodit pravidla a jejich sloˇzené varianty pro výpoˇcet urˇcitého integrálu.

Princip aproximace urˇcitého integrálu pomoc´ı sloˇzených vzorc˚u spoˇc´ıvá v rozdˇelen´ı intervalu na jednotlivé podintervaly. Pˇriˇcemˇz na kaˇzdý z podinterval˚u pouˇzijeme jednoduché Newton - Cotesovy vzorce. Hodnotu integrálu poté aproximujeme hodnotou, kterou vypoˇcteme jako souˇcet obsahu ploch na jednotlivých intervalech, ˇc´ımˇz dostáváme pˇresnˇejˇs´ı aproximaci daného integrálu.

(30)

ˇS´ıˇrku intervalu h urˇc´ıme jako h = ^b−a_m , kde m je poˇcet podinterval˚u. Dalˇs´ı informace m˚uˇzeme nal´ezt v [2] nebo [16].

Vˇeta 3.1.1. Kvadraturn´ı formule z´ıskaná integrac´ı interpolaˇcn´ıho polynomu pro uzlové body [x_i, f (x_i)] pro i = 0, . . . , n má stupeˇn pˇresnosti alespoˇn n.

D˚ukaz. Uvaˇzujme integrál I(f ), který aproximujeme Lagrangeovým interpolaˇcn´ım polynomem L_n stupnˇe n s chybou E_n(f )

I(f ) =

b

Z

a

L_n(x) + E_n(f )dx. (3.8)

Uprav´ıme pravou stranu rovnosti

I(f ) =

b

Z

a n

X

i=0

f (xi)li(x) + En(f )dx (3.9)

=

b

Z

a n

X

i=0

f (x_i)l_i(x)dx +

b

Z

a

E_n(f )dx.

=

n

X

i=0

f (x_i)

b

Z

a

l_i(x)dx +

b

Z

a

E_n(f )dx, kde

b

Z

a

l_i(x)dx = α_i.

Nyn´ı vyj´adˇr´ıme chybu En(f )

I(f ) −

n

X

i=0

f (xi)αi =

b

Z

a

En(f )dx. (3.10)

Pˇriˇcemˇz pro chybu Lagrangeovy interpolace plat´ı E_n(f ) = f⁽ⁿ⁺¹⁾(ξ)ω_n+1(x)

(n + 1)!. (3.11)

Nyn´ı m˚uˇzeme vidˇet, ˇze vyjádˇren´ı chyby En(f ) obsahuje (n + 1)-n´ı derivac´ı funkce f . Z toho pˇr´ımo plyne, ˇze chyba interpolaˇcn´ı kvadraturn´ı formule bude nulová pro vˇsechny polynomy stupnˇe nejvýˇse n a ˇrád pˇresnosti bude tedy alespoˇn n.

Pozdˇeji ukáˇzeme, ˇze interpolaˇcn´ı kvadraturn´ı formule má ˇrád pˇresnosti pro n + 1 bod˚u nejvýˇse 2n + 1. Dalˇs´ı informace v [17].

3.1.1 (Sloˇ zen´ e) obd´ eln´ıkov´ e pravidlo

Je nejjednoduˇsˇs´ım vzorcem pro kvadraturu a ˇrad´ıme jej mezi otevˇrené Newton - Cotesovy vzorce. Jeho ˇrád pˇresnosti je jedna, ovˇsem je vhodné jej vyuˇz´ıt v pˇr´ıpadˇe,

(31)

kdy nem˚uˇzeme vypoˇc´ıst integr´al sloˇzitˇejˇs´ımi metodami z d˚uvodu nemoˇznosti urˇcen´ı funkˇcn´ı hodnoty na kraj´ıch integr´alu.

Princip tohoto pravidla spoˇc´ıvá v urˇcen´ı funkˇcn´ı hodnoty stˇredu intervalu. Hod- nota urˇcitého integrálu je tak aproximována jako obsah obdéln´ıku, kdy jedna strana je urˇcena ˇs´ıˇrkou intervalu a druhá funkˇcn´ı hodnotou.Dostaneme

Z b a

f (x)dx ≈ hf a + b 2

. (3.12)

Pokud jsme ovˇsem schopni urˇcit funkˇcn´ı hodnoty i na kraj´ıch integraˇcn´ıch mez´ı, je vhodnˇejˇs´ı pouˇz´ıt d´ale popsan´e metody.

Pˇredpis pro aproximaci urˇcitého integrálu pomoc´ı sloˇzeného obdéln´ıkového pravidla:

Z b a

f (x)dx ≈ h

m

X

i=1

f x_i−1+ x_i 2

, (3.13)

kde m je poˇcet podinterval˚u. Dalˇs´ı informace lze nal´ezt v [3].

3.1.2 (Sloˇ zen´ e) lichobˇ eˇ zn´ıkov´ e pravidlo

Jedná se o jednoduchou metodu s ˇrádem pˇresnosti jedna, která je snadno im- plementovatelná pomoc´ı poˇc´ıtaˇce. Hodnota integrálu je vypoˇctena jako obsah li- chobˇeˇzn´ıku, kdy funkci f nahrazujeme Lagrangeovým interpolaˇcn´ım polynomem prvn´ıho ˇrádu. Pokud vyuˇzijeme sloˇzenou variantu tohoto pravidla je interval rozdˇelen na ˇcásti a na kaˇzdé této ˇcásti je zkonstruován lichobˇeˇzn´ık. Hodnota integrálu je urˇcena jako souˇcet obsahu jednotlivých lichobˇeˇzn´ık˚u.

Pouˇzit´ım vztahu (3.3) odvod´ıme vztah pro lichobˇeˇzn´ıkov´e pravidlo. Urˇc´ıme ˇs´ıˇrku intervalu h = x₁ − x₀, kde a = x₀ < x₁ = b. Integr´al funkce f na intervalu ha, bi aproximujeme pomoc´ı Lagrangeova interpolaˇcn´ıho polynomu prvn´ıho stupnˇe

Z b a

f (x)dx ≈

1

X

i=0

α_if (x_i). (3.14)

Pro koeficienty α_i dle vztahu (3.7) plat´ı α_i = h

Z 1 0

n

Y

j=0,j6=i

t − j i − j

dt, (3.15)

α0 = h Z 1

0

t − 1

0 − 1 dt = h Z 1

0

1 − tdt = h

t − t²

2

1 0

= h

2, (3.16) α1 = h

Z 1 0

t − 0

1 − 0 dt = h Z 1

0

tdt = h t² 2

1 0

= h

2. (3.17)

T´ımto dostáváme vztah pro lichobˇeˇzn´ıkové pravidlo Z b

a

f (x)dx ≈ h

2[f (a) + f (b)] . (3.18)