7 Výpočet určitého integrálu a numerické řešení diferenciálních

7 Výpočet určitého integrálu a numerické řešení di‐ ferenciálních rovnic Studijní cíl V sedmém bloku se seznámíme s možnostmi MATLABu pro řešení dalších dvou důležitých matematických úloh ‐ numerický výpočtem určitých integrálů a oby‐ čejných diferenciálních rovnic. Výpočet určitých integrálů (jedno‐ i vícenásob‐ ných) je důležitý při řešení mnoha úloh – např. při určování délky křivky, obsa‐ hu rovinných útvarů, při určování povrchu, objemu a hmotnosti těles, při určo‐ vání těžiště či momentu setrvačnosti atd. Na řešení diferenciálních rovnic vede mnoho úloh z elektrotechniky, mechaniky, termodynamiky a zejména jsou využívány při popisu chování dynamických systémů. O významu soustav oby‐ čejných nelineárních diferenciálních rovnic svědčí i to že řešením pouze této úlohy se zabývá specializovaná nadstavba MATLABu zvaná SIMULINK, které se budeme později věnovat v blocích 9‐12.

Doba nutná k nastudování

1 ‐ 2 hodiny

Průvodce studiem Předpokládá se, že čtenář je seznámen s prostředím MATLABu, základními příkazy, tvorbou grafů a tvorbou skriptů a funkcí v rozsahu prvních čtyř bloků. Ze šestého bloku jsou potřeba znalosti o předávání funkce jako parametru jiné funkce. Dále se předpokládají základní znalosti integrálního a diferenciálního počtu. Při studiu je vhodné mít spuštěný MATLAB a jednotlivé příkazy či funkce hned vyzkoušet a případně získat podrobnější informace využíváním nápovědy k programu. Na závěr jsou uvedeny řešené příklady na procvičení. V textu jsou použity následující typografické konvence: Calibri 11, Bold, Italic nové pojmy či informace k zapamatování Calibri 11 označení klávesy Courier New 10, Bold názvy nástrojů MATLABu Courier New 10, Bold upřesnění nápovědy (help téma) Courier New 10 názvy příkazů, funkcí a objektů Courier New 10 označení části příkazu, která se může vynechat Courier New 10, Italic názvy proměnných použitých v programu Courier New 9 příkaz příkazové řádky / výpis programu Courier New 9 výpis programu – klíčová slova Courier New 9 výpis programu – řetězec Courier New 9 výpis programu – komentář

KŘP/IMSW Modelování ve výpočtových software

7‐1 (20) 4.7.12

František Dušek KŘP FEI Univerzita Pardubice

7.1 Numerický výpočet určitého integrálu Další matematickou úlohou, pro kterou lze stanovit obecný postup řešení bez znalosti konkrétní funkce, je výpočet určitého integrálu. Problém b I1 = f ( x )dx je principiálně jednoduchý v případě jednorozměrného inte‐ grálu tj. pro jednorozměrnou funkci. V případě vícerozměrných a I 2 = f ( x, y )dxdy integrálů je problém komplikovanější v případě, že inte‐ grujeme přes jinou než obdélníkovou oblast (v případě S dvourozměrného) nebo oblast tvaru kvádru (v I3 = f ( x, y , z )dxdydz případě třírozměrného integrálu).

∫

∫∫

∫∫∫ V

Z pohledu řešení v MATLABu jde opět o funkci funkcí, kdy je řešící funkci dodán handle konkrétní funkce, která se má na zadaném intervalu či oblasti integro‐ vat. Pro výpočet jednotlivých integrálů jsou určeny následující základní řešící funkce (blíže viz help quadl, help dblquadl, help triplequad) I=quadl(f,a,b) výpočet určitého jednorozměrného integrálu z funkce f na intervalu od a do b I=dblquad(f,xmin,xmax,ymin,ymax) výpočet určitého dvou‐ rozměrného integrálu z funkce f na obdélníkové oblasti P vymezené rozsa‐ hem od xmin do xmax a od ymin do ymax I=triplequad(f,xmin,xmax,ymin,max,zmin,zmax) výpočet určitého třírozměrného integrálu z funkce f na kvadratické oblasti V vyme‐ zené rozsahem od xmin do xmax, od ymin do ymax a od zmin do zmax Defaultní přesnost určení určitého integrálu je ve všech případech 1e‐6. Uživa‐ telskou funkci f je potřeba vytvořit tak, aby akceptovala vektor pro proměnnou x (a skaláry pro případné další nezávisle proměnné y či z) a vracela odpovídající vektor funkčních hodnot.

7.1.1 Výpočet jednoduchého určitého integrálu ‐ příklady Při řešení některých úloh může být užitečná znalost geometrické interpretace určitého integrálu. Hodnota určitého integrálu nějaké funkce od bodu a do bodu b je plocha mezi osou x a průběhem funkce ohraničená kolmicemi v bodec a a b. Pozor ale, že plocha pod osou x se uvažuje jako záporná. Uveďme jako příklad výpočet obsahu kružnice o poloměru r (který dokážeme spočítat i bez určitého integrálu). Vyjdeme‐ li z Obrázku 7‐1 pak obsah kružnice je určen dvoj‐ násobkem plochy pod křivkou na obrázku tj. vztahem

y

x

+r

S = 2 ∫ r − x dx 2

2

‐r

−r

r f ( x) = r 2 − x 2

Řešení v MATLABu s ověřením pomocí vzorce pro výpočet obsahu kružnice pak může vypadat např. (pro případ r=10) jako

Obrázek 7‐1 Výpočet obsahu kružnice


7‐2 (20) 4.7.12


>> r=10; >> P1=2*quadl(@(x) sqrt(r^2-x.^2),-r,r); >> P2=pi*r^2; >> e=P1-P2 e = -1.6325e-006

Jako další příklad uveďme výpočet hodnoty určitého integrálu z funkce sinus od 0 do 2π, kde umíme určit i bez výpočtu integrálu, že má vyjít 0. Víte proč? >> P=quadl(@sin,0,2*pi) P = 1.7588e-016

Jednoduchý určitý integrál je možné použít i pro výpočet objemu rotačních těles. Vznikne‐li nějaké těleso rotací křivky f(x) okolo osy x, pak je jeho objem určen vztahem B2=[0,0.5] B1=[‐1,0.5]

b

V = π ∫ f ( x)dx

B3=[1,0.5]

y

2

a

0.5 0.6 Chceme‐li tedy určit např. objem sudu s geometrickými rozměry podle x Obrázku 7‐2, pak je řešení velmi jedno‐ duché za předpokladu znalosti funk‐ ce popisující zaoblení sudu a vhodné 2 volby souřadného systému. Pokud budeme funkci zakřivení aproximovat Obrázek7‐2 Objem sudu polynomem 2. stupně (parabolou) procházející body B1 až B3 pak můžeme řešení v MATLABu napsat např. jako >> x=[-1,0,1]; y=[0.5,0.6,0.5]; >> p=polyfit(x,y,2); >> V=quadl(@(x) pi*polyval(p,x).^2,-1,1) V = 2.0232

Poslední příklad ukazuje použití určitého integrálu pro určení délky křivky v prostoru. Máme‐li křivku v prostoru popsanou parametrickými rovnicemi ve tvaru x = f x (t ) y = f y (t ) z = f z (t ) pak délka křivky je dána vztahem t =b

L=

∫

t =a

2

2

2

⎛ dx ⎞ ⎛ dy ⎞ ⎛ dz ⎞ ⎜ ⎟ + ⎜ ⎟ + ⎜ ⎟ dt ⎝ dt ⎠ ⎝ dt ⎠ ⎝ dt ⎠

1

Mějme křivku na Obrázku 7‐3 popsanou rovnicemi

( )

x = sin t 2 y = sin(t ) cos(t ) z = cos t 2

( )

( )

x′ = 2t cos t 2 y′ = cos 2 (t ) − sin 2 (t ) z′ = −2t sin t 2

0

-0.5

( )

-1 0.5

pak vykreslení průběhu křivky pro 0 ≤ t ≤ π a výpočet délky křivky je možné vypočítat např. jako >> >> >> >>

0.5

1 0.5

0 0 -0.5 -0.5

-1

t=linspace(0,pi,1000); Obrázek 7‐3 Křivka v prostoru x=sin(t.^2); y=sin(t).*cos(t); z=cos(t.^2); plot3(x,y,z), grid f=@(t) (2*t.*cos(t.^2)).^2 + (cos(t).^2-sin(t).^2).^2 + + (2*t.*sin(t.^2)).^2;


7‐3 (20) 4.7.12


>> L=quadl(@(x) sqrt(f(x)),0,pi) L = 10.3071

7.1.2 Výpočet dvojného určitého integrálu – příklady Geometrická interpretace dvojného integrálu z funkce f(x,y) přes integrační oblast S(x,y) tj.

∫∫ f ( x, y)dxdy je objem omezený shora plochou danou funkcí S

f(x,y), zdola uzavřenou oblastí v rovině x‐y a z boků pláštěm kolmým na rovinu x‐y procházejícím hranicí oblasti S. Opět s konvencí, že objem pod rovinou x‐y je záporný. V případě pravoúhlé integrační oblasti a ≤ x ≤ b, c ≤ y ≤ d (úloha bd

∫∫ f (x, y)dydx ) je přímo použitelná funkce dblquad(f,a,b,c,d). a c

Jako příklad uveďme výpočet objemu části polokoule o poloměru r popsané funkcí f ( x, y ) = r 2 − x 2 − y 2 nad oblastí ‐ r /2 ≤ x ≤ + r /2, ‐ r /2 ≤ y ≤ +r/2 (viz Obrázek 7‐4). Matematická formulace pro uvedený příklad je r

V=

2

r

2

∫∫

r 2 − x 2 − y 2 dydx

−r 2 −r 2

Příkazy pro vykreslení vzniklého útvaru pro r=10 a výpočet jeho objemu jsou

10 8 6

>> r=10; x=-r/2:r/19:r/2; y=x; 4 >> [X,Y]=meshgrid(x,y); 2 >> Z=sqrt(r^2-X.^2-Y.^2); 0 10 >> Z1=Z*0; >> h=meshz(x,y,Z); 5 >> set(h,'LineWidth',1) 0 >> axis equal -5 -5 >> axis([-r,r,-r,r,0,r]) -10 -10 >> hold on >> surf(x,y,Z1) Obrázek 7‐4 Část polokoule >> f=@(x,y,p) sqrt(p^2-x.^2-y^2) >> V=dblquad(@(x,y) f(x,y,r),-r/2,r/2,-r/2,r/2) V = 910.9658

10 5 0

Pokud chceme spočítat objem celé polokoule o poloměru r (integrační oblastí je kruh o poloměru r) pak je matematická formulace problému

⎡ r 2 − x2 ⎤ V = ∫⎢ ∫ r 2 − x 2 − y 2 dy ⎥dx . Protože funkce dblquad má integrační ⎥⎦ −r ⎢ ⎣− r 2 − x2 r

meze konstantní a tedy volá vnitřní funkci pro body z obdélníkové1 integrační oblasti, není možné vnitřní funkci vytvořit pouhým přepisem matematického zadání jako v předchozím případě. Možné řešení je vnitřní funkci realizující matematickou funkci ze zadání rozšířit tak, aby pro hodnoty x, y mimo požado‐ 1

Tj. i pro body mimo požadovanou kruhovou integrační oblast, kde funkce není definována


7‐4 (20) 4.7.12


vanou kruhovou integrační oblast vracela nulové hodnoty, tj. aby funkční hod‐ noty pro tyto body neovlivnily výslednou hodnotu integrálu. Příkladem takové funkce může být např. PoloKoule.m ==================================================== function z=PoloKoule(x,y,r) % výpočet výšky bodu polokoule o poloměru r nad rovinou x-y % pro zadaný bod [x,y] v=r^2-x.^2-y^2; z=zeros(size(v)); for k=1:length(v), % !!! x (tj. i v) může být vektor if v(k)>0, z(k)=sqrt(v(k)); end end ==================================================== PoloKoule.m

Využijeme‐li tuto funkci i pro vykreslení útvaru (Obrázek 7‐5), jehož objem zjišťujeme, pak můžeme zapsat např. skript r=10; x=-r:r/20:r; 10 y=x; 8 nx=length(x); 6 ny=length(y); 4 Z=zeros(nx,ny); 2 for k=1:nx, 0 10 for l=1:ny, Z(k,l)=PoloKoule(x(k),y(l),r); end end h=mesh(x,y,Z); set(h,'LineWidth',2) axis equal axis([-r,r,-r,r,0,r]) V=dblquad(@(x,y) PoloKoule(x,y,r),-r,r,-r,r);

10

5 5

0 0 -5

-5 -10

-10

Obrázek 7‐5 Polokoule

Jiná varianta vnitřní funkce s „chytřejším“ způsobem řešení rozhodování zda hodnoty leží v integrační oblasti, je ukázána dále PoloKoule.m ==================================================== function z=PoloKoule(x,y,r) % výpočet výšky bodu polokoule o poloměru r nad rovinou x-y % pro zadaný bod [x,y] v=max(r^2-x.^2-y^2,0); z=sqrt(v); ==================================================== PoloKoule.m

nebo přímo v definici anonymní funkce V=dblquad(@(x,y) sqrt(max(r^2-x.^2-y^2,0)),-r,r,-r,r); V = 2.0944e+003

7.1.3 Výpočet trojného určitého integrálu – příklady Geometrická interpretace trojného integrálu z funkce f(x,y,z) přes integrační oblast O(x,y,z) tj.

∫∫∫ f ( x, y, z)dxdydz

je hmota oblasti (objemu) O

O


7‐5 (20) 4.7.12


s proměnnou hustotou popsanou funkcí f(x,y,z). Pokud je funkce f(x,y,z) v celé integrační oblasti jednotková, pak hodnota integrálu odpovídá objemu oblasti.

V případě pravoúhlé integrační oblasti a ≤ x ≤ b, c ≤ y ≤ d, e ≤ z ≤ f (úloha b d f

∫∫∫ f (x, y, z )dzdydx ) lze použít funkci triplequad(f,a,b,c,d,e,f). a c e

Použití funkce triplequad si ukážeme opět na příkladu – výpočtu hmoty tělesa s proměnnou hustotou pomocí trojného integrálu. Hmota tělesa je obecně dána vztahem

∫∫∫ ρ (x, y, z )dxdydz kde V je oblast v prostoru, přes V

kterou integrujeme, a ρ(x,y,z) je rozložení hustoty v prostoru. Integrační oblast V bude kvadratická 0 ≤ x ≤ a, 0 ≤ y ≤ b, 0 ≤ z ≤ c a rozložení hustoty v prostoru je popsáno vztahem ρ ( x, y , z ) =

1 1 + x2 + y2 + z 2

. Příkazy pro řešení pak

mohou být např. >> f=@(x,y,z) 1./sqrt(1+x.^2+y^2+z^2); >> M=triplequad(f,0,1,0,2,0,3) M = 2.8423

Využití určitých integrálů není omezeno jen pro výpočty objemu. Používají se pro výpočet plochy, povrchu, hmoty, statických momentů, polohy těžiště, mo‐ mentu setrvačnosti tělesa (plochy, křivky), pro výpočet např. práce po křivce v silovém poli atd.

7.2 Numerické řešení obyčejných diferenciálních rovnic V dalším textu se budeme zabývat pouze obyčejnými diferenciálními rovnicemi (Ordinary Differential Equation ODE) jejichž řešením je funkce jedné (nezávisle) proměnné. Parciálními diferenciálními rovnicemi, jejichž řešením je funkce více (nezávisle) proměnných se zabývat nebudeme. Diferenciální rovnice popisuje závislost mezi neznámou funkcí, jejími derivacemi a případně funkcemi nezá‐ visle proměnné a neznámé funkce. Obecným analytickým řešením je určení množiny (prostoru) všech funkcí2, které po derivování a dosazení splňují dife‐ renciální rovnici. Pro výběr jedné konkrétní funkce z prostoru řešení je potřeba zvolit bod v prostoru řešení, kterým má funkce (řešení) procházet. Protože nás obvykle zajímá řešení na určitém intervalu hodnot nezávisle proměnné tak se bod, kterým má procházet řešení nejčastěji volí buď na začátku či konci tohoto intervalu. Podle umístění vybraného bodu v prostoru řešení se pak diferenciál‐ ní rovnice označuje jako diferenciální rovnice s počátečními (začátku intervalu) či okrajovými (na konci intervalu) podmínkami. Řešení může být buď v tvaru analytickém (vyjádřeno funkcí, předpisem pro výpočet souřadnice bodu řešení 2

Prostor všech funkcí představujících řešení diferenciální rovnice je popsán (pokud v analytické podobě vůbec existuje) jednou funkcí, ve které se vyskytuje tolik blíže neurčených parametrů, kolikátého řádu daná diferenciální rovnice je. Počáteční či okrajové podmínky slouží k určení konkrétních hodnot těchto parametrů.


7‐6 (20) 4.7.12


na základě hodnoty nezávisle proměnné) nebo numerickém (vyjádřené tabul‐ kou souřadnic bodů představujících řešení).

Analytické řešení diferenciální rovnice s počátečními či okrajovými podmínka‐ mi je tedy představováno funkcí, která pro hodnoty nezávisle proměnné umožní spočítat odpovídající funkční hodnotu. Analytické řešení obyčejných nelineárních diferenciálních rovnic existuje (v uzavřeném tvaru) pouze ve spe‐ ciálních případech. Pro obyčejné lineární3 diferenciální rovnice analytické řeše‐ ní existuje ale jeho získání je pracné a výsledná funkce může být značně složitá. Jako příklad analytického řešení ukažme řešení diferenciální rovnice

⎛ x2 ⎞ dy + y = xe − x s obecným řešením y ( x) = ⎜⎜ C + ⎟⎟e − x a konkrétním řešením 2⎠ dx ⎝ ⎛ x2 ⎞ y ( x) = ⎜⎜ 2 + ⎟⎟e − x pro počáteční podmínku y(x=0)=2. 2⎠ ⎝ Pokud nám však stačí znát pro zadané hodnoty nezávislé proměnné hodnoty, které jsou řešením dané diferenciální rovnice s počátečními (či okrajovými) podmínkami, pak můžeme použít numerické řešení diferenciální rovnice, které lze nalézt vždy (pokud vůbec řešení dané diferenciální rovnice existuje). Řeše‐ ním v tomto případě je tedy tabulka, kde pro každou hodnotu nezávisle pro‐ měnné je spočítána odpovídající funkční hodnota. MATLAB obsahuje několik různých metod (algoritmů) numerického výpočtu obyčejných diferenciálních rovnic opět ve formě funkce funkcí. Aby bylo možné použít tyto funkce pro obyčejnou diferenciální rovnici libovolného stupně4, vyžadují zápis diferenciální rovnice ve tvaru soustavy diferenciálních rovnic prvního stupně. Soustavu N obyčejných diferenciálních rovnic prvního stupně včetně počátečních či okrajových podmínek

dy1 dy2 dyN

dx

= f1 ( x, y1 , y2 ,K, yN )

y1 ( x = x0 ) = y10 y2 ( x = x0 ) = y20

dx

= f 2 ( x, y1 , y2 ,K, yN )

dx

= f N ( x, y1 , y2 ,K, yN ) yN ( x = x0 ) = yN 0

M

M

lze zapsat ve formálně stejném maticovém zápisu

dY = f (x, Y) Y0 = Y(x = x0 ) dx kde

x je označení nezávisle proměnné, Y = [y1(x), y2(x), … , yN(x)]T je vektor (neznámých) funkcí řešících diferen‐ ciální rovnici,

3

Obyčejná lineární diferenciální rovnice obsahuje pouze lineární kombinace (sčítání/odečítání, násobení konstantou) neznámé funkce a jejich derivací a případně funkce nezávisle proměnné.

4

Stupeň či řád diferenciální rovnice je určen řádem nejvyšší derivace, která se v diferenciální rovnici vyskytuje.


7‐7 (20) 4.7.12


f = [f1(x, y1(x), y2(x), … , yN(x), … , fN(x, y1(x), y2(x), … , yN(x)] T je vektor funkcí popisujících výpočet hodnoty derivace příslušné funkce v závislosti na hodnotě nezávisle proměnné x a Y0 = [y1(x0), y2(x0), … , yN(x0)]T je vektor počátečních podmínek tj. hodnot neznámých funkcí pro hodnotu nezávisle proměnné x=x0.

Jestliže je počet rovnic a počet hledaných funkcí stejný a jsou známé počáteční podmínky pak lze pro takovou soustavu rovnic v MATLABu získat (poud řešení existuje) numerické řešení pro libovolné hodnoty x na intervalu x0 ≤ x ≤ xmax.

7.2.1 Převod diferenciální rovnice vyššího řádu na soustavu diferen‐ ciálních rovnic řádu prvního Velmi často máme nalézt řešení obyčejné diferenciální rovnice (ODE) vyššího stupně než prvého. Každou obyčejnou diferenciální rovnici n‐tého stupně lze převést na soustavu n diferenciálních rovnic prvního stupně, kde buď přímo jedna z funkcí, nebo kombinace několika nových funkcí představuje hledané řešení původní rovnice. Převod původní diferenciální rovnice vyššího řádu na soustavu diferenciálních rovnic řádu prvého není jednoznačný tj. k jedné ODE vyššího řádu lze nalézt (nekonečně) mnoho různých soustav ODE řádu prvého. Ukážeme si jeden z možných způsobů převodu, kdy jedna funkce je přímo ře‐ šením původní rovnice, a ostatní funkce jsou derivacemi řešení. Obecně lze každou ODE stupně N s počátečními podmínkami zapsat jako

⎛ dNy dy d 2 y d N −1 y ⎞ ⎜ ⎟ = f x y K , , , , , ⎜ dx N dx dx 2 dx N −1 ⎟⎠ ⎝ y ( x0 ) = y0

dy = y1 dx x= x0

d2y dx 2

= y2 x = x0

d N −1 y K dx N −1

= y N −1 x = x0

Zavedeme‐li substituce (nové funkce v1(x), …, vN‐1(x)) jako

v1 ( x0 ) =

dy = v1 dx dv1 = v2 dx M dv N −2 = v N −1 dx

⎛ d2y ⎞ ⎜⎜ = 2 ⎟⎟ ⎝ dx ⎠

dy dx

= y1 x = x0

d2y v2 ( x0 ) = 2 dx

= y2 x = x0

⎛ d N −1 y ⎞ d N −1 y ⎜⎜ = N −1 ⎟⎟ v N −1 ( x0 ) = N −1 dx ⎝ dx ⎠

= y N −1 x = x0

a původní rovnici s využitím substitucí přepíšeme jako

dv N −1 = f ( x , y , v1 , v 2 , K , v N −1 ) dx získáme soustavu N rovnic pro N funkcí včetně počátečních podmínek.


7‐8 (20) 4.7.12


= v1 v1 ( x0 ) = v2 ( x0 ) = dx = v2 M dvN − 2 vN −1 ( x0 ) = dx = v N −1 dvN −1 y ( x0 ) = dx = f ( x, y , v1 , K , v N −1 ) dy

dx

dv1

y1 y2

y N −1 y0

Postup ukažme na příkladu ODE 3. stupně ve tvaru

(

)

y′′′ + 1 − x 2 y′′ + ( y′) + xy = sin( x ) 2

y (0) = 1 y′(0) = 2

y′′(0) = 3

Odpovídající soustava tří diferenciálních rovnic je po substituci

dy

= v1

v1 ( 0 ) = 2

dx

= v2

v2 (0) = 3

dx

= sin( x ) − xy − v12 − 1 − x 2 v 2

dx dv1 dv 2

(

)

y (0) = 1

7.2.2 Řešení soustavy diferenciálních rovnic prvního řádu Soustavu N diferenciálních rovnic pro N funkcí odpovídající ODE N‐tého řádu můžeme zapsat v maticovém tvaru např. jako

⎡ y ⎤ ⎢ v ⎥ ⎢ 1 ⎥ Y=⎢ M ⎥ ⎢ ⎥ ⎢v N − 2 ⎥ ⎢⎣ vN −1 ⎥⎦

v1 ⎡ dy dx ⎤ ⎡ ⎤ ⎡ y0 ⎤ ⎢ dv1 ⎥ ⎢ ⎥ ⎢ y ⎥ v2 dx ⎥ ⎢ ⎢ ⎥ ⎢ 1 ⎥ dY ⎥ Y0 = ⎢ M ⎥ = ⎢ M ⎥ = f ( x, Y ) = ⎢ M dx ⎢dv ⎥ ⎢ ⎥ ⎢ ⎥ N −2 vN −1 dx ⎥ ⎢ ⎥ ⎢ ⎢ y N −2 ⎥ ⎢⎣ dvN −1 dx ⎥⎦ ⎢⎣ f ( x, y, v1 ,K, vN −1 )⎥⎦ ⎢⎣ y N −1 ⎥⎦

Pro výše uvedený příklad obyčejné diferenciální rovnice třetího stupně je pak maticový zápis – je potřeba zvolit pořadí funkcí ve sloupcovém vektoru Y a pak ho dále dodržovat – např. pro předchozí příklad

⎡ y⎤ ⎡1 ⎤ ⎢ ⎥ Y = ⎢ v1 ⎥ Y0 = ⎢⎢2⎥⎥ ⎢⎣v2 ⎥⎦ ⎢⎣3⎥⎦ ⎡ dy dx = v1 ⎤ dY ⎢ dv1 ⎥ = ⎢ dx = v2 ⎥ = f ( x, Y ) dx 2 2 dv 2 ⎢⎣ dx = sin( x) − xy − v1 − (1 − x )v2 ⎥⎦ a jeho přepis do formy m‐funkce MATLABu je Funkce.m ======================================================= function dY=Funkce(x,Y) % zápis soustavy diferenciálních rovnic % pomocné proměnné pro přehlednost zápisu


7‐9 (20) 4.7.12


y=Y(1); v1=Y(2); v2=Y(3); % vektor derivací musí být sloupcový dY=zeros(3,1); dY(1)=v1; dY(2)=v2; dY(3)=sin(x)-x*y-v1^2-(1-x^2)*v2; ======================================================= Funkce.m

Máme‐li diferenciální rovnici ve tvaru soustavy N rovnic pro N neznámých funkcí a známe‐li počáteční podmínky, je získání numerického řešení v MATLABu triviální. Použijeme některou z funkcí pro řešení obyčejných dife‐ renciálních rovnic – nejčastěji5 se používá funkce ode45 se základní syntaxí (blíže viz help ode45) [x,YN]=ode45(fce,X,Y0)

kde

dY=fce(x,Y) je funkce pro výpočet vektoru derivací dY na základě aktuální hodnoty nezávisle proměnné x a vektoru Y aktuálních hodnot hledaných funkcí, X=[x0, (x1, ..., xm,) xmax] je vektor minimálně dvou hodnot (počáteční a konečná) nezávisle proměnné, mezi kterými se má počítat řešení, případně vektor hodnot, ve kterých chceme získat řešení, Y0 je sloupcový vektor počátečních podmínek tj. hodnot jednotlivých funkcí pro hodnotu x0, x je sloupcový vektor6 hodnot nezávisle proměnné, ve kterých jsou počítány hodnoty funkcí uložené v odpovídajícím řádku matice YN (viz Tabulka 1) a YN je matice řešení, jednotlivé sloupce odpovídají hodnotám příslušných funkcí (viz Tabulka 1) Tabulka 1 Uspořádaní výstupních parametrů funkce ode45

x0 x1 ... xm xmax x

YN YN(:,1) YN(:,2) YN(:,3) y(x0) v1(x0) v2(x0) y(x1) v1(x1) v2(x1) ... ... ... y(xm) v1(xm) v2(xm) y(xmax) v1(xmax) v2(xmax)

... ... ... ... ...

YN(:,N) vN‐1(x0) vN‐1(x1) ... vN‐1(xm) vN‐1(xmax)

Význam sloupců matice řešení YN (průběhy jednotlivých funkcí představujících řešení) je dán volbou uspořádání vektoru Y při přípravě soustavy rovnic. Záleží pouze na uživateli, jak funkce uspořádá. Uspořádáním vektoru Y je dáno upo‐ 5

Funkce ode45, která implementuje jednokrokový algoritmus Runge‐Kutta 4‐5. řádu, se dopo‐ ručuje použít jako první volbu. Pokud nelze nalézt řešení nebo víme, že jde o ODE se speciálními vlastnostmi (např. tzv. stiff systémy) jsou k dispozici další funkce využívající jiné algoritmy řešení: ode23, ode113, ode15s, ode23s, ode23t a ode23tb.

6

Počet hodnot mezi počáteční a konečnou hodnotou nezávisle proměnné si určí výpočetní algo‐ ritmus sám v závislosti na průběhu řešení tak, aby bylo dosaženo zadané přesnosti a zároveň byl počet bodů minimální (tj. i minimální počet volání funkce výpočtu derivací = minimalizace výpo‐ četního času). V případě, že je zadán vektor hodnot nezávisle proměnné probíhá výpočet opět v bodech, který si algoritmus určí sám a do zadaných bodů jsou výsledky interpolovány.


7‐10 (20) 4.7.12


řádání (význam) vektoru derivací dY i vektoru počátečních podmínek Y0. První řádek matice YN představuje počáteční podmínky uložené ve vektoru Y0. První a poslední hodnota vektoru x nezávisle proměnných je dána parametrem X při volání funkce ode45. Numerické řešení výše uvedeného příkladu

(

)

y′′′ + 1 − x 2 y′′ + ( y′) + xy = sin( x ) 2

y (0) = 1 y′(0) = 2

y′′(0) = 3

pro 0 ≤ x ≤ 1.5 s využitím m‐funkce Funkce.m včetně vykreslení průběhu všech tří funkcí lze získat příkazy MATLABu >> >> >> >> >> >> >> >>

Y0=[1;2;3]; [x,Y]=ode45('Funkce',[0,1.5],Y0); h=plot(x,Y,'LineWidth',2), grid set(h(1),'LineWidth',3) h=legend('y(x)','v_1=dy/dx','v_2=d^2y/dx^2'); set(h,'Location','SouthWest') hold on plot([0;0;0],Y0,'o','LineWidth',2)

a prvních deset hodnot (ze 45 celkem) numeric‐ kého řešení je 4

x 0 0.0215 0.0431 0.0646 0.0861 0.1236 0.1611 0.1986 0.2361 0.2736

y 1.0000 1.0437 1.0888 1.1351 1.1826 1.2679 1.3562 1.4470 1.5400 1.6347

v1 2.0000 2.0630 2.1226 2.1790 2.2319 2.3160 2.3896 2.4525 2.5048 2.5465

v2 3.0000 2.8482 2.6942 2.5382 2.3806 2.1024 1.8209 1.5373 1.2527 0.9687

2

0

-2

-4

-6

y(x) v1=dy/dx Obrázek 7‐6 Průběh řešení diferenciální rovnice 2 2 v2=d y/dx

-8


0

7‐11 (20) 4.7.12

0.5

1

1.5


7.3 Příklady na procvičení V kapitole 7.3.1 jsou zadání a v kapitole 7.3.2 řešení jednoduchých příkladů na probíranou problematiku. Pro řešení příkladů jsou nutné znalosti z prvních čtyř bloků. Předpokládá se, že při řešení příkladů má čtenář k dispozici MATLAB a aktivně využívá nápovědu k programu.

7.3.1 Zadání příkladů Př01. Jednoduchý určitý integrál – nalezněte řešení 2π

∫

1

b) cos(n arccos( x ) )dx pro n=1, 2, 3, 4, 5

∫

2

a) cos ( x)dx

−1

0

Př02. Dvojný určitý integrál – nalezněte řešení π π

a)

∫π ∫π

1 + sin( x) cos( y )dxdy

− −

b)

∫∫ x sin( y) + y cos( x)dxdy kde oblast S je popsaná rovnicí x

2

+ y 2 ≤ 1 .

S

Př03. Trojný určitý integrál – nalezněte řešení π π π

a)

⎛ xyz

∫π ∫π ∫π cos⎜⎝

− − −

⎞⎟dxdydz

π3⎠

b) Pomocí trojného integrálu určete objem krychle o délce strany 2. Návod: využijte vzorec

∫∫∫ dxdydz pro výpočet objemu a integrační oblast V

tvaru krychle Př04. Využití určitých integrálů. a) Výpočet délky křivky. Vypočtěte délku spirály popsané parametrickými rov‐ nicemi x(t ) = t cos(t ) y (t ) = t sin(t ) pro 0 ≤ t ≤ 2π. Návod: využijte vztah pro výpočet délky křivky popsané parametrickými rovni‐ cemi L =

t2

∫ t1

2

2

⎛ dx ⎞ ⎛ dy ⎞ ⎜ ⎟ + ⎜ ⎟ dt . ⎝ dt ⎠ ⎝ dt ⎠

b) Výpočet objemu rotačního tělesa. Vypočtěte objem tělesa vzniklého rotací křivky y ( x) = xe − x kolem osy x pro 0 ≤ x ≤ 5. Návod: využijte vztah pro výpočet objemu rotačního tělesa, jehož plášť vznikne rotací křivky y=f(x) kolem osy x v mezích a ≤ x ≤ b V = π

b

∫f

2

( x)dx .

a


7‐12 (20) 4.7.12


c) Výpočet plochy rovinného útvaru. Vypočtěte plochu elipsoidu popsaného rovnicí

x2 y2 + = 1 pomocí dvojného integrálu a výsledek ověřte výpočtem a 2 b2

s jednoduchým integrálem pro ekvivalentní popis parametrickými rovnicemi x(t ) = a cos(t ) y (t ) = b sin(t ) 0 ≤ t ≤ 2π pro a=2 a b=3. Návod: využijte vztahy pro výpočet plochy rovinného obrazce nad danou uza‐ vřenou oblastí O(x,y) S =

t2

∫∫ dxdy a S = ∫ f (t ) g ′(t )dt pokud je oblast da‐

t1 ná parametrickými rovnicemi x=f(t) a y=g(t) O( x, y )

Př05. Diferenciální rovnice – nalezněte řešení a nakreslete jeho průběh a)

dx + x2 = t 2 dt

x(t = 0) = 1 0 ≤ t ≤ 10

d 2z dz dz + x 2 = 2 z ( x = 0) = 1 = −1 0 ≤ x ≤ 6 b) 2 + x dx dx dx x=0 y′′ + y′ − x = 2te − t c) x′′ + 2 x′ + y = te −t

y (t = 0) = +1 y′(t = 0) = 2 x(t = 0) = −1 x′(t = 0) = 0 0≤t ≤5

Př06. Dráha střely – využití diferenciální rovnice Nakreslete dráhu střely v gravitačním poli země (g=9.81 m.s‐2)vystřelené pod úhlem α v prvních 10 s po opuštění hlavně (ústí hlavně je h0=1 m nad povr‐ ‐1 h chem, úsťová rychlost v0=800 m.s ) při uvažování odporu prostředí úměrného kvadrátu rychlosti (konstanta úměrnosti je k=0.004 m‐1). Průběh dráhy nakreslete pro hodnoty α=20°, 30°, 40°, v0 50°, 60° a 70° do jednoho grafu. Časový průběh výškové α (h) i délkové (x) souřadnice je x popsán dife‐ h0 renciálními rovnicemi ve tvaru

d 2h dh dh +k + g = 0 h(t = 0) = h0 2 dt dt dt d 2x dx dx +k + g = 0 x(t = 0) = 0 2 dt dt dt


dh = v0 sin(α ) dt t =0 dx = v0 cos(α ) dt t =0

7‐13 (20) 4.7.12


Řešení příkladů Př01. Jednoduchý určitý integrál – nalezněte řešení 2π

∫

2 a) cos ( x)dx =3.1416 0

>> I=quadl(@(x) cos(x).^2,0,2*pi) I = 3.1416 1

b) cos(n arccos( x ) )dx pro n=1, 2, 3, 4, 5

∫

−1

>> ans ans ans ans ans

for n=1:5,quadl(@(x) cos(n*acos(x)),-1,1), end = 1.1126e-016 = -0.6667 = 9.6435e-017 = -0.1333 = 1.9547e-016

Př02. Dvojný určitý integrál – nalezněte řešení π π

a)

∫∫

1 + sin( x) cos( y )dxdy =37.8239

−π −π

>> I=dblquad(@(x,y) sqrt(1+sin(x)*cos(y)),-pi,pi,-pi,pi) I = 37.8239

b)

∫∫ xy cos(1 + xy)dxdy kde oblast S je popsaná rovnicí x

2

+ y2 ≤ 1

S

Dvojný integrál přes oblast S lze integrovat přes libovolnou obdélníkovou ob‐ last (zahrnující oblast S) pokud při vyhodnocení funkce pod integrálem je za‐ jištěno, že tato vrací hodnotu 0 pro hodnoty x a y mimo oblast S. Pro body v oblasti S platí ‐1 ≤ x ≤ 1 a současně y2 ≤ 1‐x2. Vytvoříme‐li tedy např. funkci Funkce.m ===================================================== function z=Funkce(x,y) % zápis funkce pod integrálem s kontrolou definičního oboru % parametr x může být vektor n=length(x); z=zeros(size(x)); for k=1:n, if (x(k)^2<=1) && (y^2<=1-x(k)^2), z(k)=x(k)*y*cos(1+x(k)*y); else z(k)=0; end end ===================================================== Funkce.m

pak se řešení nalezne příkazem (pro libovolné meze zahrnující oblast S) >> I=dblquad(@Funkce,-2,3,-3,5) I = -0.1080


7‐14 (20) 4.7.12


Př03. Trojný určitý integrál – nalezněte řešení π π π

a)

⎛ xyz

∫π ∫π ∫π cos⎜⎝

− − −

⎞⎟dxdydz =243.5384

π3⎠

>> I=triplequad(@(x,y,z) cos(x*y*z/pi^3),-pi,pi,-pi,pi,-pi,pi) I = 243.5384

b) Pomocí trojného integrálu určete objem krychle o délce strany 2. Funkce pod integrálem musí pro všechny hodnoty x, y a z vracet hodnotu 1 a zároveň musí pro vektor hodnot x vracet vektor funkčních hodnot. Řešení tedy je např. >> I=triplequad(@(x,y,z) ones(size(x)),0,2,0,2,0,2) I = 8

Př04. Využití určitých integrálů. a) Vypočtěte délku spirály x(t ) = t cos(t )

y (t ) = t sin(t ) pro 0 ≤ t ≤ 2π.

výpočet potřebných derivací

d (t cos(t ) ) = cos(t ) − t sin(t ) dt

d (t sin(t ) ) = sin(t ) + t cos(t ) dt

výpočetní vztah t2

L=∫ t1

2

2

⎛ dx ⎞ ⎛ dy ⎞ ⎜ ⎟ + ⎜ ⎟ dt = ⎝ dt ⎠ ⎝ dt ⎠

2π

∫ [cos(t ) − t sin(t )] + [sin(t ) + t cos(t )] dt 2

2

0

příkazy MATLABu >> f=@(t) sqrt((cos(t)-t.*sin(t)).^2+(sin(t)+t.*cos(t)).^2); >> L=quadl(f,0,2*pi) L = 21.2563

b) Výpočet objemu rotačního tělesa. Vypočtěte objem tělesa vzniklého rotací křivky y ( x) = xe − x kolem osy x pro 0 ≤ x ≤ 5. b

5

( )

V = π ∫ f ( x)dx = π ∫ xe− x dx 2

a

0

2

příkazy MATLABu >> V=pi*quadl(@(x) (x.*exp(-x)).^2,0,5) V = 0.7832

c) Výpočet plochy rovinného útvaru O(x,y):

x2 y2 + = 1 pomocí dvojného a2 b2

integrálu a jednoduchým integrálem pro ekvivalentní popis parametrickými rovnicemi x (t ) = a cos(t ) y (t ) = b sin(t ) 0 ≤ t ≤ 2π pro a=2 a b=3. Chceme‐li použít dvojný integrál a oblast O(x,y) není obdélníková, je nutné zajistit, aby pomocná funkce pod integrálem vracela hodnotu 1 pouze pro

(

2 2 2 2 hodnoty x a y z oblasti O (jinak hodnotu 0) tj. x ≤ a ∧ y ≤ b 1 − x

2

a2

).


7‐15 (20) 4.7.12


Pomocná funkce s parametry a,b Funkce.m ===================================================== function z=Funkce(x,y,a,b) % zápis funkce pod integrálem s kontrolou definičního oboru % parametr x může být vektor n=length(x); z=zeros(size(x)); for k=1:n, if (x(k)^2<=a^2) && (y^2<=b^2*(1-(x(k)/a)^2), z(k)=1; else z(k)=0; end end ===================================================== Funkce.m

příkazy MATLABu >> a=2; b=3; >> S1=dblquad(@(x,y) Funkce(x,y,a,b),-a,a,-b,b) S1 = 18.8495

Při použití jednoduchého integrálu a parametricky zadané oblasti je nutné určit derivaci g ′ = b cos(t ) a výpočetní vztah pak je S =

2π

∫ ab cos (t )dt 2

0

příkazy MATLABu >> S2=quadl(@(t) a*b*cos(t).^2,0,2*pi) S2 = 18.8496

Př05. Diferenciální rovnice – nalezněte řešení a nakreslete jeho průběh a)

dx + x2 = t 2 dt

x(t = 0) = 1 0 ≤ t ≤ 10

10

9

příkazy MATLABu

8

>> [t,X]=ode45(@(t,x) t^2-x^2,[0,10],1); >> plot(t,X) >> grid

7

6

5

d 2z dz + x + x2 = 2 2 dx dx dz = −1 b) z ( x = 0) = 1 dx x=0

4

3

2

1

0

0≤ x≤6

0

1

2

3

4

5

6

7

8

9

10

převod na soustavu dvou rovnic řádu prvního

⎡z⎤ ⎡ z (0) = 1 ⎤ Z = ⎢ ⎥ Z0 = ⎢ ⎥ ⎣v ⎦ ⎣v(0) = −1⎦ ⎤ dZ ⎡ z′ = v = dx ⎢⎣v′ = 2 − x 2 − xv ⎥⎦ KŘP/IMSW Modelování ve výpočtových software

7‐16 (20) 4.7.12


zápis ve formě m‐funkce Funkce.m =========================== function dZ=Funkce(x,Z) % zápis soustavy dif. rovnic

2 z v 0

% pomocné proměnné pro přehlednost z=Z(1); v=Z(2); dZ=[v; 2-x^2-x*v]; =========================== Funkce.m

-2

-4

-6

příkazy MATLABu >> >> >> >> >> >>

-8

Z0=[1;-1]; [x,Z]=ode45('Funkce',[0,6],Z0); h=plot(x,Z); grid legend('z','v') set(h(1),'LineWidth',3)

y′′ + y′ − x = 2te − t c) x′′ + 2 x′ + y = te

−t

-10

-12

0

1

2

3

4

5

6

y (t = 0) = +1 y′(t = 0) = 2 x(t = 0) = −1 x′(t = 0) = 0 0≤t ≤5

převod na soustavu dvou rovnic řádu prvního

⎡y⎤ ⎢y ⎥ Z = ⎢ 1⎥ ⎢x⎥ ⎢ ⎥ ⎣ x1 ⎦

⎡ y ( 0) = 1 ⎤ ⎢ y (0) = 2 ⎥ ⎥ Z0 = ⎢ 1 ⎢ x(0) = −1⎥ ⎥ ⎢ ⎣ x1 (0) = 0 ⎦

⎡ y′ ⎢ dZ ⎢ y1′ = dx ⎢ x′ ⎢ ⎣ x1′

= y1 ⎤ −t = 2te + x − y1 ⎥⎥ ⎥ = x1 ⎥ = te −t − y − 2 x1 ⎦

zápis ve formě m‐funkce 3

Funkce.m =========================== function dZ=Funkce(t,Z) % zápis soustavy dif. rovnic

2

1

% pomocné proměnné pro přehlednost y=Z(1); y1=Z(2); x=Z(3); x1=Z(4); dZ=zeros(4,1); dZ(1)=y1; dZ(2)=2*t*exp(-t)+x-y1; dZ(3)=x1; dZ(4)=t*exp(-t)-y-2*x1; =========================== Funkce.m

0

-1

-2

-3

příkazy MATLABu >> >> >> >> >> >>

y y1 x x1

-4

0

0.5

1

1.5

2

2.5

3

3.5

4

4.5

5

Z0=[1;2;-1;0]; [x,Z]=ode45(@Funkce,[0,5],Z0); h=plot(x,Z); grid legend('y','y1','x','x1') set(h([1,3]),'LineWidth',3)


7‐17 (20) 4.7.12


Př06. Dráha střely – využití diferenciální rovnice Nakreslete dráhu střely v gravitačním poli země (g=9.81 m.s‐2)vystřelené pod úhlem α v prvních 10 s po opuštění hlavně (ústí hlavně je h0=1 m nad povr‐ chem, úsťová rychlost v0=800 m.s‐1) při uvažování odporu prostředí úměrné‐ ho kvadrátu rychlosti (konstanta úměrnosti je k=0.004 m‐1). Průběh dráhy nakreslete pro hodnoty α=20°, 30°, 40°, 50°, 60° a 70° do jednoho grafu. Ča‐ sový průběh výškové (h) i délkové (x) souřadnice je popsán diferenciálními rovnicemi ve tvaru

d 2h dh dh +k + g = 0 h(t = 0) = h0 2 dt dt dt d 2x dx dx +k + g = 0 x(t = 0) = 0 2 dt dt dt

dh = v0 sin(α ) dt t =0 dx = v0 cos(α ) dt t =0

Převod na soustavu čtyř diferenciálních rovnic prvního stupně

⎡h⎤ ⎡h(0) = h0 ⎤ ⎢v ⎥ ⎢v (0) = v sin(α ) ⎥ 0 ⎥ Z = ⎢ h ⎥ Z0 = ⎢ h ⎢x⎥ ⎢ x(0) = 0 ⎥ ⎢ ⎥ ⎢ ⎥ ⎣v x ⎦ ⎣ x1 (0) = v0 cos(α )⎦

⎡h′ ⎢ dZ ⎢v′h = dx ⎢ x′ ⎢ ⎣v′x

= vh

⎤ = − g − k vh vh ⎥⎥ ⎥ = vx ⎥ = − k vx vx ⎦

zápis ve formě m‐funkce Funkce.m ===================================================== function dZ=Funkce(t,Z,k) % zápis soustavy diferenciálních rovnic % s parametrem k % pomocné proměnné pro přehlednost zápisu h=Z(1); vh=Z(2); x=Z(3); vx=Z(4); g=9.81; % vektor derivací musí být sloupcový dZ=zeros(4,1); dZ(1)=vh; dZ(2)=-g-k*abs(vh)*vh; dZ(3)=vx; dZ(4)=-k*abs(vx)*vx; ==================================================== Funkce.m

příkazy MATLABu (skript) skript.m ===================================================== % skript pro výpočet dráhy střely v0=800; % úsťová rychlost ko=0.004; % odpor prostředí h0=1; % výška ústí hlavně x0=0; % x-ová poloha ústí (vzdálenost) al=[20,30,40,50,60,70]; % náměr ve stupních % příprava pro vykreslení drah střely axes % prázdný graf hold on % režim překreslování te=cell(size(al)); % pro text legendy % pole buněk - řetězce pro barvy grafů co={'b','g','r','c','m','k','b','g','r','c','m','y','k'};


7‐18 (20) 4.7.12


for k=1:length(al), alr=al(k)*pi/180; % převod na radiány vh0=v0*sin(alr); % poč. rychlost v h vx0=v0*cos(alr); % poč. rychlost v x Z0=[h0;vh0;x0;vx0]; % poč. podmínky [t,Z]=ode45(@(t,Z) Funkce(t,Z,ko), [0,20],Z0); x=Z(:,3); % časový průběh dráhy (vzdálenost) y=Z(:,1)'; % časový průběh dráhy (výška) hp=plot(x,y,co{k},'LineWidth',2); te{k}=sprintf('uhel %2d ^o',al(k)); end grid legend(te,'Location','NorthWest') ==================================================== skript.m

700 uhel 20 o uhel 30 o 600

uhel 40 o uhel 50 o uhel 60 o

500

uhel 70 o 400

300

200

100

0

-100

0

200

400

600

800

1000

1200


7‐19 (20) 4.7.12


Pojmy k zapamatování Okruhy problémů: výpočet určitých integrálů (jednoduchý, dvojný a troj‐ ný), výpočet plochy a objemu, numerické řešení obyčejných diferenciál‐ ních rovnic Použité nástroje:

Příkazy a funkce:

dblquad, ode45, quadl, triplequad

Otázky na procvičení 1.

2. 3. 4. 5.

Uvádí se, že geometrická interpretace určitého integrálu nějaké funkce je plocha shora ohraničená danou funkcí a zdola osou x. V čem je tato před‐ stava nepřesná? Funkce MATLABu pro výpočet dvojného a trojného určitého integrálu nej‐ sou úplně obecné. Jaké mají omezení? Jaký je rozdíl mezi analytickým a numerickým řešením diferenciální rovni‐ ce? Proč je potřeba pro numerické řešení v MATLABu převést ODE vyššího řádu na soustavu ODE řádu prvého? Co je nutné kromě diferenciální rovnice ještě znát pro výpočet numerické‐ ho řešení?

Odkazy a další studijní prameny

•

on‐line dokumentace k programu nebo www.mathworks.com/help/ – části „Numerical Integration (Quadrature)“, „Ordinary Differential Equations“


7‐20 (20) 4.7.12


7 Výpočet určitého integrálu a numerické řešení diferenciálních

Recommend Documents