Programování pro fyziky 2022/23 – Numerické recepty č. 0

Náš průvodní text je především v poznámkách k přednášce. Pomáhat si budeme i odkazy do knih Numerical Recipes in C a Numerical Recipes in Fortran (krátce: NR), které jsou (v Second Edition) na webu volně k dispozici. Vidět je ovšem lze jen z Adobe Acrobat Readeru po instalaci pluginu FileOpen (tj. ve Firefoxu myší na Save Link As a v Downloads myší na Open in Adobe Acrobat Reader; interní PDF prohlížeče Firefoxu a Chromu nám nepomohou). Zde si připravíme ukázky skriptů pro Matlab/Octave a také Python s NumPy a SciPy. V těchto prostředích je k demonstraci numerických metod krátká cesta.

Mini-algoritmy

Řešení kvadratické rovnice, NR kap. 5.6 (ještě jednou: Save Link As & Open in Adobe Acrobat Reader). 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 odečítají blízká čísla (katastrofa!) a ztrácí se přesnost. Lepší to bude podle vzorců, 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 zjevně druhou z našich 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
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
q=-(b+numpy.sign(b)*math.sqrt(b*b-4*a*c))/2; x1=q/a; x2=c/q
print(x1,x2)     # robustní postup
a=1; b=1e7; c=1; x=numpy.roots([a,b,c])
print(x)         # černá skříňka
print([f'{x[i]:.15e}' for i in range(x.size)])
# -9999999.9999999 -9.96515154838562e-08
# -9999999.9999999 -1.00000000000001e-07
# [-1.e+07 -1.e-07]
# ['-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 fuč! 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í modulu 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
import matplotlib.pyplot as plt
x=numpy.arange(0,numpy.pi*3,.01)  # výpočetní síť
y=numpy.sin(x)                    # funkční hodnoty na ní
dydx=numpy.diff(y)/numpy.diff(x)  # numerická derivace
# for i in range(x.size-1): print(f'{x[i]:20.15f} {numpy.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 operace * 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:

% 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
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=numpy.polyval([1,0,0],numpy.arange(1,6))
print(yy)                          # vyčíslení polynomu x**2 na síti
x=numpy.linspace(0,2,21); p=[0.5,0,-1]
plt.plot(x,numpy.polyval(p,x))     # vykreslení polynomu 0.5*x**2-1
plt.show()