Ř ÍZE APLIKACE NUMERICKÝCH METOD PRO Ř EŠENÍDYNAMIKY ODVÍJENÍ P TECHNICKÁ UNIVERZITA V LIBERCI

TECHNICKÁ UNIVERZITA V LIBERCI

Fakulta mechatroniky a mezioborových inženýrských studií

Katedra: Katedra řídící techniky

Studijní program: M2612 - Elektrotechnika a informatika Obor: 3902T005 – Automatické řízení a informatika

APLIKACE NUMERICKÝCH METOD PRO ŘEŠENÍ DYNAMIKY ODVÍJENÍ PŘÍZE

(Application of Numerical Metods for Solving Dynamics of Unwinding Yarn)

Jiří Pliska

Vedoucí práce: Ing. Libor Tůma, CSc.

Konzultant: Doc.Ing. Vladimír Kracík, CSc.

Ing. Jakub Kašše

Počet

stran práce a příloh stran textu obrázků tabulek příloh

61 45 13 8 18

V Liberci dne:

Prohlášení

Byl(a) jsem seznámen(a) s tím, že na mou diplomovou práci se plně vztahuje zákon č.

121/2000 Sb. o právu autorském, zejména § 60 – školní dílo.

Beru na vědomí, že Technická univerzita v Liberci (TUL) nezasahuje do mých autorských práv užitím mé diplomové práce pro vnitřní potřebu TUL.

Užiji-li diplomovou práci nebo poskytnu-li licenci k jejímu využití, jsem si vědom(a) povinnosti informovat o této skutečnosti TUL; v tomto případě má TUL právo ode mne požadovat úhradu nákladů, které vynaložila na vytvoření díla, až do jejich skutečné výše.

Diplomovou práci jsem vypracoval(a) samostatně s použitím uvedené literatury a na základě konzultací s vedoucím diplomové práce a konzultantem.

V Liberci dne: Jiří Pliska

Poděkování

Rád bych touto cestou poděkoval vedoucímu diplomové práce Ing. Liboru Tůmovi CSc.

a také konzultantům Ing. Jakubovi Kašše a Doc. Ing. Vladimíru Kracíkovi CSc. za podnětné rady a zodpovědné vedení. Dále bych chtěl poděkovat celé své rodině a kamarádům za všestrannou podporu při studiu a při tvorbě této diplomové práce.

Abstrakt

Cílem diplomové práce je naprogramovat různé metody pro numerické řešení obyčejných diferenciálních rovnic s počátečními podmínkami a s pomocí těchto metod řešit model dynamiky odvíjení příze vytvořený na Technické univerzitě v Liberci.

V rámci této práce byla naprogramována samostatná knihovna, ve které je implementován poměrně značný počet metod pro řešení obyčejných diferenciálních rovnic.

Tato knihovna byla vytvořena tak, aby bylo umožněno její univerzální použití i pro řešení jiných úloh než je odvíjení příze. Dále byly provedeny simulační výpočty na modelu dynamiky odvíjení příze s použitím vytvořené knihovny. Pro porovnání výsledků těchto výpočtů byla vytvořena samostatná aplikace.

Abstract

The aim of this diploma thesis is to program different numerical methods for ordinary differential equations with initial conditions and to solve a model of the dynamics of unwinding yarn with these methods. The model was made in the Technical university of Liberec.

In the framework of this thesis, an independent library was programmed, in which a relatively considerable number of numerical methods for differential equations was implemented. This library was made for universal use, also to solve other problems than the model of unwinding yarn. Simulation computations of the model of unwinding yarn were performed with use of the library. A standalone application was made to allow the comparison the results of these computations.

Obsah:

1. Úvod ...8

2. Numerické řešení obyčejných diferenciálních rovnic ...9

2.1. Formulace úlohy...9

2.2. Eulerova metoda. Problematika odhadu chyby ...11

2.3. Rungovy–Kuttovy metody ...12

2.3.1. Obecná jednokroková metoda ...12

2.3.2. Metoda Taylorova rozvoje...13

2.3.3. Klasické Rungovy–Kuttovy metody ...14

2.3.4. Odhad chyby Rungových–Kuttových metod. Automatická volba kroku...16

2.4. Lineární k-krokové metody ...21

2.4.1. Metody numerické integrace ...23

2.4.2. Použití k-krokových metod, metody prediktor-korektor ...24

3. Knihovna pro řešení diferenciálních rovnic ...26

3.1. DLL – dynamicky linkovaná knihovna...26

3.2. Popis vytvořené knihovny ...27

3.2.1. Pomocné procedury ...28

3.2.2. Popis procedur pro řešení diferenciálních rovnic...30

4. Porovnání metod pro řešení diferenciálních rovnic...36

4.1. Model odvíjení příze...36

4.2. Aplikace pro porovnání výsledků simulace dynamiky odvíjení příze ...39

4.3. Simulační výpočty ...41

4.3.1. Simulační výpočty s pevným krokem ...42

4.3.2. Simulační výpočty s proměnným krokem...46

4.3.3. Shrnutí výsledků simulačních výpočtů...48

5. Závěr...52

6. Literatura ...53

7. Přílohy ...54

1. Úvod

Procesy odvíjení příze patří mezi základní technologické operace, které jsou součástí celé řady textilních technologií a současně výrazným způsobem tyto technologie ovlivňují, např. z hlediska maximální výkonnosti, spolehlivosti, dynamického chování a podobně.

Odvíjení příze z cívky je v řadě případů spojeno s rotací příze kolem pevné osy a dochází tak k fyzikálnímu jevu, který je označován jako balónování příze. Dynamice balónující příze je věnována řada studií, které kromě několika případů řeší pouze dynamiku ve stacionárním režimu, tj. po odeznění přechodové fáze děje. Na Technické univerzitě v Liberci je řešena problematika dynamiky příze z obecného hlediska. Řeší se přechodový proces při odvíjení příze z cívky z počátečních podmínek, kdy je příze v klidu, až do ustáleného - stacionárního stavu.

V rámci výzkumného záměru fakulty strojní byl vytvořen matematický model, který popisuje daný jev. Tento model je řešen jako Cauchyho úloha, tedy jako popis soustavou obyčejných diferenciálních rovnic s počátečními podmínkami. Pro numerické řešení diferenciálních rovnic je možné použít celou řadu metod a výhodnost či nevýhodnost konkrétní metody podstatně závisí na individuálním charakteru daného problému. Zvolit pro daný problém konkrétní metodu řešení tedy většinou znamená vybrat z celé řady možností. Exaktní řešení problému takovéhoto výběru je však značně komplikované a většinou se omezuje na výběr podle určitých zkušeností.

V rámci této diplomové práce je pro řešení úlohy dynamiky odvíjení příze použito různých numerických metod pro řešení obyčejných diferenciálních rovnic s počátečními podmínkami. Počítačová implementace těchto metod si navíc klade za cíl určitou nezávislost na řešené úloze. To znamená, že je vytvořena tak, aby ji bylo možné použít i pro řešení jiných úloh, než je simulace dynamiky odvíjení příze. Nejdříve je věnována pozornost teoretickému přehledu numerických metod vhodných pro řešení diferenciálních rovnic, pak následuje popis jejich počítačové realizace spolu s možností oddělit statickou a dynamickou část simulačního experimentu a na závěr jsou shrnuty výsledky simulačních výpočtů na modelu dynamiky odvíjení příze.

2. Numerické řešení obyčejných diferenciálních rovnic

V této kapitole je uveden přehled některých numerických metod pro řešení obyčejných diferenciálních rovnic s počátečními podmínkami a nastíněna problematika odhadu chyb, které vznikají při numerickém řešení diferenciálních rovnic. Jsou zde opomenuty otázky stability a konzistence, které jsou nad rámec této kapitoly. Více o této problematice je možno se dozvědět v ([Vit87], [Ral78] a [Rek00]).

2.1. Formulace úlohy

Úlohou s počátečními podmínkami pro soustavu m diferenciálních rovnic prvního řádu

) , , , (

), , , , (

), , , , (

1 1 2 2

1 1 1

y y

x f y

y y

x f y

y y

x f y

m m

m

m m

K M

K K

′=

′=

′=

(2.1.1)

je nazývána úloha nalézt m funkcí ¹ proměnné x definovaných, spojitých a spojitě diferencovatelných v nějakém intervalu

y y, K, ^m

b

a, , které vyhovují rovnicím (2.1.1) pro každé b

a

x∈ , a pro něž platí

, ) ( , , ) ( , )

( ₁ ² ₂

1

m my

y

y ξ =η ξ =η K ξ =η (2.1.2)

kde ηr =(η₁,_K,η_m)^T je daný m-dimenzionální vektor a ξ je pevný bod z intervalu a,b . Doplňující podmínky (2.1.2) kladené na řešení soustavy (2.1.1) jsou nazvány počáteční podmínky.

Soustava (2.1.1) a počáteční podmínky (2.1.2) je zapsána stručně ve tvaru

, ) (

), , (

η ξ ^r r

r r r

=

′= y

y x f

y (2.1.3)

kde

,

2 1





















= y y y y

m

M r

^



















=

) , , , (

) , , , (

) , , , ( )

, (

1 1 2

1 1

y y

x f

y y

x f

y y

x f y

x f

m m

m m

L M

L L r r

2 .

1





















=

ηm

η η

η M

r

Protože každou diferenciální rovnici m-tého řádu ) , , , ,

( ′ ^{( −1)}

= ^m

m f x y y y

y K (2.1.4)

lze psát po zavedení nových neznámých funkcí rovnicemi

) 1 ( 2 1

, ,

= −

= ′

=

m

my y

y y

y y

L (2.1.5)

jako soustavu m rovnic prvního řádu

), , , , (

, , ,

1 1

3 2

2 1

y y

x f y

y y

y y

y y

m m

m m

L L

′=

′=

′=

′=

−

(2.1.6)

je vidět, co je úloha s počátečními podmínkami pro rovnici (2.1.4).

Popis konkrétních metod bude dále většinou omezen na jednu diferenciální rovnici )

, (x y f

y′= (2.1.7)

s počáteční podmínkou

η

= ) (a

y (2.1.8)

Pravá strana f diferenciální rovnice (2.1.7) musí být definovaná a spojitá v intervalu )

, ( ,

, ∈ −∞ ∞

∈ a b y

x a v tomto pásu splňovat tzv. Lipschitzovu podmínku vzhledem k y s konstantou nezávislou na x, tj. že existuje konstanta L taková, že platí

z y L z x f y x

f( , )− ( , ) ≤ − (2.1.9)

pro každé x∈ a,b a pro libovolná y a z. Ke splnění této podmínky stačí např., aby ve zmíněném pásu měla pravá strana dané diferenciální rovnice omezenou první parciální derivaci podle y. To je splněno např. v případě lineární diferenciální rovnice. Uvedené předpoklady o pravé straně dané diferenciální rovnice zaručují existenci a jednoznačnost řešení úlohy (2.1.7), (2.1.8) v celém intervalu a,b .

Z konkrétních metod pro řešení diferenciálních rovnic s počátečními podmínkami jsou v této kapitole podrobně popsány Rungovy–Kuttovy metody a lineární k-krokové metody včetně jejich použití způsobem prediktor–korektor.

2.2. Eulerova metoda. Problematika odhadu chyby

Eulerova metoda je jakýmsi základem všech ostatních numerických metod pro řešení úloh s počátečními podmínkami a je na ní možné jednoduše ukázat různá hlediska, která je třeba brát v úvahu při užití jakékoli metody. Přibližná hodnota řešení diferenciální rovnice (2.1.7) v bodě je počítána z rekurence

yn y(x) xn

x=

b a x N

n y x hf y y y

n n

n n

n ( , ), 0,1, , 1 ,

,

1 0

∈

−

= +

=

=

+ K

η (2.2.1)

Veličina h=(b−a)/N se nazývá integrační krok.

Geometrická interpretace Eulerovy metody je taková, že diferenciální rovnice (2.1.7) udává v intervalu x∈ a,b směrové pole, a nalézt její řešení splňující danou počáteční podmínku (2.1.8) značí tedy nalézt křivku, která prochází bodem (a,η) a která má v každém novém bodě směrnici určenou tímto směrovým polem. Eulerova metoda v této interpretaci tedy značí aproximaci hledané křivky lomenou čarou procházející body (x_i,y_i), i =0,1,K, přičemž směrnice úseček, které ji vytvářejí, souhlasí se směrovým polem vždy v levém krajním bodě.

Z popsané geometrické interpretace však nelze činit žádné závěry o chování rozdílu . Tento rozdíl je nazván celkovou nebo také akumulovanou diskretizační chybou. Aby jakákoli metoda pro řešení úlohy (2.1.7), (2.1.8) byla „rozumná“, je nezbytné, aby bylo možné učinit celkovou diskretizační chybu libovolně malou, tak jak je potřeba.

Jediným volným parametrem v metodě (2.2.1) je integrační krok h. Aby bylo možno učinit si představu o přesnosti numerického řešení, je potřeba mít k dispozici nějaké informace o chování celkové diskretizační chyby jako funkce integračního kroku h. Silnou informací je odhad chyby

) ( _n

n

n y y x

e = −

)

v ≤ψ(h

e , kde ψ je známá funkce, pro niž platí ψ(h)→0 pro . Pro danou metodu bývá však velmi obtížné odhad chyby nalézt a i když se to podaří, je možnost jeho praktického použití velmi malá, neboť při jeho konstrukci je třeba se orientovat na nejnepříznivější možný případ a to má za následek, že odhad chyby je značně pesimistický, tj. ve většině případů, kde je užit, dává mnohonásobně větší hodnotu, než je skutečnost.

→0 h

Další veličinou, která dává představu o celkové diskretizační chybě je lokální diskretizační chyba dané metody. To je chyba, která vzniká při provedení jednoho kroku dané

metody za předpokladu, že všechny hodnoty, které jsou k jeho realizaci potřeba, jsou přesné.

V konkrétním případě Eulerovy metody je vztah pro lokální diskretizační chybu dán výrazem )),

( , ( ) ( ) ( ) ), (

(y x h y x h y x hf x y x

L ≡ + − − (2.2.2)

kde y je přesné řešení dané diferenciální rovnice. Lokální dikretizační chyba numerického řešení dané diferenciální rovnice se může měnit v každém kroku a celková diskretizační chyba je dána součtem lokálních diskretizačních chyb, přičemž každý krok vychází z hodnot, které jsou zatíženy chybou z předešlého průběhu. Je tedy zapotřebí, aby daná metoda měla tu vlastnost, že v ní nedochází k výrazné akumulaci lokálních diskretizačních chyb. Metodu, u níž tomu tak je, je možno nazvat stabilní metodou.

Doposud bylo předpokládáno, že přibližné řešení splňuje rovnici Eulerovy metody (nebo rovnice složitější metody) přesně. Tak tomu vlivem zaokrouhlování (aritmetické operace je možno provádět pouze s konečným počtem čísel o konečném počtu cifer) není a místo veličiny je vypočítána ve skutečnosti veličina , která splňuje rovnici (2.2.1) pouze přibližně. Celková chyba numerického řešení se tedy skládá ze dvou částí: celkové diskretizační chyby a veličiny , která se nazývá celková zaokrouhlovací chyba.

Tato chyba vznikne nakupením lokálních zaokrouhlovacích chyb yn

yn y_n^*

n n

n y y

r = ^* −

εn definovaných rovnicí

n n n n

n y hf x y

y^*₊₁ = ^* + ( , ^*)+ε (2.2.3)

Dvě části chyby, kterou je zatíženo numerické řešení, se chovají velmi odlišně. Celková diskretizační chyba se zmenšováním h u „rozumných“ metod klesá. Naopak celková zaokrouhlovací chyba při zmenšování h roste. Při výpočtu pro velké h převládne celková diskretizační chyba a při zmenšování h se bude celková chyba zmenšovat. Existuje však mezní velikost integračního kroku , při které převládne celková zaokrouhlovací chyba, takže chyba při dalším zmenšování integračního kroku roste. Pro danou metodu a pro danou velikost lokální zaokrouhlovací chyby tedy existuje určitá mezní přesnost, kterou nelze překročit.

h0

2.3. Rungovy–Kuttovy metody 2.3.1. Obecná jednokroková metoda

Obecnou jednokrokovou metodou je jakýkoli algoritmus pro řešení úlohy (2.1.7), (2.1.8), v němž se přibližné řešení y_n+1 v bodě x_n+1 počítá pouze ze znalosti veličin x_n, y_n,h.

Eulerova metoda je tedy speciálním případem obecné jednokrokové metody. Závislost přibližného řešení y_n+1 na veličinách x_n, y_n,h je užitečné zapsat ve tvaru

, (x hΦ _n +

) ( ) y x h

( (y

L −

( ) );h =O h

)

; (x h − y

)− h

(h^p+1

Φ

) ,

)(

1

( x y +

) h =

) , y

) , ,

f y

x + ∂

, (

, (

) (

) 0 (

x x

s

),

1 y y ,h

y_n+ = _{n n} (2.3.1)

kde Φ je funkce tří proměnných, která závisí na dané diferenciální rovnici.

O chování celkové diskretizační chyby obecné jednokrokové metody rozhoduje v podstatě chování lokální diskretizační chyby. Lokální diskretizační chyba obecné jednokrokové metody (2.3.1) je dána výrazem

), ), ( , ( (

)

);h y x h x y x h

x = + − Φ (2.3.2)

kde y je přesné řešení problému (2.1.7), (2.1.8). Největší celé číslo p pro než platí ),

(

(y x ^p+1

L (2.3.3)

kde číslo p se nazývá řád metody, který závisí na hladkosti řešení dané diferenciální rovnice.

Odhad celkové diskretizační chyby má stejné negativní vlastnosti jako odhad celkové diskretizační chyby Eulerovy metody a to zejména proto, že je velice pesimistický a je složité ho získat. Proto se v praxi častěji využívá odhadu chyby metodou polovičního kroku

[

( ; /2)

]

( ),

1 2 ) 2 (

;

( + ⁺¹

= _p −^p y x h O h^p

x y x

y (2.3.4)

kde je přibližné řešení v bodě x získané obecnou jednokrokovou metodou užitou s krokem h. Tento odhad chyby je zkonstruován z tvrzení o asymptotickém chování chyby.

V praxi se člen na pravé straně vztahu (2.3.4) zanedbává.

)

; ( hx y

) O

V následujícím textu jsou uvedeny některé konkrétní jednokrokové metody, které vzniknou speciální volbou funkce .

2.3.2. Metoda Taylorova rozvoje

V této speciální jednokrokové metodě je funkce Φ zapsána ve tvaru , ) ,

! ( 1 2

) 1 , ( ,

,

( h ¹f^{( 1)} x y

hf p y

x f y

x + + ^{p− p−}

Φ K (2.3.5)

kde funkce f^(s)(x je definována rekurentně vzorci

. ) , ) ( , ( ) (

) , ( )

) 1 ( )

1 (

y x y f

y x x

y f f

y x f y f

s s

∂

∂

= ∂

=

−

−

Funkce je rovna (s+1)-ní derivaci funkce y(x), která je řešením diferenciální rovnice (2.1.7). Lokální chyba metody (2.3.5) je rovna veličině , z čehož plyne, že se jedná o metodu řádu p. Metoda Taylorova rozvoje umožňuje konstruovat obecnou jednokrokovou metodu libovolně vysokého řádu. Přesto se jí v praxi užívá jen velice zřídka. Hlavní důvod je ten, že pro výpočet hodnoty funkce

)) ( ,

)(

( x y x

f ^s

)!

1 /(

)

( ¹

) 1

( ⁺ h ⁺ p+

y ^p ξ ^p

Φ je nutno derivovat pravou stranu dané diferenciální rovnice, což často vede k velice komplikovaným výrazům.

2.3.3. Klasické Rungovy–Kuttovy metody

Rungovy–Kuttovy metody jsou nejčastěji užívané obecné jednokrokové metody a hlavní nevýhoda metody Taylorova rozvoje se u nich odstraňuje tak, že funkce Φ je ve tvaru

, )

, ,

(x y h =w₁k₁ + w₂k₂ + +w_sk_s

Φ K (2.3.6)

kde

. , , 2 ,

, ) , (

1 1 1

s i

h y h x f k

y x f k

i

j ij

i

i  = K





 + +

=

=

∑

⁻

=

β α

Číselné parametry w_i,α_i,β_ij se volí tak, aby se funkce Φ lišila od pravé strany rovnice (2.3.5) až teprve ve členech řádu h^p, tj. tak, aby vzniklá metoda byla řádu p. Při výpočtu hodnoty funkce Φ je tedy nutnost počítat derivace pravé strany dané diferenciální rovnice nahrazena tím, že je třeba počítat hodnoty funkce f ve více bodech.

Maximální dosažitelný řád Rungovy–Kuttovy metody p závisí na čísle s tak, že pokud

∞pak

→

s p→∞, kromě toho platí, že p = s pro s = 1, 2, 3, 4, p = 4 pro s = 5, p = 5 pro s = 6, p = 5 pro s = 7, p = 6 pro s = 8, , a řádu p lze dosáhnout nekonečně mnoha volbami příslušných parametrů.

K

V následujících řádcích jsou uvedeny některé důležité konkrétní Rungovy–Kuttovy metody. Metoda druhého řádu



 



 + +

+

+ = ( , )

2 , 1 2 1

1 n n n n n

n y hf x h y hf x y

y (2.3.7)

a metoda

[

( , ) ( , ( , ))

]

,

1 n n n n n n n

n y hf f x y f x h y hf x y

y ₊ = + + + + (2.3.8)

která je rovněž druhého řádu. Obě se v literatuře vyskytují pod názvem modifikovaná Eulerova metoda.

Dvě velmi známé metody čtvrtého řádu jsou

, ) ,

(

2 , , 1 2 1

2 , , 1 2 1 , ) , (

, ) 2

2 6 (

1

3 4

2 3

1 2

1

4 3 2 1 1

hk y h x f k

hk y

h x

f k

hk y

h x

f k

y x f k

k k k k h y

y

n n

n n

n n

n n n n

+ +

=



 



 + +

=



 



 + +

=

=

+ + + +

+ =

(2.3.9)

a

. ) ,

(

3 , , 1 3 2

3 , , 1 3 1 , ) , (

, ) 3

3 8 (

1

3 2 1 4

2 1 3

1 2

1

4 3 2 1 1

hk hk hk y h x f k

hk hk y

h x

f k

hk y

h x f k

y x f k

k k k k h y

y

n n

n n

n n

n n n n

+

− + +

=



 



 + − +

=



 



 + +

=

=

+ + + +

+ =

(2.3.10)

Daleko nejužívanější Rungovou–Kuttovou metodou je metoda (2.3.9). Je to jakási standardní Rungova–Kuttova metoda.

Rungovy–Kuttovy metody vyšších řádů se používají zřídka, přesto budou uvedeny alespoň dva příklady takových metod. Prvním příkladem je tzv. Nyströmova metoda, která je pátého řádu, a vyžaduje výpočet šesti funkčních hodnot a je dána rovnicemi

. ) 8 10 36

6 75 ( , 1

5 4

, ) 8 50 90

6 81 ( , 1

3 2

, ) 15 12

4 ( , 1

, ) 6 4 25 ( , 1

5 2

3 , , 1 3 1 , ) , (

, ) 125 81

125 23

192 ( 1

4 3 2

1 6

4 3 2

1 5

3 2

1 4

2 1 3

1 2

1

6 5

3 1

1



 



 + + + + +

=



 



 + + + − +

=



 



 + + − +

=



 



 + + +

=



 



 + +

=

=

+

− +

+

+ =

k k k

k h y

h x

f k

k k k

k h y

h x

f k

k k

k h y

h x f k

k k h y

h x

f k

hk y

h x f k

y x f k

k k

k k

h y

y

n n

n n

n n

n n

n n

n n n n

(2.3.11)

a druhým Huťova metoda

, ) 72 9

454 834

1002 2079

716 82 ( , 1

, ) 3 80 66

472 678

183 48 (

, 1 6 5

, ) 102

867 981

221 9 ( , 1 3 2

, ) 6 24 27

5 8 ( , 1 2 1

, ) 4 3 6 (

, 1 3 1

, ) 3 24 (

, 1 6 1

9 , , 1 9 1 , ) , (

, ) 41 216

27 272

27 216

41 840 (

1

7 6

5 4

3 2

1 8

6 5 4

3 2

1 7

5 4 3

2 1

6

4 3 2

1 5

3 2 1 4

2 1 3

1 2

1

8 7

6 5

4 3

1 1



 



 + + − + + − − +

=



 



 + + − + − − + +

=



 



 + + − + − +

=



 



 + + − + − +

=



 



 + + − +

=



 



 + + +

=



 



 + +

=

=

+ +

+ +

+ +

+

+ =

k k

k k

k k

k h y

h x f k

k k k

k k

k h

y h x

f k

k k k

k k

h y

h x

f k

k k k

k h y h x

f k

k k k h y

h x f k

k k h y

h x f k

hk y

h x

f k

y x f k

k k

k k

k k

k h y

y

n n

n n

n n

n n

n n

n n

n n

n n n n

(2.3.12)

která je šestého řádu a vyžaduje v každém kroku osmkrát vypočítat hodnotu pravé strany dané diferenciální rovnice.

Obě tyto metody nejsou užívány často, ačkoliv k tomu není žádného vážného důvodu.

Jejich nevýhodou je, že k provedení jednoho kroku je potřeba více početních operací (výpočet pravé strany), avšak jsou vysokého řádu, takže je lze užít s větším integračním krokem.

2.3.4. Odhad chyby Rungových–Kuttových metod. Automatická volba kroku Pro posouzení chyby, kterou je řešení diferenciální rovnice zatíženo, je třeba mít k dispozici informace o velikosti lokální diskretizační chyby. Jako příklad je uveden odhad lokální diskretizační chyby klasické Rungovy–Kuttovy metody čtvrtého řádu (2.3.9). Pokud ve vhodném okolí bodu (x_n,y_n), v němž je lokální chyba odhadnuta, platí

p j M i

L y x M f

y x

f ⁱ_i ^j _i < ⁱ_j^j + ≤

∂

∂

< , ∂^{+ +}₋ , )

,

( ₁ (2.3.13)

pak pro lokální chybu Rungovy–Kuttovy metody čtvrtého řádu (2.3.9) platí 720 .

) 73

; (

(y x h h⁵ML⁴

L _n = (2.3.14)

Uvedený odhad je, jak je vidět, značně komplikovaný a není nikterak jednoduché jej získat.

Mimoto jsou tyto odhady většinou značně pesimistické a volba integračního kroku na nich založená vede k neefektivním algoritmům.

S problematikou odhadu lokální diskretizační chyby úzce souvisí možnost automatické volby integračního kroku. Jednokrokový charakter Rungových–Kuttových metod po technické stránce neklade žádné překážky změně velikosti integračního kroku, třeba i v každém kroku metody. Je proto výhodné řídit jeho velikost velikostí lokální diskretizační chyby. Komplikované výrazy, které byly pro její odhad uvedeny, jsou k tomuto účelu zřejmě nevhodné. Proto se častěji používají aposteriorní odhady, tj. odhady, které jsou vypočítány zároveň s přibližným řešením diferenciální rovnice.

Nejjednodušším způsobem získání odhadu lokální diskretizační chyby pro obecnou jednokrokovou metodu je metoda „zdvojování kroku“, která k odhadu lokální chyby využívá metodou polovičního kroku (2.3.4). Při výpočtu přibližného řešení diferenciální rovnice je vypočítáno řešení v každém kroku dvakrát. Jednou s integračním krokem velikosti h, tím je získáno přibližné řešení, které je označeno jako y1. Pro výpočet druhého přibližného řešení y2

jsou použity dva kroky o velikosti h/2. Takto získaná řešení y1 a y2 mají tu vlastnost, že jejich rozdílem je získán odhad lokální diskretizační chyby, podle které je pak možné přizpůsobit velikost integračního kroku pro další výpočet. Nevýhodou toho způsobu je velké množství výpočtů pravé strany diferenciální rovnice. Zatímco např. u standardní Rungovy–Kuttovy metody je pravá strana vypočítána čtyřikrát, tak u zdvojování kroku je třeba znát jedenáct hodnot pravé strany diferenciální rovnice.

Jiným a výhodnějším způsobem je užití Rungovy–Kuttovy metody s více hodnotami pravé strany, než je zapotřebí k dosažení daného řádu. To umožní počítat nejen přibližné řešení, ale také i jeho lokální chybu.

Takovou metodu objevil Fehlberg, který vytvořil metodu pátého řádu

40 , 11 4104

1859 2565

2 3544 27

, 8 2 1

4104 , 845 513

8 3680 216

, 439

, ) 7296 7200

1932 2197 (

, 1 13 12

, ) 9 3 32 ( , 1

8 3

4 , , 1 4 1 , ) , (

55 , 2 50

9 56430

28561 12825

6656 135

16

5 4

3 2

1 6

4 3

2 1 5

3 2

1 4

2 1 3

1 2

1

6 5

4 3

1 1



 



 



 



− + − + −

+ +

=



 



 



 



 − + −

+ +

=



 



 + + − +

=



 



 + + +

=



 



 + +

=

=



 



 + + − +

+

+ =

k k

k k

k h

y h x

f k

k k

k k h y h x f k

k k

k h y

h x

f k

k k h y

h x f k

hk y

h x

f k

y x f k

k k

k k

k h y y

n n

n n

n n

n n

n n

n n n n

(2.3.15)

tak, že šest hodnot k v ní užitých lze zkombinovat tak, že 5 , 1 4104 2197 2565

1408 216

25

5 4 3

1

*

1 



 



 + + −

+

+ = y h k k k k

y_{n n} (2.3.16)

je jiná aproximace řešení, která je čtvrtého řádu. Číslo lze pak pokládat za odhad lokální diskretizační chyby přibližného řešení .

* 1

1 +

+ − _n

n y

y

* +1

yn

Jiným příkladem metody podobné Rungově–Kuttově–Fehlbergově metodě je metoda čtvrtého řádu

, ) 378 54

546 125

28 625 ( , 1

5 1

, ) 10

7 27 ( , 1

3 2

, )) 2 (

, (

, ) 4 (

, 1 2 1

2 , , 1 2 1 , ) , (

, ) 4

6 ( 1

5 4

3 2

1 6

4 2 1

5

3 2 4

2 1 3

1 2

1

4 3 1 1



 



 + + − + + −

=



 



 + + + +

=

+

− + +

=



 



 + + +

=



 



 + +

=

=

+ + +

+ =

k k

k k

k h y

h x

f k

k k k

h y

h x

f k

k k h y h x f k

k k h y

h x

f k

hk y

h x

f k

y x f k

k k k h y

y

n n

n n

n n

n n

n n

n n n n

(2.3.17)

s odhadem lokální chyby

) 125 162

21 224

42 336 ( ) 1 );

(

(y x h h k₁ k₃ k₄ k₅ k₆

L = − − − + + (2.3.18)

Všimněme si, že pokud není vyžadován odhad chyby, vyžaduje užití vzorce (2.3.17) pouze čtyři hodnoty pravé strany dané diferenciální rovnice.

Další metodou je inovovaná Rungova–Kuttova-Felhbergova metoda s koeficienty od pánů Cashe a Karpa, která by měla být dokonalejší než původní Rungova–Kuttova- Fehlbergova metoda. Je to opět metoda pátého řádu.

4096 , 253 110592

44275 13824

575 512

175 55296

, 1631 8 7

27 , 35 27

70 2

5 54 , 11

, ) 12 9

3 10 ( , 1

5 3

, ) 9 3 40 ( , 1

10 3

5 , , 1 5 1 , ) , (

1771 , 512 594

125 621

250 378

37

5 4

3 2

1 6

4 3

2 1 5

3 2

1 4

2 1 3

1 2

1

6 4

3 1

1



 



 



 



 + + + +

+ +

=



 



 



 



− + − +

+ +

=



 



 + + − +

=



 



 + + +

=



 



 + +

=

=



 



 + + +

+

+ =

k k

k k

k h

y h x

f k

k k

k k h

y h x f k

k k

k h y

h x f k

k k h y

h x

f k

hk y

h x f k

y x f k

k k

k k

h y y

n n

n n

n n

n n

n n

n n n n

(2.3.19)

pomocné řešení čtvrtého řádu je získáno ze vzorce

4 . 1 14336

277 55296

13525 48384

18575 27648

2825

6 5

4 3

1

*

1 



 



 + + + +

+

+ = y h k k k k k

y_{n n} (2.3.20)

Jako další metoda je uvedena metoda třetího řádu, která nese po svých tvůrcích jméno Bogacki–Sampline.

, ) 4 3 2 9 ( , 1

4 , , 3 4 3

2 , , 1 2 1 , ) , (

, ) 4 3 2 9 ( 1

3 2 1 4

2 3

1 2

1

3 2 1 1



 



 + + + +

=



 



 + +

=



 



 + +

=

=

+ + +

+ =

k k k h y

h x f k

hk y

h x

f k

hk y

h x

f k

y x f k

k k k h y

y

n n

n n

n n

n n n n

(2.3.21)

pomocné řešení čtvrtého řádu je získáno ze vzorce 8 . 1 3 1 4 1 24

7

4 3 2 1

*

1 



 



 + + +

+

+ = y h k k k k

y_{n n} (2.3.22)

Na závěr je popsána metoda pátého řádu od pana Vernera



 



 



 



 − + − + +

+ +

=



 



 



 



− + − − +

+ +

=



 



 



 



 − + − +

+ +

=



 



 



 



− + − +

+ +

=

+

− +

+

=



 



 + + +

=



 



 + +

=

=



 



 + + + + +

+

+ =

7 5

4 3

2 1

8

5 4

3 2

1 7

5 4

3 2

1 6

4 3

2 1

5

3 2

1 4

2 1

3

1 2

1

8 7

5 4

3 1

1

26703 3850 84065

24068 2322

319 52632

297275 43

300 1720

, 3501

10625 , 2484 250

81 680

643 75

125 15000

, 8263 15

1

255 , 88 36

11 612

8 4015 5

, 12

96 , 85 64

425 6

55 64

, 165 6 5

) 23 . 3 . 2 ( ,

)) 15 16

5 6 ( , 1 3 ( 2

, ) 16 4 75 ( , 1

15 4

6 , , 1 6 1 , ) , (

616 , 43 11592

125 1955

264 72

23 2244

875 40

3

k k

k k

k k

h y h x f k

k k

k k

k h

y h x

f k

k k

k k

k h y h x f k

k k

k k

h y h x

f k

k k

k h y

h x

f k

k k h y

h x

f k

hk y

h x

f k

y x f k

k k

k k

k k

h y y

n n

n n

n n

n n

n n

n n

n n

n n n n

Pomocné řešení pátého řádu je dáno vzorcem

44 . 3 85

12 16

5 5984

2375 160

13

6 5

4 3

1

1 



 



 + + + +

+

+ = y h k k k k k

y_{n n} (2.3.24)

Na závěr je věnována pozornost samotnému určení ideální velikosti integračního kroku ze znalosti odhadu velikosti lokální diskretizační chyby. Nejčastěji se používá vzorec

,

1

1 *

p

n n n

n y y

h Tol

h 











= −

+ (2.3.25)

kde p je konstanta volená tak, že p−1 je rovno řádu pomocného řešení diferenciální rovnice, potom p je ve většině případů rovno řádu dané metody, ale existují i jiné volby tohoto parametru. Konstanta Tol určuje, s jakou přesností je výpočet numerického řešení prováděn.

V některých případech je potřeba, aby změna velikosti integračního kroku nebyla náhlá (skoková), ale aby se měnila pomaleji.

Takovou volbu integračního kroku zajišťuje následující vztah ,

* 1 1 1 *

β α











 −













= − ^{− −}

+ Tol

y y y

y h Tol

h ^{n n}

n n n

n (2.3.26)

kde konstanty α a β jsou voleny např. podle Gustaffsona

5 . 2 10 , 7

p p

=

= β α

(2.3.27)

Vztah (2.3.26) lze obohatit o bezpečnostní konstanty ,

2

*

* 2 1 1

β α











 −













= −

+ S Tol

y y y

y Tol S S

h

h ^{n n}

n n n

n (2.3.28)

kde konstanta je absolutní a volí se z intervalu (0, 1), nejčastěji se její velikost pohybuje okolo 0.8. Velikost konstanty je závislá na řádu metody a na řešené úloze. V některých případech bývá velikost konstanty volena rovna velikosti integračního kroku h. Další běžnou modifikací vztahu (2.3.25) je

S1

S2

S2

* 1

* 2 1 1

* 1

1

* 2 1

1 ,

n n p

n n n n

n n p

n n n n

y y Tol y pro

y Tol S S

h h

y y Tol y pro

y Tol S S

h h

−

 <













= −

−

 >













= −

+

+ +

(2.3.29)

Z předchozího textu je patrné, že automatická volba integračního kroku je závislá na použité metodě .

2.4. Lineární k-krokové metody

V předchozím textu bylo uvedeno, že výpočet přibližného řešení dané diferenciální rovnice jednokrokovou metodou probíhá tak, že po provedení jednoho kroku je řešen nový problém s počáteční podmínkou zadanou v bodě, do kterého řešení předtím došlo. Výhodou tohoto postupu je, že se taková metoda jednoduše programuje a změna integračního kroku je formálně jednoduchá apod. Nevýhodou může být to, že v tom okamžiku, kdy je vypočítáno přibližné řešení v jednom bodě, jsou zapomenuty všechny informace o řešení v předcházejících bodech. Existují však také metody, které k získání přibližného řešení užívají

nejen informace z bezprostředně předchozího bodu, ale i informace z několika předchozích bodů a tyto metody se nazývají k-krokové metody.

K-korkovou nebo mnohokrokovou metodou je označena metoda daná vzorcem

∑ ∑

= + = = + =

k k

n

n h f n

y

0 0

, , 1 , 0 ,

ν αν ν ν βν ν _{L (2.4.1)}

kde α₀,K,α_k a β₀,K,β_k jsou konstanty a stručný zápis funkční hodnoty . Všude v dalším výkladu bude předpokládáno, že

ν +

fn

) , _ν

ν +

+ yn

(x_n

f α_k =1 a že čísla α₀ a β₀

nejsou současně rovna nule. Z rovnice (2.4.1) má být určena hodnota za předpokladu, že hodnoty jsou už známé a tento postup opakovat pro n = 0, 1, … tak dlouho, až je dosaženo bodu, kde je přibližné řešení požadováno. Z toho je patrné, že k použití vzorce (2.4.1) je třeba znát přibližné řešení v bodech

. Mnohokroková metoda (2.4.1) tedy vyžaduje k počátečních podmínek, které musí být před započetím výpočtu podle rovnice (2.4.1) nějakým způsobem získány, to znamená, že tato metoda není „samostartující“, na rozdíl od jednokrokových metod.

K získání počátečních hodnot se často používá klasická Rungova–Kuttova metoda.

n+k

, K y

1 0, y y , 1

,K y_n+k−

+1

yn n, y

1 1,K,x_k−

,y_k−1 0, x

x

V případě, že je β_k =0, je rovnicí (2.4.1) přibližné řešení jednoznačně určeno a v tomto případě se jedná o explicitní metodu. Obecně, když

k

yn₊

≠0

βk , se jedná o implicitní metodu a rovnice (2.4.1) představuje obecně nelineární rovnici pro určení y_n+k.

Lokální diskretizační chyba je popsána rovnicí

∑ ∑

= =

′ +

− +

= ^k

i

k

i i

iy x ih h y x ih

h x y l

0 0

, ) ( )

( )

);

(

( α β _(2.4.2)

kde y je přesné řešení diferenciální rovnice. Jestliže je tato funkce dostatečně hladká, pak ji můžeme zapsat ve tvaru

K K

K K

K K

, 3 , 2 ,

) )!(

1 ( ) 1

!( 1

,

, )

( )

( )

( )

);

( (

1 1

1 0 0

) ( 1

0

= +

− +

− +

+

=

+ +

=

+ +

′ + +

=

− p

p k p k

C C

h x y C h

x y C x y C h x y l

k p k

p p

k

p p p

β β

α α

α

α (2.4.3)

Potom platí, že lineární mnohokroková metoda je řádu p, když platí .

0 ,

0 ₁

1

0 =C = =C_p = C_p+ ≠

C K

2.4.1. Metody numerické integrace

Metody numerické integrace jsou nejpoužívanějšími k-krokovými metodami. Zde jsou popsány algoritmy pro získání řešení diferenciální rovnice pomocí explicitní Adamsovy–

Bashforthovy a Nyströmovy metody a implicitní Adamsovy–Moultonovy metody. Tyto metody vycházejí z rovnice

∫

=

− ^r

s

x

x s

r y x f t y t dt

x

y( ) ( ) (, ( )) . (2.4.4.)

Explicitní Adamsova–Bashforthova metoda pro řešení úlohy (2.1.7), (2.1.8) je pro jednotlivá k (podle počtu kroků) dána rovnicemi

, 2 ,

)) , ( ) , ( 3 2 ( 1

1 1

1 = + − _{− −} =

+ y h f x y f x y k

y_{n n n n n n} (2.4.5)

, 3 ,

)) , ( 5 ) , ( 16 ) , ( 23 12 (

1

2 2 1

1

1 = + − _{− −} + _{− −} =

+ y h f x y f x y f x y k

y_{n n n n n n n n} (2.4.6)

Dále u většího počtu kroků je postupováno analogicky a jsou uvedeny jen koeficienty ,

4 ,

) 9 , 37 , 59 , 55 24 (

1 h − − k = (2.4.7)

, 5 ,

) 251 , 1274 ,

2616 , 2774 ,

1901 720 (

1 h − − k = (2.4.8)

. 6 ,

) 425 , 2627 , 6798 ,

9482 , 7673 ,

4227 1440 (

1 h − − − k = (2.4.9)

Implicitní Adamsova–Moultonova metoda pro řešení úlohy (2.1.7), (2.1.8) je pro jednotlivá k (podle počtu kroků) dána rovnicemi

, 1 ,

)) , ( ) , ( 2 ( 1

1 1

1 = + _{+ +} + =

+ y h f x y f x y k

y_{n n n n n n} (2.4.10)

, 2 ,

)) , ( ) , ( 8 ) , ( 5 12 (

1

1 1 1

1

1 = + _{+ +} + − _{− −} =

+ y h f x y f x y f x y k

y_{n n n n n n n n} (2.4.11)

Dále u většího počtu kroků je postupováno analogicky a jsou uvedeny jen koeficienty ,

3 ,

) 1 , 5 , 19 , 9 24 (

1 h − k = (2.4.12)

, 4 ,

) 19 , 106 , 264 , 646 , 251 720 (

1 h − − k = (2.4.13)

. 5 ,

) 27 , 173 , 482 , 798 , 1427 , 475 1440 (

1 h − − k = (2.4.14)