Programování pro fyziky 2023/24 – Numerické recepty č. 0

V druhé části našich přednášek se budeme věnovat numerickým metodám, tvořícím hlavní stavební kameny fyzikálních výpočtů. Zahrneme mezi ně základní metody lineární algebry, aproximaci funkcí, integrování, řešení nelineárních rovnic a soustav obyčejných diferenciálních rovnic. V numerické praxi si pro tyto metody chodíme do knihoven psaných typicky v C/C++ a Fortranu. Čerpat z těchto knihoven je ovšem nejjednodušší v prostředích interaktivních numerických systémů, k nimž se řadí Python s balíčky NumPy a SciPy, Matlab (resp. jeho volně dostupný klon GNU Octave), Mathematica a jiné.

Při výkladu si budeme často pomáhat odkazy do knih Numerical Recipes in C++, Numerical Recipes in C a Numerical Recipes in Fortran (krátce: NR) a také Numerical Computing with MATLAB, které jsou na webu volně k dispozici. Výklad budou provázet ukázky skriptů pro Python s balíčky NumPy a SciPy a také pro Matlab a GNU Octave. V těchto prostředích je k demonstraci numerických metod krátká cesta.

Mini-algoritmy

Do začátku si ukážeme, jak omezená přesnost reálné aritmetiky ovlivňuje i nejběžnější algoritmy numerické matematiky. Hlavní potíží zde bude odčítání blízkých čísel provázené katastrofickou ztrátou platných číslic. Omezený výkon reálné aritmetiky (řádově miliardy základních operací za sekundu) nás bude nutit sledovat také asymptotickou časovou složitost algoritmů, a obdobně typická velikost paměti počítačů povede k limitům na paměťové nároky algoritmů, případně k potřebě paralelizace algoritmů pro jejich provozování na výpočetních klastrech s distribuovanou pamětí.

Řešení kvadratické rovnice, NR kap. 5.6. Přímočarý pokus o \(x_{1,2}=(-b\pm\sqrt{b^2-4ac})/(2a)\)

% Matlab/Octave: řešení kvadratické rovnice podle běžného vzorce
format compact; format long
a=1; b=1e7; c=1; D=sqrt(b*b-4*a*c); x1=(-b-D)/(2*a), x2=(-b+D)/(2*a)
% x1 = -9999999.999999899
% x2 = -9.965151548385620e-08

numericky selže, v závorkách výrazu pro x2 se odčítají blízká čísla a ztrácí se přesnost. Lepší to bude postupem, kde se nic neodečítá, tedy s využitím vztahu \(x_1 x_2=c/a\):

% Matlab/Octave: řešení kvadratické rovnice robustně
q=-(b+sign(b)*sqrt(b*b-4*a*c))/2; x1=q/a, x2=c/q
% x1 = -9999999.999999899
% x2 = -1.000000000000010e-07

Funkce roots pro kořeny polynomu (zadaného v argumentu jako konstruktor pole s prvky v pořadí od nejvyšší mocniny) jde patrně druhou z těchto cest:

% Matlab/Octave: řešení kvadratické rovnice černou skříňkou
x=roots([a b c])
% x = -9.999999999999899e+06
%     -1.000000000000010e-07

Reálná aritmetika Pythonu je stejně jako v Matlabu defaultně 8bytová, tedy následujícím skriptem obdržíme stejné výsledky:

# Python: řešení kvadratické rovnice
import math
import numpy as np
a=1; b=1e7; c=1; D=math.sqrt(b*b-4*a*c); x1=(-b-D)/(2*a); x2=(-b+D)/(2*a)
print(x1,x2)     # podle běžného vzorce      ... -9999999.9999999 -9.96515154838562e-08
q=-(b+np.sign(b)*math.sqrt(b*b-4*a*c))/2; x1=q/a; x2=c/q
print(x1,x2)     # robustní postup           ... -9999999.9999999 -1.00000000000001e-07
a=1; b=1e7; c=1; x=np.roots([a,b,c])
print(x)         # černá skříňka             ... [-1.e+07 -1.e-07]
print([f'{x[i]:.15e}' for i in range(x.size)]) # ['-9.999999999999899e+06', '-1.000000000000010e-07']

Jak vidíme na funkci numpy.roots (a budeme vidět i dále), modul NumPy (a SciPy také) bere pro názvy a rozhraní svých funkcí hojně inspiraci v Matlabu.

Numerické derivování, NR kap. 5.7, obsahuje odečítání blízkých čísel už ve své definici, \(f'(x)\approx (f(x+h)-f(x))/h\). To je zapeklitá situace, dokonce lze ukázat (NR (5.7.5) a (5.7.6)), že i při optimální volbě kroku \(h\) se relativní přesnost vyčíslení funkce \(\epsilon_f\) (tedy řekněme 1e-15 v 8bytovém float) sníží na \(\sqrt{\epsilon_f}\) … polovina platných míst je pryč! Pro přesné počítání je numerické derivování oříšek.

Ale kreslení to nevadí, tam pracujeme jen s 3-4 platnými místy. V následujících skriptech vypočteme na intervalu \(\langle 0,3\pi\rangle\) s krokem 0.1 funkční hodnoty funkce sin a její derivace, získané jako pole podílů obsahujících v čitateli rozdíly sousedících prvků pole y (funkce diff) a ve jmenovateli podobně pro x, tedy pole aproximovaných hodnot první derivace funkce sin. Pole y a dydx pomocí funkce plot vykreslíme do společného obrázku, který prosím získejte sami.

% Matlab/Octave: numerické derivování
x=0:.1:pi*3; y=sin(x);     % výpočetní síť a funkční hodnoty na ní
dydx=diff(y)./diff(x);     % numerická derivace
plot(x,y,x(1:end-1),dydx)  % obrázek s dvěma čarami
saveas(gcf,'dydx-m.png')   % uložení obrázku (current figure) do PNG

V Pythonu je zvykem kreslit grafy z NumPy polí pomocí balíčku Matplotlib. Buď se rovnou replikují kroky obdobné jako v Matlabu (plot NumPy polí a uložení do PNG, PDF aj. podle přípony souboru), nebo se parametry obrázku konfigurují podrobněji a pracuje se v objektové syntaxi:

# Python: numerické derivování
import numpy as np
import matplotlib.pyplot as plt
x=np.arange(0,np.pi*3,.01)        # výpočetní síť
y=np.sin(x)                       # funkční hodnoty na ní
dydx=np.diff(y)/np.diff(x)        # numerická derivace
# for i in range(x.size-1): print(f'{x[i]:20.15f} {np.cos(x[i]):20.15f} {dydx[i]:20.15f}')
# fig=plt.figure(figsize=(16,8),layout='tight')  # volitelná konfigurace obrázku
plt.plot(x,y,x[:-1],dydx)         # obrázek s dvěma čarami
plt.savefig('dydx-py.png')        # uložení do PNG
plt.show()                        # vykreslení na obrazovku
# fig.savefig('dydx-py.png')      # alternativní metoda pro uložení PNG

Vyčíslení polynomu – Hornerovo schéma, NR kap. 5.3

Pro výpočet \(P_2(x)=c_1x^2+c_2x+c_3\) bychom pro úsporu jedné ze tří operací * místo

% Matlab/Octave: vyčíslení polynomu
n=2; c=[1 0 0]; x=0.5;
p=c(1)*x*x+c(2)*x+c(3)
% p = 0.250000000000000

měli vytýkat:

% Matlab/Octave: Hornerovo schéma ručně
p=(c(1)*x+c(2))*x+c(3)
% p = 0.250000000000000

Pro polynomy vyššího stupně se vyplatí cyklus s lineární časovou náročností (počet + a * úměrný n):

% Matlab/Octave: Hornerovo schéma v cyklu
p=c(1); for i=2:n+1, p=p*x+c(i); end; p
% p = 0.250000000000000

Jinak je ovšem pro vyčíslování polynomů připravena funkce polyval, kterou můžeme s polem x v druhém argumentu použít rovnou pro vykreslení polynomu pomocí plot:

% Matlab/Octave: vyčíslení polynomu funkcí polyval
polyval([1 0 0],1:5)                         % vyčíslení polynomu x**2 na síti
% ans =    1    4    9   16   25
x=0:.1:2; p=[0.5 0 -1]; plot(x,polyval(p,x)) % vykreslení polynomu 0.5*x**2-1

Totéž v Pythonu, včetně funkce numpy.polyval:

# Python: vyčíslení polynomu
import numpy as np
import matplotlib.pyplot as plt
n=2; c=[1,0,0]; x=0.5
p=c[0]*x*x+c[1]*x+c[2]             # doslovný zápis s mocninami
print(p)
p=(c[0]*x+c[1])*x+c[2]             # Hornerovo schéma ručně
print(p)
p=c[0]
for i in range(1,n+1): p=p*x+c[i]  # Hornerovo schéma v cyklu
print(p)
yy=np.polyval([1,0,0],np.arange(1,6))
print(yy)                          # vyčíslení polynomu x**2 na síti
x=np.linspace(0,2,21); p=[0.5,0,-1]
plt.plot(x,np.polyval(p,x))        # vykreslení polynomu 0.5*x**2-1
plt.show()