Jdi na navigaci předmětu

12: Diferenciální rovnice

Tento notebook je výukovým materiálem v předmětu BI-JUL.21 vyučovaném v zimním semestru akademického roku 2021/2022 Tomášem Kalvodou. Tvorba těchto materiálů byla podpořena NVS FIT.

Hlavní stránkou předmětu, kde jsou i další notebooky a zajímavé informace, je jeho Course Pages stránka.

versioninfo()

Julia Version 1.11.1
Commit 8f5b7ca12ad (2024-10-16 10:53 UTC)
Build Info:
  Official https://julialang.org/ release
Platform Info:
  OS: Linux (x86_64-linux-gnu)
  CPU: 8 × Intel(R) Core(TM) i5-8250U CPU @ 1.60GHz
  WORD_SIZE: 64
  LLVM: libLLVM-16.0.6 (ORCJIT, skylake)
Threads: 1 default, 0 interactive, 1 GC (on 8 virtual cores)

1. Úvod

V tomto notebooku budeme používat balíček DifferentialEquations.jl. Před použitím tohoto notebooku si ho nezapomeňte nainstalovat pomocí ] add DifferentialEquations. Instalace může chvíli trvat. Balíček z očividných důvodů patří k těm objemnějším.

Jeho dokumentace je opravdu rozsáhlá a podrobná. Zájemcům vřele doporučuji podívat se minimálně na Tutoriály.

using DifferentialEquations, PyPlot

1.1 Obyčejné diferenciální rovnice (ODE)

Diferenciálních rovnic existuje mnoho a mnoho typů a druhů. My se budeme soustředit na soustavy obyčejných diferenciálních rovnic (ODE).

Pod ODE máme na mysli problém, který lze vyjádřit ve tvaru $u' = f(u, t), \quad u(t_0) = u_0 \tag{ODE},$

kde $t_0 \in \mathbb{R}$ a $u_0 \in \mathbb{R}^n$ a $f: \mathbb{R}^n \times \mathbb{R} \to \mathbb{R}^n$ jsou zadány a $u: M \to \mathbb{R}^n$ je neznámá vektorová funkce jedné proměnné $t \in M$ . Čárka označuje derivaci po složkách. O funkci $u(t)$ mluvíme jako o řešení (ODE) pokud splňuje předepsanou počáteční podmínku v čase $t_0$ a v každém čase $t \in M$ splňuje rovnost $u'(t) = f(u(t), t)$ .

Důležitost obyčejných diferenciálních rovnic je těžké přecenit. Nacházejí využití při popisu obrovského množství přírodních jevů, ale jejich aplikace přesahují i do dalších odvětví.

Vedle obyčejných diferenciálních rovnic existují i parciální diferenciální rovnice, zpožděné diferenciální rovnice, stochastické diferenciální rovnice, algebraicko-diferenciální rovnice, integro-diferenciální rovnice,... Těmi se zde ale vůbec zabývat nebudeme. Na téma diferenciální rovnic by šlo udělat dvousemestrální kurz, který by i tak byl napnutý k prasknutí.

1.2 Základní použití balíčku `DifferentialEquation.jl`

Balíček DifferentialEquations.jl poskytuje metodu solve, která se snaží zadanou diferenciální rovnici (numericky!) řešit. K zadání obyčejné diferenciální rovnice musíme vytvořit objekt ODEProblem, jehož konstruktoru musíme předat popořadě následující:

metodu představující pravou stranu, jejíž hodnota se zapisuje do předalokovaného pole, tj. tato funkce má signaturu f(du, u, p, t), kde du je vektor pro uložení pravé strany (tj. $f(u,t)$ ), u je aktuální stav, t je čas a p jsou případné parametry modelu (viz níže).
počáteční podmínku,
časový rozsah (t_start, t_end),
případné parametry modelu (vektor, pojmenovaný tuple...),
další nepovinné keyword argumenty.

Poté tento model předáme metodě solve, která vrátí řešení (v ideálním případě). Návratová hodnota je speciální objekt solution, který obsahuje napočtená data a dokáže mezi nimi funkční hodnotu interpolovat (aproximace řešení je totiž napočtena pouze v jistých diskrétních časových krocích). Hodnotu proměnných získáme prostým funkčním voláním solution(t), kde t je požadovaný čas.

Metoda solve se snaží úlohu řešit automaticky, v některých případech bude ale vhodné předat ji metodu, která se má k řešení použít (druhý poziční argument; zde pod "metoda" máme na mysli algoritmus, ne "Julia metodu"), nebo další upřesňující keyword argumenty (viz Common Solver Options). Ze všech zmiňme aspoň tyto:

adaptive (true/false): adaptivní změna délky časového kroku.
dt (číslo): počáteční délka časového kroku.
callback (objekt Callback): možnost odchytávat různé události a reagovat na ně.

2. Příklady

Podobně jako v předchozí kapitole si ukážeme různorodé příklady modelů využívajících diferenciální rovnice a vyřešíme je pomocí DifferentialEquations.jl.

2.1 Harmonický oscilátor

Jako úvodní příklad si ukažme tzv. harmonický oscilátor. Jde o systém popsaný soustavou obyčejných diferenciálních rovnic $\begin{aligned} x' &= \omega p, \\ p' &= -\omega x, \end{aligned}$

doplněných počátečními podmínkami $(x_0, p_0)$ . Zde $\omega$ je kladný parametr hrající roli frekvence. Definujme si nejprve funkci aktualizující pravou stranu.

Poznámka: Tuto úlohu lze popsat i pomocí jedné obyčejné diferenciální rovnice druhého řádu $x'' + \omega^2 x = 0$ .

function hoRHS!(du, u, params, t)
    x, p = u
    ω    = params[:ω]
    
    du[1] =  ω * p
    du[2] = -ω * x
end

hoRHS! (generic function with 1 method)

Dále definujme ODE problém:

tspan = (0., 10.)
problem = ODEProblem(hoRHS!, [2.2, -3.3], tspan, (ω = 2,))

ODEProblem with uType Vector{Float64} and tType Float64. In-place: true
timespan: (0.0, 10.0)
u0: 2-element Vector{Float64}:
  2.2
 -3.3

Poté se můžeme pokusit získat řešení.

solution = solve(problem)

retcode: Success
Interpolation: 3rd order Hermite
t: 31-element Vector{Float64}:
  0.0
  0.07579156867605023
  0.21137977042437245
  0.3666007170571094
  0.5339867565650207
  0.7594105520505979
  1.0131515462312843
  1.2926624543863836
  1.5748314780662025
  1.9114460731464276
  2.2404156064704854
  2.5958910311136916
  2.971272748591635
  ⋮
  5.7481107700711265
  6.182229543650441
  6.587797216627079
  7.021620938026792
  7.436815033438661
  7.864980601513402
  8.289948399876323
  8.706551407009796
  9.141352308501986
  9.541248934349767
  9.985353453739242
 10.0
u: 31-element Vector{Vector{Float64}}:
 [2.2, -3.3]
 [1.6764622584871014, -3.594367020915053]
 [0.652393013093067, -3.91208166806662]
 [-0.5738529836520146, -3.9243715009903024]
 [-1.8313635923192124, -3.517968924560801]
 [-3.1812517997527254, -2.3684665331540415]
 [-3.931346012674242, -0.5239447728666089]
 [-3.6107676365382138, 1.6408447084088613]
 [-2.17326370810606, 3.317673246807072]
 [0.3695829648233682, 3.948867807248034]
 [2.7071570252195625, 2.8985414967837086]
 [3.942815278566368, 0.42970834804961067]
 [3.175978212946925, -2.375697024129057]
 ⋮
 [3.951736445399354, 0.34466393945772655]
 [2.816468923312182, -2.7935160938031838]
 [-0.08603439429780817, -3.966042933321715]
 [-3.08125435037571, -2.498795369070701]
 [-3.9232478393417587, 0.5891150631265886]
 [-2.1253876177469166, 3.350035186604909]
 [1.1141103534440433, 3.8078571881462335]
 [3.5675942548151958, 1.7361273142461158]
 [3.628124227036286, -1.6061798741384505]
 [1.3761827779102054, -3.7215333661986696]
 [-2.0199597584959306, -3.4153962950207073]
 [-2.119126383673231, -3.3547686429921315]

typeof(solution);

solution(0.5)

2-element Vector{Float64}:
 -1.5881893236551443
 -3.634234474548677

A konečně vykreslit časový průběh složek řešení (tj. $x$ a $p$ ):

ts = LinRange(tspan..., 100)

fig, ax = plt.subplots()

ax.grid()
plt.xlabel("t")
plt.ylabel("x, p")
plot(ts, solution.(ts))
ax.legend(["x", "p"], loc="upper right")

PyObject <matplotlib.legend.Legend object at 0x7b04aea6b860>

Řešení lze skutečně analyticky vyjádřit pomocí sinu a kosinu s frekvencí $\omega$ .

Ve fázovém prostoru $(x,p)$ je každá trajektorie vázána na kružnici se středem v bodě $(0,0)$ . Skutečně, celková energie se v tomto modelu zachovává: $\left( \frac{1}{2}(p^2 + x^2) \right)' = p p' + x x' = -\omega p x + \omega x p = 0.$

Tj. pro každé správné řešení musí být hodnota $\frac{1}{2}(x(t)^2 + p(t)^2)$ konstantní. Zobrazením více řešení ve fázovém prostoru se to zdá býti jako splněné:

ts = LinRange(tspan..., 200)

fig, ax = plt.subplots()

ax.grid()
ax.axis("equal")
plt.xlabel("x")
plt.ylabel("p")

p0 = 0
for x0 in LinRange(0.1, 10.1, 10)
    problem  = ODEProblem(hoRHS!, [x0, p0], tspan, (ω = 2,))
    solution = solve(problem)
    data = hcat(solution.(ts)...)
    plot(data[1,:], data[2,:])
end

Podívejme se ale na zachování energie podrobněji, v dlouhodobém měřítku. Pokud bychom měli přesné řešení, pak by následující graf měl být grafem nulové funkce:

tspan    = (0., 10_000.)
x0, p0   = [2.2, -3.3]
problem  = ODEProblem(hoRHS!, [x0, p0], tspan, (ω = 2,))
solution = solve(problem)

ts = LinRange(tspan...,1_000)
es = map(xp -> xp[1]^2 + xp[2]^2 - x0^2 - p0^2, solution.(ts))

plt.grid()
plt.xlabel("t")
plt.ylabel("chyba energie")
plot(ts, es)

1-element Vector{PyCall.PyObject}:
 PyObject <matplotlib.lines.Line2D object at 0x7b04adaa9b80>

To je velmi, velmi, velmi špatné! V tomto případě si ale můžeme trochu pomoci modifikací některých parametrů řešiče.

Pro takovéto systémy je výhodné použít implicitní řešič s pevnou délkou kroku:

tspan    = (0., 10_000.)
x0, p0   = [2.2, -3.3]
problem  = ODEProblem(hoRHS!, [x0, p0], tspan, (ω = 2,))
solution = solve(problem, ImplicitMidpoint(), adaptive=false, dt=0.1)

ts = LinRange(tspan...,1_000)
es = map(xp -> xp[1]^2 + xp[2]^2 - x0^2 - p0^2, solution.(ts))

plt.grid()
plt.xlabel("t")
plt.ylabel("chyba energie")
plot(ts, es)

1-element Vector{PyCall.PyObject}:
 PyObject <matplotlib.lines.Line2D object at 0x7b04ad72adb0>

Toto je daleko lepší.

Volba metody není triviální a vyžaduje jisté znalosti. Jak vidíte, tak výchozí volba nemusí být to nejlepší. Balíček poskytuje velmi rozsáhlou škálu metod. Například "klasickou" Runge-Kutta metodu 4. řádu s pevným časovým krokem můžeme použít takto:

tspan    = (0., 10_000.)
x0, p0   = [2.2, -3.3]
problem  = ODEProblem(hoRHS!, [x0, p0], tspan, (ω = 2,))
solution = solve(problem, RK4(), dt=0.025, adaptive=false)

ts = LinRange(tspan...,1_000)
es = map(xp -> xp[1]^2 + xp[2]^2 - x0^2 - p0^2, solution.(ts))

plt.grid()
plt.xlabel("t")
plt.ylabel("energie")
plot(ts, es)

1-element Vector{PyCall.PyObject}:
 PyObject <matplotlib.lines.Line2D object at 0x7b04ad87d220>

Pointa předchozích experimentů je zejména následující:

I když většinu ODE neumíme analyticky vyřešit, tak nám jejich analýza může dát kvalitativní vlastnosti, které by každé řešení mělo splňovat. To pro nás může být indikátorem nakolik se strojovému řešení dá věřit.
Různé metody mohou dát dost různé "výsledky"!

2.2 Lotka-Voltera

Jedním z velmi názorných a pěkných modelů využívajících obyčejné diferenciální rovnice je tzv. Lotka-Voltera model popisující populační dynamiku dvou druhů, kde jeden druh představuje predátora (např. vlci) a druhý jeho potravu (např. králíci). Jde o jeden z jednodušších modelů popisujících takovouto dynamiku. Lotka s Volterou se tomuto modelu věnovali v dvacátých letech dvacátého století, oba k němu dospěli z jiných oblastí (Lotka při popisu jistých chemických reakcí a Voltera právě při analýze populační dynamiky).

Popis modelu je poměrně jednoduchý: označme jako $x$ počet králíků (potrava) a $y$ počet vlků (predátor). Tyto hodnoty mohou nabývat reálných čísel: počty jedinců mohou být vysoké, je v nich neurčitost (z chemického pohledu jde o koncentrace jistých látek, kde celočíselnost nehraje roli). Lotka-Voltera model postuluje následující "evoluční zákon": $\begin{aligned} x' &= \alpha x - \beta x y, \\ y' &= \delta xy - \gamma y, \end{aligned}$

který musí být doplněný počátečními podmínkami $(x_0, y_0) \in \mathbb{R}_+^2$ . Model závisí na čtyřech kladných parametrech $\alpha, \beta, \gamma, \delta$ . Tyto parametry vstupují do faktorů na pravé straně rovnic a kontrolují jejich "velikost", znaménko derivací a tedy to, jestli daná proměnná roste nebo klesá:

$\alpha x$ udává přírustek počtu králíků v důsledku jejich množení,
$-\beta x y$ udává úbytek králíků v důsledku vlčího lovu (čím víc králíků nebo vlků, tím je lov snazší a úbytek větší),
$\delta xy$ udává přírustek počtu vlků, závisí přímo úměrně na počtu králíků (snazší lov),
$-\gamma y$ úbytek počtu vlků.

Definice pravé strany včetně parametrů by byla například následující (rozbalení proměnných používáme pouze pro přehlednější notaci):

function lotka_voltera!(du, u, p, t)
    x, y = u
    α, β, γ, δ = p[:α], p[:β], p[:γ], p[:δ]
    
    du[1] = α * x     - β * x * y
    du[2] = δ * x * y - γ * y
end

lotka_voltera! (generic function with 1 method)

Než se pustíme do ukázkových řešení problémů, tak je dobré si povšimnout jednoho důležitého pozorování. Z diferenciální rovnice jsme často schopni určit tzv. stacionární řešení, tedy řešení, která se nemění v čase. V kontextu Lotka-Voltera modelu by to odpovídalo situaci, kdy populace králíků a vlků jsou stabilní, systém je v "ekologické" rovnováze.

Jak takové konstantní řešení najít? Stačí, když bude nulovat pravou stranu (derivace konstant je nula). Tj. stacionární řešení jsou řešení (nelineární) soustavy rovnic $f(y, t) = 0$ . V našem konkrétním případě $\begin{aligned} \alpha x - \beta x y &= 0, \\ \delta xy - \gamma y &= 0 \end{aligned}$

Řešením této soustavy je například $(0,0)$ . Takové řešení by ale asi nebylo příliš zajímavé (na druhou stranu, je dobré, že jde o stacionární řešení; bylo by zvláštní, kdyby nám králicí vstávali z mrtvých).

Druhým řešením (po pár jednoduchých úpravách) je $x_* = \frac{\gamma}{\delta} \quad \text{a} \quad y_* = \frac{\alpha}{\beta}.$

Pojďme si nyní několik řešení sestavit a prozkoumat je.

params  = (α = 40, β = 2, γ = 50, δ = 2)
u0      = [20, 15]
tspan   = [0, 1]
problem = ODEProblem(lotka_voltera!, u0, tspan, params)

ODEProblem with uType Vector{Int64} and tType Int64. In-place: true
timespan: (0, 1)
u0: 2-element Vector{Int64}:
 20
 15

sol = solve(problem)

retcode: Success
Interpolation: 3rd order Hermite
t: 43-element Vector{Float64}:
 0.0
 0.029444011763495117
 0.0449319421065196
 0.0724435279232002
 0.08878749408015585
 0.11006532405567868
 0.13069173524598982
 0.15374001848734195
 0.17778927431123967
 0.20264535914842624
 0.225156650893443
 0.24904938064188714
 0.26941842799297877
 ⋮
 0.7318753613663668
 0.7583299034510007
 0.7893012158746888
 0.8108958927021223
 0.833168830300051
 0.8583647080362742
 0.8835544852150293
 0.909847555450319
 0.9346812545548209
 0.9575575483147452
 0.9792729670788525
 1.0
u: 43-element Vector{Vector{Float64}}:
 [20.0, 15.0]
 [28.586579707474876, 14.08541273343719]
 [33.06652121975078, 16.963418880151877]
 [30.214118742716217, 26.420587177672104]
 [23.459885008555418, 28.039191525209787]
 [18.240826433527506, 23.009743919580295]
 [18.284694300352136, 17.166090404653282]
 [22.874818231999946, 13.77809403081678]
 [30.656814745596513, 14.953302444024303]
 [33.5468769358451, 22.370517891661166]
 [25.675468479926277, 28.16575810068539]
 [18.61474239676238, 23.895288555179015]
 [18.011978667462085, 17.923882034343684]
 ⋮
 [26.219157561439133, 13.677517073641251]
 [33.623780751236794, 18.17309365153466]
 [27.220614943351986, 27.7483471786044]
 [19.773574261255188, 25.579132229740807]
 [17.85950676919653, 19.011480451673403]
 [21.566224066398394, 14.241283583443439]
 [29.357279176248586, 14.446402297021248]
 [33.72927622461636, 21.422984992746017]
 [25.755516402697193, 28.035698523248797]
 [18.879465095301814, 24.217319330600855]
 [18.10210294750075, 17.92120770065849]
 [21.508226921909987, 14.273198232542526]

Přibližnou hodnotu řešení například v čase $t=0.5$ pak získáme takto:

sol(0.5)

2-element Vector{Float64}:
 29.395420945881572
 26.87696477251207

Pojďme vizualizovat závislost populací králíků a vlků na čase:

ts = LinRange(tspan..., 1000)

fig, ax = plt.subplots()

ax.grid()
plt.xlabel("t")
plt.ylabel("populace")
plot(ts, sol.(ts))
ax.legend(["králící", "vlci"], loc="upper right")

PyObject <matplotlib.legend.Legend object at 0x7b04ad80c860>

Ačkoliv tyto křivky vypadají jako sinus/kosinus není tomu tak.

Alternativně můžeme vývoj vyzualizovat v rovině $(x,y)$ , červeným puntíkem také znázorníme stacionární řešení:

ts = LinRange(tspan..., 1000)
xy = hcat(sol.(ts)...)

2×1000 Matrix{Float64}:
 20.0  20.204   20.4155  20.6346  …  20.8011  21.0295  21.2652  21.5082
 15.0  14.8536  14.7148  14.5836     14.6068  14.488   14.3768  14.2732

fig, ax = plt.subplots()

ax.grid()
plt.xlabel("králíci")
plt.ylabel("vlci")
scatter(params[:γ] / params[:δ], params[:α] / params[:β], s=4, color="red")
plot(xy[1, :], xy[2, :])

1-element Vector{PyCall.PyObject}:
 PyObject <matplotlib.lines.Line2D object at 0x7b04ad9db3b0>

Prozkoumejme tento "fázový prostor" podrobněji, pomocí více trajektorií:

fig, ax = plt.subplots()

ax.grid()
plt.xlabel("králíci")
plt.ylabel("vlci")
scatter(params[:γ] / params[:δ], params[:α] / params[:β], s=4, color="red")

y0 = 20
tspan = [0, 0.1]
ts = LinRange(tspan..., 1000)

for x0 in LinRange(10, 24, 15)
    problem = ODEProblem(lotka_voltera!, [x0, y0], tspan, params)
    sol = solve(problem)
    xy = hcat(sol.(ts)...)
    
    plot(xy[1, :], xy[2, :])
end

Co zde vidíme? Populace opisují "bramboroid" proti směru hodinových ručiček kolem stacionárního řešení.

2.3 Balistické křivky

V úvodní přednášce jsme jako jednu z motivací pro vývoj počítačů zmínili dělostřelcovu touhu počítat balistické křivky. Pojďme se k tomuto problému vrátit, alespoň v relativně zjednodušené formě.

Uvažme trojrozměrný prostor $\mathbb{R}^3$ popsaný standardními souřadnicemi $(x,y,z)$ . V tomto prostoru mějme rovinu $z=0$ a dělo umístěné v bodě $(0,0,0)$ . Z tohoto děla v okamžiku $t=0$ vystřelíme projektil hmotnosti $m$ s počáteční rychlostí $v_0 = (v_{0x}, v_{0y}, v_{0z})$ , $v_{0z} > 0$ , gravitace působí v záporném směru osy $z$ . Předpokládejme, že

na tento projektil působí gravitační síla $F_g = (0, 0, -g)$
a dále nechť fouká vítr paralelně s rovinou $z = 0$ a působí na projektil silou $F_w = (w_x, w_y, 0)$
a konečně nechť proti pohybu projektilu působí odporová síla vzduchu, která je úměrná druhé mocnině jeho rychlosti.

Newtonovy pohybové rovnice takovéhoto projektilu jsou soustava tří obyčejných diferenciálních rovnic druhého řádu pro tři neznámé funkce $x(t), y(t), z(t)$ , kterou můžeme zapsat ve vektorovém tvaru

$m \begin{pmatrix} \ddot x \\ \ddot y \\ \ddot z \end{pmatrix} = F_g + F_w - c \sqrt{\dot x^2 + \dot y^2 + \dot z^2} \begin{pmatrix} \dot x \\ \dot y \\ \dot z \end{pmatrix} \tag{B}$

a doplnit počátečními podmínkami

$(x(0), y(0), z(0)) = (0,0,0), \quad (\dot x(0), \dot y(0), \dot z(0)) = v_0.$

Tečky zde mají význam derivací (jedna tečka -- první derivace, dvě tečky -- druhá derivace). Od dob Newtona bývá zvykem v klasické mechanice derivace takto značit.

Nejprve musíme rovnice upravit do tzv. standardní tvaru, v kterém máme na jedné straně první derivace a na druhé již pouze výrazy bez derivací. Obecně tedy do tvaru $\dot u = f(u, t)$ , kde $u: \mathbb{R} \to \mathbb{R}^n$ je vektorová funkce obsahující všechny neznámé a funkce $f: \mathbb{R}^n \times \mathbb{R} \to \mathbb{R}^n$ představuje pravou stranu diferenciální rovnice. Samozřejmě nesmíme zapomenout na počáteční podmínky $u(0) \in \mathbb{R}^n$ .

V našem případě Newtonových rovnic (B) výše to obecně není problém. Druhých derivací se snadno zbavíme zavedením pomocných funkcí, majících význam rychlosti nebo hybnosti. Zvolíme druhou možnost, abychom se nedostali do přílišné kolize ve značení. Položme prostě $(p_x, p_y, p_z) := m(\dot x, \dot y, \dot z)$ .

Úloha (B) pak přejde na soustavu česti obyčejných diferenciálních rovnic (každá obsahující pouze první derivace na levé straně) ve tvaru (první tři rovnice jsou pouze definice hybnosti)

$\begin{pmatrix} \dot x \\ \dot y \\ \dot z \\ \dot p_x \\ \dot p_y \\ \dot p_z \end{pmatrix} = \begin{pmatrix} p_x / m \\ p_y / m \\ p_z / m \\ w_x - \frac{c p_x}{m^2} \sqrt{p_x^2 + p_y^2 + p_z^2} \\ w_y - \frac{c p_y}{m^2} \sqrt{p_x^2 + p_y^2 + p_z^2} \\ -g - \frac{c p_z}{m^2} \sqrt{p_x^2 + p_y^2 + p_z^2} \end{pmatrix}$

a počáteční podmínkou

$(x(0), y(0), z(0), p_x(0), p_y(0), p_z(0)) = (0,0,0,v_{0x},v_{0y},v_{0z}).$

V Julia tedy definujeme pravou stranu takto:

function ballistic!(du, u, p, t)
    # pro zpřehlednění zápisu
    x, y, z, px, py, pz = u
    m, c, g, wx, wy     = p[:m], p[:c], p[:g], p[:wx], p[:wy]
    p                   = sqrt(px^2 + py^2 + pz^2)
    
    # in-place aktualizace pravé strany
    du[1] = px / m
    du[2] = py / m
    du[3] = pz / m
    du[4] = wx - c * px * p / m^2
    du[5] = wy - c * py * p / m^2
    du[6] = -g - c * pz * p / m^2
end

ballistic! (generic function with 1 method)

Počáteční podmínky, parametry, časový rozsah a definice problému následuje níže.

u0     = [.0, .0, .0, 16., 17., 18.]
params = (m = 2.0, c = 0.4, wx = -3, wy = 5.0, g = 9.81049)
tspan  = (0.0, 100.0)

problem = ODEProblem(ballistic!, u0, tspan, params)

ODEProblem with uType Vector{Float64} and tType Float64. In-place: true
timespan: (0.0, 100.0)
u0: 6-element Vector{Float64}:
  0.0
  0.0
  0.0
 16.0
 17.0
 18.0

Navíc použijeme tzv. callback (Event Handling and Callback Functions), pomocí kterého detekujeme okamžik dopadu projektilu (průsečík jeho trajektorie s rovinou $z=0$ ). V ten okamžik výpočet přerušíme.

condition(u, t, integrator) = u[3] <= 0       # podmínka dopadu
affect!(integrator) = terminate!(integrator)  # akce při dopadu
cb = DiscreteCallback(condition, affect!)

sol = solve(problem, callback = cb, adaptive = false, dt = 0.1)

retcode: Terminated
Interpolation: 3rd order Hermite
t: 23-element Vector{Float64}:
 0.0
 0.1
 0.2
 0.30000000000000004
 0.4
 0.5
 0.6
 0.7
 0.7999999999999999
 0.8999999999999999
 0.9999999999999999
 1.0999999999999999
 1.2
 1.3
 1.4000000000000001
 1.5000000000000002
 1.6000000000000003
 1.7000000000000004
 1.8000000000000005
 1.9000000000000006
 2.0000000000000004
 2.1000000000000005
 2.1000000000000005
u: 23-element Vector{Vector{Float64}}:
 [0.0, 0.0, 0.0, 16.0, 17.0, 18.0]
 [0.6946593193612562, 0.7569047317180045, 0.7666917147613571, 12.115394264557052, 13.599030201911495, 13.058839240907991]
 [1.23454841949801, 1.3824981670015686, 1.3332246671963957, 9.642415470084707, 11.584227985093289, 9.795117617419136]
 [1.6710435605618985, 1.9272797859245852, 1.7606102153131387, 7.908681638909297, 10.29701347387022, 7.408525403705943]
 [2.0325600382375306, 2.4192455025274358, 2.08254113296492, 6.608007587241312, 9.436896383679242, 5.535086381281603]
 [2.336366039731629, 2.875424890088702, 2.3194799995139697, 5.580760358429829, 8.84615663456045, 3.9858678168355235]
 [2.593650642928485, 3.306802026386341, 2.4847253338166246, 4.73551018140385, 8.432933546796793, 2.6538583503944473]
 [2.8120071129524034, 3.720697746734642, 2.5874067186099703, 4.016310452133416, 8.139080766258854, 1.4750986608636947]
 [2.996780180625896, 4.1220377575457166, 2.6341288997350345, 3.387422093931877, 7.925255652362837, 0.41045902529265543]
 [3.151861571125855, 4.5140879281301265, 2.629953183933119, 2.8255331579009693, 7.763557090155721, -0.563824569805977]
 [3.2801938778620925, 4.89892210474606, 2.5790173937381398, 2.315446644967804, 7.633762853410554, -1.4616454151830995]
 [3.3841055904243036, 5.277750912874637, 2.484934125924327, 1.8474426466867007, 7.521317971546094, -2.2906154801243668]
 [3.4655377057205614, 5.651172750001926, 2.3510376571938827, 1.4154909199318275, 7.416122311263543, -3.054616038292798]
 [3.5261941427889503, 6.019373454242738, 2.1805218359115597, 1.0159537481723113, 7.3116379898475214, -3.755638655447442]
 [3.5676362751896775, 6.382284786694065, 1.9765008094647094, 0.6466328385304402, 7.204117217000627, -4.3950464774073295]
 [3.5913369544728693, 6.739706118683464, 1.742019269932894, 0.3060934530321363, 7.091901560514065, -4.974348240345325]
 [3.598706647912436, 7.091393037156368, 1.4800340426245953, -0.006783186354844201, 6.9748027832111354, -5.495588948683319]
 [3.5911018840928883, 7.437117438364615, 1.1933834061258584, -0.2930682378514368, 6.853584049457119, -5.961468080325097]
 [3.569823749246869, 7.7767043967266005, 0.8847551654695279, -0.5539225941847681, 6.729550041132281, -6.37528664573538]
 [3.5361118813920567, 8.110051203123364, 0.5566599154350831, -0.7906461110299547, 6.604242108229249, -6.740804066661667]
 [3.4911374694222763, 8.437133537558031, 0.2114124659905975, -1.0046731554159314, 6.479225933880741, -7.062062552116155]
 [3.435997270078849, 8.758002994026224, -0.14887793979395353, -1.1975398472478236, 6.3559550338499395, -7.343215595858066]
 [3.435997270078849, 8.758002994026224, -0.14887793979395353, -1.1975398472478236, 6.3559550338499395, -7.343215595858066]

V atributu u objektu sol je vektor vektorů. Každý z vektorů představuje stav $(x,y,z,p_x,p_y,p_z)$ systému v čase, který odpovídá tomu v atributu t.

typeof(sol.u)

Vector{Vector{Float64}} (alias for Array{Array{Float64, 1}, 1})

Pro vykreslování si napočteme řešení v dostatečném počtu bodů:

ts   = LinRange(minimum(sol.t), maximum(sol.t), 100)
data = hcat(sol.(ts)...)

6×100 Matrix{Float64}:
  0.0   0.164282   0.318533   0.463833  …   3.46056      3.44849     3.436
  0.0   0.175452   0.341987   0.500661      8.62263      8.69045     8.758
  0.0   0.184107   0.355563   0.515651      0.00568303  -0.0712935  -0.148878
 16.0  14.998     14.1068    13.3083       -1.1182      -1.15832    -1.19754
 17.0  16.104     15.3166    14.6202        6.40795      6.38189     6.35596
 18.0  16.7403    15.6122    14.5945    …  -7.22857     -7.28672    -7.34322

A vykreslíme trajektorii i její projekce (bez vlivu vzduchu by šlo o řešitelný model a viděli bychom paraboly):

using PyPlot

fig = figure()
ax  = Axes3D(fig)

plot3D(data[1,:], data[2,:], data[3,:])

1-element Vector{PyCall.PyObject}:
 PyObject <mpl_toolkits.mplot3d.art3d.Line3D object at 0x7b049b85cdd0>

Průmět do roviny $z=0$ .

fig = figure()
plot(data[1,:], data[2,:], "gray")

1-element Vector{PyCall.PyObject}:
 PyObject <matplotlib.lines.Line2D object at 0x7b049b8ea510>

2.4 Hénonův-Heilesův model

Tento model byl motivován popisem hvězdy obíhající kolem galaktického jádra (Michel Hénon a Carl Heiles, 1962). Hvězda se pohybuje v rovině v potenciálu $V(x, y) = \frac{1}{2}(x^2 + y^2) + \lambda\left( x^2 y - \frac{y^3}{3} \right).$

Celková energie pak je $E(p_x, p_y, x, y) = \frac{1}{2}(p_x^2 + p_y^2) + V(x,y),$

zde $p_x$ a $p_y$ jsou hybnosti hvězdy. Odpovídající pohybové rovnice pak jsou $\begin{aligned} \dot x &= p_x, \\ \dot y &= p_y, \\ \dot p_x &= -x - 2\lambda xy, \\ \dot p_y &= -y - \lambda(x^2 - y^2). \end{aligned}$

Parametr $\lambda$ se většinou klade roven $1$ , tím se budeme řídit i my.

hhV(x, y) = (x^2 + y^2) / 2 + x^2 * y - y^3 / 3

function hhE(u)
    x, y, px, py = u
    
    return (px^2 + py^2) / 2 + hhV(x, y)
end

function hhRHS!(du, u, t, p)
    x, y, px, py = u
    
    du[1] = px
    du[2] = py
    du[3] = -x - 2 * x * y
    du[4] = -y - x^2 + y^2
end

hhRHS! (generic function with 1 method)

Vizualizace potenciálu $V$ .

xs = LinRange(-5, 5, 100)
ys = LinRange(-5, 5, 100)
vs = [hhV(x, y) for x in xs, y in ys]

contourf(xs, ys, vs, 20)

PyObject <matplotlib.contour.QuadContourSet object at 0x7b04ad7bd460>

plot_surface(xs, ys, vs)

PyObject <mpl_toolkits.mplot3d.art3d.Poly3DCollection object at 0x7b049bf6c140>

Pojďme si sestrojit první pokusné řešení.

tspan = (0., 100.)
u0    = [0., -0.3, 0.3, -0.1]

@info "Energy: $(hhE(u0))"

problem = ODEProblem(hhRHS!, u0, tspan)
solution = solve(problem)

[ Info: Energy: 0.10400000000000001

retcode: Success
Interpolation: 3rd order Hermite
t: 146-element Vector{Float64}:
   0.0
   0.005741209690377193
   0.06315330659414911
   0.2257443495023194
   0.44321694069609213
   0.7016479189509297
   1.0398805960603983
   1.4265203321301576
   1.860900326724969
   2.3121656670744146
   2.837664486153681
   3.3748386375482995
   3.866488786626565
   ⋮
  92.33067641938057
  93.06400277385926
  93.83122920564945
  94.56579797801598
  95.29003173414279
  96.04463492723552
  96.74565158055898
  97.46794222601206
  98.21355525283542
  98.91695260577673
  99.66539964220061
 100.0
u: 146-element Vector{Vector{Float64}}:
 [0.0, -0.3, 0.3, -0.1]
 [0.001722359127735067, -0.30056768846788734, 0.29999802606582054, -0.0977583160509226]
 [0.01894102845687605, -0.30553141113425186, 0.2997652995704978, -0.07508436020296941]
 [0.0675029388137467, -0.3124171476284529, 0.2970997792791078, -0.009427319340873397]
 [0.1313203894151157, -0.30499855081263844, 0.2888667352059383, 0.0767558338431096]
 [0.20384421245445508, -0.2731541493068877, 0.270622873143306, 0.16646557301333664]
 [0.28871241841965306, -0.20194999794899646, 0.2267159298954665, 0.24624665203882923]
 [0.3606799162176398, -0.09942642022530164, 0.13797286766448538, 0.27239966587128167]
 [0.3894396786804943, 0.011733055839921066, -0.01360163757319022, 0.22903495866031304]
 [0.34126105348099733, 0.09618232682988556, -0.20060492723545786, 0.14226591737489122]
 [0.1853502662370312, 0.1443807759919453, -0.3779320866423944, 0.04535498555263215]
 [-0.03884501154183106, 0.1482098969430036, -0.4306379134977741, -0.028316424368282862]
 [-0.23401297430666493, 0.11800087459725607, -0.3435613871551339, -0.09723144912317]
 ⋮
 [-0.15311764978206713, -0.32103697613409427, 0.19712556467545891, -0.19729986144000675]
 [-1.266474473323274e-6, -0.34218419208994916, 0.2142128026514793, 0.1471576893500995]
 [0.15644346046949645, -0.12225650991164327, 0.17857186079597176, 0.3805946627474249]
 [0.23763399534072674, 0.15590217528301062, 0.01789844089342173, 0.33942827753980886]
 [0.16398759812216285, 0.3424615654398167, -0.21734145847080724, 0.16840743140955508]
 [-0.04881441019004617, 0.3995688930356338, -0.2996677291430254, -0.01552700032375385]
 [-0.21029102752570042, 0.3270236363311624, -0.12994512771614375, -0.19445800265305832]
 [-0.21351612042388335, 0.12245990216093643, 0.11060969062918129, -0.35880711497822565]
 [-0.0782308100533857, -0.1604571497026776, 0.22517183962066092, -0.35797393963909163]
 [0.08421642049453698, -0.3407794398336221, 0.22764071181025916, -0.12225315437385839]
 [0.24282959273601368, -0.3058648742750256, 0.18817913082326806, 0.201623319439166]
 [0.29930587487676885, -0.2221538890395912, 0.14536733442398025, 0.2898509549712682]

ts = LinRange(tspan..., 10000)
xy = hcat(solution.(ts)...)

plt.xlabel("x")
plt.ylabel("y")
plt.grid()
plot(xy[1, :], xy[2, :])

1-element Vector{PyCall.PyObject}:
 PyObject <matplotlib.lines.Line2D object at 0x7b049b34e420>

ts = LinRange(tspan..., 10000)
es = [ hhE(solution(t)) for t in ts ]

plt.xlabel("t")
plt.ylabel("energy error")
plt.grid()
plot(ts, es)

1-element Vector{PyCall.PyObject}:
 PyObject <matplotlib.lines.Line2D object at 0x7b049b7d0080>

solution = solve(problem, ImplicitMidpoint(), dt=0.025)

retcode: Success
Interpolation: 3rd order Hermite
t: 4001-element Vector{Float64}:
   0.0
   0.025
   0.05
   0.07500000000000001
   0.1
   0.125
   0.15
   0.175
   0.19999999999999998
   0.22499999999999998
   0.24999999999999997
   0.27499999999999997
   0.3
   ⋮
  99.72500000000555
  99.75000000000556
  99.77500000000556
  99.80000000000557
  99.82500000000557
  99.85000000000558
  99.87500000000558
  99.90000000000559
  99.9250000000056
  99.9500000000056
  99.97500000000561
 100.0
u: 4001-element Vector{Vector{Float64}}:
 [0.0, -0.3, 0.3, -0.1]
 [0.0074995340651211615, -0.30237753456872835, 0.29996272520969297, -0.09020276549826883]
 [0.014997220355343774, -0.30450904289611674, 0.299852178008116, -0.08031790069280084]
 [0.022491245632948263, -0.3063924905411441, 0.29966984420024323, -0.07035791090938812]
 [0.0299798305085048, -0.3080261575047128, 0.2994169458442796, -0.06033544617611021]
 [0.03746122281084378, -0.30940864154712266, 0.299094438342839, -0.050263277216680544]
 [0.04493369089554918, -0.31053886090028354, 0.2987030084335926, -0.04015427103618901]
 [0.052395516914448534, -0.31141605636533304, 0.29824307307835557, -0.030021366167770548]
 [0.0598449900685336, -0.3120397927880632, 0.2977147792484499, -0.019877547650643295]
 [0.06728039986667783, -0.3124099599063402, 0.29711800460308824, -0.00973582181151498]
 [0.07470002941242217, -0.31252677256551586, 0.29645235905645817, 0.0003908090774591075]
 [0.0821021487409796, -0.3123907702996704, 0.29571718722813567, 0.010489372190179225]
 [0.08948500822846135, -0.31200281627837867, 0.2949115717704046, 0.020546949513159808]
 ⋮
 [0.22934613549873178, -0.28970449658861486, 0.18036868172901113, 0.23066885968028544]
 [0.2338244897257552, -0.28383921940925744, 0.17789965643286212, 0.23855331466830967]
 [0.2382396461904894, -0.2777804012854154, 0.17531286074587304, 0.24615213523905066]
 [0.2425886080582018, -0.2715352554673943, 0.17260408867111726, 0.2534595302026375]
 [0.24686827408531037, -0.2651111332994666, 0.16976919349756994, 0.26047024323158097]
 [0.2510754403105187, -0.2585155107961857, 0.16680410451909658, 0.26717955703089513]
 [0.2552068021656354, -0.25175597513196335, 0.16370484389024112, 0.2735832961068912]
 [0.2592589570089489, -0.24484021107872733, 0.16046754357483992, 0.27967782815198977]
 [0.26322840708290324, -0.23777598742598763, 0.1570884623415066, 0.28546006406718627]
 [0.26711156289664906, -0.2305711434170505, 0.15356400275815923, 0.290927456647784]
 [0.27090474703282624, -0.22323357523443232, 0.14989072813601503, 0.29607799796167034]
 [0.2746041983758561, -0.21577122256843909, 0.14606537937275593, 0.3009102154516957]

ts = LinRange(tspan..., 10000)
xy = hcat(solution.(ts)...)

plt.xlabel("x")
plt.ylabel("y")
plt.grid()
plot(xy[1, :], xy[2, :])

1-element Vector{PyCall.PyObject}:
 PyObject <matplotlib.lines.Line2D object at 0x7b049b20b530>

ts = LinRange(tspan..., 10000)
es = [ hhE(solution(t)) for t in ts ]

plt.xlabel("t")
plt.ylabel("energy error")
plt.grid()
plot(ts, es)

1-element Vector{PyCall.PyObject}:
 PyObject <matplotlib.lines.Line2D object at 0x7b049b236840>

Pojďme nyní vizualizovat Poincarého řez fázového prostoru $x = 0$ .

tspan = (0., 2000.)
u0    = [0., -0.3, 0.3, -0.1]

@info "Energy: $(hhE(u0))"

condition(u, t, integrator) = u[1]
affect!(integrator) = nothing

problem  = ODEProblem(hhRHS!, u0, tspan)
cb       = ContinuousCallback(condition, affect!, save_positions=(false, true))
solution = solve(problem, callback=cb, save_everystep=false)

[ Info: Energy: 0.10400000000000001

retcode: Success
Interpolation: 1st order linear
t: 646-element Vector{Float64}:
    0.0
    3.28496165858154
    6.7625682223383325
    9.536432791592647
   12.878756450259576
   15.624500869071511
   18.34167671734951
   21.734926599116044
   24.538652925539417
   27.983754484747937
   31.318342378546163
   34.21907589265747
   37.72436267717064
    ⋮
 1967.4119343845352
 1969.6478251323192
 1972.4687440656835
 1975.5206291586728
 1977.8603525397054
 1982.6665351853248
 1986.0231416116108
 1989.3928653913526
 1994.1266658926984
 1996.4863791124942
 1999.5292196788585
 2000.0
u: 646-element Vector{Vector{Float64}}:
 [0.0, -0.3, 0.3, -0.1]
 [1.0411385642507222e-15, 0.15024287516742435, -0.43290929342605255, -0.016859687318181354]
 [-1.1843167079127948e-17, -0.2692135601739056, 0.2822490131047243, 0.20712361039673732]
 [6.327792668540799e-16, 0.2951594357500116, -0.36637110048369276, -0.06281268999865522]
 [-1.7709441537000195e-16, -0.3115992315103706, 0.21822724596725376, 0.20839354271302798]
 [3.07911183817352e-16, 0.40189334720886244, -0.3003900735552479, 0.004227643338867945]
 [-1.2252347097720328e-16, -0.2974662159906258, 0.2219431278575453, -0.2309898527689653]
 [5.819741392101887e-16, 0.281290714795524, -0.3750543560243322, 0.06256191683278484]
 [-4.710894954934781e-16, -0.2768622624250179, 0.28753527506219806, -0.18870238436278153]
 [1.4283565146415631e-15, 0.14618951493472176, -0.4357533433305143, 0.009450618436039412]
 [-5.909671756225816e-17, -0.29540491223927096, 0.2977945984333608, 0.1278418827924485]
 [2.578034124931838e-15, 0.24217949293515106, -0.3962036117819447, -0.05922591967356661]
 [-1.8865027910731523e-16, -0.26408934769320636, 0.23450890684325393, 0.2699621985228413]
 ⋮
 [-5.1199363343202496e-15, 0.1185140948560184, 0.33468932056728345, 0.46554633393887335]
 [2.2257989430679385e-14, 0.5435298909784136, -0.373437060244712, -0.1188611441093199]
 [-2.537514830176611e-14, -0.4616457399202934, 0.22905881226723754, -0.10458131622176137]
 [6.357249641061809e-14, 0.4830835713908742, -0.4140419495392252, 0.11310039963846714]
 [-4.16438404788219e-14, -0.020641926896298998, 0.3583406577649276, -0.46239959032593914]
 [5.826341685677745e-14, -0.1511300012819624, -0.5559636405862047, 0.09261626679931677]
 [-4.743128893113723e-14, -0.30526335837632196, 0.48034638144575104, -0.002340600835392743]
 [1.2597221179305677e-13, -0.15523615362913357, -0.5554005086191284, -0.08906722428873336]
 [-7.982306938534966e-14, -0.048300182829756254, 0.35923689060571684, 0.4600364902575698]
 [7.676477532866244e-15, 0.47798648530826476, -0.419772499217002, -0.10730709313381258]
 [-6.74167581162664e-15, -0.4625142987022959, 0.2343897707125048, 0.09384423221927592]
 [0.10968797293670488, -0.34915400708879146, 0.22909324974301246, 0.3704457789893359]

data = hcat(solution.u...)

4×646 Matrix{Float64}:
  0.0   1.04114e-15  -1.18432e-17   6.32779e-16  …  -6.74168e-15   0.109688
 -0.3   0.150243     -0.269214      0.295159        -0.462514     -0.349154
  0.3  -0.432909      0.282249     -0.366371         0.23439       0.229093
 -0.1  -0.0168597     0.207124     -0.0628127        0.0938442     0.370446

plt.grid()
plt.xlabel(L"$y$")
plt.ylabel(L"$p_y$")
scatter(data[2, :], data[4, :], 2)

PyObject <matplotlib.collections.PathCollection object at 0x7b049b6dc170>

using ProgressMeter

function get_points(energy, y, py; tmax=100, maxiters=1_000_000_000, solver=RK4(), dt=1e-4)
    tspan = (0., tmax)
    px    = sqrt(2(energy - y^2/2 + y^3/3) - py^2)
    u0    = [0., y, px, py]
    
    condition(u, t, integrator) = u[1]
    affect!(integrator) = nothing

    problem  = ODEProblem(hhRHS!, u0, tspan)
    cb       = ContinuousCallback(condition, affect!, save_positions=(false, true))
    solution = solve(problem, solver, callback=cb, save_everystep=false, save_start=true, maxiters=maxiters, dt=dt)
    
    return hcat(solution.u...)
end

function draw_poincare(energy; n=10, tmax=100, point_size=2, solver=RK4(), dt=1e-4)
    y = 0
    
    fig, ax = plt.subplots()
    
    ax.axis("equal")
    plt.grid()
    plt.xlabel(L"$y$")
    plt.ylabel(L"$p_y$")
    
    @showprogress "Computing..." for py in LinRange(-sqrt(2energy)*0.9, sqrt(2energy)*0.9, n)
        data = get_points(energy, y, py; tmax=tmax, solver=solver, dt=dt)
        scatter(data[2, :], data[4, :], point_size)
    end
end

draw_poincare (generic function with 1 method)

draw_poincare(0.125, tmax=1000, n=60, dt=1e-8, point_size=0.5)

Computing... 100%|███████████████████████████████████████| Time: 0:00:10

draw_poincare(0.025, tmax=1000, n=60, dt=1e-8, point_size=0.5)

Computing... 100%|███████████████████████████████████████| Time: 0:00:01

draw_poincare(0.005, tmax=1000, n=60, dt=1e-8, point_size=0.5)

Computing... 100%|███████████████████████████████████████| Time: 0:00:00

Reference

Velmi obsáhlou dokumentaci balíčku DifferentialEquations.jl naleznete zde.

O obyčejných diferenciálních rovnicích existují dále celé stohy literatury. Zajímavou přehlednou monografií jsou například Hairer, Nørsett, Wanner: Solving Ordinary Differential Equations I, II, Springer, 1987.