Numerické řešení 
diferenciálních rovnic

Antonín Slavík
Matematicko-fyzikální fakulta UK

Interpolační funkce

Objekt typu InterpolatingFunction představuje funkci, která vznikla po částech polynomiální interpolací konečného počtu hodnot.
Interpolační funkci můžeme získat např. použitím Interpolation[seznam]:

hodnoty={{0,0},{1,1},{2,0},{3,1},{4,0}}

f=Interpolation[hodnoty]

Plot[f[x],{x,0,4},Epilog->{PointSize[Large],Point[hodnoty]}]

f1=Interpolation[hodnoty,InterpolationOrder->1]

Plot[f1[x],{x,0,4},Epilog->{PointSize[Large],Point[hodnoty]}]

Interpolační funkce se chová jako jakákoliv jiná funkce, lze ji např. derivovat:

Plot[f1'[x],{x,0,4},Epilog->{PointSize[Large],Point[hodnoty]}]

Table[Plot[Interpolation[hodnoty,InterpolationOrder->n][x],{x,0,4}],{n,1,5}]

(Pro interpolaci řádu n je potřeba zadat aspoň n+1 bodů.)

Použitím balíčku DifferentialEquations`InterpolatingFunctionAnatomy` lze 
z interpolační funkce získat zpět hodnoty, ze kterých je zkonstruována.

Needs["DifferentialEquations`InterpolatingFunctionAnatomy`"]

Definiční obor:

InterpolatingFunctionGrid[f]

Funkční hodnoty:

InterpolatingFunctionValuesOnGrid[f]

Clear[hodnoty,f,f1]

NDSolve a NDSolveValue

NDSolve hledá přibližné numerické řešení zadané diferenciální rovnice, výsledkem je lokální pravidlo s interpolační funkcí na pravé straně:

NDSolve[{x'[t]==t Cos[t]-t ((Cos[π t])^2) x[t],x[0]==0},x[t],{t,0,20}]

reseni=x[t]/.First[%]

Plot[reseni,{t,0,20},PlotStyle->Thick]

NDSolveValue vrací jen pravou stranu lokálního pravidla.

reseni=NDSolveValue[{x'[t]==t Cos[t]-t ((Cos[π t])^2) x[t],x[0]==0},x[t],{t,0,20}]

Plot[reseni,{t,0,20},PlotStyle->Thick]

Zkouška: Porovnání derivace nalezené interpolační funkce s pravou stranou diferenciální rovnice.

derivace=D[reseni,t]

Plot[derivace-(t Cos[t]-t ((Cos[π t])^2) reseni),{t,0,20},PlotRange->All]

Clear[reseni,derivace]

Dvě varianty použití NDSolve/NDSolveValue:

NDSolveValue[{x'[t]==t Cos[t]-t ((Cos[π t])^2) x[t],x[0]==0},x[t],{t,0,20}]

Výsledkem je hodnota funkce v bodě t, tj. výraz. Vhodné např. pro zobrazení řešení pomocí Plot.

NDSolveValue[{x'[t]==t Cos[t]-t ((Cos[π t])^2) x[t],x[0]==0},x,{t,0,20}]

Výsledkem je funkce. Vhodné např. pro počítání funkčních hodnot.

%[10]

Numerické hledání primitivní funkce

f[x_]:=Sin[Sin[x]+Cos[x]]

Integrate[f[x],x]

a=0;b=50;

Numerická integrace pomocí NIntegrate (pomalé):

F1[x_]:=NIntegrate[f[t],{t,0,x}]

Timing[Plot[F1[x],{x,a,b}]]

Převedení problému na diferenciální rovnici, numerické řešení:

F2=NDSolveValue[{g'[t]==f[t],g[a]==0},g,{t,a,b}]

Timing[Plot[F2[x],{x,a,b}]]

NDSolve: podporované typy rovnic

Obyčejné diferenciální rovnice libovolného řádu

Příklad: Logistický růstový zákon x'(t)=a x(t)(b-x(t))

NDSolveValue[{x'[t]==x[t](1-x[t]),x[0]==1/5},x[t],{t,0,10}]

Plot[%,{t,0,10}, PlotRange->All,AxesOrigin->{0,0}]

Nalezení řešení pro více počátečních hodnot:

Map[NDSolveValue[{x'[t]==x[t](1-x[t]),x[0]==#},x[t],{t,0,10}]&,Range[0,2,1/5]]

Plot[%,{t,0,10},PlotRange->All,AxesOrigin->{0,0},PlotTheme->"Marketing"]

Totéž pomocí ParametricNDSolveValue:

reseni=ParametricNDSolveValue[{x'[t]==x[t](1-x[t]),x[0]==x0},x,{t,0,10},{x0}]

reseni[0]

Plot[Evaluate[Table[reseni[x0][t],{x0,0,2,1/5}]],{t,0,10},PlotRange->All,AxesOrigin->{0,0},PlotTheme->"Marketing"]

Clear[reseni]

Soustavy obyčejných diferenciálních rovnic

Příklad: Volterrův-Lotkův model (dravec a kořist)

Clear[a,b]

eqns={x'[t]==a x[t]-b x[t] y[t],     (* kořist *)
       y'[t]==-c y[t]+d x[t] y[t],(* dravec *)
        x[0]==x0,y[0]==y0};

coeffs={a->4,b->(3/2),c->3,d->1};

NDSolveValue[eqns/.coeffs/.{x0->5,y0->3},{x[t],y[t]},{t,0,10}]

ParametricPlot[%,{t,0,10},AxesOrigin->{0,0},PlotStyle->Thick,BaseStyle->12,AxesLabel->{"kořist","dravec"}]

K vykreslení řešení soustavy dvou rovnic lze použít také StreamPlot. Po zadání vektorové funkce {f(x,y), g(x,y)} zobrazí vybraná řešení soustavy x'(t)=f(x(t),y(t)),  y'(t)=g(x(t),y(t)) vyplňující zadanou oblast.

StreamPlot[{a x-b x y,-c y+d x y}/.coeffs,{x,0,8},{y,0,8},Frame->False,Axes->True,AxesLabel->{"kořist","dravec"},BaseStyle->12]

Clear[eqns,coeffs]

Příklad: Lorenzův systém

NDSolveValue[{
x'[t]==10(y[t]-x[t]),
y'[t]==28x[t]-y[t]-z[t]x[t],
z'[t]==-8/3z[t]+x[t]y[t],
x[0]==5,y[0]==5,z[0]==5},{x[t],y[t],z[t]},{t,0,50}]

ParametricPlot3D[%,{t,0,50},PlotStyle->Thin]

Obyčejné diferenciální rovnice ve vektorovém tvaru

Příklad: Pohyb planety kolem Slunce (Newtonův gravitační a pohybový zákon)

NDSolveValue[{X''[t]==-((X[t])/((Norm[X[t]])^3)),X[0]=={0.5,0},X'[0]=={0,1}},X[t],{t,0,2 π}]

ParametricPlot[Evaluate[%],{t,0,2 π},PlotStyle->Thick]

Excentricita závisí na počáteční rychlosti:

Manipulate[sol=NDSolveValue[{X''[t]==-((X[t])/((Norm[X[t]])^3)),X[0]=={0.5,0},X'[0]=={0,v0}},X[t],{t,0,2π}];
ParametricPlot[Evaluate[sol],{t,0,2π},PlotRange->{{-1.5,1},{-1,1}},PerformanceGoal->"Quality",PlotStyle->Thick],{{v0,1},.1,(Sqrt[3])}]

Obyčejné diferenciální rovnice s okrajovými podmínkami

Příklad: Nalezení nejkratší spojnice dvou bodů na ploše (geodetické křivky)

Jednotková sféra parametrizovaná pomocí zeměpisných souřadnic:

f[u_,v_]:={Cos[u] Cos[v],Cos[v] Sin[u],Sin[v]}

ParametricPlot3D[f[u,v],{u,-π,π},{v,-(π/2),(π/2)},Boxed->False,Axes->False]

Výběr dvou bodů a soustava diferenciálních rovnic pro geodetické křivky na sféře:

ptA={-60 °,-10 °};
ptB={20°,45°};

rovnice={-2 Tan[c2[t]] c1'[t] c2'[t]+c1''[t]==0,(1/2) Sin[2 c2[t]] ((c1'[t])^2)+c2''[t]==0,
c1[0]==ptA[[1]],c2[0]==ptA[[2]],c1[1]==ptB[[1]],c2[1]==ptB[[2]]};

reseni=NDSolveValue[rovnice,{c1[t],c2[t]},{t,0,1}]

Show[{
ParametricPlot3D[f[u,v],{u,-π,π},{v,-(π/2),(π/2)},Mesh->None,PlotPoints->50],
ParametricPlot3D[{f@@reseni},{t,0,1},PlotStyle->{Thickness[.01],Red}],Graphics3D[{AbsolutePointSize[10],Point[{f@@ptA,f@@ptB}]}]},Boxed->False,Axes->False]

Clear[rovnice,reseni]

Diferenciální rovnice se zpožděním

Příklad: Člověk stojící ve sprše, z níž teče studená voda o teplotě (Subscript[T, 0])=20 °C, se pomocí kohoutku s teplou vodou snaží nastavit teplotu (Subscript[T, 1])=40 °C. Rychlost, jakou otáčí kohoutkem, je přímo úměrná rozdílu mezi požadovanou teplotou a teplotou právě vytékající vody. Voda se ohřívá v kotli a její cesta z kotle do sprchy trvá τ=5 s.

κ=0.2;(Subscript[T, 0])=20;(Subscript[T, 1])=40;τ=5;

NDSolveValue[{T'[t]==κ((Subscript[T, 1])-T[t-τ]),(* T[t] je teplota vody vytékající z kotle v čase t *)
T[t/;t≤0]==(Subscript[T, 0])},(* počáteční podmínka: teplota v intervalu [-τ,0] *)
T[t],{t,0,100}]

Plot[%,{t,0,100},PlotRange->All,AxesLabel->{"čas","teplota"},BaseStyle->12,PlotStyle->Thick]

Interaktivní verze s parametry κ, τ, (Subscript[T, 0]), (Subscript[T, 1]).

Manipulate[sol=NDSolveValue[{T'[t]==κ((Subscript[T, 1])-T[t-τ]),T[t/;t≤0]==(Subscript[T, 0])},T[t],{t,0,100}];Plot[sol,{t,0,100},PlotRange->All,AxesLabel->{"čas","teplota"},BaseStyle->12,PlotStyle->Thick],{{κ,0.2},0.01,0.3,0.01,Appearance->"Labeled"},{{(Subscript[T, 0]),25},0,30,1,Appearance->"Labeled"},{{(Subscript[T, 1]),40},30,50,1,Appearance->"Labeled"},{{τ,5},1,10,Appearance->"Labeled"}]

Clear[κ,τ,sol];
Unset[(Subscript[T, 0])];Unset[(Subscript[T, 1])];

Parciální diferenciální rovnice

Příklad (rovnice vedení tepla): Nechť u(x,t) značí teplotu homogenní tyče v bodě x v čase t. Pak vývoj teploty v čase je popsán diferenciální rovnicí

({{({{((∂u)/(∂t))(x,t)=c ((∂(u^2))/(∂(x^2)))}})(x,t),}})

kde c>0 je konstanta závislá na materiálu.

Uvažujme např. tyč jednotkové délky, jejíž počáteční teplota je nulová, v levém krajním bodě udržujeme nulovou teplotu a v pravém krajním bodě teplota periodicky osciluje.

temp[x_,t_]=NDSolveValue[{D[u[x,t],t]==D[u[x,t],{x,2}],u[x,0]==0,u[0,t]==0,u[1,t]==Sin[t]},u[x,t],{x,0,1},{t,0,20}]

Manipulate[Plot[temp[x,t],{x,0,1},PlotRange->{-1,1}],{t,0,20}]

Plot3D[temp[x,t],{x,0,1},{t,0,20}, PlotRange->All,AxesLabel->Automatic]

Jak funguje NDSolve?

Hledáme řešení obyčejné diferenciální rovnice x'(t)=f(t,x(t)) na intervalu [a,b] s počáteční podmínkou x(a)=(Subscript[x, 0]). 
NDSolve postupně konstruuje body a=(Subscript[t, 0])<(Subscript[t, 1])<(Subscript[t, 2])<⋯<(Subscript[t, n])=b a funkční hodnoty (Subscript[x, 0]),(Subscript[x, 1]),(Subscript[x, 2]),…, (Subscript[x, n]), které by měly přibližně ležet na grafu hledané funkce. 
Výsledkem NDSolve je interpolační funkce zkonstruovaná z bodů ((Subscript[t, i]),(Subscript[x, i])).
Nalezení každé další dvojice ((Subscript[t, i]),(Subscript[x, i])) se nazývá krok. Počet kroků je omezen volbou MaxSteps.

ifun=NDSolveValue[{x'[t]==t Cos[t]-t ((Cos[π t])^2) x[t],x[0]==0},x,{t,0,500}]

ifun=NDSolveValue[{x'[t]==t Cos[t]-t ((Cos[π t])^2) x[t],x[0]==0},x,{t,0,500},MaxSteps->Infinity]

Plot[ifun[t],{t,0,100},PlotRange->All]

Číslo Δ(Subscript[t, i])=(Subscript[t, i])-(Subscript[t, i-1])se nazývá délkou i-tého kroku. Každý krok má obvykle jinou délku (např. kratší kroky v místech, kde funkce osciluje).
Hodnoty ((Subscript[t, i]),(Subscript[x, i])) můžeme z interpolační funkce získat použitím balíčku DifferentialEquations`InterpolatingFunctionAnatomy`.

Needs["DifferentialEquations`InterpolatingFunctionAnatomy`"]

tHodnoty=InterpolatingFunctionGrid[ifun]

Length[%]

xHodnoty=InterpolatingFunctionValuesOnGrid[ifun]

Plot[ifun[t],{t,3,4.5},Epilog->{PointSize[Medium],Point[Transpose[{Flatten[tHodnoty],xHodnoty}]]}]

Jinou možností je použití StepMonitor, Sow a Reap:

hodnoty=Reap[NDSolve[{x'[t]==t Cos[t]-t ((Cos[π t])^2) x[t],x[0]==0},x[t],{t,0,100},MaxSteps->Infinity,StepMonitor:>Sow[{t,x[t]}]]][[2,1]];

Plot[ifun[t],{t,3,4.5},Epilog->{PointSize[Medium],Point[hodnoty]}]

Clear[reseni,ifun,tHodnoty,xHodnoty,hodnoty]

Jak poznat nesprávné řešení?

NDSolve volí délky kroků Δ(Subscript[t, i]) tak, aby chyba v každém kroku vyhovovala kritériím AccuracyGoal (absolutní chyba) a PrecisionGoal (relativní chyba – viz dokumentaci).
Výchozí nastavení je voleno tak, že maximální přípustná absolutní i relativní chyba v jednom kroku je polovina WorkingPrecision (tj. při počítání s přesností MachinePrecision je přípustná chyba přibližně (10^(-8))).

Některé rovnice jsou velmi citlivé na počáteční podmínky:
Vzdálenost mezi řešeními s různými počátečními podmínkami může s časem růst až exponenciálně (viz např. rovnici x'(t)=x(t), x(0)=(Subscript[x, 0])).
Vyplatí se proto otestovat, zda změnou počátečních podmínek v řádu (10^(-8)) nedostaneme zcela jiné řešení.

Příklad: Duffingova rovnice s nulovou pravou stranou

duffing=(x''[t]-0.5 x[t]+0.5 ((x[t])^3)+0.015x'[t]==0);

sol1=NDSolveValue[{duffing,x[0]==1,x'[0]==0.54},x[t],{t,0,100}]

Plot[sol1,{t,0,100},PlotStyle->{Blue,Thickness[0.010]}]

ϵ=(10^(-8));

sol2=NDSolveValue[{duffing,x[0]==1+ϵ,x'[0]==0.54+ϵ},x[t],{t,0,100}]

Plot[{sol1,sol2},{t,0,100},PlotStyle->{{Blue,Thickness[0.015]},{Orange,Thickness[0.005]}},ImageSize->Large]

{sol1,sol2}/.t->100

%//FullForm

Příklad: Duffingova rovnice s jinými koeficienty a oscilující pravou stranou

duffing=(x''[t]-x[t]+x[t]^3+0.05x'[t]==0.3Cos[t]);

{sol1,sol2}=Map[NDSolveValue[{duffing,x[0]==#[[1]],x'[0]==#[[2]]},x[t],{t,0,100}]&,{{-1,1},{-1+ϵ,1+ϵ}}]

Plot[{sol1,sol2},{t,0,100},PlotStyle->{{Blue,Thickness[0.015]},{Orange,Thickness[0.005]}},ImageSize->Large]

Řešení se shodují přibližně do t=60. Dále se rozcházejí, ani jedno z nich pravděpodobně nebude správné.

Zkusme zvýšit WorkingPrecision na 30. Tím se AccuracyGoal a PrecisionGoal automaticky zvýší na 15.

V rovnici se nesmějí vyskytovat desetinná čísla – ta mají přesnost MachinePrecision a zvýšení WorkingPrecision by pak nemělo žádný účinek (NDSolve vypíše varování).

duffing=(x''[t]-x[t]+x[t]^3+5/100x'[t]==3/10Cos[t]);

ϵ=10^(-15);

{sol1,sol2}=Map[NDSolveValue[{duffing,x[0]==#[[1]],x'[0]==#[[2]]},x[t],{t,0,100},WorkingPrecision->30]&,{{-1,1},{-1+ϵ,1+ϵ}}]

Plot[{sol1,sol2},{t,0,100},PlotStyle->{{Blue,Thickness[0.015]},{Orange,Thickness[0.005]}},ImageSize->Large]

Clear[duffing,ϵ,sol1,sol2]

Detekce událostí

WhenEvent[událost,akce] v průběhu řešení rovnice sleduje, zda nedošlo k výskytu události; pokud ano, provede předepsanou akci.

Příklad: Volný pád předmětu z výšky h.
Chceme najít řešení pouze do okamžiku, kdy předmět dopadne na zem (y=0).

g=9.81;h=5;
NDSolveValue[{y''[t]==-g,y[0]==h,y'[0]==0,WhenEvent[y[t]==0,tmax=t;"StopIntegration"]},y[t],{t,0,20}]

Plot[%,{t,0,tmax}]

Místo “StopIntegration” (přerušení řešení) lze zadat libovolnou jinou akci, která se provede při výskytu události.
Příklad: Odraz předmětu od země, snížení rychlosti na 95 %. Vypisujeme všechny okamžiky, kdy došlo k odrazu.

NDSolveValue[{y''[t]==-g,y[0]==5,y'[0]==0,WhenEvent[y[t]==0,Print[t];y'[t]->-0.95y'[t]]},y[t],{t,0,20}]

Plot[%,{t,0,20}]

Příklad: Cesta na Měsíc
Francis Godwin, The Man in the Moone (1638):
Astronaut je přivázán k hejnu labutí, které jej dopraví na Měsíc. Labutě se vznesou ze Země a míří stále k Měsíci. Jak dlouho jim potrvá cesta k cíli, jestliže se pohybují konstantní rychlostí?

Lunární jednotka: 1 LU = 384 000 km (vzdálenost Země od Měsíce)

M[t_?NumericQ]:={Cos[((2 π t)/p)],Sin[((2 π t)/p)]}; (* poloha Měsíce měřená v LU (Země v počátku soustavy souřadnic) *)
p=27.3; (* perioda oběhu Měsíce *)
Rz=0.0167; (* poloměr Země v LU *)
Rm=0.0045; (* poloměr Měsíce v LU *)
v=0.25; (* rychlost v LU za den *)

Hledáme křivku t↦X(t), po které se pohybují labutě s astronautem. Úloha vede na vektorovou diferenciální rovnici

X'(t)=v·((M(t)-X(t))/(‖M(t)-X(t)‖)), X(0)=((Subscript[R, Z]),0).

Pravá strana rovnice není definována pro X(t)=M(t). Řešení diferenciální rovnice zastavíme v okamžiku, kdy labutě dorazí na Měsíc, tj. pro ||M(t)-X(t)||=(Subscript[R, M]).

sol=NDSolveValue[{X'[t]==((v (M[t]-X[t]))/(Norm[M[t]-X[t]])),X[0]=={Rz,0},WhenEvent[Norm[M[t]-X[t]]==Rm,tmaxfound=t;"StopIntegration"]},X[t],{t,0,100}]

Manipulate[ParametricPlot[Evaluate[sol],{t,0,tmax},Epilog->{Dashed,Circle[],Blue,PointSize[Large],Point[M[tmax]],Red,Point[sol/.t->tmax]},PlotRange->{{-1,1},{-1,1}},PlotRangePadding->0.1,Axes->False],{tmax,0.1,tmaxfound}]

Vylepšená verze z Wolfram Demonstrations:

http://demonstrations.wolfram.com/FlyingToTheMoon/
Autor: Stan Wagon, Macalester College

Module[{t,v, tmax, dirQ,p=27.3,Rz=0.0167,Rm=0.0045},
Manipulate[Module[{ tmaxfound=∞},
M[t_?NumericQ]:={Cos[((2 π t)/p)],Sin[((2 π t)/p)]};
 p=27.3;
Show[sol=NDSolveValue[{X'[t]==((v (M[t]-X[t]))/(Norm[M[t]-X[t]])),X[0]=={Rz,0},WhenEvent[Norm[M[t]-X[t]]==Rm,tmaxfound=t;"StopIntegration"]},X[t],{t,0,tmax}];
LastPoint=sol /. t ->If[NumericQ[ tmaxfound], tmaxfound, tmax];
newtmax=Min[tmaxfound, tmax];
ParametricPlot[Evaluate[sol], {t, 0, newtmax}, PlotStyle -> Thickness[0.004]],
Graphics[{Thickness[0.004],Circle[],
{Blue,Disk[{0,0},Rz]},
{PointSize[0.03],{Darker[Green], Point[M[newtmax]]}, PointSize[0.02],Point[LastPoint]},
If[dirQ && Norm[LastPoint-M[newtmax]]> 0.015,{{Thickness[0.001], Line[{LastPoint,M[newtmax]}]},
Blue,Arrow[{LastPoint,LastPoint + v Normalize[(M[newtmax]-LastPoint)]}]},{}],
Text["Moon orbit", {0.85,0.8}], Text["Earth",{0.07,-0.07}], Text["Moon",M[newtmax]1.12]}],
Frame->True,FrameTicks->{{{-1,0,1}, {}},{{-1,0,1}, {}}}, BaseStyle->12,Axes->False, ImageSize->400, PlotRange-> 1.25,ImagePadding->15]], 
{{v, 0.15, "flying speed of swans [LU/day]"},0.01,0.3, Appearance->"Labeled"}, 
{{tmax, 0.1, "flight time [days]"}, 0.1, 100, Appearance->"Labeled"},
{{dirQ, True, "swan velocity vector"}, {True, False}}, TrackedSymbols->Manipulate]]

Cvičení

Uvažujme dvojité kyvadlo tvořené dvěma deskami ve tvaru čtverců o straně délky L. Jeden roh první desky je uchycen v počátku soustavy souřadnic, v sousedním rohu první desky je připevněn roh druhé desky. Obě desky se mohou otáčet kolem bodů, ve kterých jsou připevněny. Kyvadlo má tedy dva stupně volnosti a jeho poloha je popsána dvojicí úhlů (Subscript[θ, 1]), (Subscript[θ, 2]), které odpovídají odchylkám úhlopříček čtverců od svislého směru tak, jak ukazuje následující obrázek.

Pokud se kyvadlo pohybuje, jsou úhly (Subscript[θ, 1]), (Subscript[θ, 2]) funkcemi času t, které splňují následující soustavu diferenciálních rovnic:

Přitom g=9.81 je gravitační zrychlení, ϕ=(π/4)-(Subscript[θ, 1])+(Subscript[θ, 2]) a tečky značí derivaci vzhledem k proměnné t.

Vyřešte následující úlohy:

Naprogramujte funkci, která pro zadané hodnoty L, (Subscript[θ, 1]), (Subscript[θ, 2]) vykreslí obrázek kyvadla. (Otočené čtverce lze získat pomocí Rotate.)

Naprogramujte funkci, která pro zadané hodnoty L, (Subscript[x, 1]), (Subscript[x, 2]), (Subscript[v, 1]), (Subscript[v, 2]) najde numerické řešení výše uvedené soustavy rovnic pro t∈[0, 60] s počátečními podmínkami (Subscript[θ, 1])(0)=(Subscript[x, 1]), (Subscript[θ, 2])(0)=(Subscript[x, 2]), (Subscript[θ, 1]^('))(0)=(Subscript[v, 1]), (Subscript[θ, 2]^('))(0)=(Subscript[v, 2]).

Pomocí Animate nebo Manipulate vyrobte animaci znázorňující pohyb kyvadla s počátečními podmínkami (Subscript[θ, 1])(0)=π, (Subscript[θ, 2])(0)=π, (Subscript[θ, 1]^('))(0)=0, (Subscript[θ, 2]^('))(0)=0, přičemž L=15.

Pomocí Animate nebo Manipulate vyrobte animaci znázorňující pohyb čtyř kyvadel s počátečními podmínkami (Subscript[θ, 1])(0)=(π/2)+i/10, (Subscript[θ, 2])(0)=0, (Subscript[θ, 1]^('))(0)=0, (Subscript[θ, 2]^('))(0)=0, kde i∈{1,2,3,4}, přičemž L=15. (K zobrazení čtyř obrázků vedle sebe lze použít Row nebo GraphicsRow, první varianta je rychlejší).