Demonstrace k NDSolve

Tomáš Kalvoda, KAM FIT, ZS 2013/2014

Popis problému

Úkolem je vyřešit soustavu obyčejných diferenciálních rovnic

y'=f(t,y), y(0)=(Subscript[y, 0]),

kde y=((Subscript[y, 1]),...,(Subscript[y, n])) je vektor n neznámých funkcí, y'=((Subscript[y, 1]^(')),...,(Subscript[y, n]^('))) je vektor jejich derivací a f:ℝ×(ℝ^n)->(ℝ^n) je vektorová funkce pravých stran. Nezávisle proměnná t má význam času a (Subscript[y, 0])=((Subscript[y, 10]),...,(Subscript[y, n0]))∈(ℝ^n) je vektor předepsaných počátečních podmínek v čase 0. Tento problém se z očividných důvodů také někdy nazývá “počáteční úloha”. Protože se v celém problému vyskytují funkce pouze jedné proměnné, času t, mluvíme také o tomto problému jako o soustavě obyčejných diferenciálních rovnic (ordinary differential equations, ODE).

V tomto notebooku se pokusím ukázat několik nejjednodušších numerizkých způsobů řešení zadaného problému.

Eulerova metoda

Abychom mohli aproximovat analytické (přesné) řešení původního problému je nejprve nutné diskretizovat nezávisle proměnnou t. Přibližné řešení y proto hledáme v diskrétních časových bodech (Subscript[t, k])=Δ k, k=0,1,2,..., kde Δ představuje délku časového kroku. Sofistikovanější metody mohou délku tohoto kroku pružně měnit. My se zde zaměříme na nejjednodušší metody s konstantním časovým krokem. Hledané hodnoty závisle proměnné (funkce) y v diskrétních časech (Subscript[t, k]) si označíme odpovídajícím způsobem jako (Subscript[y, k])=y((Subscript[t, k])), k=0,1,2,.... Všimněme si, že první hodnota, (Subscript[y, 0])=y(0) je známa, je to předepsaná počáteční podmínka.

Jakým způsobem se zbavit derivace v diferenciální rovnici? Místo tečny (derivaci) přejdeme k sečně (přímka procházející dvěma sousedními body). Nahradíme proto derivaci (dopřednou) diferencí a tu pak hodnotou pravé strany v předcházejícím okamžiku,

(((Subscript[y, k+1])-(Subscript[y, k]))/Δ)≈y'((Subscript[t, k]))≈f((Subscript[t, k]),(Subscript[y, k])), k=0,1,2,....

nebo jinak zapsáno

(Subscript[y, k+1])=(Subscript[y, k])+Δ f((Subscript[t, k]),(Subscript[y, k])), k=0,1,2,...,
(Subscript[y, 0])=y(0).

Jedná se o rekurentně zadanou posloupnost, kterou můžeme libovolně dlouho iterovat. Přesnost aproximace zřejmě závisí na velikosti časového kroku Δ. Čím bude tento krok menší, tím bude aproximace lepší (aspoň na krátkém časovém intervalu). Eulerova metoda se však v praxi nepoužívá. Nedává příliš přesné výsledky, jak bude patrné z následujících příkladů.

Jako první příklad uvažme soustavu popisující jednorozměrný harmonický oscilátor (těleso na pružině),

y'=-(ω^2) x,
x'=y,
x(0)=(Subscript[x, 0]), y(0)=(Subscript[y, 0]).

Zde x představuje vzdálenost od rovnovážné polohy, y rychlost tělesa a ω>0 frekvenci kmitání.

Manipulate[
data=NestList[{First[#]+delta,Last[#]+delta{{0,1},{-(omg^2),0}}.Last[#]}&,{0,{x0,y0}},Floor[tmax/delta]];
err={First[#],(1/2)(y0^2)+(1/2)(omg^2)(x0^2)-((1/2)((#[[2,2]])^2)+(1/2)(omg^2)((#[[2,1]])^2))}&/@data;
Column[{
Show[
ListPlot[{First[#],#[[2,1]]}&/@data,PlotRange->{-5,5},ImageSize->500,AxesLabel->{t,x}],
Plot[x0 Cos[omg t]+(y0/omg)Sin[omg t],{t,0,tmax},PlotStyle->Lighter[Red]]
],
ListPlot[err,PlotStyle->Lighter[Red],ImageSize->500,AxesLabel->{t,Style["E",Italic]},PlotRange->All]
}],
Style["Eulerova metoda",14,Bold],
"",
Style["Počáteční podmínky",12],
{{x0,1,"(Subscript[x, 0])"},-3,3},
{{y0,1,"(Subscript[y, 0])"},-3,3},
Style["Parametry modelu",12],
{{omg,1,"ω"},.1,4},
Style["Parametry Eulera",12],
{{delta,.1,"Δ"},.001,1},
{{tmax,10,"(Subscript[t, max])"},1,20},
"",
Style["Body představují numerickou aproximaci,\nčervená křivka představuje přesné řešení.\n",10],
Style["Na spodním grafu je změna energie,\nod počáteční hodnoty. Energie musí\nbýt v realitě konstantní (neuvažujeme-li\ntření). Výsledek je tedy špatný.",10],
ControlPlacement->Left
]

Centrální Eulerova metoda

Centrální Eulerova metoda se od eulerovy metody liší pouze tím, kdy vyhodnocuje pravé strany. Máme

(((Subscript[y, k+1])-(Subscript[y, k]))/Δ)≈y'((Subscript[t, k]))≈f((((Subscript[t, k+1])+(Subscript[t, k]))/2),(((Subscript[y, k+1])+(Subscript[y, k]))/2)), k=0,1,2,....

Nyní ihned nezískáváme explicitní vyjádření (Subscript[y, k+1]), tato rovnice ho pouze implicitně zadává. Uvážíme-li stejný příklad jako dříve, pak řešíme

Solve[(1/delta){x[k+1]-x[k],y[k+1]-y[k]}=={{0,1},{-(omg^2),0}}.{(x[k+1]+x[k])/2,(y[k+1]+y[k])/2},{x[k+1],y[k+1]}]

Podívejme se, jaké dává toto numerické schéma výsledky.

Manipulate[
data=NestList[{First[#]+delta,
{
(((4-(delta^2)(omg^2))#[[2,1]]+4delta#[[2,2]])/(4+(delta^2)(omg^2))),
(((4-(delta^2)(omg^2))#[[2,2]]-4delta (omg^2)#[[2,1]])/(4+(delta^2)(omg^2)))
}
}&,{0,{x0,y0}},Floor[tmax/delta]];
err={First[#],(1/2)(y0^2)+(1/2)(omg^2)(x0^2)-((1/2)((#[[2,2]])^2)+(1/2)(omg^2)((#[[2,1]])^2))}&/@data;
Column[{
Show[
ListPlot[{First[#],#[[2,1]]}&/@data,PlotRange->{-5,5},ImageSize->500,AxesLabel->{t,x}],
Plot[x0 Cos[omg t]+(y0/omg)Sin[omg t],{t,0,tmax},PlotStyle->Lighter[Red]]
],
ListPlot[err,PlotStyle->Lighter[Red],ImageSize->500,AxesLabel->{t,"E"},PlotRange->All]
}],
Style["Centrální Euler",14,Bold],
"",
Style["Počáteční podmínky",12],
{{x0,1,"(Subscript[x, 0])"},-3,3},
{{y0,1,"(Subscript[y, 0])"},-3,3},
Style["Parametry modelu",12],
{{omg,1,"ω"},.1,4},
Style["Parametry Centrálního Eulera",12],
{{delta,.1,"Δ"},.001,1},
{{tmax,10,"(Subscript[t, max])"},1,20},
"",
Style["Body představují numerickou aproximaci,\nčervená křivka představuje přesné řešení.",10],
Style["Na spodním grafu je změna energie,\nod počáteční hodnoty. Energie musí\nbýt konstantní. V tomto případě to\nprakticky platí.",10],
ControlPlacement->Left
]

Runge-Kutta metoda (RK4)

Runge-Kuttove metoda (RK4) je explicitní metoda 4. řádu. Příslušné numerické schéma nebudeme podrobně odvozovat, pouze zmíníme jeho tvar:

(Subscript[k, 1])=f((Subscript[t, k]),(Subscript[y, k])),
(Subscript[k, 2])=f((Subscript[t, k])+(Δ/2),(Subscript[y, k])+(Δ/2)(Subscript[k, 1])),
(Subscript[k, 3])=f((Subscript[t, k])+(Δ/2),(Subscript[y, k])+(h/2)(Subscript[k, 2])),
(Subscript[k, 4])=f((Subscript[t, k])+Δ,(Subscript[y, k])+Δ (Subscript[k, 3])),
(Subscript[y, k+1])=(Subscript[y, k])+(1/6)Δ((Subscript[k, 1])+2(Subscript[k, 2])+2(Subscript[k, 3])+(Subscript[k, 4])).

Manipulate[
f={{0,1},{-(omg^2),0}};
rungeStep[in_]:=Module[{tk=First[in],yk=Last[in],k1,k2,k3,k4},
k1=f.yk;
k2=f.(yk+delta/2*k1);
k3=f.(yk+delta/2*k2);
k4=f.(yk+delta*k3);
Return[{tk+delta,yk+delta/6*(k1+2k2+2k3+k4)}]
];
data=NestList[rungeStep,{0,{x0,y0}},Floor[tmax/delta]];
err={First[#],(1/2)(y0^2)+(1/2)(omg^2)(x0^2)-((1/2)((#[[2,2]])^2)+(1/2)(omg^2)((#[[2,1]])^2))}&/@data;
Column[{
Show[
ListPlot[{First[#],#[[2,1]]}&/@data,PlotRange->{-5,5},ImageSize->500,AxesLabel->{t,x}],
Plot[x0 Cos[omg t]+(y0/omg)Sin[omg t],{t,0,tmax},PlotStyle->Lighter[Red]]
],
ListPlot[err,PlotStyle->Lighter[Red],ImageSize->500,AxesLabel->{t,"E"},PlotRange->All]
}],
Style["Runge-Kutta RK4",14,Bold],
"",
Style["Počáteční podmínky",12],
{{x0,1,"(Subscript[x, 0])"},-3,3},
{{y0,1,"(Subscript[y, 0])"},-3,3},
Style["Parametry modelu",12],
{{omg,1,"ω"},.1,4},
Style["Parametry RK4",12],
{{delta,.1,"Δ"},.001,1},
{{tmax,10,"(Subscript[t, max])"},1,20},
"",
Style["Body představují numerickou aproximaci,\nčervená křivka představuje přesné řešení.",10],
Style["Na spodním grafu je změna energie,\nod počáteční hodnoty. Energie musí\nbýt konstantní.",10],
ControlPlacement->Left
]

NDSolve Mathematica (defaultní chování)

Defaultní NDSolve.

Manipulate[
sol=NDSolve[{x'[t]==y[t],y'[t]==-(omg^2)x[t],x[0]==x0,y[0]==y0},{x,y},{t,0,tmax}];
Column[{
Show[
Plot[x[t]/.sol,{t,0,tmax},PlotRange->{-5,5},ImageSize->500,AxesLabel->{t,x}],
Plot[x0 Cos[omg t]+(y0/omg)Sin[omg t],{t,0,tmax},PlotStyle->Lighter[Red]]
],
Plot[(1/2)(y0^2)-(1/2)(omg^2)(x0^2)-((1/2)((y[t])^2)+(1/2)(omg^2)((x[t])^2))/.sol,{t,0,tmax},PlotStyle->Lighter[Red],ImageSize->500,AxesLabel->{t,"E"},PlotRange->All]
}],
Style["NDSolve",14,Bold],
"",
Style["Počáteční podmínky",12],
{{x0,1,"(Subscript[x, 0])"},-3,3},
{{y0,1,"(Subscript[y, 0])"},-3,3},
Style["Parametry modelu",12],
{{omg,1,"ω"},.1,4},
Style["Parametry NDSolve",12],
{{tmax,10,"(Subscript[t, max])"},1,20},
"",
Style["Modrá křivka představuje numerické řešení,\nčervená křivka představuje přesné řešení.",10],
Style["Na spodním grafu je změna energie,\nod počáteční hodnoty. Energie musí\nbýt konstantní.",10],
ControlPlacement->Left
]

Závěr

Na jednoduchém fyzikálním problému jsme porovnali výsledky napočítané pomocí Eulerovy, Centrální Eulerovy a Runge-Kutta 4 metody s výsledky Mathematica napočtených pomocí NDSolve. Celkem lze říci, že Eulerova metoda byla vyloženě špatná. (Defaultní) NDSolve a Runge-Kutta 4 dávají dobré výsledky pro většinu voleb parametrů, ale pro některé rozsahy už nestačí a chyby v předpovědi vývoje energie (musí být konstantní - Zákon zachování energie) jsou podstatné. Nejlepší výsledky v tomto případě dává Centrální Eulerova metoda. To není náhoda, jedná se totiž o symplektickou metodu, která zachovává některé fundamentální vlastnosti původního klasického mechanického systému. Poznamenejme, že příkaz NDSolve je velmi mocný a umožňuje modifikovat používanou metodu podle přání uživatele (vizte rozsáhlou dokumentaci této funkce v Mathematica).

Raketa v homogenní gravitačním poli pod vlivem větru a odporu vzduchu se spotřebou paliva

Needs["DifferentialEquations`InterpolatingFunctionAnatomy`"];

V této sekci si ukážeme jak zábavně popisovat raketu v homogenním gravitačním poli pod vlivem větru (fouká a započítáváme i odpor vzduchu). Tento systém je popsán soustavou Newtonových pohybových rovnic (tj. předpokládáme, že rychlost rakety se neblíží rychlosti světla a nepohybujeme se v blízkosti gravitační singularity, vystačíme tak bez speciální i obecné teorie relativity).

m(t)(({{x''}, {y''}, {z''}}))=(({{0}, {0}, {-g}}))+(({{(Subscript[w, x])}, {(Subscript[w, y])}, {0}}))-k(({{x'}, {y'}, {z'}}))+((M(t))/(v(t)))(({{x'}, {y'}, {z'}})),

kde (x,y,z) představuje polohu rakety v trojrozměrném prostoru, m(t) měnící se váhu rakety (v důsledku úbytku paliva), první člen na pravé straně představuje gravitační zrychlení, druhý člen vliv (konstatní) větru, třetí odpor vzduchu s konstantou k a čtvrtý zrychlení v důsledku motorů (funkce M(t) je kladná a nenulová pouze po dobu kdy máme palivo, v(t) je velikost okamžité rychlosti).

Experimentováním s následujícím kódem můžeme docílit mnoha zajímavých výsledků. Nejprve definujeme parametry modelu.

(* palivo stačí na ... *)
tpalivo=3;
(* ubývající hmotnost *)
m[t_]:=Which[t<tpalivo,1-.25t/tpalivo,True,.75];
(* motor jede jen do té doby kdy má palivo *)
motor[t_]:=Which[t<tpalivo,10,True,0];
(* tření vzduchu *)
k=2;
(* gravitace *)
g=2;
(* větry *)
wx=5;
wy=1;
(* počáteční rychlost *)
vx0=.2;
vy0=.2;
vz0=motor[0];

Dále zaveďme členy v pohybové rovnici.

gravitace={0,0,-g};
vitr={wx,wy,0};
odpor=-k m[t]{x'[t],y'[t],z'[t]};
pohon=motor[t](({x'[t],y'[t],z'[t]})/(Sqrt[((x'[t])^2)+((y'[t])^2)+((z'[t])^2)]));

A konečně vlastní Newtonova rovnice s počátečními podmínkami (start je z bodu (0,0,0)).

newton=m[t]{x''[t],y''[t],z''[t]}==gravitace+vitr+odpor+pohon;
start={x[0]==0,y[0]==0,z[0]==0,x'[0]==vx0,y'[0]==vy0,z'[0]==vz0};

Řešení pomocí NDSolve. Řešení přerušíme v okamžik kdy raketa spadne na zem (rovina (x,y,0)).

sol=NDSolve[
{newton,start},
{x,y,z},{t,0,100},
Method->{"EventLocator","Event"->z[t]}
];
(* čas dopadu *)
dopad=InterpolatingFunctionDomain[First[z/.sol]][[1,2]];

Interaktivní analýza letu. Pro správné funguvání je potřebné spustit první příkaz v této sekci.

Manipulate[
GraphicsGrid[{
{ParametricPlot3D[
{x[t],y[t],z[t]}/.sol,{t,0,dopad},
PlotRange->All,ImageSize->300,
AxesLabel->{"x","y","z"},
PlotStyle->{Red,Thick}
],
Plot[z[t]/.sol,{t,0,dopad},
PlotStyle->{Blue,Thick},
Epilog->{Red,PointSize[Large],Point[Flatten[{tt,z[tt]/.sol}]]},
Frame->True,
FrameLabel->{"čas","výška"},
PlotLabel->"Výškový profil",
ImageSize->300
]
},
{
Plot[(Sqrt[((x'[t])^2)+((y'[t])^2)+((z'[t])^2)])/.sol,{t,0,dopad},
PlotStyle->{Green,Thick},
Epilog->{Red,PointSize[Large],Point[Flatten[{tt,(Sqrt[((x'[tt])^2)+((y'[tt])^2)+((z'[tt])^2)])/.sol}]]},
Frame->True,
FrameLabel->{"čas","rychlost"},
PlotLabel->"Rychlostní průběh",
ImageSize->300
],
Plot[motor[t],{t,0,dopad},
PlotStyle->{Brown,Thick},
Epilog->{Red,PointSize[Large],Point[Flatten[{tt,motor[tt]}]]},
Frame->True,
FrameLabel->{"čas","motor"},
ImageSize->300]
}
}
],
{{tt,0,"Čas"},0,dopad,ContinuousAction->False}
]