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.9.4
Commit 8e5136fa297 (2023-11-14 08:46 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
  LIBM: libopenlibm
  LLVM: libLLVM-14.0.6 (ORCJIT, skylake)
  Threads: 2 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,))

[38;2;86;182;194mODEProblem with uType [38;2;86;182;194mVector{Float64} and tType [38;2;86;182;194mFloat64. In-place: [38;2;86;182;194mtrue
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: specialized 4th order "free" interpolation, specialized 2nd order "free" stiffness-aware interpolation
t: 31-element Vector{Float64}:
  0.0
  0.07579156867605023
  0.2113797760612301
  0.3666006994986526
  0.5339867672948962
  0.7594105570395719
  1.013151553353126
  1.292662455541687
  1.5748315053446973
  1.9114460701164906
  2.2404155934064187
  2.595891039770255
  2.9712727077929757
  ⋮
  5.74811068267924
  6.182229465097288
  6.587797087473309
  7.021620776712173
  7.436814811717523
  7.864980354860273
  8.289948074882714
  8.706551102699885
  9.141351974082045
  9.54124865230331
  9.985353129203363
 10.0
u: 31-element Vector{Vector{Float64}}:
 [2.2, -3.3]
 [1.6764622584871014, -3.5943670209150524]
 [0.6523929689893416, -3.9120816754215073]
 [-0.5738528458399711, -3.9243715211423273]
 [-1.831363667814173, -3.5179688852600397]
 [-3.1812518233849123, -2.3684665014118895]
 [-3.931346020137157, -0.5239447168696503]
 [-3.6107676327472085, 1.6408447167498865]
 [-2.1732635270984, 3.317673365379704]
 [0.36958294087007976, 3.9488678094785783]
 [2.707156949452848, 2.8985415675273356]
 [3.9428152860054717, 0.42970827978921966]
 [3.175978406811751, -2.3756967648732417]
 ⋮
 [3.9517363850603133, 0.3446646303105456]
 [2.816469362240124, -2.79351565118836]
 [-0.0860333695698528, -3.966042955400602]
 [-3.0812535437518065, -2.498796363321305]
 [-3.9232481002892596, 0.5891133227626457]
 [-2.125389270665472, 3.350034137398218]
 [1.114107877293287, 3.807857911940257]
 [3.567593197240925, 1.7361294860802405]
 [3.6281253010761803, -1.60617744630582]
 [1.3761848777997723, -3.7215325890267907]
 [-2.0199575403610166, -3.415397605994239]
 [-2.119126382372928, -3.354768642905599]

typeof(solution)

ODESolution{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, SciMLBase.AutoSpecialize, FunctionWrappersWrappers.FunctionWrappersWrapper{Tuple{FunctionWrappers.FunctionWrapper{Nothing, Tuple{Vector{Float64}, Vector{Float64}, NamedTuple{(:ω,), Tuple{Int64}}, Float64}}, FunctionWrappers.FunctionWrapper{Nothing, Tuple{Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, NamedTuple{(:ω,), Tuple{Int64}}, Float64}}, FunctionWrappers.FunctionWrapper{Nothing, Tuple{Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{Float64}, NamedTuple{(:ω,), Tuple{Int64}}, ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}}, FunctionWrappers.FunctionWrapper{Nothing, Tuple{Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, NamedTuple{(:ω,), Tuple{Int64}}, ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}}}, false}, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, CompositeAlgorithm{Tuple{Tsit5{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}, Rosenbrock23{1, false, Nothing, typeof(OrdinaryDiffEq.DEFAULT_PRECS), Val{:forward}, true, nothing}}, OrdinaryDiffEq.AutoSwitchCache{Tsit5{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}, Rosenbrock23{0, false, Nothing, typeof(OrdinaryDiffEq.DEFAULT_PRECS), Val{:forward}, true, nothing}, Rational{Int64}, Int64}}, OrdinaryDiffEq.CompositeInterpolationData{ODEFunction{true, SciMLBase.AutoSpecialize, FunctionWrappersWrappers.FunctionWrappersWrapper{Tuple{FunctionWrappers.FunctionWrapper{Nothing, Tuple{Vector{Float64}, Vector{Float64}, NamedTuple{(:ω,), Tuple{Int64}}, Float64}}, FunctionWrappers.FunctionWrapper{Nothing, Tuple{Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, NamedTuple{(:ω,), Tuple{Int64}}, Float64}}, FunctionWrappers.FunctionWrapper{Nothing, Tuple{Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{Float64}, NamedTuple{(:ω,), Tuple{Int64}}, ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}}, FunctionWrappers.FunctionWrapper{Nothing, Tuple{Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, NamedTuple{(:ω,), Tuple{Int64}}, ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}}}, false}, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, Vector{Vector{Float64}}, Vector{Float64}, Vector{Vector{Vector{Float64}}}, OrdinaryDiffEq.CompositeCache{Tuple{OrdinaryDiffEq.Tsit5Cache{Vector{Float64}, Vector{Float64}, Vector{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{true, ODEFunction{true, SciMLBase.AutoSpecialize, FunctionWrappersWrappers.FunctionWrappersWrapper{Tuple{FunctionWrappers.FunctionWrapper{Nothing, Tuple{Vector{Float64}, Vector{Float64}, NamedTuple{(:ω,), Tuple{Int64}}, Float64}}, FunctionWrappers.FunctionWrapper{Nothing, Tuple{Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, NamedTuple{(:ω,), Tuple{Int64}}, Float64}}, FunctionWrappers.FunctionWrapper{Nothing, Tuple{Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{Float64}, NamedTuple{(:ω,), Tuple{Int64}}, ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}}, FunctionWrappers.FunctionWrapper{Nothing, Tuple{Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, NamedTuple{(:ω,), Tuple{Int64}}, ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}}}, false}, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, Vector{Float64}, NamedTuple{(:ω,), Tuple{Int64}}}, SciMLBase.UJacobianWrapper{true, ODEFunction{true, SciMLBase.AutoSpecialize, FunctionWrappersWrappers.FunctionWrappersWrapper{Tuple{FunctionWrappers.FunctionWrapper{Nothing, Tuple{Vector{Float64}, Vector{Float64}, NamedTuple{(:ω,), Tuple{Int64}}, Float64}}, FunctionWrappers.FunctionWrapper{Nothing, Tuple{Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, NamedTuple{(:ω,), Tuple{Int64}}, Float64}}, FunctionWrappers.FunctionWrapper{Nothing, Tuple{Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{Float64}, NamedTuple{(:ω,), Tuple{Int64}}, ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}}, FunctionWrappers.FunctionWrapper{Nothing, Tuple{Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, NamedTuple{(:ω,), Tuple{Int64}}, ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}}}, false}, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, Float64, NamedTuple{(:ω,), Tuple{Int64}}}, LinearSolve.LinearCache{Matrix{Float64}, Vector{Float64}, Vector{Float64}, SciMLBase.NullParameters, LinearSolve.DefaultLinearSolver, LinearSolve.DefaultLinearSolverInit{LinearAlgebra.LU{Float64, Matrix{Float64}, Vector{Int64}}, LinearAlgebra.QRCompactWY{Float64, Matrix{Float64}, Matrix{Float64}}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, LinearAlgebra.LU{Float64, Matrix{Float64}, Vector{Int64}}, Tuple{LinearAlgebra.LU{Float64, Matrix{Float64}, Vector{Int64}}, Vector{Int64}}, Nothing, Nothing, Nothing, LinearAlgebra.SVD{Float64, Float64, Matrix{Float64}, Vector{Float64}}, LinearAlgebra.Cholesky{Float64, Matrix{Float64}}, LinearAlgebra.Cholesky{Float64, Matrix{Float64}}, Tuple{LinearAlgebra.LU{Float64, Matrix{Float64}, Vector{Int32}}, Base.RefValue{Int32}}, Tuple{LinearAlgebra.LU{Float64, Matrix{Float64}, Vector{Int64}}, Base.RefValue{Int64}}, LinearAlgebra.QRPivoted{Float64, Matrix{Float64}, Vector{Float64}, Vector{Int64}}}, LinearSolve.InvPreconditioner{LinearAlgebra.Diagonal{Float64, Vector{Float64}}}, LinearAlgebra.Diagonal{Float64, Vector{Float64}}, Float64, Bool}, FiniteDiff.JacobianCache{Vector{Float64}, 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}()}, Float64, Rosenbrock23{1, false, Nothing, typeof(OrdinaryDiffEq.DEFAULT_PRECS), Val{:forward}, true, nothing}, Nothing}}, OrdinaryDiffEq.AutoSwitchCache{Tsit5{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}, Rosenbrock23{0, false, Nothing, typeof(OrdinaryDiffEq.DEFAULT_PRECS), Val{:forward}, true, nothing}, Rational{Int64}, Int64}}}, SciMLBase.DEStats, Vector{Int64}}

solution(0.4)

2-element Vector{Float64}:
 -0.8345222876484746
 -3.8773158506498473

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

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

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

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

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)

[38;2;86;182;194mODEProblem with uType [38;2;86;182;194mVector{Int64} and tType [38;2;86;182;194mInt64. In-place: [38;2;86;182;194mtrue
timespan: (0, 1)
u0: 2-element Vector{Int64}:
 20
 15

sol = solve(problem)

retcode: Success
Interpolation: specialized 4th order "free" interpolation, specialized 2nd order "free" stiffness-aware interpolation
t: 43-element Vector{Float64}:
 0.0
 0.029444011763495117
 0.04493194238433268
 0.07244352789696233
 0.08878749431187169
 0.11006532388099227
 0.13069173509284515
 0.15374001793211753
 0.17778927287865376
 0.20264535597257305
 0.2251566489768519
 0.24904937448464598
 0.26941842311088804
 ⋮
 0.7318753294354083
 0.7583298738268807
 0.7893011337166145
 0.8108958033660345
 0.8331687214205961
 0.8583646038783384
 0.8835543628509601
 0.9098474323294037
 0.9346811953727427
 0.957557431519584
 0.9792728605938944
 1.0
u: 43-element Vector{Vector{Float64}}:
 [20.0, 15.0]
 [28.586579707474897, 14.085412733437183]
 [33.06652127552501, 16.963418956192186]
 [30.214118752458123, 26.42058717322761]
 [23.459884918475737, 28.039191506987397]
 [18.240826449311992, 23.009743973833178]
 [18.284694282077226, 17.166090437795987]
 [22.874818073647365, 13.778094060636139]
 [30.65681430601017, 14.953302199636772]
 [33.5468774474866, 22.370516680577154]
 [25.675469299351352, 28.16575803077007]
 [18.614743280907067, 23.895290470093673]
 [18.01197827680752, 17.923883270586877]
 ⋮
 [26.21914712450837, 13.677516015606994]
 [33.623777185606876, 18.173084488677343]
 [27.220646891120822, 27.748338815370143]
 [19.773591638538, 25.579154900186285]
 [17.859502354923713, 19.011507859322954]
 [21.56619948305448, 14.241292432323876]
 [29.35724201323928, 14.446387061344435]
 [33.729288633108126, 21.422941638059523]
 [25.755539423494593, 28.035697196428483]
 [18.879482239825144, 24.217353554013325]
 [18.102094270515966, 17.92123287524174]
 [21.50822714483569, 14.273197185421274]

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.395421182161744
 26.876964922808096

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

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

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)

[38;2;86;182;194mODEProblem with uType [38;2;86;182;194mVector{Float64} and tType [38;2;86;182;194mFloat64. In-place: [38;2;86;182;194mtrue
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: specialized 4th order "free" interpolation, specialized 2nd order "free" stiffness-aware interpolation
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.694715907465591, 0.754666451814098, 0.7667541716983192, 12.118485315197486, 13.513693316977092, 13.062224994963582]
 [1.2349052965163516, 1.3742493981277597, 1.333609845772295, 9.651615456513705, 11.43111865181711, 9.80486107954095]
 [1.6720392806208801, 1.909879676299767, 1.7616570625440786, 7.925172666176938, 10.085224634373276, 7.42526693630583]
 [2.0345769969616616, 2.3899137794742655, 2.084596851584516, 6.632456148289781, 9.172241837677344, 5.558616286118252]
 [2.339814109908635, 2.8316335793886624, 2.322867063767838, 5.613615851140362, 8.53311784459356, 4.015418306726688]
 [2.5989572859098775, 3.2462329323844106, 2.489713802490473, 4.777017438306389, 8.075512510152375, 2.6880884020548623]
 [2.8196051748773123, 3.6412291533694576, 2.59418360042679, 4.066421416239204, 7.741169429902883, 1.5120158223887001]
 [3.007090407614397, 4.021742167990888, 2.6427645357760543, 3.4456831586612746, 7.49072924318579, 0.44739356642250805]
 [3.1652701003071466, 4.391229768102873, 2.640369255521907, 2.891005873990715, 7.296194278705291, -0.530120091190858]
 [3.2970275963877627, 4.751948998697828, 2.590963645593964, 2.386707453843273, 7.137084115711886, -1.4347488283992553]
 [3.4046098717680398, 5.105276731277106, 2.4979803640036757, 1.9226775828677947, 6.998431539402986, -2.274083796822538]
 [3.4898608491461927, 5.451954511822144, 2.364584048964602, 1.4926624411138456, 6.869638906184121, -3.051621172233515]
 [3.5543809909188973, 5.792286094287444, 2.1938258209090877, 1.0929925455252831, 6.743669558361979, -3.768679594567921]
 [3.5996310640288987, 6.126297573585438, 1.988715354222696, 0.7216089825635632, 6.6163390731556895, -4.425787437399276]
 [3.6269935908271536, 6.453863083901447, 1.752235922898612, 0.3773396168342227, 6.485645378896944, -5.023585844525951]
 [3.6378037685089293, 6.774798237939459, 1.4873249366594494, 0.05939652669399107, 6.351153338918994, -5.56332119448248]
 [3.63336007168259, 7.088924793363241, 1.1968380548694775, -0.23294077100577595, 6.213463588805524, -6.047027243022824]
 [3.6149227856178285, 7.396111276823617, 0.883509767431333, -0.5004985632748605, 6.073783612208934, -6.477500785067149]
 [3.583706579944861, 7.69629479127645, 0.5499184496370302, -0.7442807400346245, 5.933603221732246, -6.858160167960367]
 [3.54087125421271, 7.989489028724744, 0.19845997247331676, -0.9654755407941197, 5.794465072454237, -7.192853823936723]
 [3.4875131655354363, 8.275782869020622, -0.16866884184775094, -1.1654255339008455, 5.657814758372536, -7.485663717512558]
 [3.4875131655354363, 8.275782869020622, -0.16866884184775094, -1.1654255339008455, 5.657814758372536, -7.485663717512558]

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

2-element Vector{PyCall.PyObject}:
 PyObject <matplotlib.lines.Line2D object at 0x7f267ae4ecd0>
 PyObject <matplotlib.lines.Line2D object at 0x7f267ae4f190>

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

plot_surface(xs, ys, vs)

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

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: specialized 4th order "free" interpolation, specialized 2nd order "free" stiffness-aware interpolation
t: 146-element Vector{Float64}:
   0.0
   0.005741209690377193
   0.06315330659414911
   0.2257443459017004
   0.4432169809223693
   0.7016479771086079
   1.0398806400661986
   1.426520416085224
   1.860900456196591
   2.3121657783112908
   2.837664627783342
   3.3748387755374267
   3.8664888778483824
   ⋮
  92.33067632295064
  93.0640025487916
  93.83122888576897
  94.56579751391101
  95.2900310978671
  96.04463424770525
  96.74565092814949
  97.46794137213709
  98.21355465976676
  98.91695191029633
  99.66539907985606
 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.3055314111342519, 0.2997652995704978, -0.07508436020296941]
 [0.06750293774400351, -0.31241714759450884, 0.29709977937029275, -0.00942732080079851]
 [0.13132040103515416, -0.30499854772502927, 0.288866733145736, 0.07675584916039081]
 [0.20384422819326714, -0.27315413962561275, 0.270622867764747, 0.1664655908221013]
 [0.2887124283964392, -0.20194998711273354, 0.22671592232199175, 0.24624665905241475]
 [0.3606799278015389, -0.09942639735566645, 0.13797284340567495, 0.27239966412641375]
 [0.3894396769209354, 0.011733085494153087, -0.013601689176277179, 0.2290349375215154]
 [0.3412610311676431, 0.09618234265482715, -0.20060497249851325, 0.14226589474873508]
 [0.18535021271130175, 0.14438078241491686, -0.377932120475736, 0.0453549631883803]
 [-0.038845070965751856, 0.1482098930340257, -0.43063790654983813, -0.028316441995897126]
 [-0.23401300564686664, 0.11800086572432625, -0.34356136077225263, -0.09723146361510786]
 ⋮
 [-0.15311766972809868, -0.3210369574190884, 0.19712555924202063, -0.19729989879661924]
 [-1.3156542873976223e-6, -0.34218422455686126, 0.21421280274148263, 0.14715758700080037]
 [0.15644340250389166, -0.12225663058485754, 0.17857189889049135, 0.3805946267632359]
 [0.2376339864037652, 0.15590201863250572, 0.017898585893331462, 0.3394283644588199]
 [0.1639877361117294, 0.34246145882034157, -0.21734128199610403, 0.1684075913287333]
 [-0.04881420628671565, 0.399568903716536, -0.2996677881720466, -0.015526836342412935]
 [-0.21029094221076505, 0.32702376294918145, -0.12994535450207084, -0.19445783073048065]
 [-0.21351621445918834, 0.12246020805192959, 0.11060946304369675, -0.35880698440001546]
 [-0.07823094376709035, -0.16045693786168208, 0.22517180738534984, -0.3579740461322317]
 [0.08421626150502108, -0.34077935487828276, 0.2276407298568833, -0.12225346645351692]
 [0.2428294858307684, -0.30586498727835754, 0.18817918339315165, 0.20162312864658288]
 [0.29930587365890166, -0.2221538884446682, 0.14536733406743907, 0.28985095559190444]

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

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