Diferenciální rovnice 1 Úvod Diferenciální rovnice jsou matematické rovnice, ve kterých jako proměnné vystupují funkce a jejich derivace. Každý z vás se už s diferenciální rovnicí setkal, i když ve vyřešeném tvaru. Vzpomínáte si na hodiny fyziky a výpočet rychlosti tělesa za pomocí dráhy a času? Zde jste používali již vyřešenou diferenciální rovnici. Diferenciální rovnice nejsou ale pouze součástí fyziky, setkáme se s nimi téměř ve všech oblastech lidského bádání. Při řešení diferenciálních rovnic narazíme na pojem diferenciál1, se kterým budeme pracovat. 1.1 Základní rozdělení diferenciálních rovnic Základní rozdělení diferenciálních rovnic můžeme popsat následovně: • Obyčejné diferenciální rovnice jsou matematické rovnice, které obsahují neznámou funkci jedné nezávislé proměnné a její derivace. Název „obyčejné" se používá jako protiklad k termínu parciální diferenciální rovnice, ve kterých se vyskytuje více než jedna nezávislá proměnná. Nejjednodušší třídou obyčejných diferenciálních rovnic jsou lineární diferenciální rovnice, které mají řešení, jež lze sčítat a násobit koeficienty. Řešení těchto rovnic lze obvykle zapsat v analytickém tvaru pomocí elementárních funkcí. • Parciální diferenciální rovnice je v matematice rovnice obsahující neznámou funkci několika nezávisle proměnných a její parciální derivace dle těchto proměnných. Nej-vyšší z řádů parciálních derivací vyskytujících se v rovnici se nazývá řád parciální diferenciální rovnice. Parciální diferenciální rovnice jsou zobecněním obyčejných diferenciálních rovnic, které obsahují neznámou funkci jedné proměnné a její derivace. Každá obyčejná diferenciální rovnice je současně i parciální diferenciální rovnicí. Parciální diferenciální rovnice může být používána k popisu jevů, jako jsou zvuk, teplo, elektrostatika, elektrodynamika, proudění tekutin, pružnost, nebo kvantová mechanika. diferenciál v matematice vyjadřuje závislost změny hodnoty funkce na malé změně jejího argumentu. Tuto závislost aproximuje jako přímou úměrnost v okolí zvoleného bodu. Pro funkce více proměnných se používá totální diferenciál. 1 1.1.1 Řád diferenciální rovnice Rád diferenciální rovnice je řád nejvyšší derivace, která je v ní obsažená. Za řád soustavy diferenciálních rovnic považujeme hodnotu nejvyšší derivace, která se v soustavě vyskytuje. Podle řádu bývají diferenciální rovnice děleny na diferenciální rovnice prvního řádu a diferenciální rovnice vyšších řádů. Diferenciální rovnice, v nichž se hledaná funkce vyskytuje pouze lineárně, přičemž se nikde nevyskytují ani součiny hledané funkce s jejími derivacemi, ani součiny derivací této funkce, označujeme jako lineární diferenciální rovnice. Pokud jedna z uvedených podmínek není splněna, hovoříme o nelineárních diferenciálních rovnicích. 1.1.2 Řešení diferenciální rovnice Za řešení (integrál) diferenciální rovnice (v daném oboru) považujeme každou funkci, která má příslušné derivace a vyhovuje dané diferenciální rovnici. Řešením (integrálem) soustavy diferenciálních rovnic je množina funkcí s derivacemi potřebného řádu, které vyhovují všem rovnicím dané soustavy. Řešení diferenciálních rovnic dělíme na • obecné - Jako obecné řešení označujeme takové řešení diferenciální rovnice, které obsahuje libovolnou integrační konstantu. Přiřadíme-li každé konstantě obecného řešení určitou číselnou hodnotu, získáme řešení partikulární. • partikulární (částečné) - Partikulární (částečné) řešení je řešení diferenciální rovnice, které získáme přiřazením určité číselné hodnoty každé integrační konstantě obecného řešení. • singulární (výjimečné) - Některá řešení nelze získat z obecného řešení. Taková řešení, která se vyskytují pouze u některých rovnic, popř. v některých bodech oboru, označujeme jako singulární nebo výjimečná. • homogenní Partikulární řešení můžeme v případě jednoduchých diferenciálních rovnic vypočítat analyticky. Nicméně ve velkém množství případů je analytické řešení příliš obtížné a diferenciální rovnice se řeší numericky. Diferenciální rovnice y' = f(x,y) přiřazuje bodu v rovině [x,y] právě jednu hodnotu y'(x), tedy hodnotu derivace hledané funkce. Z dřívějších lekcí víme, že derivace v daném bodě je směrnicí tečny ke grafu v tomto bodě. Můžeme si ji představit jako krátkou úsečku se středem v bodě [x,y]. Tato úsečka se nazývá lineární element. Množina všech lineárních e lementů diferenciální rovnice se nazývá směrové pole. Směrové pole nám tak pomáhá zobrazit tvar hledaných integrálních křivek. 2 2 Diferenciální rovnice 1. řádu Doposud jsme při řešení rovnic, které byly ve tvaru fix) = 0 (funkce jedné proměnné), hledali reálné číslo (kořen rovnice), které splňovalo podmínku. Nyní budeme řešit zcela jiný typ rovnic a jejich řešením nebude číslo, ale funkce. Příkladem může být zadání: Hledejme takovou funkci, která je rovna své derivaci, tedy y = y'. Rovnici například vyhovuje naše velmi dobře známá funkce ex. Této rovnici ale vyhovují násobky funkce ex. Řešení y = k.ex, kde k je libovolná konstanta, se nazývá obecné řešení. Pokud budeme mít stanovenou podmínku y(0) = 2, dosazením do obecného řešení dostáváme 2 = k.e°, odtud k = 2. Řešení y = 2ex je partikulární řešení, které splňuje počáteční podmínku pro x = 0 je y = 2. 2.1 Diferenciální rovnice se separovanými proměnnými Jsou to rovnice ve tvaru: V = f(x)g(y) (1) kde / a g jsou spojité funkce. Derivace funkce y' je rovna podílu diferenciálů dy &dx, tedy y' = dosazením do rovnice 1 dostáváme: 5r = MM (2) dx Za předpokladu, že funkce g(y) je nenulová, separujeme proměnné, tedy upravíme rovnici tak, že výrazy s proměnnou y převedeme na jednu stranu rovnice a výrazy s proměnnou x převedeme na druhou stranu rovnice. Získáme rovnici ve tvaru: ^ = fix)dx (3) 9(y) Tuto rovnici můžeme zintegrovat: dy = / f(x)dx (4) g{y) J Nesmíme zapomenout, že primitivní funkce se liší o konstantu. Přičtením konstanty k jedné straně rovnice obdržíme obecné řešení diferenciální rovnice. Pokud známe počáteční podmínky pro řešení, můžeme vypočítat hodnotu přičtené konstanty a tím získat partikulární řešení diferenciální rovnice, které splňuje zadané počáteční podmínky. Příklad 1 Nyní si ukázkově vyřešíme dříve zmíněnou diferenciální rovnici y = y'. = dy dx dx = — 3 x + C = lny y = ex+c y = ex.ec protože, ec je konstanta, můžeme ji nahradit konstantou K = ec a získáme hledanou rovnici: y = Kex Příklad 2 Nalezněte obecné řešení diferenciální rovnice y' = 2xy. Řešení: Rovnici si přepíšeme do diferenciálového tvaru. dy 0 — = Zxy odtud za předpokladu, že y ^ 0: dx — = 2x dx y = 2x dx Obě strany rovnice zintegrujeme: dy y Vyřešením integrálů na obou stranách rovnice dostaneme tvar: ln \y\ = x2 + K Aplikujeme pravidla pro počítání s logaritmy: ln \y | = x2 + K = ln ex2 + ln eK = ln ex\K Po odlogaritmování dostáváme řešení: \y\ = ex2eK Konstantu eK můžeme přepsat do tvaru C = eK, kde C G IR a rovnici převedeme na tvar: y = Cex2 V případě, že C = 0 bude zde zahrnuto i řešení y = 0. 4 2.2 Lineární diferenciální rovnice 1. řádu Jedná se o diferenciální rovnice tvaru: y' +p{x)y = f{x) (5) kde pix) a fix) jsou spojité funkce. Rovnici y' + p(x)y = 0 (6) nazýváme homogenní lineární diferenciální rovnicí - HLDR. Rovnici y +p(x)y = f{x) (7) nazýváme nehomogenní lineární diferenciální rovnicí - NLDR. Vyřešení HLDR je snadné - rovnici lze převést na rovnici se separovanými proměnnými. Řešení NLDR získáme tak, že k obecnému řešení příslušné HLDR přičteme jedno partikulární řešení NLDR. Příklad 3 Nalezněte obecné řešení rovnice y' = 2y + x. Řešení: 1. krok: Nejprve vyřešíme homogenní rovnici: y' = 2y. Protože se jedná o rovnici se separovanými proměnnými, její řešení lehce zvládneme. — = 2y dx d-y=Í2dx y J Řešením rovnice ln |y| = 2x + K a po odlogaritmování obdržíme obecné řešení homogenní rovnice Vo = Ce2x kde C G IR. 2. krok: Pro hledání partikulárního řešení nehomogenní rovnice použijeme metodu variace konstanty, kdy konstantu C nahradíme funkcí C(x) proměnné x a vypočteme její první derivaci podle x. yp = C{x)e2x (8) Derivací podle proměnné x dostáváme rovnici (jedná se o derivaci součinu) y' = C'(x)e2x + 2C(x)e2x 5 Dosazením za y a y' do původní rovnice získáme rovnici: C'{x)e2x + 2C{x)e2x = 2C{x)e2x + x úpravou této rovnice (odečtením výrazu 2C(x)e2x a vydělením výrazem e2x dostáváme výraz: C'(x) = xe -2x Tuto rovnici zintegrujeme (použijeme metodu per partes) C(x) = Jxe~2x = -X-e~2x - l-e~2x Dosazením Cix) do partikulárního řešení 8 obdržíme: yv = C(x)e2x = (--e~2x - -e~2x)e2x = -- - -yP v > y 2 4 i 2 A 3. krok: Obecné řešení nehomogenní rovnice je dáno součtem obecného řešení nehomogenní rovnice a partikulárního řešení nehomogenní rovnice. V našem případě: y = Ce2x ----y 2 4 Příklad 4 Nalezněte řešení diferenciální rovnice y' — 2^ = 2x3, které vyhovuje počáteční podmínce y(l) = |. Řešení: 1. krok Nalezneme řešení homogenní rovnice y' — 2^ = 0, tedy šl = 2a. dy ľ 2dx y řešení vede na logaritmy: ln \y | = 2 ln x + ln C => ln \y \ = \n{Cx2 odtud \y\ = Cx2 protože hledáme řešení, vzhledem k počáteční podmínce, v intervalu (0; oo), můžeme danou rovnici přepsat do tvaru: y = Cx2 6 krok Vytvoříme a vyřešíme příslušnou nehomogenní lineární rovnici. yp(x) = C(x)x2. y'p{x) = C'{x)x2 + C{x)2x dosazením do původní rovnice obdržíme rovnici: C'(x)x2 + C(x)2x - 2C(x)x2 = 2x3 x po úpravách dostáváme: C'{x) = 2x^ C{x) = J 2xdx = x2 partikulární řešení má tvar: krok Součet obecného řešení HLDR a partikulárního řešení NLDR je obecným řešením zadané diferenciální rovnice. y{x) = Cx2 + x4 Hledané partikulární řešení obdržíme dosazením počátečních podmínek do obecného řešení a výpočtem konkrétní hodnoty konstanty C. yil) = l + C = \ odtud Hledané partikulární řešení je tedy: y{x) = x + —x kde x E (0; oo). 7 Příklad 5 Kromě metody variace konstanty můžeme použít i metodu integračního faktoru. Máme-li NLDR tvaru: y' +p{x)y = f(x) můžeme ji řešit tak, že obě strany rovnice vynásobíme výrazem I{x) = e^p^dx, který nazýváme integračním faktorem, a poté rovnici z integrujeme. Při řešení NLDR y + 3x2y = 6x2 je integračním faktorem výraz: I{x) = eJ'3x2dx = ex?> po vynásobení obou stran rovnice výrazem ex?' dostáváme rovnici: exV + 3xV3y = 6x2ex3 Při pozorném zkoumání předchozí rovnice zjistíme, že na levé straně máme rozepsanou derivaci součinu. Rovnici proto můžeme upravit do tvaru: (ex3yy = 6xV3 nyní můžeme obě strany rovnice zintegrovat: ex3y = J6x2ex3dx Integrál na pravé straně rovnice vyřešíme substitucí t = x3, odtud dx = ^ dostaneme tak výraz: ex3y = 2ex3 + C úpravou rovnice dostáváme: y = 2 + Ce~x3 Tento výraz je obecným řešením zadané lineární diferenciální rovnice. 3 Aplikace diferenciálních rovnic prvního řádu v chemii 3.1 Radioaktivita a radioaktivní rozpad Radioaktivita je jev, při němž dochází k vnitřní přeměně složení nebo energetického stavu atomových jader, přičemž je zpravidla emitováno vysokoenergetické ionizující záření. Radioaktivitu objevil v roce 1896 Henri Becquerel u solí uranu. K objasnění podstaty radioaktivity zásadním způsobem přispěli francouzští fyzikové Pierre Curie a jeho žena Marie 8 Curie-Sklodowská, která byla polského původu. Další z vědců zabývajících se radioaktivitou byl Ernest Rutherford, který dokázal, že rychlost radioaktivního rozpadu je přímo úměrná počtu atomů příslušného prvku N. Tento jev je závislý na čase N(t). Radioaktivní rozklad můžeme popsat diferenciální rovnicí: N' = -XN kde A > 0 je konstanta přeměny. Jedná se o rovnici se separovanými proměnnými, u které můžeme zavést počáteční podmínku ÍV(0) = No, tedy v čase počátku měření máme No atomů. Změnu počtu atomů můžeme zapsat pomocí rovnice: dN = -XN Po separaci proměnných obdržíme: Tuto rovnici zintegrujeme dN -= -Xdt N Odtud dostáváme řešení dN Iv N(t) = Ke -Xdt -xt Aby byla splněna počáteční podmínka musí platit N0 = Ke0 a řešení počáteční úlohy má tvar: N{t) = N0e -xt Pro úlohy tohoto typu je důležitý pojem poločas rozpadu, což je doba, za kterou se přemění právě polovina atomů. V praxi se tohoto jevu využívá k určování stáří organických látek, kde sledujeme izotop uhlíku 14C, jehož poločas rozpadu je 5568 let. Čím delší je poločas rozpadu, tím stabilnější je daný izotop. V lékařské praxi se využívají naopak izotopy s malým poločasem rozpadu. Radioaktivní prvky s dlouhým poločasem rozpadu mohou být při kontaminaci ve velké koncentraci nebezpečné pro zdraví živočichů. Příklad 6 Pro izotop uhlíku 14C známe poločas rozpadu, který je 5568 let. Za jak dlouho se přemění | jader izotopu uhlíku 14C. Dosadíme-li do vztahu za N(t) = |iVo a t = 5568 dostaneme: l-N0 = iVoe-5568A 9 vydělením rovnice výrazem No dostáváme: 1 _ e-5568A 2 Rovnici zlogaritmujeme a získáme výraz: In 2 -ln 2 = -5568A => A = 5568 Pro přeměnu | jader izotopů platí: 3 „ „ ln 2 Odtud pro t platí: ^iV0 = N0e 5568* 5568 ln | t =----i = 2310 ln2 Čtvrtina všech jader izotopu uhlíku C se přemění za 2310 let. 3.2 Kinetika chemických reakcí Předmětem zájmu chemické kinetiky je průběh chemických reakcí před ustavením chemické rovnováhy. Kinetika studuje zejména rychlost reakcí, katalytické působení některých látek a mechanizmy reakcí, tj. jejich skutečné průběhy na molekulární úrovni. 3.2.1 Základní pojmy Rychlost chemické reakce je definována analogicky jako ve fyzice rychlost tělesa. Rychlost chemické reakce bude bude závislá na čase stejně jako rychlost tělesa ve fyzice. Při chemické reakci ale dochází ke změně koncentrace jak výchozích látek - jejich spotřebě, tak produktů - nárůst jejech koncentrace. Rychlost chemické reakce tedy vyjádříme jako podíl změny koncentrace za časový úsek. d(t) Pro snazší zápis budeme pro zápis koncentrace látky A místo c(A) používat nám již dobře známý zápis [A]. Protože při chemické reakci dochází ke spotřebě výchozích látek budeme jejich koncentraci uvádět se záporným znaménkem, naopak produktů při reakci přibývá a proto budeme změnu jejich koncentrace vyjadřovat kladným znaménkem. v = —^jj pro výchozí látky a v = pro produkty. 10 Jak ale víme z praxe, nemusí vždy rychlost chemické reakce záviset pouze na koncentraci jedné látky a v případě některých reakcí může dokonce záviset na mocnině koncentrace chemické látky, nebo na součinu koncentrací několik látek. Uvažujme reakci: aA + I3B ->■ 7C Při vyjadřování rychlosti chemické reakce tedy využíváme stechiometrické koeficienty. [Ar v = d(t) Dalším důležitým pojmem, který budeme používat je řád chemické reakce. Budeme rozlišovat dílčí řád chemické reakce, který bude vyjádřen exponentem (stechiometrickým koeficientem) u koncentrace v rovnici vyjadřující rychlost chemické reakce. Součet dílčích řádů chemické reakce vyjadřuje celkový řád chemické reakce. Uvažujeme-li kinetickou rovnici: v = k \A]a.[BY hovoříme o dílčím řádu reakce a vůči koncentraci látky A a dílčím řádu reakce /3 vůči koncentraci látky B. Celkový řád reakce je roven součtu a + /3. K určení rychlostního zákona a tím i řádu dané reakce nestačí znát jen chemickou rovnici, musí se v každém jednotlivém případě určit speciálními metodami. Dalším důležitým pojmem v kinetice chemických reakcí je rychlostní konstanta - k. Rychlostní konstanta je konstanta vyskytující se ve vztahu mezi rychlostí reakce a okamžitými koncentracemi reagujících látek (v rychlostním zákonu). Položíme-li všechny koncentrace v kinetické rovnici rovny jedné, pak tato rovnice přejde na tvar: v = k. Z toho vidíme, že rychlostní konstanta z kinetické rovnice nám číselně udává, jak rychle by reakce probíhala při jednotkových koncentracích všech výchozích látek v reakci. Hodnota rychlostní konstanty je závislá na tom, jaké látky reagují a na podmínkách pokusu. Největší vliv na změnu rychlostní konstanty má teplota - se vzrůstající teplotou se hodnota k zvětšuje. Dalšími faktory, které mají vliv na rychlost reakce, jsou: tlak, rozpouštědlo a přítomnost některých látek, jež se reakcí nezmění. Takové látky nazýváme katalyzátory a o procesu urychlení reakce mluvíme jako o katalýze. V případě chemické reakce probíhající v plynné fázi budeme místo koncentrací látek používat jejich parciální tlaky. 3.2.2 Reakce nultého řádu Pro tento druh chemických reakcí je charakteristické, že reakce probíhá konstantní rychlostí a není závislá na koncentraci reagujících látek. Matematicky tento vztah vyjádříme: 11 v = k [A}°[B}° tedy v = k .1.1 = k Rovnici rychlosti chemické reakce můžeme zapsat následovně: dc —r = k dt Tuto diferenciální rovnici můžeme řešit metodou separace proměnných. —dc = k .dt dc = —k.dt Před integrací musíme ještě stanovit meze, protože řešíme reálný případ, který probíhá v konečném časovém úseku a za určité počáteční koncentrace. Příklad řešíme pro koncentraci reaktantu, který se v průběhu chemické reakce spotřebovává, proto je před změnou koncentrace záporné znaménko. Integraci budeme provádět pro časový úsek v čase od 0 do času t. Na počátku bude koncentrace reaktantu cq a po čase t bude koncentrace reaktantu rovna hodnotě c. Matematicky zapsáno: ŕ dc = —k / dt c0 Jo [c]cC0 = -k[tf0 c — cq = —kit — 0) co — c = k t Jak zjistíme dále, bude mít rychlostní konstanta různé jednotky pro různé řády chemických reakcí. V našem případě určíme jednotku rychlostní konstanty následovně. k =^ t koncentrace je udávána v jednotkách mol.dm-3, čas je v sekundách s. Rychlostní konstanta má tedy jednotku: r r -3 -1 moí mol.dm .s = dm3.s Dalším důležitým ukazatelem kinetiky chemické reakce je poločas chemické rekace. Podobně jako v případě radioktivního rozpadu hovoříme o poločasu rozpadu látky, je 12 poločas chemické reakce definován jako čas, kdy sledovaná koncentrace reaktantu poklesne na polovinu. V čase t je tedy koncentrace rovna ^cq. Kinetická rovnice má tedy tvar: 1 co - t:co k .t 1 k .t £o 2k Poslední otázkou je, jak vypadá graf rychlosti chemické reakce. Na osu x budeme vynášet průběh sledovaného času a na osu y koncentraci sledované látky. Grafem bude přímka, jak ukazuje následující obrázek: t [min] Obrázek 1: Graf rychlosti chemické reakce pro reakci nultého řádu Směrnice přímky je rovna rychlostní konstantě k. 3.2.3 Reakce prvního řádu V tomto případě je rychlost chemické reakce závislá na první mocnině koncentrace jedné sledované látky v = k [A]. Rychlost rovnice reakce prvního řádu můžeme zapsat následovně: I v tomto případě řešíme diferenciální rovnici metodou separace proměnných a opět integrujeme pro časový úsek od 0 do ŕ a koncentrace od cq do c. c [mol/dm3] dc ďt k .c dc = k .c.dt dc kdt c 13 [\ncfc=-k.t ln c — ln co = —kt Protože rozdíl logaritmů hodnot je roven logaritmu podílu daných hodnot, můžeme rovnici přepsat do tvaru: ln — = -kt co ln^ = kt c Již z tohoto výsledku je zřejmé, že rychlostní konstanta bude mít pro reakci 1. řádu jinou jednotku, než pro reakci nultého řádu. Proto si ji odvodíme. ln^ k = —c-t V čitateli zlomku máme přirozený logaritmus podílu koncentrací, tedy bezrozměrnou veličinu. Ve jmenovateli je jednotka času, tedy sekunda s. Rychlostní konstanta má tedy jednotku s_1. Podobně jako v předchozím případě si odvodíme poločas chemické reakce 1. řádu. 1 , , co x 1 , , 1 x 1 , ři = í-ta = *ln2 Výsledná rovnice má tvar: ln2 2 k Je vidět, že poločas chemické reakce 1. řádu nezávisí na počáteční koncentraci výchozí látky. 3.2.4 Reakce druhého řádu pro případ, kdy je rychlost chemické rakce závislá na koncentraci pouze jedné látky V tomto případě lze rychlost chemické reakce vyjádřit rovnicí: v = k .[A]2 nebo také: 4 = k -c2 dt 14 Tuto diferenciální rovnici řešíme stejně jako předchozí pomocí metody separace proměnných. Postup řešení vypadá následovně: dc = -k .dt c -k.[t][ = -k.t c c0 -- — = k.t c c0 Stejně jako v předchozích případech, i zde se podíváme na jednotky rychlostní konstanty chemické reakce. Opět vidíme, že i v tomto případě bude mít rychlostní konstanta opět jiný rozměr. Osamostatněním rychlostní konstanty získáme: 1,1 1 , k = v t c Co' jednotka rychlostní konstanty tedy bude: i „-i dm .mol .s Pro vytvoření grafu bude výhodnější na osu y vynést hodnoty ~c než hodnoty koncentrace c. 1/cO t [min] Obrázek 2: Graf rychlosti chemické reakce pro reakci 2. řádu Případ, kdy rychlost chemické reakce závisí na koncentraci dvou látek je nad rámec tohoto kurzu. Pro řešení uvedeného případu je třeba umět řešit integrály racionálních lomených funkcí. Pro zájemce o tuto oblast doporučuji učebnice matematické analýzy. 15 3.2.5 Řešené příklady Esterifikace Zadání: Esterifikace CH3COOH + C2H5OH ->■ CH3COOC2H5 + H20 byla sledována měřením koncentrace kyseliny octové po dobu 31 minut a výsledky byly zapisovány do tabulky. Objem systému a teplota reakční směsi se během reakce nemění. Určete rychlost chemické reakce v čase t = SOmín, dále určete rychlost úbytku koncentrace ethanolu (koncentrace na začátku reakce byla 1 mol.dm-3 a také rychlost přírůstku koncentrace esteru pro případ, kdy jeho počáteční koncentrace byla nulová. čas reakce [min] 0 . 29 30 31 koncentrace kys. octové [mol.dm~3] 0.5 . . 0.101 0.100 0.099 Řešení: Pro případ, kdy jsou rozdíly koncentrací malé lze rychlost vypočítat jako poměr diferencí koncentrace a času. Rozdíl koncentrací mezi 29. a 30. minutou je 0.001 mol.dm~3 pro časovou diferenci 1 min. rychlost chemické reakce je tedy 2^21 _ z tabulky plyne, že počáteční koncentrace kyseliny octové je 0.5 mol.dm~3'. Z látkové bilance je zřejmé, že úbytek koncentrace reaktantů a přírůstek koncentrace produktů mají stejné hodnoty. Úbytek koncentrace kyseliny octové je po 30 minutách roven 0.5 - 0.1 = 0.4 mol.dm~3, úbytek koncentrace ethanolu je 1 - 0.4 = 0.6 mol.dm~3'. Přírůstek koncentrace esteru je roven 0.4 mol.dm~3'. Reakce 2NO + Cl2 ->■ 2NOCI Zadání: Uvažujme chemickou reakci 2NO + Cl2 —> 2NOCI, pro jejíž rychlost platí vztah v = k [NO]2[Cl] a rychlostní konstanta k = 3.84 • \Q~Adm&mol~2s_1. Vypočítejte rychlost reakce pro reakční nádobu o objemu 2.50 dm3, která obsahuje 4.00 mol NO a 3.50 mol Cl2. Jak se změní rychlost reakce, použijeme-li poloviční množství NO. Řešení: Nejdříve musíme vypočítat koncentrace reaktantů. [NO] = |§ = ímmol.dm-3 [Cl2] = |f = lÁOmol.dm-3 Dosazením do rychlostní rovnice získáme hodnotu: v = 3.84.10"4 • (1.60)2.(1.40) = l.?>%.lQ-3mol.dm-3.s'1 Protože rychlost reakce závisí na druhé mocnině koncentrace NO, při použití polovičního množství bude rychlost reakce čtvrtinová. 16 4 Numerické řešení diferenciálních rovnic V numerické matematice je numerické řešení obyčejných diferenciálních rovnic postup, kterým můžeme získat přibližné řešení obyčejných diferenciálních rovnic (ne parciálních). Používá se v případech, kdy by bylo nalezení přesného (analytického) řešení náročné nebo v případech, kdy analytické řešení nelze najít. 4.1 Eulerova metoda Základní myšlenkou této metody je aproximace řešení lomenou čarou. Začneme v bodě zadaném počáteční podmínkou a v okolí tohoto bodu nahradíme integrální křivku její tečnou. Tím se dostaneme do dalšího bodu, odkud opět integrální křivku aproximujeme tečnou. Směrnici tečny zjistíme z diferenciální rovnice přímo z derivace. Uvažujme počáteční úlohu: V = f(x,y), y(x0) = y0 Potřebujeme nalézt přibližné řešení y(x) pro x E [xo; xq + a]. Postupujeme tak, že interval rozdělíme na n podintervalů délky ht, takto dostaneme dělící body: xi = x0 + hi, x2 = xi + h2, ■ ■ ■, xn = + hn = x0 + a kde h\ + hi + . .. + hn = a. Za pomoci funkčních hodnot vypočteme y'0 = f(xo,yo) a položíme yi = vo + hy'o = yo + hif{x0, y0) Podobně určíme y2 = y± + h2f(xi, yi) a stejně budeme postupovat dále. Dostaneme tak přibližné řešení: y{x) =y% + f{xu y%){x - Xi) pro x E [xí, xl+i], (i = 0,1,... , n — 1). Nejjednodušším způsobem dělení intervalu je použití stejně vzdálených dělících bodů. V tomto případě bude Eulerova metoda popsána následovně: 3^1+1 = xi + h yl+1 =y% + hf{x%, yi), i = (0,1, 2,.. ., n) Modelový příklad Pomocí Eulerovy metody určete přibližné řešení počáteční úlohy y' = x + y pro y(0) = 1 s krokem h = 0,1. Získaný výsledek porovnejte s analytickým řešením. 17 Řešení: Máme dáno: h = 0,1, xq = 0, yo = 1 a f(x, y) = x + y. Podle výše popsaného postupu budeme postupně počítat: Ví = yo + hf{x0,y0) = 1 + 0,1(0 + 1) = 1,1 V2 = yi + hf(Xl, yi) = 1,1 + 0,1(0,1 + 1,1) = 1, 22 V3 = V2 + hf{x2, y2) = 1, 22 + 0,1(0, 2 + 1, 22) = 1, 362 Podobně pak postupujeme dále. Výsledky si zapíšeme do tabulky: i Ví i Ví 1 0,1 1,100000 6 0,6 1,943122 2 0,2 1,220000 7 0,7 2,197434 3 0,3 1,362000 8 0,8 2,487178 4 0,4 1,528200 9 0,9 2,815895 5 0,5 1,721020 10 1,0 3,187485 Jedná se o lineární diferenciální rovnici, jejíž řešením je y(x) = 2ex — x — 1 a pro y(l) = 2e — 2 ~ 3, 436564. Chyba našeho numerického řešení od analytického řešení je 0, 249079. 4.2 Runge - Kuttovy metody Zatímco Eulerova metoda používá jako směrnici k = f(xl,yl) hodnotu směrového pole v bodě, ze kterého vychází, nejjednodušší Runge - Kuttova metoda sleduje hodnoty směrového pole v bodě, kam bychom došli v případě polovičního kroku od výchozího bodu - Metoda RK druhého řádu. Ještě sofistikovanější metoda je RK4 - metoda čtvrtého řádu. Zde podobně jako u metody druhého řádu uděláme fiktivně půl kroku směrem k\ podle Eulerovy metody a ze směrového pole v bodě do kterého se dostaneme získáme směr k2. Poté podobně provedeme opět fiktivně půl kroku směrem k2 a ze směrového pole v bodě, do kterého se dostaneme, získáme směr k^. Konečně, směrem k% provedeme fiktivně celý krok a získáme směr k±. Ze všech těchto směrů vypočítáme vážený průměr ve kterém jsou k2 a &3 zastoupeny dvojnásobnou vahou oproti k\ a k± a získáme směr pro provedení dalšího kroku metody. Matematicky vyjádříme Eulerovu metodu takto: xl+1 = xl + h; yl+1 =yl + h.k; k = kx = f(xt, y{) Runge - Kuttovu metodu druhého řádu vyjádříme takto: h h xl+1 = xl + h; yl+1 = yl + h.k] k = k2 = fix, + -,y, + kx-) kde k\ = f(xi, yij jako v Eulerově metodě. 18 Runge - Kuttovu metodu čtvrtého řádu vyjádříme takto: xl+1 = xl + h; yi+1 = yl + h.k; k = \{k1 + 2k2 + 2k3 + k4) o kde kx = f(xu y%), k2 = /(xi + f, yl + k1^), k3 = /(xj + §, ylJrk2\) ak4 = f(xi + h, yl + h). Popsaná metoda RK4 dosahuje v jednom kroku chyby odpovídající h5, kumulativní chyba je pak v řádu h4. Neuvažujeme-li vliv zaokrouhlovacích chyb, tak menší krok obvykle vede k přesnějšímu odhadu, avšak za cenu více počítání. 19