Tekst źródłowy strony TI/Wybrane zagadnienia numeryczne

==Rozwiązywanie równań liniowych==
; Przykład<nowiki>:</nowiki>

Rozwiązanie układu równań liniowych:
::<math>
\left\{ \begin{array}{l}
x_1+3x_3=3\\
2x_1+x_2+5x_3=2\\
4x_1+8x_2+8x_3=1
\end{array} \right.
</math><br>
<source lang= python>
import numpy as np

A = np.array([[1, 0, 3],
              [2, 1, 5],
              [4, 8, 8]])
b = np.array([3, 2, 1])

x = np.linalg.solve(A, b)
print x
</source>

; Rezultat

<source lang='shell_example'>
[-12.75   1.25   5.25]
</source>

; Jak to działa?

Nasz wyjściowy układ równań liniowych możemy zapisać w postaci macierzowej <math>Ax=b</math>:<br><br>
::<math>
\begin{array}{cccc}
\left( \begin{array}{ccc}1 & 0 & 3\\2 & 1 & 5\\4 & 8 & 8\end{array} \right)&
\left( \begin{array}{c}x_1\\x_2\\x_3\end{array} \right)&
=&
\left( \begin{array}{c}3\\2\\1\end{array} \right)\\
\ &\ &\ &\ \\
A&x&\ &b \end{array}
</math>. <br><br>
Macierz <math>A</math> konstruujemy ze współczynników mnożących zmienne <math>x_1</math>, <math>x_2</math> i <math>x_3</math> w poszczególnych równaniach &mdash; każdy wiersz macierzy <math>A</math> opisuje jedno równanie, każda jej kolumna opisuje kolejno jedną z niewiadomych zmiennych <math>x_i</math>. Po pomnożeniu <math>A</math> przez wektor <math>x</math>, którego każdy wiersz zawiera odpowiednią niewiadomą, otrzymamy lewą stronę naszego układu równań. Prawą stronę tworzy wektor <math>b</math>, którego wiersze utworzone są odpowiednio z prawych stron poszczególnych równań układu.

Funkcja <tt>solve</tt> modułu Numpy zwraca nam listę zawierającą poszukiwane przez nas rozwiązanie układu równań: ''x''<sub>1</sub> = &minus;12,75, ''x''<sub>2</sub> = 1,25, ''x''<sub>3</sub> = 5,25.

; Inne podejście

Możemy zauważyć, że rozwiązanie układu równań można formalnie uzyskać mnożąc obie strony równania (lewostronnie) przez macierz odwrotną do ''A'':<br>
::<math>
\begin{array}{l}
Ax=b\quad\quad\quad\quad\quad\quad \mid A^{-1}\cdot \\
\underbrace{A^{-1}\cdot A}_{=1} \cdot x=A^{-1} \cdot b \\
x=A^{-1}b
\end{array}
</math>
<br><br>
W poniższym przykładzie do odwrócenia macierzy ''A'' stosujemy funkcję <tt>inv</tt> z modułu <tt>numpy.linalg</tt>. Mnożenie macierzy realizuje zaś funkcja <tt>dot</tt>.

<source lang= python>
# -*- coding: utf-8 -*-
import numpy as np

A = np.array([[1, 0, 3],
              [2, 1, 5],
              [4, 8, 8]])
b = np.array([3, 2, 1])

A1 = np.linalg.inv(A)
print "Macierz odwrotna do A:"
print A1
print "Sprawdzenie macierzy odwrotnej:"
print np.dot(A1, A)
x = np.dot(A1, b)
print "Rozwiązanie:"
print x
</source>

; Rezultat

<source lang='shell_example'>
Macierz odwrotna do A:
[[-8.    6.   -0.75]
 [ 1.   -1.    0.25]
 [ 3.   -2.    0.25]]
Sprawdzenie macierzy odwrotnej:
[[  1.00000000e+00   1.77635684e-15   1.77635684e-15]
 [ -1.11022302e-16   1.00000000e+00  -2.22044605e-16]
 [ -2.22044605e-16  -4.44089210e-16   1.00000000e+00]]
Rozwiązanie:
[-12.75   1.25   5.25]
</source>

<br>

==Wyrażenia lambda==
W takich obliczeniach jak te opisywane w niniejszym rozdziale, występuje potrzeba definiowania króciusieńkich funkcji.
Okazuje się, że można je definiować na dwa sposoby:
# zwyczajnie, z wykorzystaniem słowa kluczowego <tt>def</tt>,
# jako tzw. ''wyrażenie &lambda;'', z wykorzystaniem słowa kluczowego <tt>lambda</tt>.

Po słowie kluczowym <tt>lambda</tt> podajemy nazwy, które będą pełnić rolę zmiennych w naszej funkcji, a po dwukropku wypisujemy wyrażenie (używając podanych przed dwukropkiem nazw jako zmiennych).
<source lang='shell_example'>
>>> def f(x): return x**3
... 
>>> print f
<function f at 0x1c285f0>
>>> f1 = f
>>> print f1
<function f at 0x1c285f0>

>>> g = lambda x: x**3
>>> print g
<function <lambda> at 0x1c286e0>
</source>
W tym przykładzie, funkcje <tt>f</tt>, <tt>f1</tt> i <tt>g</tt> działają tak samo &mdash; biorą jeden argument i&nbsp;zwracają jego sześcian.
Różnica jest taka, że funkcja <tt>f</tt> &bdquo;wie&rdquo; że nazywa się &bdquo;f&rdquo;, co widać w napisie wypisywanym po <tt>print f</tt> czy <tt>print f1</tt>. Natomiast funkcja <tt>g</tt> nazywa się &bdquo;<lambda>&rdquo;, tak samo jak wszystkie inne funkcje zdefiniowane jako wyrażenie &lambda;, czyli nie ma swojej nazwy. Zmienna <tt>g</tt> jest tylko dowiązaniem do obiektu funkcji, tak samo jak zmienna <tt>f1</tt>.

Funkcja zdefiniowana jako wyrażenie &lambda; jest równoważna funkcji definiowanej zwyczajnie zwracającej to samo wyrażnie.
Niemniej, większości rzeczy które można zrobić w funkcji, nie można zrobić w wyrażeniu &lambda; ze względu na ograniczenia
składni: 
# Wyrażenie &lambda; musi być pojedynczym wyrażeniem.
# Wyrażenie &lambda; nie może zawierać ''docstringa''.

Wyrażenia &lambda; są wygodne ze względu na to, że wymagają mniej stukania w klawisze.
==Rozwiązywanie równań dowolnych==
<source lang= python>
import scipy.optimize as so

def func(x):
    return x**3 + 3. * x - 0.3

print so.fsolve(func, 0.5)
print so.fsolve(lambda x: x**3 + 3*x - 0.3, 0.5)
</source>

; Rezultat<nowiki>:</nowiki>

<source lang='shell_example'>
0.099669956223525716
0.099669956223525716
</source>

; Jak to działa?

Funkcja <tt>fsolve</tt> z modułu <tt>scipy.optimize</tt> poszukuje miejsc zerowych funkcji, czyli takich wartości ''x'', dla których <math>f(x)=0</math>. Załóżmy, że chcemy rozwiązać równanie:<br>
<math>x^3+3x=0,3</math>.<br>
Po przepisaniu tego równania do postaci<br>
<math>x^3+3x-0,3=0</math><br>
widzimy, że rozwiązanie tego równania jest równoważne poszukiwaniu miejsc zerowych funkcji:<br>
<math>f(x)=x^3+3x-0,3</math>.<br>
Funkcja <tt>fsolve</tt> oczekuje badanej funkcji jako pierwszego argumentu. Możemy tu przekazać nazwę funkcji osobno zdefiniowanej, jak na przykład <tt>func</tt>. Jeśli chcemy szukać miejsc zerowych funkcji, którą da się zapisać prostym wyrażeniem złożonym z operacji na dostępnych wielkościach i funkcjach możemy użyć jako argumentu ''wyrażenia lambda''. Wyrażenie to zastępuje &bdquo;pełną&rdquo; definicję funkcji w momencie, gdy tak naprawdę taka pełna definicja nie jest nam potrzebna. W podanym przykładzie szukaną funkcję przekazaliśmy dwoma sposobami &mdash; raz z użyciem nazwy uprzednio zdefiniowanej funkcji <tt>func</tt>, a raz z użyciem wyrażenia lambda.<br>
Drugim argumentem tej funkcji jest punkt startowy do poszukiwania miejsca zerowego. Powinniśmy najpierw zgrubnie oszacować położenie poszukiwanego miejsca zerowego i wpisać je tutaj. W zależności od funkcji parametr ten jest mniej lub bardziej istotny dla znalezienia poprawnego wyniku. Jest on zaś szczególnie istotny, jeśli badana funkcja ma więcej niż jedno miejsce zerowe (w tym przypadku wskazane może być wykreślenie funkcji).

<br><br>

==Numeryczne obliczanie całek oznaczonych==
; Przykład<nowiki>:</nowiki>

<source lang= python>
# -*- coding: utf-8 -*-
import scipy.integrate as si
import numpy as np

def func(x):
    return x**3/(np.exp(x)-1)

print 'Całka od',0,'do',1,'wynosi',si.quad(func,0,1)
print 'Całka od',0,'do','+∞','wynosi',si.quad(func,0,np.inf)
print 'Całka od',0,'do','π','wynosi',si.quad(lambda x: np.sin(x),0,np.pi)
print 'Całka od',0,'do',1,'wynosi',si.quad(lambda x: x,0,1)
</source>

; Rezultat<nowiki>:</nowiki>

<source lang='shell_example'>
Całka od 0 do 1 wynosi (0.22480518802593821, 2.4958389580158414e-15)
Całka od 0 do +∞ wynosi (6.4939394022668298, 2.6284700289248249e-09)
Całka od 0 do π wynosi (2.0, 2.2204460492503131e-14)
Całka od 0 do 1 wynosi (0.5, 5.5511151231257827e-15)
</source>

; Jak to działa?

Do numerycznego obliczania całek oznaczonych używamy funkcji <tt>quad</tt> z modułu <tt>scipy.integrate</tt>. Funkcja ta jako pierwszego argumentu wywołania oczekuje funkcji, którą będziemy całkować. Podobnie jak w poprzednim przykładzie, możemy tu przekazać nazwę funkcji osobno zdefiniowanej (<tt>func</tt>) albo wyrażenia lambda.<br>
Drugi i trzeci argument funkcji <tt>quad</tt> to początek i koniec zakresu całkowania. Zauważ, że oprócz &bdquo;normalnych&rdquo; liczb możliwe są zakresy nieskończone, z użyciem stałej <tt>inf</tt> z modułu <tt>numpy</tt>. Rezultatem funkcji <tt>quad</tt> jest krotka dwuelementowa. Pierwszy jej element to wynik całkowania, natomiast drugi określa dokładność, z jaką udało się wynik ten policzyć.
<br>
<br>

==Rozwiązywanie równań różniczkowych==
; Przykład<nowiki>:</nowiki>

Rozwiązanie równania różniczkowego <math>\frac{dy(t)}{dt}=y</math> z warunkiem początkowym <math>y(t=0)=1</math>.

<source lang= python>
import scipy.integrate as si
import numpy as np
import pylab as p

df = lambda y, t: y
X = np.linspace(0, 5, 51)        # 51 punktów od 0 do 5 włącznie
wynik = si.odeint(df, 1, X)

p.plot(X, np.zeros_like(X), 'ro',
       X, wynik, 'go',
       X, numpy.exp(X), '--')
p.show()
</source>

; Rezultat

[[Grafika:ODE1.png]]

; Jak to działa?
Poszukujemy funkcji <math>y(t)</math> spełniającej nasze równanie różniczkowe <math>y^\prime=y</math>. Funkcja <tt>odeint</tt> (skrót ODE pochodzi z angielskiego określenia ''ordinary differential equation'' czyli równanie różniczkowe zwyczajne) oblicza wartości poszukiwanej przez nas funkcji w punktach podanych jako jej trzeci argument. 
Z teorii wiemy, że równanie to spełnia funkcja <math>y(t)=\exp(t)</math> (spełnia je również funkcja <math>y(t) \equiv0 </math>, ale ze względu na przyjęty warunek początkowy ta funkcja nie może być naszym rozwiązaniem).
<br>
<br>

== Równania wyższych stopni i układy równań różniczkowych ==
Poniższy przykład ilustruje jak można rozwiązać równanie różniczkowe zwyczajne wyższego rzędu. Jednocześnie będzie to przykład na rozwiązywanie układów równań różniczkowych.

Rozważmy masę <math>m</math> umieszczoną na sprężynie o stałej sprężystości <math>k</math>. Siła działająca na masę zależna jest od jej wychylenia <math>x</math> z położenia równowagi i wynosi &minus;''kx''. Równanie ruchu tej masy to (kropka nad symbolem zmiennej oznacza jej różniczkowanie po czasie: jedna kropka &mdash; pierwszą pochodną, dwie kropki &mdash; drugą pochodną):
:<math>
m\ddot{x} = - k x
</math>

Dzieląc to równanie stronami przez <math>m</math> otrzymujemy standardowe równanie oscylatora harmonicznego o częstości <math>\omega = \sqrt{\frac{k}{m}} </math>:
:<math>
\ddot{x} = - \omega^2 x
</math>

Funkcje całkujące równania różniczkowe w SciPy radzą sobie tylko z równaniami i układami równań pierwszego rzędu. Musimy zatem przepisać nasze równanie na układ równań pierwszego rzędu. Można to zrobić wprowadzając dodatkową zmienną. Tą zmienną jest <math>v = \dot{x}</math> (prędkość masy).
Teraz nasz układ wygląda tak:
:<math>
\left\{ \begin{array}{l}
\dot{x} = v \\
\dot{v} = - \omega^2 x
\end{array}
\right.
</math>

<source lang= python>
# -*- coding: utf-8 -*-

from scipy.integrate import odeint
import numpy as np
import pylab as p

def prawa_strona_rownania(w, t, params):
    ''' Argumenty:
         w: wektor stanu (u nas x, v)
         t: czas
         params: wektor parametrów (u nas omega_kwadrat)

     W wektorze f funkcja zwraca
              obliczone dla danego wektora stanu
              wartości prawej strony równania
    '''
    x, v = w                  # dla czytelności równania wypakowuję zmienne z wektora "w"
    omega_kwadrat, = params   # i podobnie z parametrami "params"
    # poniżej do tablicy f w kolejnych wierszach wpisuję
    # kolejne prawe strony równań stanowiących układ
    f = [v,                   # wartość pochodnej dx/dt
         -omega_kwadrat * x]  # wartość pochodnej dv/dt
    return f

t = np.linspace(0, 5, 51)
params = [2]
w = [0, 3]                    # warunek początkowy (t=0) dla x i v
print "Wektor stanu w chwili początkowej: ",
print prawa_strona_rownania(w, t[0], params)

# argumentami odeint są:
# - nazwa funkcji,
# - wektor stanu początkowego,
# - wektor zawierający chwile czasu, dla których ma być zwrócony stan układu
# - krotka zawierająca dodatkowe parametry, które mają być przekazane do funkcji
#           opisującej prawe strony równań

wynik = odeint(prawa_strona_rownania, w, t, args=(params,) ) 

x = wynik[:, 0]
v = wynik[:, 1]

p.plot(t,x, t,v)
p.legend(('x', 'v'))
p.grid(True)
p.show()
</source>

Funkcja <tt>odeint</tt> oczekuje przynajmniej trzech parametrów: funkcji obliczającej pochodną poszukiwanej funkcji (albo pochodne poszukiwanych funkcji w przypadku układu równań), warunku początkowego dla szukanych funkcji oraz sekwencji punktów czasu, w których będzie obliczone rozwiązanie. W naszym przypadku będą to:
# nazwa funkcji obliczającej pochodne &mdash; <tt>prawa_strona_rownania</tt>, zwraca wektor <math>[\dot{x},\ \dot{v}]</math>, czyli <math>[v,\ -\omega ^2 x]</math>;
# warunek początkowy &mdash; lista <tt>w = [0, 3]</tt>, czyli [''x''(0), ''v''(0)];
# sekwencja punktów czasu &mdash; wektor <tt>t = np.linspace(0,5,51)</tt>, czyli [0, 0,1, 0,2,... 4,9, 5].
U nas jest jeszcze czwarty parametr (nazwany, o nazwie <tt>args</tt>). Powinna być to krotka zawierająca dodatkowe parametry, jakich może potrzebować funkcja obliczająca wartości pochodnych. U nas jest taki jeden dodatkowy parametr: &omega;<sup>2</sup>, która w naszym programie jest przechowywana w zmiennej <tt>params</tt>. Tworzymy więc z niej krotkę jednoelementową i przekazujemy ją funkcji <tt>odeint</tt> jako czwarty parametr. Zostanie ona użyta podczas liczenia prawej strony równań rozwiązywanego układu, czyli w funkcji <tt>prawa_strona_rownania</tt>.


; Rezultat

Otrzymujemy wykres zmiennych ''x'' i ''v'' czyli położenia i prędkości masy zaczepionej na sprężynie. Ponieważ zastosowaliśmy  warunek początkowy ''x''(0) = 0, ''v''(0) = 3, więc w chwili początkowej masa mija punkt równowagi ''x''=0 z prędkością 3 jednostek. Otrzymujemy drgania z amplitudą teoretycznie 3/&omega; = 2,12 jednostek.

[[Grafika:ODE2.png]]

==Transformata Fouriera==
<source lang= python>
import pylab as p
import numpy as np
import scipy.fftpack as sf

x = np.loadtxt('c4spin.txt')

p.subplot(2,1,1)
p.plot(x)
p.subplot(2,1,2)
z=sf.fft(x)
widmo = np.abs(z[:len(z)/2])
p.plot(widmo)
p.show()
</source>
 http://brain.fuw.edu.pl/~jarek/SYGNALY/TF/c4spin.txt

; Rezultat

[[Grafika:FFT.gif]]

===Jak to działa ===
Moduł <tt>scipy.fftpack</tt> dostarcza narzędzi do liczenia transformaty Fouriera.
Transformata Fouriera pozwala wyznaczyć skład częstotliwościowy sygnału.
W powyższym przykładzie wczytujemy dane z pliku tekstowego poleceniem <tt>np.loadtxt('c4spin.txt')</tt>. Dane te interpretujemy jako kolejne próbki pewnego sygnału.

Na górnym panelu rysunku wykreślamy przebieg sygnału w czasie. 
Funkcja <tt>sf.fft(x)</tt> oblicza transformatę Fouriera sygnału <tt>x</tt>. 
Zmienna <tt>z</tt> zawiera pełną informację o transformacie Fouriera sygnału <tt>x</tt> &mdash; jest to wektor liczb zespolonych. Aby uzyskać informację o częstościach zawartych w sygnale <tt>x</tt> trzeba wziąć wartość bezwzględną liczb <tt>z</tt>. Jeśli ponadto sygnał <tt>x</tt> składa się z liczb rzeczywistych (a tak jest w naszym przypadku) to wynikowy wektor <tt>z</tt> zawiera pełną informację o składzie częstotliwościowym w swojej pierwszej połowie. Dlatego do zmiennej <tt>widmo</tt> przepisujemy tylko wycinek <tt>z[:len(z)/2]</tt>.
Na dolnym panelu wyrysowujemy zmienną <tt>widmo</tt>.