Programowanie Wykład 5 - Przegląd bibliotek naukowych: SciPy

33  Download (0)

Pełen tekst

(1)

Programowanie

Wykład 5 - Przegląd bibliotek naukowych: SciPy

Plan wykładu:

Wprowadzenie Funkcje specjalne Wielomiany

Całkowanie numeryczne

Równania różniczkowe zwyczajne Transformaty Fouriera

Algebra liniowa Optymalizacja Interpolacja Statystyka Materiały do wykładu:

wykład J.R. Johanssona "Scientific Python" dostępny pod adresem

http://github.com/jrjohansson/scientific-python-lectures (http://github.com/jrjohansson/scientific-python- lectures)

domumentacja SciPy: http://docs.scipy.org/doc/scipy/reference/

(http://docs.scipy.org/doc/scipy/reference/)

Wprowadzenie

In [1]:

Pakiet SciPy to zestaw algorytmów numerycznych i przydatnych funkcji, dzięki którym Python może zastąpić dedykowane środowiska do obliczeń numerycznych (Matlab, GNU Octave, Scilab). Pakiet wykorzystuje NumPy do tworzenia wielowymiarowych struktur danych i operacji na nich. Pakiet podzielony jest na podmoduły, z których większość odpowiada tematycznie wybranej dziedzinie metod numerycznych. Oto wybrane z nich:

Moduł Funkcjonalność

scipy.special funkcje specjalne scipy.integrate całkowanie numeryczne scipy.optimize optymalizacja

scipy.interpolate interpolacja

scipy.fftpack transformaty Fouriera scipy.signal przetwarzanie sygnałów scipy.linalg algebra liniowa

#więcej na kolejnych wykładach

%matplotlib inline

(2)

Moduł Funkcjonalność

scipy.sparse zagadnienia własne macierzy rzadkich scipy.stats statystyka

scipy.ndimage przetwarzanie obrazów scipy.io operacje wejścia/wyjścia Moduł importujemy poleceniem:

In [2]:

Jeśli nie potrzebujemy całej funkcjonalności SciPy, możemy importować pojedyncze podmoduły, np.:

In [3]:

Z poziomu interpretera mamy dostęp do dokumentacji poszczególnych podmodułów:

In [4]:

lub pojedynczych funkcji:

====================================

Linear algebra (:mod:`scipy.linalg`)

====================================

.. currentmodule:: scipy.linalg Linear algebra functions.

.. seealso::

`numpy.linalg` for more linear algebra functions. Note that although `scipy.linalg` imports most of them, identically named functions from `scipy.linalg` may offer more or slightly differi ng

functionality.

Basics

======

import scipy as sp import numpy as np

import matplotlib.pyplot as plt

import scipy.linalg as la

sp.info(la)

(3)

Możemy również podglądnąć kod źródłowy interesującej nas funkcji z poziomu interpretera:

lu(a, permute_l=False, overwrite_a=False, check_finite=True) Compute pivoted LU decomposition of a matrix.

The decomposition is::

A = P L U

where P is a permutation matrix, L lower triangular with unit diagonal elements, and U upper triangular.

Parameters ---

a : (M, N) array_like Array to decompose permute_l : bool, optional

Perform the multiplication P*L (Default: do not permute) overwrite_a : bool, optional

Whether to overwrite data in a (may improve performance) check_finite : bool, optional

Whether to check that the input matrix contains only finite numbe rs.

Disabling may give a performance gain, but may result in problems (crashes, non-termination) if the inputs do contain infinities or NaNs.

Returns ---

**(If permute_l == False)**

p : (M, M) ndarray Permutation matrix l : (M, K) ndarray

Lower triangular or trapezoidal matrix with unit diagonal.

K = min(M, N) u : (K, N) ndarray

Upper triangular or trapezoidal matrix

**(If permute_l == True)**

pl : (M, K) ndarray Permuted L matrix.

K = min(M, N) u : (K, N) ndarray

Upper triangular or trapezoidal matrix Notes

---

This is a LU factorization routine written for Scipy.

sp.info(la.lu)

(4)

In [6]:

In file: /usr/local/lib/python3.6/dist-packages/scipy/linalg/decomp_l u.py

def lu(a, permute_l=False, overwrite_a=False, check_finite=True):

"""

Compute pivoted LU decomposition of a matrix.

The decomposition is::

A = P L U

where P is a permutation matrix, L lower triangular with unit diagonal elements, and U upper triangular.

Parameters ---

a : (M, N) array_like Array to decompose permute_l : bool, optional

Perform the multiplication P*L (Default: do not permute) overwrite_a : bool, optional

Whether to overwrite data in a (may improve performance) check_finite : bool, optional

Whether to check that the input matrix contains only finite n umbers.

Disabling may give a performance gain, but may result in prob lems

(crashes, non-termination) if the inputs do contain infinitie s or NaNs.

Returns ---

**(If permute_l == False)**

p : (M, M) ndarray Permutation matrix l : (M, K) ndarray

Lower triangular or trapezoidal matrix with unit diagonal.

K = min(M, N) u : (K, N) ndarray

Upper triangular or trapezoidal matrix **(If permute_l == True)**

pl : (M, K) ndarray Permuted L matrix.

K = min(M, N) u : (K, N) ndarray

Upper triangular or trapezoidal matrix Notes

---

This is a LU factorization routine written for Scipy.

"""

sp.source(la.lu)

(5)

Oczywiście, nadal działa polecenie dir:

In [7]:

SciPy kontra NumPy

SciPy wykorzystuje NumPy jako silnik do szybkiego manipulowania wielowymiarowymi strukturami danych.

Modyfikuje przy tym niektóre operacje:

if len(a1.shape) != 2:

raise ValueError('expected matrix')

overwrite_a = overwrite_a or (_datacopied(a1, a)) flu, = get_flinalg_funcs(('lu',), (a1,))

p, l, u, info = flu(a1, permute_l=permute_l, overwrite_a=overwrit e_a)

if info < 0:

raise ValueError('illegal value in %d-th argument of '

'internal lu.getrf' % -in fo)

if permute_l:

return l, u return p, l, u

['LinAlgError', '__all__', '__builtins__', '__cached__', '__doc__', '__file__', '__loader__', '__name__', '__package__', '__path__', '__s pec__', '__version__', '_decomp_polar', '_decomp_qz', '_decomp_updat e', '_expm_frechet', '_fblas', '_flapack', '_flinalg', '_matfuncs_sqr tm', '_procrustes', '_sketches', '_solve_toeplitz', '_solvers', 'abso lute_import', 'basic', 'blas', 'block_diag', 'cho_factor', 'cho_solv e', 'cho_solve_banded', 'cholesky', 'cholesky_banded', 'circulant', 'clarkson_woodruff_transform', 'companion', 'coshm', 'cosm', 'cython_

blas', 'cython_lapack', 'decomp', 'decomp_cholesky', 'decomp_lu', 'de comp_qr', 'decomp_schur', 'decomp_svd', 'det', 'dft', 'diagsvd', 'div ision', 'eig', 'eig_banded', 'eigh', 'eigh_tridiagonal', 'eigvals', 'eigvals_banded', 'eigvalsh', 'eigvalsh_tridiagonal', 'expm', 'expm_c ond', 'expm_frechet', 'find_best_blas_type', 'flinalg', 'fractional_m atrix_power', 'funm', 'get_blas_funcs', 'get_lapack_funcs', 'hadamar d', 'hankel', 'helmert', 'hessenberg', 'hilbert', 'inv', 'invhilber t', 'invpascal', 'kron', 'lapack', 'leslie', 'linalg_version', 'log m', 'lstsq', 'lu', 'lu_factor', 'lu_solve', 'matfuncs', 'matrix_balan ce', 'misc', 'norm', 'ordqz', 'orth', 'orthogonal_procrustes', 'pasca l', 'pinv', 'pinv2', 'pinvh', 'polar', 'print_function', 'qr', 'qr_de lete', 'qr_insert', 'qr_multiply', 'qr_update', 'qz', 'rq', 'rsf2cs f', 'schur', 'signm', 'sinhm', 'sinm', 'solve', 'solve_banded', 'solv e_circulant', 'solve_continuous_are', 'solve_continuous_lyapunov', 's olve_discrete_are', 'solve_discrete_lyapunov', 'solve_sylvester', 'so lve_toeplitz', 'solve_triangular', 'solveh_banded', 'special_matrice s', 'sqrtm', 'subspace_angles', 'svd', 'svdvals', 'tanhm', 'tanm', 't est', 'toeplitz', 'tri', 'tril', 'triu']

print(dir(la))

(6)

In [8]:

In [9]:

Tworzenie macierzy w SciPy bardziej przypomina składnię Matlaba:

In [10]:

Funkcje specjalne

Pełną listę funkcji specjalnych można znaleźć pod adresem

http://docs.scipy.org/doc/scipy/reference/special.html#module-scipy.special (http://docs.scipy.org/doc/scipy/reference/special.html#module-scipy.special).

In [11]:

/usr/local/lib/python3.6/dist-packages/numpy/lib/scimath.py:262: Runt imeWarning: divide by zero encountered in log

return nx.log(x) Out[8]:

array([ 0.+0.j , -inf+0.j , 0.+3.14159265j])

Out[9]:

array([ 0., -inf, nan])

/usr/local/lib/python3.6/dist-packages/ipykernel_launcher.py:2: Runti meWarning: divide by zero encountered in log

/usr/local/lib/python3.6/dist-packages/ipykernel_launcher.py:2: Runti meWarning: invalid value encountered in log

Out[10]:

matrix([[1, 3, 5], [2, 5, 1], [2, 3, 8]])

#wersja scipy

sp.log(np.array((1.0,0.0,-1.0)))

#wersja numpy

np.log(np.array((1.0,0.0,-1.0)))

M = sp.mat('[1 3 5; 2 5 1; 2 3 8]') M

# jn i yn - funkcje Bessela pierwszego i drugiego rodzaju

# jn_zeros i yn_zeros - miejsca zerowe powyższych funkcji from scipy.special import jn, yn, jn_zeros, yn_zeros

(7)

In [13]:

In [14]:

Wielomiany

J_0(0.000000) = 1.000000 Y_0(1.000000) = 0.088257

Out[14]:

array([ 2.40482556, 5.52007811, 8.65372791, 11.79153444]) n = 0 # rząd

x = 0.0

# funkcja Bessela pierwszego rodzaju print("J_%d(%f) = %f" % (n, x, jn(n, x))) x = 1.0

# funkcja Bessela drugiego rodzaju

print("Y_%d(%f) = %f" % (n, x, yn(n, x)))

x = np.linspace(0, 10, 100) fig, ax = plt.subplots() for n in range(4):

ax.plot(x, jn(n, x), label=r"$J_%d(x)$" % n) ax.legend();

# miejsca zerowe funkcji Bessela n = 0 # rząd

m = 4 # liczba miejsc zerowych jn_zeros(n, m)

(8)

In [15]:

In [16]:

In [17]:

In [18]:

In [19]:

In [20]:

In [21]:

2

3 x + 4 x + 5

4 3 2

9 x + 24 x + 46 x + 40 x + 25

3 2

1 x + 2 x + 5 x + 6

6 x + 4

Out[19]:

12

Out[20]:

array([ 5, 12, 25])

Out[21]:

array([-0.66666667+1.1055416j, -0.66666667-1.1055416j])

#wielomian i jego reprezentacja p = sp.poly1d([3,4,5])

print(p)

#mnożenie wielomianów print(p*p)

#całkowanie wielomianów

#stała całkowania = 6 print(p.integ(k=6))

#różniczkowanie wielomianów print(p.deriv())

#wartość wielomianu w punkcie...

p(1)

#i w kilku punktach naraz p([0,1,2])

#pierwiastki wielomianu p.r

(9)

In [23]:

Całkowanie numeryczne

Do numerycznego wyliczenia całki

mamy do dyspozycji kilka róznych funkcji, m.in. quad, dblquad i tplquad do wyznaczania całek pojedynczych, podwójnych i potrójnych.

In [24]:

Obliczmy na początek całkę:

In [25]:

In [26]:

Drugi przykład, czyli całka

f (x)dx

b a

xdx =

1 0

1 2

Out[22]:

array([ 0. +8.88178420e-16j, 0. -8.88178420e-16j])

Out[23]:

array([3, 4, 5])

Wartość całki = 0.5 , błąd = 5.551115123125783e-15

#sprawdzenie p(p.r)

#tablica współczynników p.c

from scipy.integrate import quad, dblquad, tplquad

# definiujemy funkcję podcałkową def f(x):

return x

x_lower = 0 # dolna granica całkowania x_upper = 1 # górna granica całkowania val, abserr = quad(f, x_lower, x_upper)

print("Wartość całki =", val, ", błąd =", abserr)

(10)

ilustruje użycie polecenia quad w przypadku nieskończonych granic całkowania:

In [27]:

In [28]:

In [29]:

W przypadku tak prostych funkcji, jak te powyżej, zamiast definiować funkcję podcałkową, możemy użyć funkcji nienazwanej lambda:

In [30]:

Rozważmy teraz funkcję podcałkową:

In [31]:

Funkcja ta, oprócz zmiennej niezależnej x oczekuje również podania parametru n określającego rząd funkcji Bessela. Parametr ten możemy przekazać do funkcji podcałkowej przy wywołaniu polecenia quad:

dx =

−∞

e

−x2

π⎯⎯

Wartość całki = 1.7724538509055159 , błąd = 1.4202636780944923e-08

0.0

wynik numeryczny = 1.7724538509055159 błąd = 1.4202636780944923e-08 wynik analityczny = 1.77245385091

def fun(x):

return np.exp(-x**2)

val, abserr = quad(fun,-np.Inf,np.Inf)

print("Wartość całki =", val, ", błąd =", abserr)

#sprawdźmy wynik

print(val-sp.sqrt(sp.pi))

val, abserr = quad(lambda x: np.exp(-x ** 2), -np.Inf, np.Inf) print("wynik numeryczny =", val, "błąd =",abserr)

analytical = sp.sqrt(sp.pi)

print("wynik analityczny =", analytical)

def integrand(x, n):

"""

Funkcja Bessela pierwszego rodzaju o rzędzie n """

return jn(n, x)

(11)

Rozważmy teraz całkę podwójną:

Do jej obliczenia służy polecenie dblquad:

In [33]:

Wywołanie dblquad jest bardziej skomplikowane niż quad, ponieważ w przypadku ogólnym granice całkowania po y mogą być funkcjami zmiennej x. Dlatego granice całkowania po y definiujemy za pomocą funkcji użytkownika lub funkcji anonimowych, jak w powyższym przykładzie.

Równania różniczkowe zwyczajne (RRZ)

SciPy pozwala rozwiązywać równania różniczkowe zwyczajne na dwa sposoby: albo przy pomocy funkcji odeint, albo przy pomocy klasy ode. Ten pierwszy sposób jest prostszy (szybciej prowadzi do rozwiązania problemu), ten drugi oferuje większą kontrolę nad rozwiązywanym zagadnieniem.

In [34]:

Przed rozwiązaniem układu równań różniczkowych zwyczajnych należy sprowadzić go do postaci

I =dxdy = π

−∞

−∞

e

− −x2 y2

Wynik całkowania funkcji Bessela rzędu 3 0.7366751370811073 9.389126882496403e-13

3.141592653589777 +/- 2.5173086737433208e-08 x_lower = 0 # dolna granica całkowania

x_upper = 10 # górna granica całkowania n = 3 #rząd funckji Bessela

val, abserr = quad(integrand, x_lower, x_upper, args=(n,)) print("Wynik całkowania funkcji Bessela rzędu", n)

print(val, abserr)

def integrand(x, y):

return np.exp(-x**2-y**2) x_lower = -np.Inf

x_upper = np.Inf y_lower = -np.Inf y_upper = np.Inf

#val, abserr = dblquad(integrand, x_lower, x_upper, lambda x : y_lower, lambda x: y val, abserr = dblquad(integrand, x_lower, x_upper, lambda x : y_lower, lambda x: y_

print(val,"+/-", abserr)

from scipy.integrate import odeint

(12)

gdzie

a jest funkcją, która pozwala wyliczyć pochodne funkcji . Aby rozwiązać taki układ, należy znać funkcję oraz warunek początkowy .

Warto zauważyć, że równania wyższych rzędów również można sprowadzić do układu równań pierwszego rzędu. Wystarczy wprowadzić nowe zmienne dla każdej pochodnej i potraktować ich definicje jako kolejne równania.

Przykład: wahadło

Rozważmy prosty przykład z fizyki (więcej szczegółów tutaj:

http://en.wikipedia.org/wiki/Pendulum_(mathematics (http://en.wikipedia.org/wiki/Pendulum_(mathematics)):

In [35]:

Równanie ruchu wahadła (przy odpowiednio dobranych i ) ma postać:

Przekształcamy najpierw to równanie do układu równań różniczkowych pierwszego rzędu:

Chcemy rozwiązać to równanie dla warunku początkowego i :

= f (y, t) y

y = [ (t), (t), . . . , (t)] y

1

y

2

y

n

f y

i

(t) f

y(0)

l g + sin θ = 0

θ¨

= v dt

= −sinθ dv

dt

θ = 0, 3 v = 0

Out[35]:

from IPython.display import Image

Image(filename='Pendulum.jpg', width=150)

(13)

Przykład: wahadło podwójne

Kolejny, trochę bardziej skomplikowany przykład to wahadło podwójne

(http://en.wikipedia.org/wiki/Double_pendulum (http://en.wikipedia.org/wiki/Double_pendulum)):

Out[36]:

<matplotlib.legend.Legend at 0x7fa53c943f28>

#wartości początkowe th0 = 0.3

v0 = 0.0

init = [th0,v0]

#prawa strona układu równań: f(theta,v) def rhs(y,t):

return [y[1],-np.sin(y[0])]

#czas ruchu

t = np.linspace(0, 10, 250)

#rozwiązanie równania y = odeint(rhs,init,t)

plt.plot(t,y[:,0],label="theta") plt.plot(t,y[:,1],label='v') plt.legend(loc="upper right")

(14)

In [37]:

Równania ruchu wahadła mają następującą postać:

Wprowadźmy zmienną wektorową . Wówczas:

= θ˙

1 6

mℓ2

2pθ1−3 cos( − )θ1 θ2 pθ2 16−9cos2( − )θ1 θ2

= .

θ˙

2 6 mℓ2

8pθ2−3 cos( − )θ1 θ2 pθ1 16−9cos2( − )θ1 θ2

= − m [ sin( − ) + 3 sin ]

θ1 12

2

θ˙

1

θ˙

2

θ

1

θ

2 g

θ

1

= − m [− sin( − ) + sin ]

θ2 12

2

θ˙

1

θ˙

2

θ

1

θ

2 g

θ

2

x = [ , , θ

1

θ

2

p

θ1

, p

θ2

]

=

1 6

mℓ2

2 −3 cos( − )x3 x1 x2 x4 16−9cos2( − )x1 x2

=

2 6

mℓ2

8 −3 cos( − )x4 x1 x2 x3

16−9cos2( − )x1 x2

= − m [ sin( − ) + 3 sin ]

3 1

2

2

1

2

x

1

x

2 g

x

1

= − m [− sin( − ) + sin ]

4 1

2

2

1

2

x

1

x

2 g

x

2 Out[37]:

(0, 0)

mg (x1,y1)

mg (x2,y2) θ1

θ2

Image(url='http://upload.wikimedia.org/wikipedia/commons/c/c9/Double-compound-pendu

(15)

In [39]:

In [40]:

In [41]:

In [42]:

Out[42]:

<matplotlib.legend.Legend at 0x7fa53c924358>

g = 9.82 L = 0.5 m = 0.1

def dx(x, t):

"""

prawa strona układu RRZ dla wahadła podwójnego """

x1, x2, x3, x4 = x[0], x[1], x[2], x[3]

dx1 = 6.0/(m*L**2) * (2 * x3 - 3 * np.cos(x1-x2) * x4)/(16 - 9 * np.cos(x1-x2)*

dx2 = 6.0/(m*L**2) * (8 * x4 - 3 * np.cos(x1-x2) * x3)/(16 - 9 * np.cos(x1-x2)*

dx3 = -0.5 * m * L**2 * ( dx1 * dx2 * np.sin(x1-x2) + 3 * (g/L) * np.sin(x1)) dx4 = -0.5 * m * L**2 * (-dx1 * dx2 * np.sin(x1-x2) + (g/L) * np.sin(x2))

return [dx1, dx2, dx3, dx4]

# stan początkowy

x0 = [np.pi/4, np.pi/2, 0, 0]

# time coodinate to solve the ODE for: from 0 to 10 seconds t = np.linspace(0, 10, 250)

# rozwiązanie układu x = odeint(dx, x0, t)

# kąty w funkcji czasu

plt.plot(t, x[:, 0], 'r', label="theta1") plt.plot(t, x[:, 1], 'b', label="theta2") plt.legend(loc="upper right")

(16)

In [43]:

A teraz prosta animacja ruchu wahadła:

In [44]:

Out[43]:

<matplotlib.legend.Legend at 0x7fa53c890ef0>

# wykres w przestrzeni położeń x1 = + L * np.sin(x[:, 0]) y1 = - L * np.cos(x[:, 0]) x2 = x1 + L * np.sin(x[:, 1]) y2 = y1 - L * np.cos(x[:, 1])

plt.plot(x1, y1, 'r', label="pendulum1") plt.plot(x2, y2, 'b', label="pendulum2") plt.legend(loc="upper right")

from IPython.display import clear_output, display import time

(17)

Przykład: wahadło tłumione

Rozważmy jeszcze jeden przykład - tłumiony oscylator harmoniczny (http://en.wikipedia.org/wiki/Damping (http://en.wikipedia.org/wiki/Damping)). Równanie ruchu ma postać:

gdzie to położenie oscylatora, to częstość kołowa oscylatora nietłumionego, a to współczynnik tłumienia. Niech . Wówczas:

Ten przykład rozwiążemy, traktując współczynnik tłumienia jako parametr przekazywany do funkcji definiującej układ w momencie wywołania polecenia odeint:

+ 2ζ + x = 0

x d

2

dt

2

ω

0

dx

dt ω

20

x ω

0

ζ

p =

dxdt

= −2ζ p − x dp

dt ω

0

ω

20

dx = p

dt

fig, ax = plt.subplots(figsize=(4,4)) plt.show()

for t_idx, tt in enumerate(t[:200]):

x1 = + L * np.sin(x[t_idx, 0]) y1 = - L * np.cos(x[t_idx, 0]) x2 = x1 + L * np.sin(x[t_idx, 1]) y2 = y1 - L * np.cos(x[t_idx, 1])

ax.cla()

ax.plot([0, x1], [0, y1], 'r.-') ax.plot([x1, x2], [y1, y2], 'b.-') ax.set_ylim([-1.5, 0.5])

ax.set_xlim([1, -1]) display(fig)

clear_output(wait=True)

(18)

In [46]:

In [47]:

In [48]:

In [49]:

In [50]:

def dy(y, t, zeta, w0):

"""

Prawa strona układu RRZ oscylatora tłumionego """

x, p = y[0], y[1]

dx = p

dp = -2 * zeta * w0 * p - w0**2 * x return [dx, dp]

# stan początkowy y0 = [1.0, 0.0]

# czas ruchu

t = np.linspace(0, 10, 1000)

#częstość kołowa oscylatora nietłumionego w0 = 2*np.pi*1.0

# rozwiązujemy układ dla różnych współczynników tłumienia y1 = odeint(dy, y0, t, args=(0.0, w0)) # nietłumiony y2 = odeint(dy, y0, t, args=(0.2, w0)) # tłumiony słabo y3 = odeint(dy, y0, t, args=(1.0, w0)) # krytycznie tłumiony y4 = odeint(dy, y0, t, args=(5.0, w0)) # silnie tłumiony

#i oczywiście wykres

plt.plot(t, y1[:,0], 'k', label=u"brak tłumienia", linewidth=0.25) plt.plot(t, y2[:,0], 'r', label=u"słabe tłumienie")

plt.plot(t, y3[:,0], 'b', label=u"krytyczne tłumienie") plt.plot(t, y4[:,0], 'g', label=u"silne tłumienie") plt.legend();

(19)

Do obliczania transformat Fouriera SciPy wykorzystuje klasyczną, zaimplementowaną w Fortranie bibliotekę FFTPACK (http://www.netlib.org/fftpack/).

In [51]:

Policzmy transformatę Fouriera dla rozwiązania równania ruchu tłumionego oscylatora harmonicznego:

In [52]:

In [53]:

Ponieważ sygnał jest rzeczywisty, spektrum jest symetryczne, dlatego wystarczy, jeśli na wykresie przedstawimy tylko część odpowiadającą dodatnim częstotliwościom:

In [54]:

from scipy.fftpack import *

N = len(t) dt = t[1]-t[0]

# transformata rozwiązania oscylatora słabo tłumionego F = fft(y2[:,0])

# częstotliwości odpowiadające składowym transformaty w = fftfreq(N, dt)

plt.plot(w, abs(F));

indices = np.where(w > 0) # można inaczej, ale przy okazji ćwiczymy numpy w_pos = w[indices]

F_pos = F[indices]

(20)

In [55]:

Algebra liniowa

In [56]:

Układy równań liniowych Układ równań liniowych w postaci

gdzie to macierz układu, to wektor wyrazów wolnych, natomiast jest szukanym wektorem, rozwiążemy w następujący sposób:

In [57]:

In [58]:

Ax = b

A b x

Out[58]:

array([-0.33333333, 0.66666667, -0. ]) plt.plot(w_pos, abs(F_pos))

plt.xlim(0, 5);

from scipy import linalg

A = np.array([[1,2,310], [4,5,6], [7,8,9]]) b = np.array([1,2,3])

x = sp.linalg.solve(A, b) x

(21)

W podobny sposób możemy rozwiązać układ uogólniony

gdzie są macierzami. Układ taki odpowiada wielu układom równań liniowych, które mają tę samą macierz, natomiast różnią się wektorami wyrazów wolnych.

In [60]:

In [61]:

In [62]:

In [63]:

Nadokreślone układy równań liniowych

Rozważmy układ równań liniowych, w którym równań jest więcej niż zmiennych, np.:

AX = B A, B, X

x + y = 1.98 2.05x − y = 0.95 3.06x + y = 3.98

−1.02x + 2y = 0.92 4.08x − y = 2.90

Out[59]:

array([ 0.00000000e+00, -2.22044605e-16, 0.00000000e+00])

Out[62]:

array([[ 0.63775679, -0.66946063, -1.75448344], [ 0.50854862, 2.19176126, 1.95184606], [-0.54314485, -1.21683073, 0.24381299]])

Out[63]:

2.4047046924773687e-16

# sprawdzenie sp.dot(A, x) - b

A = sp.rand(3,3) B = sp.rand(3,3)

X = sp.linalg.solve(A, B)

X

# sprawdzenie

sp.linalg.norm(sp.dot(A, X) - B)

(22)

Układ taki możemy rozwiązać albo poprzez macierz pseudoodwrotną:

In [64]:

lub bezpośrednio metodą najmniejszych kwadratów:

In [65]:

Niedookreślone układy równań liniowych

In [66]:

Wartości i wektory własne

Zagadnienie na wartości własne macierzy ma postać:

gdzie jest -tym wektorem własnym, a - -tą wartością własną.

Jeśli chcemy obliczyć tylko wartości własne, mamy do dyspozycji funkcję eigvals. Do wyliczenia zarówno wartości jak i wektorów własnych służy funkcja eig:

In [67]:

x + 2y = 3

A A = v

n

λ

n

v

n

v

n

n λ

n

n

Out[64]:

matrix([[ 0.9631014 ], [ 0.98854334]])

[[ 0.9631014 ] [ 0.98854334]]

Out[66]:

matrix([[ 0.6], [ 1.2]])

c = sp.mat('[1.98; 0.95; 3.98; 0.92;2.90]')

B = sp.mat('[1 1;2.05 -1;3.06 1; -1.02 2;4.08 -1]') BI = sp.linalg.pinv2(B)

BI*c

x, res, rank, s = sp.linalg.lstsq(B,c) print(x)

b = sp.mat('[3]') A = sp.mat('[1 2]') AI = sp.linalg.pinv2(A) AI*b

A = sp.rand(3,3)

evals = sp.linalg.eigvals(A)

(23)

In [69]:

In [70]:

In [71]:

Wektor własny odpowiadający -tej wartości własnej stanowi -tą kolumnę w macierzy .

In [72]:

Operacje na macierzach In [73]:

n n evecs

Out[68]:

array([ 1.69982529+0.j , 0.03473551+0.13277785j, 0.03473551-0.13277785j])

Out[70]:

array([ 1.69982529+0.j , 0.03473551+0.13277785j, 0.03473551-0.13277785j])

Out[71]:

array([[ 0.12124268+0.j , 0.16207709-0.17778954j, 0.16207709+0.17778954j],

[ 0.72527626+0.j , 0.59666244+0.1042882j , 0.59666244-0.1042882j ],

[ 0.67769798+0.j , -0.75844565+0.j , -0.75844565 -0.j ]])

Out[72]:

5.840146961801262e-16

Out[73]:

array([[ 7.03351142, -2.10803848, 1.10296187], [ 20.94089876, -2.78521232, -0.13605678], [-27.10572388, 5.05464455, 0.02810072]]) evals

evals, evecs = sp.linalg.eig(A)

evals

evecs

#sprawdzenie n = 1

sp.linalg.norm(sp.dot(A, evecs[:,n]) - evals[n] * evecs[:,n])

# macierz odwrotna sp.linalg.inv(A)

(24)

In [74]:

In [75]:

In [76]:

In [77]:

Optymalizacja

... czyli znajdowanie miejsc zerowych, maksimów bądź minimów funkcji

In [46]:

Znajdowanie minimów:

In [47]:

Out[75]:

matrix([[ 7.03351142, -2.10803848, 1.10296187], [ 20.94089876, -2.78521232, -0.13605678], [-27.10572388, 5.05464455, 0.02810072]])

Out[76]:

0.03201878222125203

Out[77]:

(1.7927978073770681, 2.449309091933058) AM = sp.matrix(A)

AM.I

# wyznacznik sp.linalg.det(A)

# normy różnych typów

sp.linalg.norm(A, ord=2), sp.linalg.norm(A, ord=sp.Inf)

from scipy import optimize

#funkcja, której minimum będziemy szukac def f(x):

return 4*x**3 + (x-2)**2 + x**4

(25)

Do znalezienia minimów tej funkcji możemy użyć polecenia fmin_bfgs:

In [49]:

In [50]:

Sprawdźmy, jakiego algorytmu używa fmin_bfgs:

Optimization terminated successfully.

Current function value: -3.506641 Iterations: 5

Function evaluations: 24 Gradient evaluations: 8 Out[49]:

array([-2.67298151])

Optimization terminated successfully.

Current function value: 2.804988 Iterations: 3

Function evaluations: 15 Gradient evaluations: 5 Out[50]:

array([ 0.46961745])

#poglądowy wykres

x = np.linspace(-5, 3, 100) plt.plot(x, f(x));

x_min = optimize.fmin_bfgs(f, -2) x_min

optimize.fmin_bfgs(f, 0.5)

(26)

In [52]:

fmin_bfgs(f, x0, fprime=None, args=(), gtol=1e-05, norm=inf,

epsilon=1.4901161193847656e-08, maxiter=None, full_output=

0,

disp=1, retall=0, callback=None) Minimize a function using the BFGS algorithm.

Parameters ---

f : callable f(x,*args)

Objective function to be minimized.

x0 : ndarray

Initial guess.

fprime : callable f'(x,*args), optional Gradient of f.

args : tuple, optional

Extra arguments passed to f and fprime.

gtol : float, optional

Gradient norm must be less than gtol before successful terminatio n.

norm : float, optional

Order of norm (Inf is max, -Inf is min) epsilon : int or ndarray, optional

If fprime is approximated, use this value for the step size.

callback : callable, optional

An optional user-supplied function to call after each iteration. Called as callback(xk), where xk is the current parameter vector.

maxiter : int, optional

Maximum number of iterations to perform.

full_output : bool, optional

If True,return fopt, func_calls, grad_calls, and warnflag in addition to xopt.

disp : bool, optional

Print convergence message if True.

retall : bool, optional

Return a list of results at each iteration if True.

Returns ---

xopt : ndarray

Parameters which minimize f, i.e. f(xopt) == fopt.

fopt : float

Minimum value.

gopt : ndarray

Value of gradient at minimum, f'(xopt), which should be near 0.

Bopt : ndarray

Value of 1/f''(xopt), i.e. the inverse hessian matrix.

func_calls : int

Number of function_calls made.

grad_calls : int

Number of gradient calls made.

warnflag : integer

1 : Maximum number of iterations exceeded.

2 : Gradient and/or function calls not changing.

sp.info(optimize.fmin_bfgs)

(27)

Możemy również skorzystać z funkcji brent lub fminbound, które używają innych algorytmów:

In [84]:

In [85]:

Miejsca zerowe funkcji

Do znajdowania miejsc zerowych funkcji, tzn. rozwiązania zagadnienia , służy polecenie fsolve:

In [86]:

f (x) = 0

See also ---

minimize: Interface to minimization algorithms for multivariate functions. See the 'BFGS' `method` in particular.

Notes ---

Optimize the function, f, whose gradient is given by fprime using the quasi-Newton method of Broyden, Fletcher, Goldfarb, and Shanno (BFGS)

References ---

Wright, and Nocedal 'Numerical Optimization', 1999, pg. 198.

Out[84]:

0.46961743402759754

Out[85]:

-2.6729822917513886 optimize.brent(f)

optimize.fminbound(f, -4, 2)

def f(x):

return sp.sin(x)

(28)

In [87]:

In [88]:

In [89]:

In [90]:

Układy równań nieliniowych

−2 + 3xy + 4 sin y = 6 x

2

3 − 3x + 3 cos x = −4 x

2

y

2

Out[88]:

array([ 0.])

Out[89]:

array([ 3.14159265])

Out[90]:

array([ 6.28318531])

x = np.linspace(0, 10, 1000) y = f(x)

plt.plot(x, y)

plt.plot([0, 10], [0, 0], 'k') plt.ylim(-1.5,1.5);

optimize.fsolve(f, 0.1)

optimize.fsolve(f, 2.0)

optimize.fsolve(f, 5.0)

(29)

In [92]:

In [93]:

In [94]:

Interpolacja

In [95]:

In [96]:

In [97]:

Out[93]:

array([ 0.57982909, 2.54620921])

Out[94]:

[-3.602984577355528e-11, 9.8171692997084392e-11]

#funkcja definiująca układ def func(x):

out = [-2*x[0]**2 + 3*x[0]*x[1] + 4*sp.sin(x[1]) - 6]

out.append(3*x[0]**2 - 2*x[0]*x[1]**2 + 3*sp.cos(x[0]) + 4) return out

#przybliżenie początkowe x0 = [1.0,1.0]

#rozwiązanie

optimize.fsolve(func,x0)

#sprawdzenie func(_)

from scipy.interpolate import *

def f(x):

return np.sin(x)

n = np.arange(0, 10) #węzły interpolacji x = np.linspace(0, 9, 100)

y_meas = f(n) + 0.1 * np.random.randn(len(n)) # symulacja pomiarów z szumem y_real = f(x)

#interpolacja przedziałami liniowa

linear_interpolation = interp1d(n, y_meas) y_interp1 = linear_interpolation(x)

#interpolacja funkcjami sklejanymi 3 rzędu

cubic_interpolation = interp1d(n, y_meas, kind='cubic') y_interp2 = cubic_interpolation(x)

(30)

In [98]:

Statystyka

Moduł scipy.stats zawiera wiele przydatnych rozkładów, funkcji i testów statystycznych.

In [99]:

In [100]:

Out[100]:

<scipy.stats._distn_infrastructure.rv_frozen at 0x7fa83a15ea90>

#i poglądowy wykres

plt.plot(n, y_meas, 'bs', label='dane zaszumione') plt.plot(x, y_real, 'k', lw=2, label='f(x)')

plt.plot(x, y_interp1, 'r', label='interp. liniowa') plt.plot(x, y_interp2, 'g', label='funkcje sklejane') plt.legend(loc="upper center");

from scipy import stats

# dyskretna zmienna losowa z rozkładem Poissona X = stats.poisson(3.5)

X

(31)

In [102]:

n = np.arange(0,15)

fig, axes = plt.subplots(3,1, sharex=True)

# funkcja masy prawdopodobieństwa (PMF) axes[0].step(n, X.pmf(n))

# dystrybuanta (CDF) axes[1].step(n, X.cdf(n))

# histogram 1000 losowych realizacji zmiennej X axes[2].hist(X.rvs(size=1000));

# zmienna losowa ciągła o rozkładzie normalnym Y = stats.norm()

(32)

In [103]:

Statystyki:

In [104]:

In [105]:

Testy statystyczne

Chcemy sprawdzić, czy dwa zbiory niezależnych danych losowych są z tego samego rozkładu:

Out[104]:

(3.5, 1.8708286933869707, 3.5)

Out[105]:

(0.0, 1.0, 1.0)

x = np.linspace(-5,5,100)

fig, axes = plt.subplots(3,1, sharex=True)

# funkcja rozkładu prawdopodobieństwa (PDF) axes[0].plot(x, Y.pdf(x))

# dystrybuanta (CDF)

axes[1].plot(x, Y.cdf(x));

# histogram 1000 realizacji zmiennej losowej Y axes[2].hist(Y.rvs(size=1000), bins=50);

X.mean(), X.std(), X.var() # rozkład Poissona

Y.mean(), Y.std(), Y.var() # rozkład normalny

(33)

Ponieważ wartość jest bardzo duża, nie możemy odrzucić hipotezy, że zbory danych mają różne średnie.

Sprawdzamy, czy zbiór danych ma średnią 0.1 (prawdziwa średnia wynosi 0.0):

In [107]:

Mała wartośc pozwala odrzucić hipotezę, że ma średnią 0.1.

In [108]:

In [114]:

In [ ]:

p

p Y

t-statistic = -1.31524994103 p-value = 0.188576833898

Out[107]:

Ttest_1sampResult(statistic=-4.6326468776347074, pvalue=4.08715190886 94328e-06)

Out[108]:

0.0

Out[114]:

Ttest_1sampResult(statistic=0.60081835576404308, pvalue=0.54809726240 705092)

t_statistic, p_value = stats.ttest_ind(X.rvs(size=1000), X.rvs(size=1000)) print("t-statistic =", t_statistic)

print("p-value =", p_value)

stats.ttest_1samp(Y.rvs(size=1000), 0.1)

Y.mean()

stats.ttest_1samp(Y.rvs(size=1000), Y.mean())

Obraz

Updating...

Cytaty

Powiązane tematy :