Martin Isoz, st.sk 204
Matematické metody v inženýrství
Projekt M3
8.5.2010
Řešení diferenciálních rovnic
1. Zadání A. Stanovte řešení dané diferenciální rovnice popřípadě soustavy rovnic. i) Pro úlohy M3.1 až M3.12: •
uveďte matematický popis použité metody
•
sestavte vlastní program pro řešení dané metody numericky Eulerovou metodou pro různé kroky výpočtu včetně specifikace rovnice ve formě zvláštního podprogramu
•
podle možností proveďte rovněž symbolické řešení
•
porovnejte vlastní řešení s výsledky určenými pomocí programů ODE23 a ODE45
ii) Pro úlohy M3.21 až M3.32 •
uveďte matematický popis použité metody
•
podle možností proveďte rovněž symbolické řešení
•
uveďte schema a výsledky řešení dané úlohy v SIMULINKu
Při řešení použijte v návaznosti na zadání některé z funkcí ODE23, ODE45, PLOT a dále DSOLVE a EZPLOT. B. Pomocí vlastního programu stanovte symbolicky a numericky hodnotu derivace a integrálu dané funkce f x v mezích 〈 a , b 〉 při dělení na N dílků. V rámci řešení: •
uveďte matematický popis použitých metod
•
pro numerickou integraci zvolte některé ze základních pravidel (obdélníkové, lichoběžníkové, Simpsonovo) a porovnejte vlastní výsledky s řešením pomocí programů QUAD a QUAD8
•
pro numerickou derivaci zvolte vhodnou diferenční formuli graficky znázorněte průběh integrálu a derivace dané zvolené funkce v daném intervalu Při řešení použijte některé z funkcí QUAD, QUAD8 a dále INT, DIFF a EZPLOT
•
Zadání M3.25: y ' ' 4 y=cos t y 0=0
•
y ' 0 =1
(2)
Zadání M3.44 f x =3 x 25
1
(1)
(3)
Martin Isoz, st.sk 204
Matematické metody v inženýrství
8.5.2010
2. Popis řešení Zadaný problém byl řešen ve výpočetním prostředí MATLAB® 7.6.0 (R2008a). Při řešení problému M3.25 byla daná ODR s počátečními podmínkamu řešena symbolicky a následně s užitím modelu SIMULINKU. Hodnoty řešení nalezené pomocí sestrojeného modelu byly exportovány do pracovního prostředí a následně vyneseny do grafu společně se závislostí nalezenou řešením symbolickým. Nakonec byla spočtena a graficky znázorněna absolutní a relativní chyba numericky nalezeného řešení. Problém M3.44 byl také řešen nejprve symbolicky a následně numericky. Vypočtené hodnoty řešení byly vypsány do COMMAND WINDOW a graficky znázorněny spolu s absolutními a relativními odchylkami jednotlivých způsobů řešení.
3. Teoretická část 3.1 Numerické metody řešení ODR s počátečními podmínkami K řešení ODR bylo odvozeno velké množství metod, které se liší náročností i chybou generovanou během výpočtu. Nejjednodušší numerickou metodou řešení diferenciálních rovnice je metoda Eulerova, která je ale zatížena velkou chybou (chyba této metody je řádu O h ). Přesto na ní lze ukázat základní principy jednokrokových metod řešení ODR. Pro jednoduchost bude nastíněna varianta této metody pro řešení jedné diferenciální rovnice 1. řádu. Budeme tedy řešit rovnici y' = f x , y
(4)
y a=c
(5)
s počáteční podmínkou
Funkci y x budeme hledat jako diskrétně zadanou funkci, tj na množině bodů (uzlů sítě) a=x 0 x 1... x N =b budeme hledat čísla y 0=c , y 1, ... y N , která aproximují hodnoty y x 0 , y x 1 ,... y x N přesného řešení v bodech x 0 , x1 ,... , x N . Zvolíme ekvidistantní síť, tj. x n1− x n=h; n=0,1 , ... , N −1 , kde h je krok metody. Nahrazením derivace y' v bodě x = xn pomocí dvoubodové dopředné diferenční formule, y ' x n ≃
y n 1− y n h
(6)
získáme Eulerovu metodu y n1= y n h⋅ f x , y
(7)
Eulerova metoda patří do skupiny metod zvaných Runge-Kuttovy. To jsou metody, u kterých lze přírůstek hledané funkce pro y n~ y x n napsat ve tvaru y n1= y n h⋅ x n , y n ; h (8) pro výše zmíněnou Eulerovu metodu je zřejmě x n , y n ; h= f x , y . Různými úpravami funkce x n , y n ; h získáváme jednotlivé metody. Metody Runge-Kutta mají tu vlastnotst, že pro dosažení řádu chyby O hm je pro m≤4 třeba právě m dosazení do pravé strany diferenciální 2
Martin Isoz, st.sk 204
Matematické metody v inženýrství
8.5.2010
rovnice (8) a pro m4 alespoň m+1 dosazení. Z tohoto důvodu jsou metody řádů vyšších než 4 obvykle málo používány, jelikož jejich přednosti by se projevily až při extrémních nárocích na přesnost výpočtu. V prostředí programu Matlab jsou k dispozici funkce ODE23 a ODE45, reprezentující algoritmy řešení ODR s počátečními podmínkami pomocí metody Runge-Kutta 2.-3., respektive 4.5. řádu. Oba tyto programy jsou vcelku komplexní – disponují například automatickou volbou kroku. Ovšem vzhledem k faktu, že jednokrokové metody obecně nejsou A-stabilní nejsou tyto algoritmy vhodné k řešení stiff systémů ODR.
3.2 Symbolické řešení diferenciálních rovnic Rovnici (1) je možno symbolicky řešit několika způsoby. Jedná se o nehomogenní LDR 2. řádu s konstantními koeficienty a speciální pravou stranou – lze ji zapsat do tvaru: ax
(9)
0 y ' '1 y ' 2 y=e P x cos bxQ x sin bx
kde 0 ,1 , 2 , a , b∈ℝ a P,Q jsou dané polynomy. Rovnici (1) je také možno zapsat ve tvaru: 2
2
y ' '2 a y ' a b y= f t
(10)
s a , b∈ℝ Řešení takovéto rovnice je možné metodou variace konstant, ale výhodnější je hledat řešení pomocí metody odhadu či Laplaceovy transformace. 3.2.1 Metoda odhadu
Obecné řešení nehomogení LDR (NLDR) je možno zapsat jako součet obecného řešení přiřazené HLDR a partikulárního řešení NLDR. Toto budeme dále značit jako y ON = y OH y PN
(11)
v případě metody odhadu hledáme y PN ve tvaru k ax y PN = x e R x cos bxS x sin bx ,
(12)
kde čísla a, b jsou dána pravou stranou rovnice (9) a k násobnost =aib , kořene charakteristické rovnice přiřazené HLDR. V konkrétním případě rovnice (1) je tedy v rovnici (9) 0=1, 1=0, 2 =4, a=0, b=1 , P(t) = 1 a Q(t) = 0. Charakteristická rovnice přiřazené HLDR má tvar 24=0
(13)
a kořeny 1,2=±2 i . Obecné řešení přiřazené HLDR je tedy: y OH =C 1 cos 2 tC 2 sin 2 t , C 1 , C 2 ,t ∈ℝ .
3
(14)
Martin Isoz, st.sk 204
Matematické metody v inženýrství
8.5.2010
Partikulární řešení rovnice (1) budeme ze vztahu (12) hledat ve tvaru y PN =t 0 e 0t A⋅cos tB⋅sin t .
(15)
Koeficienty A, B určíme metodou neurčitých koeficientů dosazením do rovnice (1). Po úpravách 1 získáme A= , B=0 . Obecné řešení rovnice (1) je tedy: 3 1 y ON =C 1 cos 2 tC 2 sin 2 t cos t , C 1 , C 2 ,t ∈ℝ 3
(16)
1 1 Dosazením do počátečních podmínek (2) získáme hodnoty konstant C 1=− ,C 2= . 3 2 Řešením rovnice (1) s podmínkami (2) tedy je funkce: 1 1 1 y=− cos 2 t sin 2 t cos t 3 2 3
t∈ℝ
(17)
3.2.2 Laplaceova transformace
Diferenciální rovnici ve tvaru (10) je možno s výhodou řešit užitím Laplaceovy transformace. Vztahy mezi vzorem f(t) a obrazem F(p) této transformace jsou: ∞
F p=∫ f t e− pt dt , t∈ℝ , p∈ℂ 0
1 2 i
(18)
x i ∞
∫
F pe pt dp
(19)
{dtd y t }= pY p− y 0
(20)
f t=
x− i ∞
Jelikož platí: L
{
}
d2 L 2 y t = p 2 Y p− py 0− y ' 0 dt
(21)
můžeme rovnici (1) a poč. podmínku (2) s užitím vztahů (10), (18), (20) a (21) přepsat do tvaru p p 2 Y p−14 Y p= 2 (22) p 1 vyjádřením Y(p) z rovnice (22) získáme po úpravě vztah p 1 Y p= p−i pi p−2 i p2 i p−2 i p2 i
(23)
a inverzní transformací podle vzorce (19) získáme konečné řešení rovnice (1) s podmínkami (2): 1 1 1 y=− cos 2 t sin 2 t cos t , t∈ℝ (24) 3 2 3 V praxi je transformace podle vztahů (18) a (19) zdlouhavá a proto existují tabulky s předpočtenými transformacemi nejběžnějších funkcí. 4
Martin Isoz, st.sk 204
Matematické metody v inženýrství
8.5.2010
3.3 Numerická integrace a derivace 3.3.1 Diferenční formule
V případě, že je předpis funkce příliš složitý, či vůbec neexistuje, je vhodné nahradit analytickou derivaci nějakou aproximací. Základ myšlenky diferenčních formulí tkví v nahrazení funkčního předpisu interpolačním polynomem x a následné derivaci tohoto polynomu. V případě aproximace funkce Lagrangeovým polynomem získáme takto vztah n
L ' n x=∑ f i=0
x i ' n x n x − . ' n x x−x i x− x i 2
(25)
Tento vztah je však pro praktické výpočty, zvláště derivací vyšších řádů, nevhodný. Ovšem při zavedení ekvidistantní sítě uzlů s diferencí h, x i1− x i=h , i=0,1 , ... , n−1 lze v případě, že chceme derivaci vyčíslit v bodě x j ∈〈 x 0 , x n 〉 vztah (25) zapsat zjednodušeně jako n
L ' n x j =
1 ∑C f , h i=0 ji i
(26)
kde f i= f x i . Pro derivace vyšších řádů lze postupovat analogicky. Například pro druhou derivaci f(x) získáme po úpravách vztah n
1 L ' ' n x j = 2 ∑ D ji f i . h i=0
(27)
Při studiu chyby zjistíme, že nejmenší odchylky od přesných hodnot při zachování jednoduchosti dosahují symetrické formule. Tyto jsou také v praxi nejčastěji používány. 3.3.2 Kvadraturní formule
K užití aproximativních metod výpočtu určitého integrálu se uchylujeme ze stejných důvodů jako u zavádění diferenčních formulí. I princip aproximace je velmi podobný. Nahradíme funkci f x interpolujícím polynomem x a tedy d
∫ c
d
f x dx≃∫ x dx
(28)
c
zapišme x= f x 0 0 x f x 1 1 x ... f x n n x . Předpokládejme, že integrály d
∫ i x dx=ci
i=0,1 , ... , n
(29)
c
umíme analyticky spočítat. Koeficienty c i nezávisí na volbě funkce f(x) a můžeme tedy po dosazení do rovnice (29) zapsat rovnost (28) ve tvaru d
∫ f x dx≃c 0 f x 0 c 1 f x 1...c n f x n
.
c
po zavedení ekvidistantní sítě s diferencemi h lze pro aproximaci funkce Lagrangeovým 5
(30)
Martin Isoz, st.sk 204
Matematické metody v inženýrství
8.5.2010
polynomem metodou neurčitých koeficientů odvodit různé kvadraturní formule. Mezi nejběžnější patří lichoběžníkové pravidlo: d
∫ f x dx= h2 f x 0 2 f x 12 f x 2 ...2 f x n−1 f x n
(31)
c
a Simpsonovo pravidlo: d
∫ f x dx= h6 f x 0 4 f x 12 f x 2 4 f x 3...4 f x n−1 f x n
(32)
c
3.3 Symbolický výpočet derivace a primitivní funkce Funkce (3) je jako polynom 2. stupně velmi snadno derivovatelná i integrovatelná. Budou zde tedy pouze vypsány výsledky. f ' x =6 x 3
F x =x 5 x
(33) (34)
x ∈ℝ .
4. Výpočet a prezentace výsledků 4.1 Algoritmus výpočtu Zadaná diferenciální rovnice 2. řádu (1) byla řešena symbolicky s užitím funkce DSOLVE a numericky pomocí modelu sestaveného v prostředí SIMULINK. Pro snazší porovnání dosažených výsledků byl zvolen jednotný grafický výstup obou řešení do podokna typu FIGURE. Společně s řešením byly pro snadné porovnání přesnosti použité metody zobrazeny absolutní a relativní odchylky numericky dopočtených hodnot od přesného řešení. Zároveň byl do COMMAND WINDOW vypsán analytický předpis hledané funkce nalezený programem DSOLVE. Integrace a derivace zadané funkce (3) byla provedena symbolicky pomocí příkazů DIFF a INT a následně numericky. K derivaci byla užita dvoubodová symetrická formule a k integraci funkce QUADGK a Simpsonovo pravidlo. Textový výstup do COMMAND WINDOW obsahuje analytické předpisy derivace a primitivní funkce a symbolicky a numericky zjištěné hodnoty určitého integrálu na zvoleném intervalu. Grafický výstup byl rozdělen do dvou oken typu FIGURE, z nichž v prvním jsou do dvou různých grafů zobrazeny jednotlivé nalezené funkční závislosti (vzhledem k velmi jemnému zvolenému kroku byly pro přehlednost i numericky zjištěné hodnoty zobrazeny hladkou křivkou) a ve druhém jsou zobrazeny relativní a absolutní chyby vypočtených hodnot. K práci jsou přiloženy dva komentované skripty – simulink_M3_skript.m a int_der.m a model simulink_M3.mdl . Schéma modelu je položka [1] v příloze.
4.2 Prezentace výsledků Textové výstupy programů a grafické zobrazení chyb numerické integrace a derivace jsou položky [2] – [4] v příloze. Zde budou prezentovány hlavní grafické výstupy obou skriptů.
6
4.2.1 Numerické řešení ODR s počátečními podmínkami
4.2.2 Numerická integrace a derivace
Martin Isoz, st.sk 204
Matematické metody v inženýrství
8.5.2010
5. Závěr 5.1 Zadání M3.25 Smyslem práce bylo sestavit schéma v SIMULINKu s jehož pomocí by bylo možno numericky řešit zadanou rovnici (1) s počátečními podmínkami (2). V nastavení simulace byla jako řešitel zvolena funkce ODE45 používající metodu Runge-Kutta 4. - 5. řádu s relativní tolerancí 10-3 a proměnným krokem. Ze studia chyby plyne, že relativní odchylka nepřesáhla ±0.05 %, což tomuto nastevení odpovídá. Problém by mohl nastat v případě, kdybychom se tímto způsobem pokoušeli řešit stiff systém rovnic. Zde by totiž, vzhledem k tomu že metody Runge-Kutta nejsou A-stabilní docházelo k extrémnímu zmenšení kroku a výpočet by s velkou pravděpodobností selhal. Pro řešení stiff systémů MATLAB nabízí funkce ODE12s a ODE23s založené na jiných algoritmech.
5.2 Zadání M3.44 Zde bylo úkolem numericky a symbolicky nalézt hodnotu derivace a integrálu v jistých mezích. Vzhledem k tomu, že zadaná funkce (3) byla polynomem 2. stupně byly symbolické výpočty velmi dobře proveditelné. Jelikož byl interval rozdělen na 1000 stejných dílků, jsou numericky nalezené hodnoty relativně přesné. S rezervou je ovšem nutné brát grafické zobrazení hodnot numerické integrace, jelikož není matematicky ani geometricky podloženo. Chyba integrace Simpsonovým pravidlem je na zvoleném intervalu menší než 1 promile.
9
Martin Isoz, st.sk 204
Matematické metody v inženýrství
8.5.2010
6. Seznam použité literatury [1] Doc. RNDr. Daniel Turzík, CSc. a kolektiv, Matematika II ve strukturovaném studiu, Skripta VŠCHT Praha, Praha 2005 [2] Prof. RNDr Milan Kubíček, CSc., RNDr Miroslava Dubcová, Ph.D., Doc. RNDr Drahoslava Janovská, CSc., Numerické metody a algoritmy, Skripta VŠCHT Praha, Praha 2008 [3] MathWorksTM, MATLAB® Documentation http://www.mathworks.com/access/helpdesk/help/techdoc/ [4] Doc. RNDr. František Bubeník, CSc., RNDr. Milan Pultar, CSc., RNDr. Ivana Pultarová, Matematické vzorce a metody, Vydavatelství ČVUT, Praha 1994
10
Martin Isoz, st.sk 204
Matematické metody v inženýrství
8.5.2010
7. Obsah Projekt M3 Řešení diferenciálních rovnic............................................................................................1 1. Zadání..........................................................................................................................................1 2. Popis řešení..................................................................................................................................2 3. Teoretická část.............................................................................................................................2 3.1 Numerické metody řešení ODR s počátečními podmínkami................................................2 3.2 Symbolické řešení diferenciálních rovnic.............................................................................3 3.2.1 Metoda odhadu..............................................................................................................3 3.2.2 Laplaceova transformace...............................................................................................4 3.3 Numerická integrace a derivace............................................................................................5 3.3.1 Diferenční formule........................................................................................................5 3.3.2 Kvadraturní formule......................................................................................................5 3.3 Symbolický výpočet derivace a primitivní funkce...............................................................6 4. Výpočet a prezentace výsledků....................................................................................................6 4.1 Algoritmus výpočtu...............................................................................................................6 4.2 Prezentace výsledků..............................................................................................................6 4.2.1 Numerické řešení ODR s počátečními podmínkami.....................................................7 ...............................................................................................................................................7 4.2.2 Numerická integrace a derivace....................................................................................8 5. Závěr............................................................................................................................................9 5.1 Zadání M3.25........................................................................................................................9 5.2 Zadání M3.44........................................................................................................................9 6. Seznam použité literatury..........................................................................................................10 A. Přílohy.......................................................................................................................................12 [1] schéma použitého modelu simulink_M3.mdl.................................................................12 [2] textový výstup programu simulink_M3_skript.m..........................................................13 [3] textový výstup programu int_der.m................................................................................13 [4] studuim chyb numerické integrace a derivace................................................................14
11
A. Přílohy [1] schéma použitého modelu simulink_M3.mdl
[2] textový výstup programu simulink_M3_skript.m
[3] textový výstup programu int_der.m
[4] studuim chyb numerické integrace a derivace
