Seznámení se systémem Mathematica 7

Základy numerických metod

© Jaroslav Reichl, SPŠST Panská, Praha
reichl@panska.cz

Obsah

Prezentace obsahuje tato témata:

Metody řešení nelineárních rovnic
25. úkol
Metody dynamického modelování
26. úkol

Jaroslav Reichl, SPŠST Panská, Praha, © 2011

Metody řešení nelineárních rovnic

Metoda půlení intervalu (bisekce)
	1) zadanou rovnici upravíme do tvaru f(x)=0;
	2) zvolíme interval s krajními body a a b, na kterém hledáme kořen (Subscript[x, 0]) rovnice f(x)=0 (např. na základě vykreslení grafu funkce f(x));
	3) interval musí být zvolen tak, aby v něm kořen (Subscript[x, 0]) rovnice f(x)=0 určitě ležel, tj. musí platit f(a).f(b)<0;
	4) určíme střed c intervalu <a; b>: c=((a+b)/2);
	5) dále postupujeme podle podle bodů 3) a 4):
		je-li f(a).f(c)<0, pak zvolíme b=c
		je-li f(c).f(b)<0, pak zvolíme a=c a pokračujeme ve výpočtu analogicky dále (dle bodu 4) určíme nový střed intervalu, najdeme nový interval, kde bude ležet kořen, ...);
	6) celý postup popsaný v bodech 3), 4) a 5) opakujeme do té doby, až bude |f(c)|<ε (kde ε je zvolená přesnost hledání kořene); pak (Subscript[x, 0])=c.

Metoda sečen
	1) zadanou rovnici upravíme do tvaru f(x)=0;
	2) zvolíme interval s krajními body a a b, na kterém hledáme kořen (Subscript[x, 0]) rovnice f(x)=0 (např. na základě vykreslení grafu funkce f(x));
	3) interval musí být zvolen tak, aby v něm kořen (Subscript[x, 0]) rovnice f(x)=0 určitě ležel, tj. musí platit f(a).f(b)<0;
	4) napíšeme rovnici sečny s grafu funkce f, která prochází body A=[a; f(a)] a B=[b; f(b)]; směrnice této přímky tedy je k=((f(b)-f(a))/(b-a)), a tedy její rovnice je s:y=k x+q=((f(b)-f(a))/(b-a))x+((f(a).b-f(b).a)/(b-a)) (což je možné nechat vypočítat systém Mathematica);
	5) najdeme průsečík c sečny a osy x, tj. řešíme lineární rovnici s(c)=0;
	6) dále postupujeme podle podle bodů 3), 4) a 5):
		je-li f(a).f(c)<0, pak zvolíme b=c
		je-li f(c).f(b)<0, pak zvolíme a=c a pokračujeme ve výpočtu analogicky dále (dle bodu 4) napíšeme rovnici nové sečny, najdeme nový průsečík c, určíme nový interval, kde bude ležet kořen, ...);
	7) celý postup popsaný v bodech 3), 4), 5) a 6) opakujeme do té doby, až bude |f(c)|<ε (kde ε je zvolená přesnost hledání kořene); pak (Subscript[x, 0])=c.

Metoda tečen
	1) zadanou rovnici upravíme do tvaru f(x)=0;
	2) zvolíme bod a, který je blízko hledaného kořene (Subscript[x, 0]) rovnice f(x)=0 - např. na základě vykreslení grafu funkce f(x));
	3) napíšeme rovnici tečny t grafu funkce f v jejím bodě A=[a; f(a)]; směrnice této přímky tedy je k=f'(a), a tedy její rovnice je t:y=k x+q=f'(a)x+f(a)-f'(a).a (což je možné nechat vypočítat systém Mathematica);
	4) najdeme průsečík c tečny a osy x, tj. řešíme lineární rovnici t(c)=0;
	5) dále postupujeme podle podle bodů 3) a 4) - tj. dosadíme a=c a napíšeme rovnici nové tečny, najdeme její průsečík s osou x, ...;
	6) celý postup popsaný v bodech 3), 4) a 5) opakujeme do té doby, až bude |f(c)|<ε (kde ε je zvolená přesnost hledání kořene); pak (Subscript[x, 0])=c.

Komentáře
	Metoda půlení intervalu může konvergovat k hledanému řešení (Subscript[x, 0]) velmi pomalu (závisí na typu funkce f(x). K řešení ale dokonverguje.
	Metoda sečen konverguje rychleji, metoda tečen konverguje nejrychleji. Tyto metody se musí ale používat velmi opatrně, neboť se může stát, že najdou úplně jiný kořen, než který jsme najít chtěli. Závisí na typu funkce f(x) a zvoleném počátečním intervalu (resp. bodu), od kterého začínáme kořen hledat.

Jaroslav Reichl, SPŠST Panská, Praha, © 2011

25. úkol: Metody řešení nelineárních rovnic

Najděte kladný kořen rovnice (x^4)+3(x^3)-13(x^2)-15x=20 metodou půlení intervalu.

Řešení

ClearAll["Global`*"];
 barva=Black;
barvafce=Red;
barvabodu=Blue;
tloustka=0.007;
velikostbodu=0.025;
presnost=(10^(-4));

funkce[x_]:=(x^4)+3(x^3)-13(x^2)-15x-20;
graffce=Plot[funkce[x],{x,-6,5},PlotStyle->{barvafce,Thickness[tloustka]}];
pocet=0;
a=0;
b=4;
c=((a+b)/2);
While[Abs[funkce[c]]>presnost,
pocet=pocet+1;
If[funkce[a]funkce[c]<0,
b=c,
If[funkce[b]funkce[c]<0,
a=c,
Break[]]];
c=((a+b)/2)];
bod=Point[{c,0}];
Show[
graffce,
Graphics[{
{barvabodu,PointSize[velikostbodu],bod}
}],
AxesOrigin->{0,0},AxesLabel->{Style["x",15,barva,FontFamily->Times],Style["y",15,barva,FontFamily->Times]}]
Print["Hledaný kořen je ",N[c,Round[Abs[Log[10,presnost]]]]," a byl nalezen po ",pocet," iteracích."]

Najděte kladný kořen rovnice (x^4)+3(x^3)-13(x^2)-15x=20 metodou sečen.

Řešení

ClearAll["Global`*"];
 barva=Black;
barvafce=Red;
barvasecny=Orange;
barvabodu=Blue;
tloustka=0.007;
velikostbodu=0.025;
presnost=(10^(-4));

funkce[x_]:=1.(x^4)+3(x^3)-13(x^2)-15x-20;
graffce=Plot[funkce[x],{x,-6,5},PlotStyle->{barvafce,Thickness[tloustka]}];
pocet=0;
a=0;
b=4;
c=a;
While[Abs[funkce[c]]>presnost,
pocet=pocet+1;
If[funkce[a]funkce[c]<0,
b=c,
If[funkce[b]funkce[c]<0,
a=c,
Break[]]];
secnapom[x_]:=((funkce[b]-funkce[a])/(b-a))x+qpom;
reseni=Solve[secnapom[a]==funkce[a],qpom];
q=qpom/.reseni[[1,1]];
secna[x_]:=secnapom[x]/.qpom->q;
reseni=Solve[secna[cpom]==0,cpom];
c=cpom/.reseni[[1,1]]
];
bod=Point[{c,0}];
Show[
graffce,
Plot[secna[x],{x,-6,5},PlotStyle->{barvasecny,Thickness[tloustka]}],
Graphics[{
{barvabodu,PointSize[velikostbodu],bod}
}],
AxesOrigin->{0,0},AxesLabel->{Style["x",15,barva,FontFamily->Times],Style["y",15,barva,FontFamily->Times]}]
Print["Hledaný kořen je ",Round[c,1.presnost]," a byl nalezen po ",pocet," iteracích."]

Najděte kladný kořen rovnice (x^4)+3(x^3)-13(x^2)-15x=20 metodou tečen.

Řešení

ClearAll["Global`*"];
 barva=Black;
barvafce=Red;
barvatecny=Magenta;
barvabodu=Blue;
tloustka=0.007;
velikostbodu=0.025;
presnost=(10^(-4));

funkce[x_]:=1.(x^4)+3(x^3)-13(x^2)-15x-20;
graffce=Plot[funkce[x],{x,-6,5},PlotStyle->{barvafce,Thickness[tloustka]}];
pocet=0;
a=4;
While[Abs[funkce[a]]>presnost,
pocet=pocet+1;
tecnapom[x_]:=funkce'[a]x+qpom;
reseni=Solve[tecnapom[a]==funkce[a],qpom];
q=qpom/.reseni[[1,1]];
tecna[x_]:=tecnapom[x]/.qpom->q;
reseni=Solve[tecna[apom]==0,apom];
a=apom/.reseni[[1,1]]
];
bod=Point[{a,0}];
Show[
graffce,
Plot[tecna[x],{x,-6,5},PlotStyle->{barvatecny,Thickness[tloustka]}],
Graphics[{
{barvabodu,PointSize[velikostbodu],bod}
}],
AxesOrigin->{0,0},AxesLabel->{Style["x",15,barva,FontFamily->Times],Style["y",15,barva,FontFamily->Times]}]
Print["Hledaný kořen je ",Round[a,1.presnost]," a byl nalezen po ",pocet," iteracích."]

Jaroslav Reichl, SPŠST Panská, Praha, © 2011

Metody dynamického modelování

Použití metody
Metody dynamického modelování slouží k numerickému řešení diferenciálních rovnic. Jsou vhodné k řešení řady úloh z fyziky, elektrotechniky, ... Nejjednodušší z těchto metod je Eulerova metoda, jejímž autorem je švýcarský matematik Leonhard Euler (1707 - 1783).

Popis metody
Aby bylo možné vytvořit dynamický model pohybu tělesa (nebo hmotného bodu), je nutné znát:
	pohybovou rovnici - ta vychází ze druhého Newtonova zákona, který lze psát ve tvaru (Overscript[F, →])=m(Overscript[a, →])=m((d(Overscript[v, →]))/(d t))=m(((d^2)(Overscript[r, →]))/(d (t^2))), přičemž za sílu (Overscript[F, →]) dosazujeme tu sílu, která způsobuje změnu zrychlení tělesa (tíhová, gravitační, odporovou, elektrostatickou, ...);
	počáteční podmínky - tj. polohu, velikost a směr rychlosti pohybu tělesa (resp. hmotného bodu) na začátku vyšetřování pohybu tělesa.

V obecném případě je nutné řešit rovnice ve třech rozměrech (neboť polohu i rychlost pohybu tělesa resp. hmotného bodu vyšetřujeme ve trojrozměrném prostoru) , ale pro naše úvahy a pro seznámení se s metodou použijeme pouze jednorozměrné nebo dvourozměrné případy (pád tělesa v odporujícím prostředí, vrhy těles v homogenním tíhovém poli, ...).
Na základě matematického modelu je určována poloha hmotného bodu a jeho rychlost v závislosti na čase. Čas přitom v modelu neplyne rovnoměrně, ale "skokově" - dvě po sobě jdoucí polohy hmotného bodu jsou vypočteny v časech, které odděluje časový interval Δt. Pro přesnější výpočty je pochopitelně vhodné volit Δt pokud možno co nejmenší. To ovšem na druhé stranu způsobuje zpomalení výpočtu, větší nároky na paměť, ...

Princip metody
Princip metody vysvětlíme na jednorozměrném případě (např. pád tělesa v tíhovém poli Země po úsečce). Metodu lze aplikovat v několika krocích:
	1) výpočet velikosti zrychlení tělesa na základě znalosti síly, tj. (Subscript[a, i+1])=((Subscript[F, i])/m) (velikost síly a tedy i velikost zrychlení může záviset na velikosti rychlosti nebo na poloze - proto se obě tyto veličiny mohou obecně měnit);
	2) výpočet velikosti rychlosti tělesa na základě zrychlení: (Subscript[v, i+1])=(Subscript[v, i])+a Δt (na základě předchozí velikosti rychlosti (Subscript[v, i]) se vypočítá "nová" rychlosti (Subscript[v, i+1]))
	3) výpočet polohy tělesa na základě rychlosti: (Subscript[x, i+1])=(Subscript[x, i])+(Subscript[v, i+1]) Δt;
	4) opakování kroků 1) až 3) do skončení pohybu (dopad na Zem, srážka dvou těles, zastavení tělesa, ...).

Analogicky bychom postupovali i v případě vícedimenzionálního pohybu. V tom případě by bylo nutné počítat více složek rychlosti a více souřadnic pohybujícího se tělesa.

Jaroslav Reichl, SPŠST Panská, Praha, © 2011

26. úkol: Metody dynamického modelování

Pomocí metody dynamického modelování zobrazte graf závislosti výšky na čase a graf závislosti velikosti rychlosti na čase pro volný pád hmotného bodu o hmotnosti 200 g z výšky 20 metrů na vodorovnou položku.

Řešení

ClearAll["Global`*"];
hmax=20;
g=9.81; 
m=0.200;
v0=0;
deltat=0.1;
barva=Black;
barvabodu1=Red;
barvabodu2=Blue;
tloustka=0.007;
velikostbodu=0.02;
velikostpisma=12;

rychlost=Table[Null,{i,500}];
vyska=Table[Null,{i,500}];
vstara=v0;hstara=hmax;
F=m g;
pocet=1;
rychlost[[pocet]]=v0;
vyska[[pocet]]=hmax;
While[vyska[[pocet]]>0,
a=(F/m);
rychlost[[pocet+1]]=rychlost[[pocet]]+a deltat;
vyska[[pocet+1]]=vyska[[pocet]]-rychlost[[pocet+1]] deltat;
pocet=pocet+1];

dataprograf1=Table[{deltat (i-1),vyska[[i]]},{i,1,pocet-1,1}];
dataprograf2=Table[{deltat (i-1),rychlost[[i]]},{i,1,pocet-1,1}];
grafst=ListPlot[dataprograf1,PlotStyle->Directive[barvabodu1,PointSize[velikostbodu]],AxesOrigin->{0,0},AxesLabel->{Style["(t/s)",velikostpisma,barva,FontFamily->Times],Style["(s/m)",velikostpisma,barva,FontFamily->Times]},FormatType->StandardForm]
grafvt=ListPlot[dataprograf2,PlotStyle->Directive[barvabodu2,PointSize[velikostbodu]],AxesOrigin->{0,0},AxesLabel->{Style["(t/s)",velikostpisma,barva,FontFamily->Times],Style["(v/(m.(s^(-1))))",velikostpisma,barva,FontFamily->Times]},FormatType->StandardForm]

Pomocí metody dynamického modelování zobrazte graf závislosti výšky na čase a graf závislosti velikosti rychlosti na čase pro pád tělesa ve tvaru koule o hmotnosti 200 g a poloměru 5 cm z výšky 20 metrů na vodorovnou položku. Předpokládejte, že pro velikost odporové síly působící proti pohybu tělesa, platí Newtonův vztah (Subscript[F, odporová])=(1/2)C S ρ (v^2), kde S je plocha příčného řezu tělesa, C je součinitel odporu (v případě koule je to 0,48), ρ=1,3 kg.(m^(-3)) je hustota vzduchu a v je velikost aktuální rychlosti tělesa.

Řešení

ClearAll["Global`*"];
hmax=20;
g=9.81; 
m=0.200;
r=0.05;
koefC=0.48;
hustota=1.3;
v0=0;
koeficient=(1/2) koefC 4 π (r^2) hustota;
deltat=0.1;
barva=Black;
barvabodu1=Red;
barvabodu2=Blue;
tloustka=0.007;
velikostbodu=0.02;
velikostpisma=12;

rychlost=Table[Null,{i,500}];
vyska=Table[Null,{i,500}];
pocet=1;
rychlost[[pocet]]=v0;
vyska[[pocet]]=hmax;
While[vyska[[pocet]]>0,
F=m g-koeficient ((rychlost[[pocet]])^2);
a=(F/m);
rychlost[[pocet+1]]=rychlost[[pocet]]+a deltat;
vyska[[pocet+1]]=vyska[[pocet]]-rychlost[[pocet+1]] deltat;
pocet=pocet+1];

dataprograf1=Table[{deltat (i-1),vyska[[i]]},{i,1,pocet-1,1}];
dataprograf2=Table[{deltat (i-1),rychlost[[i]]},{i,1,pocet-1,1}];
grafst=ListPlot[dataprograf1,PlotStyle->Directive[barvabodu1,PointSize[velikostbodu]],AxesOrigin->{0,0},AxesLabel->{Style["(t/s)",velikostpisma,barva,FontFamily->Times],Style["(s/m)",velikostpisma,barva,FontFamily->Times]},FormatType->StandardForm]
grafvt=ListPlot[dataprograf2,PlotStyle->Directive[barvabodu2,PointSize[velikostbodu]],AxesOrigin->{0,0},AxesLabel->{Style["(t/s)",velikostpisma,barva,FontFamily->Times],Style["(v/(m.(s^(-1))))",velikostpisma,barva,FontFamily->Times]},FormatType->StandardForm]

Pomocí metody dynamického modelování zobrazte graf závislosti y-ové souřadnice pohybujícího se tělesa na jeho x-ové souřadnici a graf závislosti velikosti rychlosti na čase pro šikmý vrh tělesa ve tvaru koule o hmotnosti 500 g a poloměru 10 cm pod úhlem 35° počáteční rychlostí o velikosti 10 m.(s^(-1)). Předpokládejte, že pro velikost odporové síly, působící proti pohybu tělesa platí Newtonův vztah (Subscript[F, odporová])=(1/2)C S ρ (v^2), kde S je plocha příčného řezu tělesa, C je součinitel odporu (v případě koule je to 0,48), ρ=1,3 kg.(m^(-3)) je hustota vzduchu a v je velikost aktuální rychlosti tělesa.

Řešení

ClearAll["Global`*"];
h0=0;
g=9.81; 
m=0.500;
r=0.1;
koefC=0.48;
hustota=1.3;
v0=10;
alfa=35 Degree;
koeficient=(1/2) koefC 4 π (r^2) hustota;

deltat=0.02;
barva=Black;
barvabodu1=Red;
barvabodu2=Blue;
tloustka=0.007;
velikostbodu=0.02;
velikostpisma=12;

rychlostx=Table[Null,{i,500}];
rychlosty=Table[Null,{i,500}];
rychlost=Table[Null,{i,500}];
x=Table[Null,{i,500}];
y=Table[Null,{i,500}];
pocet=1;
rychlostx[[pocet]]=v0 Cos[alfa];
rychlosty[[pocet]]=v0 Sin[alfa];
rychlost[[pocet]]=(Sqrt[((rychlostx[[pocet]])^2)+((rychlosty[[pocet]])^2)]);
x[[pocet]]=h0;
y[[pocet]]=h0+(10^(-6));
While[y[[pocet]]>h0,
Fx=-koeficient rychlostx[[pocet]] rychlost[[pocet]];
Fy=-m g-koeficient rychlosty[[pocet]] rychlost[[pocet]];
ax=(Fx/m);
ay=(Fy/m);
rychlostx[[pocet+1]]=rychlostx[[pocet]]+ax deltat;
rychlosty[[pocet+1]]=rychlosty[[pocet]]+ay deltat;
x[[pocet+1]]=x[[pocet]]+rychlostx[[pocet+1]] deltat;
y[[pocet+1]]=y[[pocet]]+rychlosty[[pocet+1]] deltat;
rychlost[[pocet+1]]=(Sqrt[((rychlostx[[pocet+1]])^2)+((rychlosty[[pocet+1]])^2)]);
pocet=pocet+1];

dataprograf1=Table[{x[[i]],y[[i]]},{i,1,pocet-1,1}];
dataprograf2=Table[{deltat (i-1),rychlost[[i]]},{i,1,pocet-1,1}];
grafst=ListPlot[dataprograf1,PlotStyle->Directive[barvabodu1,PointSize[velikostbodu]],AxesOrigin->{0,0},AxesLabel->{Style["(x/m)",velikostpisma,barva,FontFamily->Times],Style["(y/m)",velikostpisma,barva,FontFamily->Times]},FormatType->StandardForm]
grafvt=ListPlot[dataprograf2,PlotStyle->Directive[barvabodu2,PointSize[velikostbodu]],AxesOrigin->{0,0},AxesLabel->{Style["(t/s)",velikostpisma,barva,FontFamily->Times],Style["(v/(m.(s^(-1))))",velikostpisma,barva,FontFamily->Times]},FormatType->StandardForm]

Pomocí metody dynamického modelování zobrazte graf závislosti okamžité výchylky na čase a graf závislosti velikosti rychlosti na čase pro těleso o hmotnosti 200 g kmitající na pružině o tuhosti 100 N.(m^(-1)). V okmažiku začátku pohybu mělo těleso okamžitou výchylku 15 cm a nacházelo se v klidu.

Řešení

ClearAll["Global`*"];
g=9.81; 
m=0.200;
k=100;
y0=0.15;
v0=0;

deltat=0.01;
barva=Black;
barvabodu1=Red;
barvabodu2=Blue;
tloustka=0.007;
velikostbodu=0.02;
velikostpisma=12;

rychlost=Table[Null,{i,500}];
y=Table[Null,{i,500}];
pocet=1;
rychlost[[pocet]]=v0;
y[[pocet]]=y0;
While[pocet deltat<0.5,
F=-k y[[pocet]];
a=(F/m);
rychlost[[pocet+1]]=rychlost[[pocet]]+a deltat;
y[[pocet+1]]=y[[pocet]]+rychlost[[pocet+1]] deltat;
pocet=pocet+1];

dataprograf1=Table[{deltat (i-1),y[[i]]},{i,1,pocet-1,1}];
dataprograf2=Table[{deltat (i-1),rychlost[[i]]},{i,1,pocet-1,1}];
grafst=ListPlot[dataprograf1,PlotStyle->Directive[barvabodu1,PointSize[velikostbodu]],AxesOrigin->{0,0},AxesLabel->{Style["(t/s)",velikostpisma,barva,FontFamily->Times],Style["(y/m)",velikostpisma,barva,FontFamily->Times]},FormatType->StandardForm]
grafvt=ListPlot[dataprograf2,PlotStyle->Directive[barvabodu2,PointSize[velikostbodu]],AxesOrigin->{0,0},AxesLabel->{Style["(t/s)",velikostpisma,barva,FontFamily->Times],Style["(v/(m.(s^(-1))))",velikostpisma,barva,FontFamily->Times]},FormatType->StandardForm]

Jaroslav Reichl, SPŠST Panská, Praha, © 2011