Programování pro fyziky 2019/20: ukázky MATLAB a Python

Následující stránka shrnuje ukázky skriptů pro MATLAB a Python oživujících látku z části Numerické recepty přednášky Programování pro fyziky. Průvodní poznámky k této části přednášky jsou odkazovány z webové stránky přednášky (napřímo zde). Podrobnější výklad nabízí stejnojmenná učebnice Numerical Recipes, volně dostupná ve verzích se zdrojovými texty v C (kniha v C) a Fortranu (kniha ve Fortranu).

Pro skripty demonstrující vybrané numerické recepty volíme dva programovací jazyky vyšší úrovně. Prvním je jazyk komerčního systému MATLAB s variantou jeho volně dostupného klonu GNU Octave (home). Primárním cílem MATLABu byla lineární algebra (Matrix Laboratory), čemuž je dokonale přizpůsobena syntaxe jazyka. Přímočarost tahu na branku MATLABu a Octavu v numerických aplikacích je úžasná.

Druhým jazykem je zjevně uhrančivý Python (home), neboť komunita jeho uživatelů je obrovská, stejně jako záběr rozšiřujících balíčků. Pro numerické aplikace je třeba Python doplnit alespoň trojicí NumPy (Numerical Python), SciPy (Scientific Python) a Matplotlib (plotting library inspired by MATLAB). Balíčky si lze opatřovat po jednom, jednodušší je zvolit některou distribuci Pythonu, kde jich bývají sbaleny stovky, např. Anaconda, Intel nebo jiné distribuce. Současný Python existuje ve dvou částečně nekompatibilních verzích, Python 2 a 3, naše ukázky pracují shodně v obou.

Jistě je možné a vhodné vytěžit grafické nadstavby MATLABu a Pythonu, ale pro spouštění ukázek stačí jen jejich přenesení přes schránku do souborů, např. a.m a a.py, a spouštění z příkazového řádku:

matlab a.m
octave a.m
python a.py

Předložené skripty vlastně jen chystají vstupní data, volají hotové numerické procedury (řešiče) a vracejí výstupy, včetně grafiky, lze-li. Dovnitř řešičů uživatel obvykle nevidí nebo se nedívá, použité postupy spíše jen tuší. V kontrastu k tomu jsme proto zpracovali několik numerických řešičů v klasickém programovacím jazyku užitím postupů strukturovaného programování a komentujeme je na stránkách cvičení k přednášce (2018, 2019). Nesrovnatelně ucelenější sbírka zdrojových souborů numerických metod je ovšem k dispozici v knihách zmíněných v prvním odstavci.

Lineární algebra

MATLAB: operátor \, Python a SciPy: funkce solve

V lineární algebře jsou MATLAB a Octave doma, tam je v úspornosti syntaxe nic nepředběhne. Soustavy lineárních algebraických rovnic řeší operátorem \, pro zkoušku poslouží operátor * pro maticové násobení (umí matice krát matice, matice krát sloupcový vektor, řádkový vektor krát matice). Pro vytvoření pole stačí konstruktor [] bez konverzních funkcí, řádky v konstruktoru se oddělují středníkem, přiřazovací příkazy nebo samotné výrazy nezakončené středníkem výsledek i vypisují:

A=[0 1;-1 0];           # dvouradkova matice
b=[1;0];                # sloupcovy vektor
x=A\b                   # reseni soustavy A.x=b s matici A a pravou stranou b
A*x                     # zkouska nasobenim matice soustavy a vektoru reseni

NumPy a SciPy hledají v MATLABu inspiraci, ale syntaxe je pythonovská. Pole podle NumPy vytváří konverzní funkce array, v konstruktoru matice jsou vnořeny konstruktory řádků, řádkové a sloupcové vektory lze úsporněji inicializovat pomocí triků r_ a c_, funkci solve raději kvalifikujeme jménem zdrojového balíčku, pro výpis ze skriptu potřebujeme volat funkci print a maticové násobení provedeme buď dvouargumentovou funkcí dot nebo použitím metody dot na objekt A. Funkce array a dot jsou ze jmenného prostoru balíčku Numpy, který jsme importovali celý:

from numpy import *
from scipy import linalg as LA
A=array([[0,1],[-1,0]]) # numpy matice
b=c_[[1,0]]             # numpy sloupcovy vektor
x=LA.solve(A,b)         # resic ze scipy.linalg
print(x)
print(dot(A,x))         # proceduralni syntaxe
print(A.dot(x))         # objektova syntaxe

Algebraické soustavy s pásovými a trojúhelníkovými maticemi

K jednoduchým úlohám lineární algebry řadíme řešení soustav s diagonální, trojúhelníkovou a třídiagonální maticí. Matice níže připravujeme z náhodných čísel funkcí rand pro výsledek 8bytového reálného typu v intervalu \(<0,1)\). Pro větší stabilitu řešení posilujeme diagonálu přičtením násobku identické matice získané funkcí eye. Vektory pravých stran vytváříme z jedniček funkcí ones. Mají-li všechny tyto maticové funkce jediný argument, vracejí čtvercovou matici, mají-li dva argumenty, je výsledkem obecně obdélníková matice; je-li jedním z těchto argumentů 1, jde o vektor; skalár je pro MATLAB totéž co matice 1x1. Horní a dolní trojúhelníkovou matici tvoříme nulováním příslušných prvků čtvercové matice funkcemi triu a tril (triangular upper/lower). Dvě tiknutí funkcemi tic a toc odměřují dobu běhu řešiče. Měříme na běžném desktopu s Linuxem a verzemi MATLABu a Pythonu z balíčku linuxové distribuce Ubuntu. S jiným Linuxem nebo ve Windows vycházejí časy jinak, někdy ve prospěch MATLABu, jindy Pythonu a SciPy.

# Obecny resic pro diagonalni, trojuhelnikove a tridiagonalni matice
n=5000; D=diag(rand(1,n)); b=ones(n,1); tic; x=D\b; toc;                 # 0.000 sec
n=5000; L=tril(rand(n))+n*eye(n); b=ones(n,1); tic; x=L\b; toc;          # 0.071 sec
n=5000; U=triu(rand(n))+n*eye(n); b=ones(n,1); tic; x=U\b; toc;          # 0.071 sec
n=5000; A=triu(tril(rand(n),1),-1)+eye(n); b=ones(n,1); tic; x=A\b; toc; # 3.684 sec

Zatím co časová složitost (počet aditivních a multiplikativních operací) řešení soustav s plnou maticí roste kubicky s počtem rovnic, u soustav s diagonální a třídiagonální (obecně s pásovou) maticí je růst lineární a pro trojúhelníkovou matici složitost roste kvadraticky. Z doby běhu chápeme, že operátor \ MATLABu nepochopil speciální tvar třídiagonální matice. Nezbývá než jít mu naproti explicitním převodem řídké třídiagonální matice A (řídkost = přítomnost mnoha nul) na příslušný datový typ funkcí sparse. Pro plné a trojúhelníkové matice je operátor \ MATLABu spojen s procedurami knihovny LAPACK, pro řídké matice operátor \ volá knihovnu UMFPACK z balíčku SuiteSparse. Funkce max pak potvrdí, že největší odchylka vektorů řešení získaných cestou LAPACKu a UMFPACKu není velká ani pro použité matice o tisících řádků:

# Obecny resic pro ridke matice
As=sparse(A); tic; xs=As\b; toc; max(abs(x-xs))                          # 0.000 sec

Podobné efekty v Pythonu získáme pomocí následujícího skriptu. Nejprve importujeme potřebné balíčky, tentokrát ponecháme jmenný prostor NumPy oddělený, abychom vnímali, které funkce odtud přicházejí. Balíček time poslouží pro měření času v duchu tic-toc MATLABu. Obecný řešič scipy.linalg.solve diagonální matici nerozezná, vhodnější (zde i v případě třídiagonální matice) je volat řešič solve_banded. Ten na vstupu očekává jako první argument dvouprvkový tvar pásu (počty l subdiagonál a u superdiagonál) a jako druhý argument pás matice skloněný proti směru hodinových ručiček do obdélníka o l+1+u řádcích. Pro soustavy s trojúhelníkovou maticí je určen řešič solve_triangular, obecný řešič ji nerozezná a jede s kubickou složitostí. Při ověřovacím výpisu více argumentů (zde dvou zřetězených polí) jsme potlačili nadbytečnou informaci na výstupu trikem numpy.r_[ ] neboli numpy.concatenate([ ]).

import numpy as np                                 # import numpy pod synonymem np
from scipy import linalg as LA                     # import scipy.linalg pod synonymem LA
import time

# Soustavy s diagonalni matici
print('solve for diagonal')
n=5000
D=np.diag(np.random.rand(n))                       # diagonalni matice z vektoru nahodnych cisel
b=np.ones(n)                                       # vektor jednicek
tic=time.time(); x=LA.solve(D,b); toc=time.time()  # obecny resic a diagonalni matice
print('Elapsed time is '+str(toc-tic)+' seconds.') # 3.651 sec
print(np.r_[np.dot(D[0:2,:],x),np.dot(D[-2:,:],x)]) # zkouska: prvni a posledni dva prvky prave strany
# efektivnejsi resic
print('solve_banded for diagonal')                 # resic pro pasovou matici
tic=time.time(); x=LA.solve_banded((0,0),np.random.rand(1,n),b); toc=time.time()
print('Elapsed time is '+str(toc-tic)+' seconds.') # 0.000 sec

# Soustavy s trojuhelnikovou matici
print('solve for lower triangular')
n=5000
L=np.tril(np.random.rand(n,n))+n*np.eye(n)         # dolni trojuhelnikova diagonalne dominantni matice
b=np.ones(n)
tic=time.time(); x=LA.solve(L,b); toc=time.time()  # obecny resic a trojuhelnikova matice
print('Elapsed time is '+str(toc-tic)+' seconds.') # 3.843 sec
print('solve for upper triangular')
n=5000
U=np.triu(np.random.rand(n,n))+n*np.eye(n)         # horni trojuhelnikova diagonalne dominantni matice
b=np.ones(n)
tic=time.time(); x=LA.solve(U,b); toc=time.time()  # obecny resic a trojuhelnikova matice
print('Elapsed time is '+str(toc-tic)+' seconds.') # 3.887 sec
# efektivnejsi resic
print('solve_triangular for triu')                 # resic pro trojuhelnikovou matici
tic=time.time(); x=LA.solve_triangular(U,b); toc=time.time() # 0.414 sec
print('Elapsed time is '+str(toc-tic)+' seconds.')

# Soustavy s tridiagonalni matici
print('solve for tridiagonal')
n=5000
A=np.triu(np.tril(np.random.rand(n,n),1),-1)+np.eye(n)
b=np.ones(n)
tic=time.time(); x=LA.solve(A,b); toc=time.time()  # obecny resic a tridiagonalni matice
print('Elapsed time is '+str(toc-tic)+' seconds.') # 3.789 sec
# efektivnejsi resic
print('solve_banded for tridiagonal')              # pasovy resic
tic=time.time(); x=LA.solve_banded((1,1),np.random.rand(3,n),b); toc=time.time()
print('Elapsed time is '+str(toc-tic)+' seconds.') # 0.003 sec

Algebraické soustavy s plnými maticemi

Zatím jen skripty, komentáře později. MATLAB:

# Obecny resic pro plnou matici, LU rozklad, Choleskeho rozklad
n=5000; A=rand(n); b=ones(n,1); tic; x=A\b; toc;            # 3.870 sec
n=5000; tic; [L,U,P]=lu(A); y=L\(P*b); x=U\y; toc;          # 3.881 sec
n=5000; A=A'*A; tic; R=chol(A); y=R'\b; x=R\y; toc;         # 2.011 sec
                tic; L=chol(A,'lower'); y=L\b; x=L'\y; toc; # 2.292 sec
[A(1:2,:)*x;A(end-1:end,:)*x]'

Python, NumPy a SciPy:

import numpy as np
from scipy import linalg as LA
import time

# Soustava s plnou matici - obecny solver
print('solve for general matrix')
n=5000
A=np.random.rand(n,n)+n*np.eye(n)
b=np.ones(n)
tic=time.time(); x=LA.solve(A,b); toc=time.time()
print('Elapsed time is '+str(toc-tic)+' seconds.')          # 3.550 sec

# Soustava s plnou matici - LU solver
print('lu for general matrix')
n=5000
A=np.random.rand(n,n)+n*np.eye(n)
b=np.ones(n)
tic=time.time(); lu,piv=LA.lu_factor(A); x=LA.lu_solve([lu,piv],b); toc=time.time()
print('Elapsed time is '+str(toc-tic)+' seconds.')          # 3.562 sec

# Soustava se symetrickou pozitivne definitni matici - Choleskeho solver
print('cholesky for sym posdev matrix')
n=5000
A=np.random.rand(n,n)+n*np.eye(n)
A=np.dot(np.transpose(A),A)
b=np.ones(n)
tic=time.time(); c,lower=LA.cho_factor(A); x=LA.cho_solve([c,lower],b); toc=time.time()
print('Elapsed time is '+str(toc-tic)+' seconds.')          # 1.961 sec
print(np.r_[np.dot(A[0:2,:],x),np.dot(A[-2:,:],x)]) # zkouska: prvni a posledni dva prvky prave strany