Obyčejné diferenciální rovnice – cvičení¶
Obrátíme se k soustavám obyčejných diferenciálních rovnic s počáteční podmínkou. Použijeme Eulerovu metodu 1. řádu přesnosti a Rungovu-Kuttovu metodu 4. řádu přesnosti pro numerické řešení dvou systémů s dostupným analytickým řešením a systému, jehož řešením je zase jednou fraktál.
Rungovy-Kuttovy metody¶
Úlohou je nalézt numerickou aproximaci funkcí, pro které je předepsána soustava \(R\) obyčejných diferenciálních rovnic prvního řádu,
neboli vektorově
doplněná počáteční podmínkou,
Numerickým řešením této tzv. počáteční úlohy na síti \(x_n,\ n=0,1,\dots,\) je pro každé \(r=1,\dots,R\) posloupnost \(y_{nr}\), která lépe či hůře aproximuje přesné řešení \(y_r(x_n)\).
Na přednášce je řeč o explicitní a implicitní Eulerově metodě 1. řádu, dvou Rungových-Kuttových metodách 2. řádu a jedné R.-K. metodě 4. řádu přesnosti. Zde si z těchto metod připravíme první a poslední. Pro explicitní Eulerovu metodu jsme si předvedli vzorec
vyjadřující posun z bodu \((x_n,\mathbf{y}_n)\) o horizontální krok \(h\) ve směru tečny k funkci \(\mathbf{y}(x)\), jejíž směrnice je podle zadání rovna \(\mathbf{f}(x,\mathbf{y}(x))\). To odpovídá nultému a prvnímu členu Taylorova rozvoje \(\mathbf{y}(x)\), proto mluvíme o 1. řádu přesnosti. Taylorův rozvoj do 4. členu vystihuje Rungova-Kuttova metoda 4. řádu (RK4) v následujícím tvaru:
Pro jeden krok z \(x_n\) do \(x_{n+1}\) vyčísluje metoda pravou stranu soustavy 4krát. Krok metody lze interpretovat jako posun ve směru získaném jako vážený aritmetický průměr 4 tečen sestrojených na intervalu \(\langle x_n,x_{n+1}\rangle\).
Oběma metodami nejdříve vyřešíme následující rovnici s počáteční podmínkou:
jejímž analytickým řešením je zjevně \(y(x)=e^{-x}\), a soustavu:
s analytickým řešením \(y_1(x)=\sin x,\ y_2(x)=\cos x\).
Typické rozhraní řešičů poskytovaných numerickými knihovnami pro počáteční úlohy zahrnuje mezi vstupními argumenty meze nezávisle proměnné \(x\), vektor \(\mathbf{y}_0\) s počáteční podmínkou v \(x_0\) a odkaz na funkci, která vyčísluje vektor pravé strany pro argumenty \(x,\mathbf{y}_0\). Pythonský řešič scipy.integrate.odeint(func,y0,t,...,tfirst) nazývá nezávisle proměnnou \(t\); konkrétně je t pole, v jehož nultém prvku řešič startuje počáteční podmínkou a jeho ostatní monotónně rostoucí nebo klesající prvky jsou uzly, v kterých řešič vyčíslí (interpoluje) řešení y. Vnitřní krokování řešiče může být odlišné, dané nastavenou přesností integrace. Řešení je v podobě matice y tvaru (len(t),len(y0)), tedy nultý index je indexem kroku a první index indexuje prvky vektoru řešení. Funkce vracející vektor pravé strany má rozhraní func(t,y), pokud je argument řešiče nastaven explicitně na tfirst=True, jinak je func(y,t).
Následuje program, který řeší naše dvě počáteční úlohy pomocí SciPy řešiče a také našich dvou řešičů, implementujících Eulerovu metodu 1. řádu a Rungeovu-Kuttovu metodu 4. řádu přesnosti:
# Python: řešení soustav obyčejných diferenciálních rovnic s počáteční podmínkou
import numpy
import scipy
from numba import njit
# Symbolické konstanty pro použití v konstrukci match-case
class cEquation: EXP=1; GON=2
class cSolver: EUL=1; RK4=2; SCI=0
# Pravé strany EXP a GON
@njit
def fexp(t,y):
return -y
@njit
def fgon(t,Y):
(y1,y2)=Y # přeznačení
return numpy.array((y2,-y1))
# Eulerova metoda (solver EUL)
@njit
def Euler(func,y0,t):
y=numpy.empty((len(t),len(y0)))
y[0,:]=y0
for n in range(len(t)-1):
h=t[n+1]-t[n]
y[n+1,:]=y[n,:]+h*func(t[n],y[n,:])
return y
# Rungova-Kuttova metoda 4. řádu (solver RK4)
@njit
def RungeKutta4(func,y0,t):
y=numpy.empty((len(t),len(y0)))
y[0,:]=y0
for n in range(len(t)-1):
h=t[n+1]-t[n]; h2=h*0.5
k1=func(t[n], y[n,:])
k2=func(t[n]+h2,y[n,:]+h2*k1)
k3=func(t[n]+h2,y[n,:]+h2*k2)
k4=func(t[n]+h, y[n,:]+h *k3)
y[n+1,:]=y[n,:]+h/6*(k1+2*k2+2*k3+k4)
return y
# Vytvoření PNG v Matplotlibu
def savefig(t,y,file):
import matplotlib.pyplot as plt
match Equation:
case cEquation.EXP:
fanal=lambda t: numpy.exp(-t)
plt.plot(t,fanal(t),t,y[:,0],'o')
case cEquation.GON:
fanal1=lambda x: numpy.sin(x)
fanal2=lambda x: numpy.cos(x)
plt.plot(t,fanal1(t),t,y[:,0],'o',t,fanal2(t),t,y[:,1],'x')
plt.savefig(file)
plt.show()
# Vytvoření DAT pro gnuplot
def savedat(t,y,file):
nmax=len(y)
with open(file,'w') as f:
for n in range(nmax):
print(f'{t[n]:7.3f}',end=' ',file=f)
print(*[f'{y[n,i]:7.3f}' for i in range(numpy.size(y,1))],file=f)
# Main
Equation=cEquation.GON # EXP nebo GON
Solver=cSolver.RK4 # EUL nebo RK4 nebo SCI
# Příprava vstupních dat pro vybranou soustavu
match Equation:
case cEquation.EXP: # soustava EXP: y'=-y; y(0)=[1]; y(x)=exp(-x)
func=fexp; tmin=0; tmax=2; y0=numpy.array([1])
nmax=20 if Solver==cSolver.EUL else 5
case cEquation.GON: # soustava GON: y'=(y[1],-y[0]); y(0)=[0,1]; y(x)=[sin(x),cos(x)]
func=fgon; tmin=0; tmax=2*numpy.pi; y0=numpy.array([0,1])
nmax=100 if Solver==cSolver.EUL else 25
# Interpolační síť
t=numpy.linspace(tmin,tmax,nmax+1)
# Řešení soustavy
match Solver:
case cSolver.EUL: y=Euler(func,y0,t)
case cSolver.RK4: y=RungeKutta4(func,y0,t)
case cSolver.SCI: y=scipy.integrate.odeint(func,y0,t,tfirst=True)
# Výstupní data a obrázek
if True: savedat(t,y,'odeint.dat')
if True: savefig(t,y,'odeint.png')
Pro exponenciálu na intervalu <0,2> jsme volili 20 kroků Eulerovy metody a 5 kroků RK4 metody, aby počet vyčíslení pravé strany soustavy (míra časové náročnosti těchto metod) byl v obou případech totožný. Na obrázku porovnávajícím analytické řešení s oběma numerickými pak vidíme, jak výpočet metodou RK4 lépe přiléhá k analytickému řešení:
To platí i pro numerický výpočet goniometrické funkce – všimněme si, že Eulerově metodě zde nestačí pro jednu periodu ani 100 kroků, aby uspokojila oko.
Použité skripty pro gnuplot:
set terminal pngc size 800,600
set output 'ode-exp.png'
plot [-.05:2.05][0:1.05] exp(-x) linewidth 4,\
'ode1E.dat' u 1:2:(.03) with circles lw 2 title 'Euler',\
'ode1RK.dat' u 1:2:(.05) with circles lw 2 title 'RK4'
set output
unset terminal
set terminal pngc size 800,600
set output 'ode-sin.png'
plot [-.25:pi*2+.25][-1.25:1.25] sin(x) linewidth 4,\
'ode2E.dat' u 1:2:(.06) with circles lw 2 title 'Euler',\
'ode2RK.dat' u 1:2:(.10) with circles lw 2 title 'RK4'
set output
unset terminal
Ekvivalent našeho programu uvedeme i v podobě skriptů pro Matlab a Octave, které pro řešení soustav obyčejných diferenciálních rovnic s počáteční podmínkou nabízejí solvery (Matlab) ode45 a (Octave) lsode:
% Matlab: řešení soustav obyčejných diferenciálních rovnic s počáteční podmínkou (funkce ode45)
% nastavení výstupní sítě a počáteční podmínky, volání solveru, vykreslení řešení
xspan=[0;2]; y0=[1]; [xvec y]=ode45(@oderhs1,xspan,y0); plot(xvec,y(:,1))
disp('Pokracovani za 5 sekund...'); pause(5)
xspan=[0;pi*2]; y0=[0;1]; [xvec y]=ode45(@oderhs2,xspan,y0); plot(xvec,y(:,1:2))
disp('Pokracovani za 5 sekund...'); pause(5); disp('Hotovo.')
% pravá strana soustavy EXP
function result = oderhs1(x,y), result=[-y(1)]; end
% pravá strana soustavy GON
function result = oderhs2(x,y), result=[y(2);-y(1)]; end
% Octave: řešení soustav obyčejných diferenciálních rovnic s počáteční podmínkou (funkce lsode)
0; % uvození skriptu obsahujícího definice funkcí
% pravá strana soustavy EXP
function result = oderhs1(y,x), result=[-y(1)]; end
% pravá strana soustavy GON
function result = oderhs2(y,x), result=[y(2);-y(1)]; end
% nastavení výstupní sítě a počáteční podmínky, volání solveru, vykreslení řešení
xvec=linspace(0,2,100); y0=[1]; y=lsode(@oderhs1,y0,xvec); plot(xvec,y(:,1))
disp('Pokracovani za 5 sekund...'); pause(5)
xvec=linspace(0,pi*2,100); y0=[0;1]; y=lsode(@oderhs2,y0,xvec); plot(xvec,y(:,1:2))
disp('Pokracovani za 5 sekund...'); pause(5); disp('Hotovo.')
Fortranská verze Eulerovy a Rungovy-Kuttovy metody:
! Fortran: řešení soustav obyčejných diferenciálních rovnic s počáteční podmínkou
module mOdeint
implicit none
integer,parameter :: wp=8
integer,parameter :: EXP=1,GON=2, EUL=1,RK4=2
interface ! blok rozhraní k užití procedurovým ukazatelem
function rhs(t,y) result (r)
import wp
real(wp),intent(in) :: t,y(:)
real(wp) r(size(y))
end function
end interface
contains
! Pravé strany EXP a GON
function fexp(t,y) result (r)
real(wp),intent(in) :: t,y(:)
real(wp) r(size(y))
r=-y
end function
function fgon(t,y) result (r)
real(wp),intent(in) :: t,y(:)
real(wp) r(size(y))
r(1)=y(2)
r(2)=-y(1)
end function
! Eulerova metoda (solver EUL)
function Euler(func,y0,t) result (y)
procedure(rhs),pointer :: func
real(wp),intent(in) :: y0(:),t(0:)
real(wp),allocatable :: y(:,:)
real(wp) h
integer n
allocate (y(0:size(t)-1,size(y0)))
y(0,:)=y0
do n=0,size(t)-2
h=t(n+1)-t(n)
y(n+1,:)=y(n,:)+h*func(t(n),y(n,:))
enddo
end function
! Rungova-Kuttova metoda 4. řádu (solver RK4)
function RungeKutta4(func,y0,t) result (y)
procedure(rhs),pointer :: func
real(wp),intent(in) :: y0(:),t(0:)
real(wp),allocatable :: y(:,:),k1(:),k2(:),k3(:),k4(:)
real(wp) h,h2
integer n
allocate (y(0:size(t)-1,size(y0)))
y(0,:)=y0
do n=0,size(t)-2
h=t(n+1)-t(n); h2=h*0.5
k1=func(t(n), y(n,:))
k2=func(t(n)+h2,y(n,:)+h2*k1)
k3=func(t(n)+h2,y(n,:)+h2*k2)
k4=func(t(n)+h, y(n,:)+h *k3)
y(n+1,:)=y(n,:)+h/6*(k1+2*k2+2*k3+k4)
enddo
end function
! Vytvoření DAT pro gnuplot
subroutine savedat(t,y,file)
real(wp),intent(in) :: t(0:),y(0:,:)
character(*),intent(in) :: file
integer f,n
open (newunit=f,file=file)
do n=lbound(t,1),ubound(t,1)
write (f,'(*(f7.3x))') t(n),y(n,:)
enddo
close (f)
end subroutine
end module
! Main
program testOdeint
use mOdeint
implicit none
real(wp),parameter :: pi=atan(1._wp)*4
real(wp),allocatable :: y0(:),t(:),y(:,:)
real(wp) tmin,tmax
integer nmax,n,Equation,Solver
procedure(rhs),pointer :: func
! Výběr soustavy a řešiče
Equation=GON ! EXP nebo GON
Solver=RK4 ! EUL nebo RK4
! Příprava vstupních dat pro vybranou soustavu
select case (Equation)
case (EXP) ! soustava EXP: y'=-y; y(0)=[1]; y(x)=exp(-x)
func=>fexp; tmin=0; tmax=2; y0=[1]
nmax=merge(20,5,Solver==EUL)
case (GON) ! soustava GON: y'=(y[1],-y[0]); y(0)=[0,1]; y(x)=[sin(x),cos(x)]
func=>fgon; tmin=0; tmax=2*pi; y0=[0,1]
nmax=merge(100,25,Solver==EUL)
end select
! Interpolační síť
t=[(tmin+(tmax-tmin)*n/nmax,n=0,nmax)]
! Řešení soustavy
select case(Solver)
case (EUL); y=Euler(func,y0,t)
case (RK4); y=RungeKutta4(func,y0,t)
end select
! Výstupní data
if (.true.) call savedat(t,y,'odeint.dat')
if (.true.) call execute_command_line('gnuplot odeint.gp')
end program
DÚ: Lorenzův atraktor¶
Připravené programy jsou schopny větších věcí než jen reprodukovat analytické funkce. Jako domácí úkol spočítejte a vizualizujte Lorenzův atraktor, populární exemplář systému popsaného deterministicky a přitom vykazujícího chaotické chování. Jde o trajektorii bodu (v 3D stavovém prostoru) popsanou třemi obyčejnými diferenciálními rovnicemi:
Počáteční podmínku lze zadat libovolnou nenulovou. Hledaná trajektorie je atraktorem (podivným atraktorem), takže si k ní bod najde cestu odkudkoliv už sám. Je třeba vhodně volit interval integrace a počet kroků, aby se výsledný obrázek podobal obrázku na Wikipedii nebo třeba tady. Kreslit funkcí plot Matplotlibu lze časový vývoj jednotlivých souřadnic nebo průměty atraktoru do rovin xy, xz a yz. Pokud bychom chtěli 3D atraktor vidět prostorově, hodilo by se do funkce savefig doplnit následující větev cEquation.LOR s vytvořením 3D projekce. Pokud bychom chtěli atraktor probarvit, použila by se místo funkce plot funkce scatter (nápověda Matplotlibu) s polem hodnot pro barvení v argumentu c (volíme c=t) a barevnou škálou (nabídka Matplotlibu) v argumentu cmap:
# Python: větev příkazu match ve funkci savefig pro 3D vizualizaci pole y tvaru (len(t),3)
case cEquation.LOR:
fig=plt.figure(layout='tight')
ax=plt.axes(projection='3d')
# ax.plot(*y.T) # totéž, co ax.plot(y[:,0],y[:,1],y[:,2])
ax.scatter(*y.T,s=10,c=t,cmap='Spectral')
Nezapomeňte, že 3D Lissajousovy obrazce jsme kreslili v ParaView – to by i zde vytvořilo silný dojem, silnější než gnuplot nebo Matplotlib. Pro vyšší komfort následují pythonská funkce a fortranský podprogram pro vytvoření VTK souboru pro ParaView, obsahujícího 3D kartézské souřadnice čáry včetně indexu kroku v proměnné MyScalar, aby jeho sdružením s barevnou škálou bylo možno integrační cestu atraktorem obarvit:
# Python: vytvoření VTK pro ParaView z pole y tvaru (nmax,3)
def savevtk(y,file):
nmax=len(y)
with open(file,'w') as f:
print('# vtk DataFile Version 3.0',file=f)
print('comment (max. 1023 chars)',file=f)
print('ASCII',file=f)
print('DATASET POLYDATA',file=f)
print('POINTS',nmax,'float',file=f)
for n in range(nmax): print(*[f'{y[n,i]:7.3f}' for i in range(3)],file=f)
print('VERTICES 1',nmax+1,file=f)
print(nmax,file=f)
print(*[n for n in range(nmax)],file=f)
print('LINES 1',nmax+1,file=f)
print(nmax,file=f)
print(*[n for n in range(nmax)],file=f)
print('POINT_DATA',nmax,file=f)
print('SCALARS MyScalar float 1',file=f)
print('LOOKUP_TABLE default',file=f)
print(*[n for n in range(nmax)],file=f)
! Fortran: vytvoření VTK pro ParaView z pole y tvaru (nmax,3)
subroutine savevtk(y,file)
real(wp),intent(in) :: y(:,:)
character(*),intent(in) :: file
integer f,nmax,n
nmax=size(y,1)
open (newunit=f,file=file)
write (f,'(a)') '# vtk DataFile Version 3.0'
write (f,'(a)') 'comment (max. 1023 chars)'
write (f,'(a)') 'ASCII'
write (f,'(a)') 'DATASET POLYDATA'
write (f,'(a,i0,a)') 'POINTS ',nmax,' float'
do n=1,nmax; write (f,'(3(f7.3x))') y(n,1:3); enddo
write (f,'(a,i0)') 'VERTICES 1 ',nmax+1
write (f,'(i0)') nmax
write (f,'(*(i0,x))') (n,n=0,nmax-1)
write (f,'(a,i0)') 'LINES 1 ',nmax+1
write (f,'(i0)') nmax
write (f,'(*(i0,x))') (n,n=0,nmax-1)
write (f,'(a,i0)') 'POINT_DATA ',nmax
write (f,'(a)') 'SCALARS MyScalar float 1'
write (f,'(a)') 'LOOKUP_TABLE default'
write (f,'(*(i0,x))') (n,n=0,nmax-1)
close (f)
end subroutine
Několik pojmenovaných podivných atraktorů – Hénonův, Rösslerův a Lorenzův – je popsáno a zobrazeno v našem PDF. Pro okamžité vyzkoušení ParaView přikládáme Rösslerův atraktor uložený naším programem do VTK souboru. Po jeho otevření můžeme zopakovat kroky, které jsme zkusili na třetím cvičení s Lissajousovým obrazcem, a přidat další pro vytvoření průhledné hadice s barevnými kuličkami uvnitř:
File–Open–
odeint-rossler.vtk, v panelu Properties vlevo uprostřed tlačítko Apply(pro symboly v uzlech) Filters–Alphabetical–Glyph (nebo nahoře ve třetím řádku ikonek sedmá ikonka zleva), Apply
a dál v panelu Properties… Glyph Source–Glyph Type: Sphere, Scale–Scale Array: No scale array, upravit Scale Factor, Apply
(pro lepší rozlišení koulí) Glyph Source–Theta/Phi Resolution, Apply (pokud nutno, přepnout předtím Toggle advanced properties ozubeným kolečkem pod otazníkem vpravo nahoře v Properties)
(pro pravidelný rozestup koulí) Masking–Glyph Mode–All Points (nebo Every Nth point), Apply
(pro dojem z plochy) Display–Representation–Surface
(pro obarvení koulí) Coloring–Edit a třeba vpravo u Mapping Data kliknout na ikonce se srdíčkem pro Choose Preset, (pro barvu pozadí) Background v Properties dole
(Apply automaticky) v prvním řádku ikonek kliknout na ikonce shodné s ikonkou Apply pro Apply changes to parameters automatically
(pro hadici místo čáry) vlevo v Pipeline Browser vypnout kliknutím na oko zobrazení Glyph, vrátit se kliknutím na vtk soubor a vybrat Filters–Alphabetical–Tube
a dál v panelu Properties… zvětšit Radius táhlem nebo zadáním čísla
(pro průhlednost) snížit Styling–Opacity
(pro lesk) zvýšit Lighting–Specular
(pro plný dojem při hýbání s objektem) zvýšit Level of Detail pomocí Edit–Settings–Render View–LOD Threshold (stojí za pokus, ale umí to někdy vše zhatit)
(pro uložení obrázku) File–Save Screenshot
Užitečný je i pohled na časový vývoj souřadnic:
vlevo v Pipeline Browser kliknout zpět na vtk soubor
vybrat Filter–Alphabetical–Plot Data
zobrazit vývoj souřadnice
xvýběrem Series Parameters–Points_X v panelu Properties.
Kdybychom filtru Plot Data předložili více VTK souborů podivného atraktoru, startovaných z mírně posunutých počátečních podmínek, shledali bychom po nějaké době trajektorie zcela odchylné, ačkoliv stále v rámci atraktoru – jev známý jako efekt motýlích křídel.
Zdrojové soubory .py, .f90, .m, .gp a Rösslerův atraktor ve .vtk: zip.
