Algoritmusok Python nyelven

10. Előadás: SciPy Numerikus módszerek

2020. április 23.

Tudományos csomagok (Python Scientific stack)

csomag
NumPy Hatékony N-dimenziós tömb
SciPy Numerikus számítások
Matplotlib Grafikonok és rajzok
IPython (Jupyter) Interaktív notebook
SymPy Szimbolikus számítások
Pandas Adatbányászat
In [2]:
import sympy as sym
import math
import numpy as np
import matplotlib.pyplot as plt
import scipy

Numerikus módszerek

A numerikus számítások célja, hogy a matematikai számításokat végezzünk hatékonyan.

Előnyök:

  • Gyors
  • Hasznos
  • Sok helyzetben alakalmazható

Hátrányok:

  • Csak közelítő eredményt ad
  • Emiatt figyelni kell a hiábkra
  • A matematikai intuíciónkat hátráltatja elméleti kérdésekben (pl felismernéd, hogy ez melyik algebrai szám közelítése? 2.732050807568877)

SciPy

Részletes tutorial: https://docs.scipy.org/doc/scipy/reference/tutorial/

A háttérben sokszor Fortran fut.

Alcsomagok:
scipy.cluster Vector quantization / Kmeans
scipy.constants Physical and mathematical constants
scipy.fftpack Fourier transform
scipy.integrate Integration routines
scipy.interpolate Interpolation
scipy.io Data input and output
scipy.linalg Linear algebra routines
scipy.ndimage n-dimensional image package
scipy.odr Orthogonal distance regression
scipy.optimize Optimization
scipy.signal Signal processing
scipy.sparse Sparse matrices
scipy.spatial Spatial data structures and algorithms
scipy.special Any special mathematical functions
scipy.stats Statistics
In [3]:
from scipy import linalg, optimize

Lineáris algebra

LU felbontás

In [4]:
from scipy.linalg import lu
A = np.array([[2, 5, 8, 7], [5, 2, 2, 8], [7, 5, 6, 6], [5, 4, 4, 8]])
p, l, u = lu(A)
print(A,p,l,u, sep="\n")
np.allclose(A - p @ l @ u, np.zeros((4, 4)))

[[2 5 8 7]
 [5 2 2 8]
 [7 5 6 6]
 [5 4 4 8]]
[[0. 1. 0. 0.]
 [0. 0. 0. 1.]
 [1. 0. 0. 0.]
 [0. 0. 1. 0.]]
[[ 1.          0.          0.          0.        ]
 [ 0.28571429  1.          0.          0.        ]
 [ 0.71428571  0.12        1.          0.        ]
 [ 0.71428571 -0.44       -0.46153846  1.        ]]
[[ 7.          5.          6.          6.        ]
 [ 0.          3.57142857  6.28571429  5.28571429]
 [ 0.          0.         -1.04        3.08      ]
 [ 0.          0.          0.          7.46153846]]
Out[4]:
True

QR felbontás Q ortogonális, R felső háromszög

In [5]:
A = np.random.randn(4, 4)
print(A)
q, r = linalg.qr(A)
print(q,r,sep='\n')

[[ 0.26353773  0.86973603 -0.18084589 -0.9486353 ]
 [ 1.07413653  1.3322248  -1.56374712  0.14971633]
 [ 0.04096525 -0.49581744  0.30208855  0.87311366]
 [-1.44991248 -0.46878984  1.11244664 -1.83082867]]
[[-0.14447972 -0.58402801  0.79582087  0.06860053]
 [-0.58887562 -0.49664755 -0.50494596  0.38935369]
 [-0.02245845  0.44642683  0.24949115  0.85904183]
 [ 0.79488788 -0.4614711  -0.22238005  0.32518426]]
[[-1.82404654 -1.27167399  1.82446702 -1.42601803]
 [ 0.         -1.17460964  0.50374871  1.71432922]
 [ 0.          0.          0.47366933 -0.20556853]
 [ 0.          0.          0.          0.14790021]]
In [6]:
q @ q.T
Out[6]:
array([[ 1.00000000e+00, -3.97395674e-16, -1.86832824e-17,
         2.58480705e-16],
       [-3.97395674e-16,  1.00000000e+00,  2.42618836e-17,
         1.38003999e-16],
       [-1.86832824e-17,  2.42618836e-17,  1.00000000e+00,
         1.45979610e-17],
       [ 2.58480705e-16,  1.38003999e-16,  1.45979610e-17,
         1.00000000e+00]])

Integrálás

A scipy.integrate.quad(f, a, b) integrálja a-tól b-ig az f függvényt. Két értékkel tér vissza, az inregrálás eredményével és egy hibahatárral.

In [7]:
import scipy.integrate
f= lambda x:2*x
i = scipy.integrate.quad(f, 0, 1)
print(i)

(1.0, 1.1102230246251565e-14)
In [8]:
f= lambda x:np.exp(-x**2)
i = scipy.integrate.quad(f, 0, 1)
print(i)

(0.7468241328124271, 8.291413475940725e-15)
In [9]:
## szimbolikus számítással
xval = sym.Symbol('x')
sym.integrate(sym.exp(-xval ** 2),  xval)
Out[9]:
$\displaystyle \frac{\sqrt{\pi} \operatorname{erf}{\left(x \right)}}{2}$
In [10]:
((sym.erf(1)-sym.erf(0))*np.pi**(1/2)/2).evalf()
Out[10]:
$\displaystyle 0.746824132812427$

Próbáljunk ki egy ronda függvényt

In [11]:
#sym.integrate(sym.exp(sym.sin(-sym.log(xval ** 2+xval)-xval)),  xval)
In [12]:
f= lambda x:np.exp(np.sin(-np.log(x ** 2+x)-x))
i = scipy.integrate.quad(f, 0, 1)
print(i)

(1.1008683826459191, 9.950209101106111e-09)

Több változóban is tudunk integrálni a scipy.integrate.dblquad(func, a, b, gfun, hfun) paranccsal. A gfun és hfun parancsok adják meg a belső integrál határait.

$$\int_{0}^{1/2} \big( \int_{0}^{\sqrt{1-4y^2}} 16xy \:dx\big)dy$$
In [13]:
f = lambda x, y : 16*x*y
g = lambda x : 0
h = lambda y : math.sqrt(1-4*y**2)
i = scipy.integrate.dblquad(f, 0, 0.5, g, h)
print(i)

(0.5, 1.7092350012594845e-14)

Interpoláció

In [14]:
from scipy import interpolate

x = np.linspace(0, 4, 12)
y = np.cos(x**2/3+4)
print (x,y)

[0.         0.36363636 0.72727273 1.09090909 1.45454545 1.81818182
 2.18181818 2.54545455 2.90909091 3.27272727 3.63636364 4.        ] [-0.65364362 -0.61966189 -0.51077021 -0.31047698 -0.00715476  0.37976236
  0.76715099  0.99239518  0.85886263  0.27994201 -0.52586509 -0.99582185]
In [15]:
figure, axis = plt.subplots(1, 1,figsize=(5,5))
axis.plot(x, y,'o')
plt.show()
In [16]:
f1 = interpolate.interp1d(x, y,kind = 'linear')

f2 = interpolate.interp1d(x, y, kind = 'cubic')
In [17]:
xnew = np.linspace(0, 4,30)

#plt.plot(x, y, 'o', xnew, f(xnew), '-', xnew, f2(xnew), '--')
figure, axis = plt.subplots(1, 1,figsize=(5,5))
axis.plot(x, y, 'o')
axis.plot( xnew, f1(xnew), '-')
axis.plot( xnew, f2(xnew), '--')

axis.legend(['data', 'linear', 'cubic'], loc = 'best')

plt.show()
In [18]:
from scipy.interpolate import UnivariateSpline
x = np.linspace(-3, 3, 50)
y = np.exp(-x**2) + 0.1 * np.random.randn(50)
figure, axis = plt.subplots(1, 1,figsize=(5,5))
axis.plot(x, y, 'ro', ms = 5)
plt.show()
In [19]:
f2 = interpolate.interp1d(x, y, kind = 'cubic')
xnew = np.linspace(-3, 3,30)

#plt.plot(x, y, 'o', xnew, f(xnew), '-', xnew, f2(xnew), '--')
figure, axis = plt.subplots(1, 1,figsize=(5,5))
axis.plot(x, y, 'o')
axis.plot( xnew, f2(xnew), '--')

axis.legend(['data', 'cubic' ], loc = 'best')

plt.show()
In [20]:
spl = UnivariateSpline(x, y)
xs = np.linspace(-3, 3, 1000)
figure, axis = plt.subplots(1, 1,figsize=(5,5))
axis.plot(xs, spl(xs), 'g', lw = 3)
axis.plot(x, y, 'ro', ms = 5)
plt.show()
In [24]:
spl.set_smoothing_factor(0.5)
figure, axis = plt.subplots(1, 1,figsize=(5,5))
axis.plot(x, y, 'ro') 
axis.plot(xs, spl(xs)) 
axis.plot(x,np.exp(-x**2), ms = 5)
axis.legend(['data', 'spline', 'original'], loc = 'best')
plt.show()

B-spline reprezentáció

In [25]:
t=np.linspace(0,1.75*2*np.pi,100)

x = np.sin(t)
y = np.cos(t)
z = t

# add noise
x+= np.random.normal(scale=0.1, size=x.shape)
y+= np.random.normal(scale=0.1, size=y.shape)
z+= np.random.normal(scale=0.1, size=z.shape)

# spline parameters
s=3.0 # smoothness parameter
k=2 # spline order
nest=-1 # estimate of number of knots needed (-1 = maximal)

# find the knot points
tckp,u = scipy.interpolate.splprep([x,y,z],s=s,k=k,nest=-1)

# evaluate spline, including interpolated points
xnew,ynew,znew = scipy.interpolate.splev(np.linspace(0,1,400),tckp)

figure, axis = plt.subplots(2, 2,figsize=(6,6))

axis[0,0].plot(x,y,'bo-',label='data')
axis[0,0].plot(xnew,ynew,'r-',label='fit')

axis[0,1].plot(x,z,'bo-',label='data')
axis[0,1].plot(xnew,znew,'r-',label='fit')

axis[1,0].plot(y,z,'bo-',label='data')
axis[1,0].plot(ynew,znew,'r-',label='fit')



plt.show()

Legkisebb négyzetek módszere

In [26]:
f = np.poly1d([5, 1])

x = np.linspace(0, 10, 30)
y = f(x) + 6*np.random.normal(size=len(x))
xn = np.linspace(0, 10, 200)

figure, axis = plt.subplots(1, 1,figsize=(6,6))
axis.plot(x, y, 'or')
plt.show()
In [28]:
a = np.vstack([x, np.ones(len(x))]).T
np.dot(np.linalg.inv(np.dot(a.T, a)), np.dot(a.T, y))
Out[28]:
array([5.03380086, 2.03710585])
In [29]:
np.linalg.lstsq(a, y,rcond=None)[0]
Out[29]:
array([5.03380086, 2.03710585])
In [30]:
np.polyfit(x, y, 1)
Out[30]:
array([5.03380086, 2.03710585])
In [31]:
m, c = np.polyfit(x, y, 1)
yn = np.polyval([m, c], xn)

figure, axis = plt.subplots(1, 1,figsize=(6,6))
axis.plot(x, y, 'or')
axis.plot(xn, yn)
plt.show()

Tutorial: https://mmas.github.io/least-squares-fitting-numpy-scipy

Gyök keresés

A scipy.optimize.root(fn,x0) függvényt használhatjuk. A hátterben különböző iteratív módszerek futnak, amik folyamatosan javítanak a megoldáson. A második paraméter a kezdeti tipp a gyökre.

In [32]:
from scipy.optimize import root
def func(x):
   return x*2 + 2 * np.cos(x)
sol = root(func,0)
print(sol)

    fjac: array([[-1.]])
     fun: array([0.])
 message: 'The solution converged.'
    nfev: 9
     qtf: array([-5.33573186e-13])
       r: array([-3.34722404])
  status: 1
 success: True
       x: array([-0.73908513])
In [33]:
func(-0.73908513)
Out[33]:
1.0761863178387898e-08
In [34]:
sym.solveset( xval*2 + 2 * sym.cos(xval), xval) # a szimolikus megoldó csak széttárja a kezét.
Out[34]:
$\displaystyle \left\{x \mid x \in \mathbb{C} \wedge 2 x + 2 \cos{\left(x \right)} = 0 \right\}$

Geometriai algoritmusok

In [37]:
from scipy.spatial import Delaunay
points = np.random.rand(30, 2) # 30 random points in 2-D
tri = Delaunay(points)
import matplotlib.pyplot as plt
plt.triplot(points[:,0], points[:,1], tri.simplices.copy())
plt.plot(points[:,0], points[:,1], 'o')
plt.show()
In [38]:
from scipy.spatial import ConvexHull
points = np.random.rand(30, 2) # 30 random points in 2-D
hull = ConvexHull(points)
plt.plot(points[:,0], points[:,1], 'o')
for simplex in hull.simplices:
    plt.plot(points[simplex,0], points[simplex,1], 'k-')
plt.show()

Jelfeldolgozás

In [40]:
# How many time points are needed i,e., Sampling Frequency
samplingFrequency   = 100;

# At what intervals time points are sampled
samplingInterval       = 1 / samplingFrequency;

# Begin time period of the signals
beginTime           = 0;

# End time period of the signals
endTime             = 10; 

# Frequency of the signals
signal1Frequency     = 11;
signal2Frequency     = 7;

# Time points
time        = np.arange(beginTime, endTime, samplingInterval);

# Create two sine waves
amplitude1 = np.sin(2*np.pi*signal1Frequency*time)
amplitude2 = np.sin(2*np.pi*signal2Frequency*time)

# Create subplot
figure, axis = plt.subplots(4, 1,figsize=(10,10))
plt.subplots_adjust(hspace=1)

# Time domain representation for sine wave 1
axis[0].set_title('Sine wave with a frequency of 4 Hz')
axis[0].plot(time, amplitude1)
axis[0].set_xlabel('Time')
axis[0].set_ylabel('Amplitude')

# Time domain representation for sine wave 2
axis[1].set_title('Sine wave with a frequency of 7 Hz')
axis[1].plot(time, amplitude2)
axis[1].set_xlabel('Time')
axis[1].set_ylabel('Amplitude')

# Add the sine waves
amplitude = amplitude1 + amplitude2

# Time domain representation of the resultant sine wave
axis[2].set_title('Sine wave with multiple frequencies')
axis[2].plot(time, amplitude)
axis[2].set_xlabel('Time')
axis[2].set_ylabel('Amplitude')

# Frequency domain representation
fourierTransform = np.fft.fft(amplitude)/len(amplitude)           # Normalize amplitude
fourierTransform = fourierTransform[range(int(len(amplitude)/2))] # Exclude sampling frequency
 
tpCount     = len(amplitude)
values      = np.arange(int(tpCount/2))
timePeriod  = tpCount/samplingFrequency
frequencies = values/timePeriod
# Frequency domain representation

axis[3].set_title('Fourier transform depicting the frequency components')
axis[3].plot(frequencies, abs(fourierTransform))
axis[3].set_xlabel('Frequency')
axis[3].set_ylabel('Amplitude')

plt.show()
In [ ]: