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.7.0
Commit 3bf9d17731 (2021-11-30 12:12 UTC)
Platform Info:
  OS: Linux (x86_64-pc-linux-gnu)
  CPU: Intel(R) Core(TM) i5-5200U CPU @ 2.20GHz
  WORD_SIZE: 64
  LIBM: libopenlibm
  LLVM: libLLVM-12.0.1 (ORCJIT, broadwell)

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 algoritmu, 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.

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,))

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

solution = solve(problem)

typeof(solution)

OrdinaryDiffEq.ODECompositeSolution{Float64, 2, Vector{Vector{Float64}}, Nothing, Nothing, Vector{Float64}, Vector{Vector{Vector{Float64}}}, ODEProblem{Vector{Float64}, Tuple{Float64, Float64}, true, NamedTuple{(:ω,), Tuple{Int64}}, ODEFunction{true, typeof(hoRHS!), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, CompositeAlgorithm{Tuple{Tsit5{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}, Rosenbrock23{2, false, DefaultLinSolve, Val{:forward}}}, OrdinaryDiffEq.AutoSwitchCache{Tsit5{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}, Rosenbrock23{0, false, DefaultLinSolve, Val{:forward}}, Rational{Int64}, Int64}}, OrdinaryDiffEq.CompositeInterpolationData{ODEFunction{true, typeof(hoRHS!), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Vector{Vector{Float64}}, Vector{Float64}, Vector{Vector{Vector{Float64}}}, OrdinaryDiffEq.CompositeCache{Tuple{OrdinaryDiffEq.Tsit5Cache{Vector{Float64}, Vector{Float64}, Vector{Float64}, OrdinaryDiffEq.Tsit5ConstantCache{Float64, Float64}, typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}, OrdinaryDiffEq.Rosenbrock23Cache{Vector{Float64}, Vector{Float64}, Vector{Float64}, Matrix{Float64}, Matrix{Float64}, OrdinaryDiffEq.Rosenbrock23Tableau{Float64}, SciMLBase.TimeGradientWrapper{ODEFunction{true, typeof(hoRHS!), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Vector{Float64}, NamedTuple{(:ω,), Tuple{Int64}}}, SciMLBase.UJacobianWrapper{ODEFunction{true, typeof(hoRHS!), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Float64, NamedTuple{(:ω,), Tuple{Int64}}}, DefaultLinSolve, FiniteDiff.JacobianCache{Vector{Float64}, Vector{Float64}, Vector{Float64}, UnitRange{Int64}, Nothing, Val{:forward}(), Float64}, FiniteDiff.GradientCache{Nothing, Vector{Float64}, Vector{Float64}, Float64, Val{:forward}(), Float64, Val{true}()}}}, OrdinaryDiffEq.AutoSwitchCache{Tsit5{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}, Rosenbrock23{0, false, DefaultLinSolve, Val{:forward}}, Rational{Int64}, Int64}}}, DiffEqBase.DEStats}

solution(0.4)

2-element Vector{Float64}:
 -0.8345222876475872
 -3.8773158506497607

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")

┌ Warning: `vendor()` is deprecated, use `BLAS.get_config()` and inspect the output instead
│   caller = npyinitialize() at numpy.jl:67
└ @ PyCall /home/kalvin/.julia/packages/PyCall/3fwVL/src/numpy.jl:67

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

Ř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 0x7fcb6f543df0>

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 0x7fcb6f4b0b80>

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 0x7fcbdc73f610>

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)

sol = solve(problem)

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.395421249133918
 26.876965062806708

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 0x7fcb6f507e20>

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 0x7fcbdc5dd9d0>

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 = 4.0, g = 9.81049)
tspan  = (0.0, 100.0)

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

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)

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.318538   0.463849  …   3.51137     3.49966     3.48751
  0.0   0.175343   0.341559   0.499721      8.15516     8.21562     8.27578
  0.0   0.184108   0.355568   0.515669     -0.0111358  -0.0895858  -0.168669
 16.0  14.9982    14.1075    13.3098       -1.0831     -1.12472    -1.16543
 17.0  16.0836    15.2772    14.5631        5.7154      5.68653     5.65781
 18.0  16.7405    15.613     14.5961    …  -7.36629    -7.42685    -7.48566

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

using PyPlot

fig = plt.figure()
ax  = Axes3D(fig)

plot(data[1,:], data[2,:], zeros(length(data[3,:])), "gray")
plot(data[1,:], zeros(length(data[2,:])), data[3,:], "gray")
plot(zeros(length(data[1,:])), data[2,:], data[3,:], "gray")
plot(data[1,:], data[2,:], data[3,:])

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

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 0x7fcb5eaeb940>

plot_surface(xs, ys, vs)