Cvičení č. 10: Obyčejné diferenciální rovnice

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 – za domácí úkol – 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,

\[\frac{dy_r(x)}{dx}=f_r(x,y_1(x),\dots,y_R(x)),\ r=1,\dots,R,\]

neboli vektorově

\[\frac{d\mathbf{y}(x)}{dx}=\mathbf{f}(x,\mathbf{y}(x)),\]

doplněná počáteční podmínkou,

\[\mathbf{y}(x_0)=\mathbf{y}_0.\]

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

\[\mathbf{y}_{n+1}=\mathbf{y}_n+h\mathbf{f}(x_n,\mathbf{y}_n),\quad n=0,1,\dots,\]

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:

\[\begin{split}\mathbf{y}_{n+1}=\mathbf{y}_n+h\left(\frac16\mathbf{k}_1+\frac13\mathbf{k}_2+\frac13\mathbf{k}_3+\frac16\mathbf{k}_4\right),\quad n=0,1,\dots,\\ \mathbf{k}_1=\mathbf{f}(x_n,\mathbf{y}_n), \\ \mathbf{k}_2=\mathbf{f}(x_n+\frac12h,\mathbf{y}_n+\frac12h\mathbf{k}_1), \\ \mathbf{k}_3=\mathbf{f}(x_n+\frac12h,\mathbf{y}_n+\frac12h\mathbf{k}_2), \\ \mathbf{k}_4=\mathbf{f}(x_n+h,\mathbf{y}_n+h\mathbf{k}_3).\end{split}\]

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:

\[dy(x)/dx=-y(x),\quad y(0)=1,\]

jejímž analytickým řešením je zjevně \(y(x)=e^{-x}\), a soustavu:

\[\begin{split}\begin{array}{ll} dy_1(x)/dx=y_2(x), & y_1(0)=0, \\ dy_2(x)/dx=-y_1(x),& y_2(0)=1, \end{array}\end{split}\]

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í:

_images/ode-exp.png

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.

_images/ode-sin.png

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:

\[\begin{split}\frac{d}{dt}\left(\begin{array}{r} x \\ y \\ z \\ \end{array}\right) = \left(\begin{array}{l} \sigma(y-x) \\ \rho x-y-xz \\ xy-\beta z \\ \end{array}\right), \ \sigma=10,\ \beta=8/3,\ \rho=28.\end{split}\]

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 PropertiesGlyph 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 x vý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.

_images/ode1-rossler.png _images/ode2-lorenz.png _images/ode3-lorenz.png