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
kde a a jsou zadány a je neznámá vektorová funkce jedné proměnné . Čárka označuje derivaci po složkách. O funkci mluvíme jako o řešení (ODE) pokud splňuje předepsanou počáteční podmínku v čase a v každém čase splňuje rovnost .
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í.
DifferentialEquation.jl
1.2 Základní použití balíčku 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)
, kdedu
je vektor pro uložení pravé strany (tj. ),u
je aktuální stav,t
je čas ap
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
(objektCallback
): 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
doplněných počátečními podmínkami . Zde 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. a ):
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í .
Ve fázovém prostoru je každá trajektorie vázána na kružnici se středem v bodě . Skutečně, celková energie se v tomto modelu zachovává:
Tj. pro každé správné řešení musí být hodnota 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 počet králíků (potrava) a 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":
který musí být doplněný počátečními podmínkami . Model závisí na čtyřech kladných parametrech . 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á:
- udává přírustek počtu králíků v důsledku jejich množení,
- 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ší),
- udává přírustek počtu vlků, závisí přímo úměrně na počtu králíků (snazší lov),
- ú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 . V našem konkrétním případě
Řešením této soustavy je například . 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
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 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ě , č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 popsaný standardními souřadnicemi . V tomto prostoru mějme rovinu a dělo umístěné v bodě . Z tohoto děla v okamžiku vystřelíme projektil hmotnosti s počáteční rychlostí , , gravitace působí v záporném směru osy . Předpokládejme, že
- na tento projektil působí gravitační síla
- a dále nechť fouká vítr paralelně s rovinou a působí na projektil silou
- 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 , kterou můžeme zapsat ve vektorovém tvaru
a doplnit počátečními podmínkami
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 , kde je vektorová funkce obsahující všechny neznámé a funkce představuje pravou stranu diferenciální rovnice. Samozřejmě nesmíme zapomenout na počáteční podmínky .
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ě .
Ú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)
a počáteční podmínkou
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 ). 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 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
Celková energie pak je
zde a jsou hybnosti hvězdy. Odpovídající pohybové rovnice pak jsou
Parametr se většinou klade roven , 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 .
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)