GNU Octave

Z Wikipedii

Masz nowe wiadomości (różnica z poprzednią wersją).

GNU Octave z GUI

GNU Octave - środowisko obliczeń oraz język programowania przeznaczony dla prowadzania obliczeń numerycznych.

Octave dostępny jest na większość systemów uniksowych. Rozprowadzany jest na zasadach licencji GNU GPL.

Octave jest wolnym odpowiednikem środowiska MATLAB.

[edytuj] Historia

Prace nad Octave rozpoczeły się w 1988 roku. Pełną parą ruszyły one wiosną 1992 za sprawą Johna W. Eatona. Pierwsza wersja alpha ukazała się 4 stycznia 1993. Wersja 1.0 została wydana 17 lutego 1994. Najnowsza testowa wersja środowiska nosi numer 2.9.9.

W tej chwili Octave rozprowadzany jest w postaci binarnej wraz z wieloma dystrybucjami Linuksa. Dostępna jest także wersja dla Microsoft Windows.

[edytuj] Podstawowe przykłady

Podstawową strukturą danych w Octavie jest macierz.

Wykres funkcji exp(cos(x)) wykonany za pomocą Octave'a.

Wykres funkcji sin(1/x) wykonany za pomocą Octave'a.

octave:2> a=[3 1; 2 1] 
a =

  3  1
  2  1

octave:3> det(a)
ans = 1

octave:4> inv(a)
ans =

   1  -1
  -2   3

octave:5> cond(a)
ans = 14.933

octave:7> [a,a] 
ans =

  3  1  3  1
  2  1  2  1

octave:8> [a;a]
ans =

  3  1
  2  1
  3  1
  2  1

octave:10> [v,e]=eig(a)
v =

   0.80690  -0.34372
   0.59069   0.93907

e =

  3.73205  0.00000
  0.00000  0.26795

Język dostarcza również większośś konstrukcji imperatywnych, m.in. if, while, for, function, np.:

a = [1,2,3];
function ret = kw(x,a)
  ret = a(1)*x.^2+a(2)*x+a(3);
endfunction

Octave używa (zwykle) gnuplota do rysowania wykresów, np exp(cos(x)):

octave:45> x=[-10:0.1:10];
octave:46> y=exp(cos(x))-1;
octave:47> plot(x,y);

lub $s i n (1 / x)$ :

x=[0.001:0.0001:0.5];
y=sin(1./x);
axis([0, 0.5, -1.1, 1.1]);
grid on;
plot(x,y,"-r;sin(1/x);");

[edytuj] Przykładowe zadania z rozwiązaniami

[edytuj] Generowanie danych

Wygenerować wektor długości 10 zawierający losowe liczby rzeczywiste z przedziału [-3, 5]

octave:1> rand(1,10).*8-3
ans =
   1.17900   1.20344   0.58057   -2.59145   3.03283   2.24702   4.99211   1.49550  -1.05363  -1.32788

Wygenerować macierz 3x5 (3 wiersze, 5 kolumn) zawierającą losowe liczby naturalne z przedziału [1, 99]

octave:2> floor(rand(3,5).*99).+1
ans =
  80  63  74  19  55
  81  15  81  26   3
  81  80  70  16  95

[edytuj] Sprawdzanie czasu wykonania

Obliczyć czas wykonania komendy lub podprogramu. Komenda tic włącza stoper, a toc wyłącza i podaje czas w sekundach, który upłynął od ostatniego tic'a. Na przykład:

octave:85> tic; eig(rand(500,500)); toc
ans = 1.6109
octave:86> tic; toc
ans = 0.00042700

[edytuj] Rozwiązanie układu równań liniowych

Rozwiązać układ równań $\begin{cases} 2x-y=1 \\ x+3y=11\end{cases}$

A = [2 -1; 1  3];
f = [1 11];
u = A\f'

lub alternatywnie:

u=inv(A)*f'

Wynikiem jest:

u = [2 3]

[edytuj] Rozwiązanie równania nieliniowego

Rozwiązać numerycznie równanie $(6 x 4 - x 3 - 7 x 2 + x + 1) + sin(x) = 0$ .

Szukamy miejsc zerowych funkcji $P (x) = (6 x 4 - x 3 - 7 x 2 + x + 1) + s i n (x)$ . Utwórzmy funkcję $P$ . W tym celu tworzymy osobny plik o nazwie P.m i treści:

function [y]=P(x)
   y=(6*x^4-x^3-7*x^2+x+1)+sin(x);
endfunction;

Rozwiązujemy równanie za pomocją funkcji fsolve, która implementuje metodę podobną do metody Newtona. Jako parametry zadajemy funkcję "P" oraz pierwsze przybliżenie miejsca zerowego, np. 0:

octave:5> x0=fsolve("P",0.0)
x0 = -0.27090

Uwaga: Octave musi widzieć plik P.m. Najłatwiej to osiągnąć, uruchamiając octave z tego samego katalogu, w którym znajduje się ten plik.

Miejscem zerowym jest $x_0\approx{}-0.2709$ . Sprawdźmy to:

octave:7> P(x0)   
ans = 0

[edytuj] Zadanie tolerancji

Szukanie miejsc zerowych przy zadanej tolerancji:

fsolve_options("tolerance", 1e-5);
fsolve("P",0.0);

[edytuj] Zadanie funkcji pochodnej

Podawanie ręczne funkcji pochodnej (normalnie pochodna zostanie obliczona w sposób przybliżony):

function [y]=Pprim(x)
   y=(24*x^3-3*x^2-14*x+1) + cos(x);
endfunction;

octave:13> mz=fsolve(["P"; "Pprim"], 0.0)
mz = -0.27090

[edytuj] Szczegóły fsolve

Funkcja fsolve zwraca również dodatkowe informacje:

x - przybliżone miejsce zerowe
info - informacja o wyniku wykonania, wartość 1 oznacza sukces,
msg - komunikat dla użytkownika.

octave:19> [x,info,msg]=fsolve("sin", pi/2)
x = 40.841
info = 3
msg = iteration is not making good progress

octave:20> [x,info,msg]=fsolve("sin", pi/3)
x = 0
info = 1
msg = solution converged within specified tolerance

[edytuj] Rozwiązanie uwikłanego równania nieliniowego

Graficzne rozwiązanie uwikłanego układu równań

Rozwiązać numerycznie uwikłany układ równań:

$\begin{cases} x^2+y^2=1 \\ x+y=0 \end{cases}$

Zauważmy, że układ ten opisuje przecięcie okręgu o środku 0 i promieniu 1 z prostą $y = - x$ . Istnieją zatem dwa rozwiązania i można je obliczyć dokładnie:

$\begin{cases} x= \pm \sqrt{2}/2 \approx \pm 0.707107\\ y= \mp \sqrt{2}/2 \approx \mp 0.707107 \end{cases}$

Numerycznie rozwiążemy równanie za pomocą funkcjji fsolve.

Zdefiniujmy pewną funkcję wektorową $h$ :

$h\begin{pmatrix}x \\ y\end{pmatrix} = \begin{pmatrix}x^2+y^2-1\\x+y\end{pmatrix}$

Będziemy szukać miejsc zerowych dla h, czyli wektorów $\begin{pmatrix}x_0\\y_0\end{pmatrix}$ takich, że $h\begin{pmatrix}x_0\\y_0\end{pmatrix}=\begin{pmatrix}0\\0\end{pmatrix}$

Zdefiniujmy funkcję h w pliku h.m:

function z=h(x)
   a=x(1)*x(1)+x(2)*x(2)-1;
   b=x(1)+x(2);
   z=[a,b];
endfunction;

Zdefiniujmy macierz pochodnej h w pliku hprim.m

function [y]=hprim(x)
   y = [2*x(1), 2*x(2); 1, 1 ];
endfunction;

Rozwiązujemy równanie wektorowe funkcją fsolve:

octave:5> fsolve(["h"; "hprim"], [1.0, -1.0])
ans =
  0.70711
 -0.70711
octave:6> fsolve(["h"; "hprim"], [-1.0, 1.0])
ans =
 -0.70711
  0.70711

Widać, że Octave rozwiązał równania poprawnie.

[edytuj] Rysowanie funkcji zadanej w sposób uwikłany

Rysowanie okregu danego za pomocą funkcji uwikłanej

x 2 + y 2 - 1 = 0

Funkcja zadana jest w sposób uwikłany, należy narysować jej wykres w otoczeniu zadanego punktu. Na przykład: y(x), gdzie $x 2 + y 2 - 1 = 0$ na przedziale $x\in [-1,1]$ (górny półokrąg). Użyć funkcji fsolve.

Stwórzmy funkcję okrag w pliku okrag.m.

function [z]=okrag(y)
   global x;
   z=x*x+y*y-1;
endfunction;

Zanim wywołamy fsolve(okrag) zmienna 'x' będzie odpowiednio ustawiana na zmiennej globalnej.

Zadajmy dane do obliczeń:

#liczba punktow, w ktorych przybliżamy
N=200;
#przedział [A,B], w którym przybliżamy
A=-1;
B=1;
#warunek brzegowy
v=0.0;
#funkcja ktora rozwiklujemy
fun="okrag"
#podziałka
h=(B-A)/N;
#zainicjowanie wyniku
wx=A.+h.*(0:1:N);
wy=zeros(1,N+1);
wy(1)=v;
#zmienna przechowujaca globalne, aktualnie badane x
global x

Rysowanie wykresu funkcji danej w sposób uwikłany:

y + x sin(y) + sin(x) = 0

Obliczamy wektor wy. W każdym kroku obliczamy wartość funkcji w punkcie $x=A+k \cdot h \in [A,B]$ . Jako przybliżoną początkową wartość dla $w x (k + 1)$ bierzemy wartość z poprzedniego kroku $w y (k)$ (to ważne!).

for (k=1:N)
   x=wx(k+1);
   y=fsolve(fun, wy(k));
   wy(k+1)=y;
endfor;

Rysujemy obliczony wykres funkcji:

plot(wx, wy, "-r", wx, -wy, "-b");

Inny przykład narysowany tą metodą to wykres funkcji y(x), gdzie $y + x sin(y) + s i n (x) = 0$

A=-100; B=100;
v=0.005

dla funkcji:

function [z]=u(y)
   global x;
   z=y+x*sin(y)+sin(x);
endfunction;

[edytuj] Rysowanie wykresu funkcji odwrotnej

Szczególnym przypadkiem rysowania funkcji danej w sposób uwikłany jest rysowanie funkcji odwrotnej, np $s i n - 1$ (czyli $a r c s i n$ ).

Rysowanie wykresu funkcji odwrotnej do funkcji sinus.

Stwórzmy plik mysin.m

function [z]=mysin(x)
   global y;
   z=sin(x)-y; 
endfunction;

Zadajmy warunki brzegowe:

A=-1; B=1;
v=-pi/2;
fun="mysin";

Obliczamy funkcję odwrotną:

global y;
for (k=1:N)
   y=wx(k+1);
   wy(k+1)=fsolve(fun,wy(k));
endfor;

Rysujemy:

plot(wx,wy,"-r;sin(x)-y=0;");

[edytuj] Całkowanie numeryczne

Narysować wykres dystrybuanty 1-wymiarowego rozkładu normalnego, tj. wykres funkcji:

$F(x)=\frac{1}{2\pi} \int_{-\infty}^{x}{e^{t^2/2}dt}$

Całki w Octave oblicza się za pomocą quad, na przykład, żeby obliczyć $\int_{0}^{\pi}{\sin(x)dx}$ możemy wykonać:

octave:2> quad("sin",0,pi) 
ans = 2

Dystrybuanta rozkładu normalnego $\mathcal{N}(0,1)$ .

Żeby obliczyć wartość funkcji F zdefiniujmy funkcję podcałkową:

function [y]=myexp(t)
   y=exp(-t*t/2)/sqrt(2*pi);
endfunction;

Następnie zdefiniujmy dystrybuantę (czyli funkcję F):

function [y]=dystrybuanta(x)
   y=quad("myexp", -inf, x);
endfunction;

Rysujemy wykres dystrybuanty F:

x=(-10:0.1:10);
N=length(x);
y=zeros(1,N);
for (k=1:N)
   y(k)=dystrybuanta(x(k));
endfor;
plot(x,y,"-r;Dystrybuanta N(0,1);");

[edytuj] Numeryczne rozwiązywanie równań różniczkowych zwyczajnych pierwszego rzędu

Rozwiązanie równania $y'=-ax,\,y(0)=1,\,a=0.1$ , na przedziale [0,10]

Za pomocją funkcji lsode rozwiązać równania różniczkowe, tj. znaleźć funkcję y=y(x), taką, że:

$\begin{cases}y'=-ay, x\in[0,1], a=\pm{}0.1,\pm{}1,\pm{}100 \\ y(0)=1\end{cases}$
$\begin{cases}y'=y^2, x\in[0,3] \\ y(0)=1\end{cases}$
$\begin{cases}y'=y\cdot{}x\sin(ax), x\in[0,20] \\ y(0)=1\end{cases}$

Rozwiązania powyższych równań można wyznaczyć symbolicznie, są to:

$y=e^{ax}\,$
$y=\frac{1}{1-x}\,$
$y=e^{ (\sin(ax)-ax\cos(ax))/a^2}$ , gdzie $a = 100$

Rozwiązanie równania $y'=y^2,\,y(0)=1$ , na przedziale [0, 0.99]

Zdefiniujmy funkcje, które występują w równaniach:

function y=linfun(x)
   global a;
   y=a.*x;
endfunction;

function y=kwadfun(x)
   y=x.*x;
endfunction;

function y=mixfun(x,t)
   y=x.*t.*sin(t*100);
endfunction;

Rozwiązanie równania $y'=y\cdot{}x\sin(x)$ na przedziale [0,20]. Dosalismy fale o bardzo dużym nachyleniu i coraz większych amplitudach.

Rozwiążmy pierwsze równanie: Zadajemy wartość dla współczynnika a:

global a;
a=0.1;

Zadajemy przedział, na którym będziemy rozwiązywać równanie:

t=[0:0.01:1];

Zadajemy warunek brzegowy w pierwszym punkcie czasowym t(1), czyli w t=0:

x0=1;

Rozwiązujemy równanie:

Z=lsode("linfun", x0, t);

Rysujemy rozwiązanie:

plot(t,Z,"-r")

Podobnie rozwiązujemy drugie równanie:

t=[0:0.01:0.99];
x0=1;
Z=lsode("kwadfun", x0, t);

Funkcja $1 / (1 - x)$ wybucha w $x = 1$ .

Rozwiązujemy trzecie równanie:

t=[0:0.01:20];
x0=1;
Z=lsode("mixfun", x0, t);
axis([0, 20, 0, 2]);
plot(t, Z, "-r");

[edytuj] Numeryczne rozwiązywanie równań różniczkowych zwyczajnych wyższych rzędów

Rozwiązać równanie różniczkowe 2-ego rzędu:

$\begin{cases} x''=f(x)=-x \\ x(0)=0 \\ x'(0)=1 \\ \end{cases}$

Równanie to można scałkować symbolicznie, a jego rozwiązaniem jest $x (t) = s i n (t)$ .

Aby rozwiązać równanie numerycznie, przekształcimy je do równania pierwszego rzędu, w którym zmienna jest wektorem i rozwiążemy za pomocą funkcji lsode. Przekształcone równanie dla zmiennej $\overline{x}=[x_1, x_2]=[x(t), x'(t)]$ ma postać:

$\begin{cases} x_1' = x_2 \\ x_2' = f(x_1) = -x_1 \\ x_1(0) = 0 \\ x_2(0) = 1 \end{cases}$

Wykres rozwiązania

[x, x']

równania

x'' = - x

z warunkami początkowymi

x (0) = 0

x'(0) = 1

Zdefiniujmy funkcję:

function y=ujfun(x)
   y(1)= x(2);
   y(2)=-x(1);
endfunction;

Rozwiązujemy równanie $\overline{x}'=\texttt{ujfun}(\overline{x})$ :

x0=[1;0];
t=[0:0.005:100];
Z=lsode("ujfun", x0, t);

Zauważmy, że warunek początkowy jest pionowym 2-elementowym pionowym wektorem $[1;0]$ . Narysujmy wykres $[x 1, x 2]$ :

plot(Z(:,2),Z(:,1),"-r");

Wynikiem jest okrąg, czyli zbiór $\{\,[\sin(t), \cos(t)]\,: t\in [0,2\pi]\}$ , co zgadza się z rozwiązaniem teoretycznym.

[edytuj] Numeryczne rozwiązywanie równań różniczkowych zwyczajnych z "odwróconym" warunkiem początkowym

Rozpatrzmy zagadnienie Cauchy'ego dla równania różniczkowego:

$\begin{cases} x' = f(x,t), x \in [0,1] \\ x(0) = s \end{cases}$

Równanie jest sparametryzowane parametrem $s\in[0,1]$ i dla zadanej funkcji $f$ . Oznaczmy $y s (t)$ rozwiązanie równania dla ustalonego $s$ . Dla jakiego $s$ dostaniemy $y_s(1)=\frac{1}{2}$ ? Innymi słowy, jeśli oznaczymy $F (s) = y s (1)$ należy rozwiązać równanie:

$F(s)=\frac{1}{2}$

Jak wygląda potok $y s (t)$ i wykres $F (s)$ ? Rozważmy funkcję $f (y, t) = s i n (y 2 + 1) * y * (1 - y)$ .

Wykres potoku dla zagadnienia Cauchy'ego dla równania różniczkowego

x' = sin(x 2 + 1) x (1 - x)

z warunkiem początkowym

x (0) = s

na przedziale $t \in [0,1]$ parametryzowanego parametrem $s \in [0, 1]$ . Linia, która "kończy" się w

(1,0.5)

zaczyna się w przybliżeniu w

(0,0.3)

Najpierw rozwiążemy szukane równanie. Zdefiniujmy funkcję występującą w równaniu:

function y=dfun(x)
   y=sin(x.^2.+1).*x.*(1.0.-x);
endfunction;

Rozwiązujemy równanie dla zadanego parametru $s$ na przedziale $[0,1]$ :

function y=dfunSolveODE(s)
   t=[0:0.01:1];
   y=lsode("dfun", s, t);
endfunction;

Rysujemy wykres funkcji $F (s) = y s (1) = d f u n s o l v e (s)$ :

s=[0:0.01:1];
y=dfunSolveODE(s);
plot(s,y);

Rozwiązujemy równanie $F(s)-\frac{1}{2}=0$ . Zdefiniujmy funkcję:

function y=dfunMinusPol(s)
   z=dfunSolveODE(s);
   y=z(length(z))-0.5;
endfunction;

i znajdźmy jej miejsce zerowe:

fsolve("dfunMinusPol",0.3)
ans = 0.28628

[edytuj] Numeryczne rozwiązywanie równań różniczkowych zwyczajnych z warunkami brzegowymi

Znaleźć funkcję $u = u (t)$ spełniającą równanie różniczkowe z warunkami brzegowymi:

$\begin{cases} u''+u^2=f(t) \\ u(0)=0, u(1)=0 \end{cases}$

Dla funkcji $f (t) = s i n (t)$ .

Aby rozwiązać to zagadnienie trzeba znaleźć taki parametr $s 0$ , żeby rozwiązanie zagadnienia

$\begin{cases} u''+u^2=f(t) \\ u(0)=0, u'(0)=s_0 \end{cases}$

spełniało zadany warunek brzegowy w drugim końcu przedziału, tj. $u (1) = 0$ . Można to zrobić metodą z "odwróconym" warunkiem początkowym, opisaną wcześniej. Przekształćmy zagadnienie na równanie pierwszego rzędu zmiennej $[x 1, x 2] = [u (t), u'(t)]$ :

$\begin{cases} x_1' = x_2 \\ x_2' = -x_1^2 + f(t) \\ x_1(0)=0 \\ x_2(0)=s_0 \end{cases}$

Potok fazowy dla równania

u'' + u 2 = s i n (t)

z warunkami początkowymi

u (0) = 0

u'(0) = s

dla parametrów $s \in [-1,1]$ .

Zdefiniujmy funkcje występujące w równaniu:

function y=d1fun(x,t)
   y(1)=x(2);
   y(2)=-x(1)^2+d1funF(t);
endfunction;
function y=d1funF(t)
   y=sin(t);
endfunction;

Zdefiniujmy funkcję, która rozwiązuje równanie dla zadanego parametru:

function y=d1funSolveODE(s)
   t=[0:0.01:1];
   x0=[0;s];
   y=lsode("d1fun", x0, t);
endfunction;

Narysujmy potok fazowy dla tego równania dla wartości $s_0 \in [-1, 1]$ . (Dokładniej: Funkcja d1funSolveODE działa poprawnie dla pojedynczej wartości parametru s. Aby wywołać rysowanie tak, jak poniżej, trzeba lekko zmienić jej treść).

s=[-1:0.03:1];
y=d1funSolveODE(s);
plot(t,y);

Widzimy, że szukana krzywa przechodząca przez punkty $(0,0)$ oraz $(1,0)$ ma w $t = 0$ nachylenie bliskie 0.

Rozwiązanie równania

u'' + u 2 = s i n (t)

z warunkami brzegowymi

u (0) = 0

u (1) = 0

Wyznaczymy teraz tę wartość dokładnie. Zdefiniujmy funkcję $u s (1)$

function y=d1funFunIn1(s)
   z=d1funSolveODE(s);
   w=z(:,1);
   y=w(length(w));
endfunction;

i rozwiążmy równanie:

s0=fsolve("d1funFunIn1", 0)

Dostajemy wartość:

s0 = -0.15769

Rysujemy rozwiązanie:

z=d1funSolveODE(s0);
plot(t, z(:,1);

Rozwiązanie rzeczywiście przechodzi przez punkty $(0,0)$ i $(1,0)$ .

[edytuj] Wartości własne

Obliczyć wartości własne, wektory własne i wielomian charakterystyczny dla zadanej macierzy

$A=\begin{pmatrix} 3 & 0 & 3 \\ 2 & 1 & -1 \\ 0 & 0 & 2 \end{pmatrix}$

Zadajemy macierz i używamy funkcji eig:

octave:138> A=[3 0 3; 2 1 -1; 0 0 2]
A =
   3   0   3
   2   1  -1
   0   0   2

octave:139> [v,m]=eig(A)
v =
   0.00000   0.70711  -0.39057
   1.00000   0.70711  -0.91132
   0.00000   0.00000   0.13019

m =
  1  0  0
  0  3  0
  0  0  2

Obliczyliśmy wektory własne (kolumny macierzy v) i wartości własne (wyrazy na diagonali macierzy m). Sprawdzamy: jeśli $X$ -wektor własny o wartości własnej $k$ to spełnione musi być równanie $A X - k I = 0$ , czyli

octave:171> X=v(:,2);
octave:172> k=m(2,2);
octave:173> A*X-k*X
ans =
  0
  0
  0

Wielomian charakterystyczny można obliczyć za pomocą funkcji poly:

octave:176> poly(A)
ans =
    1   -6   11   -6

Obliczając ręcznie mamy $χ(λ) = (λ - 1)(λ - 2)(λ - 3) = λ 3 - 6λ 2 + 11λ - 6$ .

[edytuj] Wartość wielomianu

Obliczyć wartość wielomianu $7 x 3 - x 2 + 5$ w punkcie $11$ .

P=[7 -1 0 1];
polyval(P, 11);

Wynikiem jest:

[edytuj] Macierz rzadka

Najprostszym sposobem przechowania macierzy w pamięci komputera jest zapamiętanie jej rozmiarów i ciągu wszystkich elementów macierzy, wiersz po wierszu. Jeśli jednak macierz jest szczególnej postaci można optymalizować sposób przechowywania macierzy w pamięci komputera. Macierz rzadka to macierz, w której występuje dużo zer, wystarczy zatem pamiętać pola niezerowe, tj. zbiór trójek: (numer_wiersza, numer_kolumny, wartość), a pozostałe elementy macierzy są zerami (Uwaga. Faktyczna implementacja macierzy rzadkiej może być zupełnie inna). Dzięki macierzom rzadkim można oszczędzić pamięć, wykonać niektóre algorytmy szybkiej, ale niektóre też wolniej, zwłaszcza jeśli wymagają konwersji na postać normalną.

Utworzyć macierz rzadką, przekształcić ją na normalną i spowrotem na rzadką, dla macierzy:

$\begin{pmatrix} 0 & 13 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 11 \\ 0 & 0 & 0 & 0 & 0 \\ -1 & 0 & 0 & 0 & 0 \end{pmatrix}$

Zadajemy macierz w formacie sparse: wektor współrzędnych x-owych, wektor współrzędnych y-owych, wektor wartości i rozmiary macierzy:

octave:6> w=[6,1,4];
octave:7> y=[1,2,5];
octave:8> v=[-1,13,11];
octave:9>A=sparse(x,y,v,6,5,0)
ans =
  (6 , 1) -> -1
  (1 , 2) -> 13
  (4 , 5) -> 11

Konwersja do trybu normalnego full i spowrotem sparse:

octave:14> B=full(A)
B =
    0   13    0    0    0
    0    0    0    0    0
    0    0    0    0    0
    0    0    0    0   11
    0    0    0    0    0
   -1    0    0    0    0
octave:15> sparse(B)
C =
  (6 , 1) -> -1
  (1 , 2) -> 13
  (4 , 5) -> 11

Zauważmy, że zamist 36 liczb zmiennoprzecinowych składających się na macierz A wystarczy zapamiętać tylko 6 małych liczb całkowitych i 3 zmiennoprzecinokowe w formacie rzadkim.

[edytuj] Regresja liniowa

Regresja liniowa, stopień dla polyfit=3

Zaimplementować regresję liniową. Zadajemy punkty $(x 1, x 2,.., x N)$ i losujemy wartości $(y 1, y 2,.., y N)$ :

X=[-10,-5,-3,-1,-0.5,0,1,2,3,3.1,3.2,4,5];
N=length(X);
B=10.0
Y=B.*rand(1,N)

Używamy gotowej funkcji polyfit, która dopasowuje wielomian zadanego stopnia do punktów $(x k, y k)$

deg=1
O=polyfit(X,Y,deg)

Wykonujemy regresję metodą najmniejszych kwadratów:

Sx = sum(X);
Sy = sum(Y);
Sxx = X*X';
Syy = Y*Y';
Sxy = X*Y';
Sxx2 = norm(X, 2);
Syy2 = norm(Y, 2);
wsp_beta = (N*Sxy-Sx*Sy)/(N*Sxx-Sxx2);
wsp_alfa = (Sy-wsp_beta*Sx)/N;
P = [wsp_beta, wsp_alfa];

Rysujemy punkty i wynik:

axis([-11,11,0,B+1]);
x=[-10:0.25:10];
Z=polyval(O,x.+0.5);
W=polyval(P, x);
plot(x,W,"-r;MNK;",(x.+0.5),Z,"xg;Polyfit;",X,Y,"*b;Dane;");

[edytuj] Interpolacja wielomianowa

Interpolacja wielomianowa dla 21 węzłów

Zaimplementować interpolację wielomianową. Znaleźć wielomian interpolacyjny dla węzłów rółnoodległych na odcinku $[ - 5,5]$ dla funkcji $f (x) = x 2 + 1$ .

Zadajemy wielomian i liczbę węzłów.

P=[1, 0, 1];
N=21;

Zadajemy punkty na podstawie danego wielomianu:

X=[-5:(10.0/N):5];
Y=polyval(P,X);
Y=Y';

Rozwiazanie macierzą Vandermonde'a, z użyciem funkcji vander:

V=vander(X);
A=inv(V)*Y;
B=V\Y;
blad_Vander_inv=norm(polyval(A,X)-Y',2)
blad_Vander_LU=norm(polyval(B,X)-Y',2)

Rozwiazanie Octave'a z użyciem polyfit

R=polyfit(X,Y',3);
blad_polyfit=norm(polyval(R,X)-Y',2)

Rysowanie danych i wyników:

_X=[-5.5:0.05:5.5];
_Y=polyval(R,_X);
_Z=polyval(A,_X);
_W=polyval(B,_X);
axis([-5.5,5.5,-2,30]);
plot(_X,_Y,"-g;Polyfit;",_X,_Z,"-b;Vander inv;",_X,_W,"-c;Vander LU;",X,Y,"*r;Dane;")

Wynik na komputerze autora (typowy PC):

blad_Vander_inv =  4.2279e-08
blad_Vander_LU = 0.024093
blad_polyfit =  7.3644e-16

[edytuj] Dyskretyzacja laplasjanu

Dyskretyzacją laplasjanu jest macierz:

$L_N=\begin{bmatrix} 2 & -1 & 0 & \cdots & 0 & 0 \\ -1 & 2 & -1 & \cdots & 0 & 0 \\ 0 & -1 & 2 & \cdots & 0 & 0 \\ \vdots & \vdots & \ddots & \ddots & \vdots & \vdots \\ 0 & \cdots & \cdots & 2 & -1 & 0 \\ 0 & \cdots & \cdots & -1 & 2 & -1 \\ 0 & \cdots & \cdots & 0 & -1 & 2 \end{bmatrix}$

(a) Utworzyć różnymi metodami macierz $L N$ i porównać czasy wykonania.

Funkcja tworząca macierz $L N$ z użyciem pętli for. Pamięć na macierz alokujemy przed wykonaniem funkcji za pomocą zeros.

function [L]=laplasjan(N)
#Funkcja tworzaca macierz -laplasjanu
L=zeros(N,N);
for k=1:(N-1)
   L(k:(k+1),k:(k+1)) = [2 -1; -1 2];
endfor;
endfunction;

Funkcja tworząca macierz $L N$ z użyciem funkcji diag.

function [L]=laplasjan_diag(N)
#Funkcja tworzaca macierz
v = -1.*ones(1,N-1);
w = 2.*ones(1,N);
L = zeros(N,N) + diag(w, 0) + diag(v, 1) + diag(v, -1);
endfunction;

Przykładowe wykonanie:

octave:178> laplasjan(5)
ans =
   2  -1   0   0   0
  -1   2  -1   0   0
   0  -1   2  -1   0
   0   0  -1   2  -1
   0   0   0  -1   2

Sprawdzamy czas wykonania dla dużych $N$ :

octave:3> tic; laplasjan(2000); toc
ans = 0.22324
octave:4> tic; laplasjan_diag(2000); toc
ans = 0.58647

Wniosek: użycie diag jest wolniejsze niż tworzenie macierzy iteracyjnie.

(b) Rozwiązać dyskretne równanie Laplace'a za pomocą odwrotności tej macierzy, policzonej na dwa różne sposoby. Patrz także zagadnienie własne dla operatora Laplace'a. Równanie ma postać:

$\begin{cases} -u''(x) = \sin(\pi{}x) \\ u(0)=u(1)=0 \end{cases}$

Jego rozwiązaniem jest funkcja $u(x)=\frac{1}{\pi^2}\sin(\pi x)$ . Dyskretyzacja równania ma postać:

$\begin{cases} -u_{n-1}+2 u_{n} - u_{n-1} = h^2 \sin(nh) \\ u_0=u_N=0 \\ \end{cases}$

gdzie $h=\frac{1}{n}$ .

Zatem musimy rozwiązać równanie: $L N u = f$ , gdzie $u=[u_1, u_2, \ldots, u_n], f = h^2 sin([h, 2h, \ldots, Nh])$ .

N=800;
h=1/N;
f=(h*h)*sin((pi*h).*(1:N));
Z=laplasjan(N);
U1=Z\f';
U2=inv(Z)*f';

Obliczamy błędy:

rozw = f*N*N/(pi*pi);
blad_LU  = norm(U1'-rozw, 2)
blad_inv = norm(U2'-rozw, 2)

co na komputerze autora dało rezultat:

blad_LU = 0.0064975
blad_inv =  9.6543e+08

[edytuj] Metoda kolejnych nadrelaksacji

Metoda kolejnych nadrelaksacji to iteracyjna metoda rozwiązywania przybliżonego równania Laplace na stacjonarny rozkład temperatury w ograniczonym obszarze.

Daną mamy kratę rozmiaru $N \times N$ , przykładamy znaną, stałą temperaturę na brzegach tej kraty. Jak będzie wyglądał rozkład temperatury wewnątrz kraty po długim czasie? Oznaczmy punkt kratowy w wierszu i-tym i kolumnie j-tej przez $u i, j$ . W metodzie kolejnych nadrelaksacji obliczamy iteracyjnie dla wewnętrznych punktów kraty:

$u_{i,j}' = (1-\omega) u_{i,j} + \omega \frac{u_{i-1,j}+u_{i+1,j}+u_{i,j-1}+u_{i,j+1}}{4}$

gdzie $\omega \in (0,2)$ jest pewnym współczynnikiem zależnym od obszaru. Dobre dobranie $ω$ ma wpływ na szybkość zbieżności algorytmu.

Uwaga. Metoda zaprezentowana poniżej na pewno nie jest najszybszą, ale jest nietrudna w implementacji.

Utożsamimy kratę $(u_{i,j})_{i,j=1}^{N}$ z wektorem długości $N 2$ utworzonym "wiersz po wierszu". Jedna iteracja algorytmu będzie polegała na mnożeniu przez pewną macierz. Każdemu punktowi kraty odpowiada jeden wiersz tej macierzy:

W wierszu odpowiadającym punktom punktom brzegowym kraty: 1 na przekątnej, pozostałe: 0
W wierszach odpowiadających punktom wewnętrzmym: $(1 - ω)$ na przekątnej oraz $ω / 4$ w czterech innych

kolumnach. Na przykład dla $N = 3$ ma postać:

$\begin{pmatrix} 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 \\ 0 & \omega & 0 & \omega & 1-\omega & \omega & 0 & \omega & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 \end{pmatrix}$

Wykonanie Metody Kolejnych Nadrelakcacji dla kraty 150x150 z temperaturą zadaną na brzegu następująco: 20.0 na jednym brzegu (czerwony), 0.0 na przeciwmym (niebieski), 6.0-8.0 na pozostałych. Obliczenie trwało na komputerze autora (zwykły PC) kilka minut. Wizualizacja zostła wykonana za pomocą IBM OpenDx.

Funkcja tworząca tę macierz dla $N > 2$ :

function [Y]=laplasjan2d_sparse(N)
   omega=1/(1+sin(pi/N))
   M=(N^2-4*N+3)*4;
   x=zeros(1,M);
   y=zeros(1,M);
   v=zeros(1,M);
   k=1;
   for (i=1:N)
      for (j=1:N)
         w=N*(i-1)+j;
         if ((i==1) || (i==N) || (j==1) || (j==N))
            x(k)=w; y(k)=w;   v(k)=1;  k++;
         else
            x(k)=w; y(k)=w;   v(k)=1-omega;  k++;
            x(k)=w; y(k)=w-1; v(k)=omega/4; k++;
            x(k)=w; y(k)=w+1; v(k)=omega/4; k++;
            x(k)=w; y(k)=w-N; v(k)=omega/4; k++;
            x(k)=w; y(k)=w+N; v(k)=omega/4; k++;
          endif;
       endfor;
   endfor;
   Y=sparse(x,y,v,N*N,N*N,0);
endfunction;

Przygotowujemy obliczenia:

#liczba punktów kratowych
N=150;
#tolerancja
epsilon=0.01;
#kolejna relaksacja
Y=laplasjan2d_sparse(N);
#początkowy stan kraty f
f = ...
f=reshape(f,1,N*N);
f=f';

Wykonujemy algorytm kolejnych relaksacji:

do
   g=Y*f;
   err=norm(f-g,2);
   f=g;
until (err<epsilon);
f=reshape(f',N,N);

Przykładowy wynik umieszczony jest na ilustracji obok.

[edytuj] Aproksymacja wielomianowa w przestrzeni Hilberta

Rozpatrzmy przestrzeń Hilberta funkcji całkowalnych z drugą potęgą $L 2 (0,1)$ , z iloczynem skalarnym $<f,g>=\int_{0}^{1}{fg}$ . Dana jest $f\in L^2(0,1)$ . Znaleźć wielomian $\mathcal{W}=\sum_{i=0}^{N}{\alpha_{i}x^i}$ stopnia $N > 1$ najlepiej przybliżający $f$ , czyli taki, który minimalizuje normę $\|f-\mathcal{W}\|_{L^2}$ .

W tym celu należy rozwiązać układ równań:

$G(1,x,x^2,...,x^N) \cdot \overline{\alpha} = \overline{f},$

gdzie $\overline{\alpha}$ jest szukanym wektorem współczyników $\mathcal{W}$ macierz G jest macierzą Grama, a $\overline{f}$ jest wektorem:

$f_i=\int_{0}^{1}{f(x) x^i dx}$

W naszym przypadku macierz Grama G sprowadza się do macierzy Hilberta, gdyż

$G_{i,j}=<x^i,x^j>=\int_{0}^{1}{x^{i+j}}=\frac{1}{i+j+1}$

Algorytm ten można wykonać w środowisku Octave następująco. Zdefiniujmy funkcję charakterystyczną $\chi_{[0,1/2]}\!$ :

function [y]=fcharakt(x)
   y=(x<=0.5);
endfunction;

Do obliczenia wektora $\overline{f}$ potrzebujemy funkcji $f\cdot{}x^n$ :

function [y]=f_razy_xn(x)
   global n;
   global fun;
   z=feval(fun,x);
   y=(z).*(x^n);
endfunction;

Dla zadanego stopnia wielomianu N tworzymy macierz Grama:

N=15;
G=hilb(N+1);

Tworzymy wektor $\overline{f}$ :

f=zeros(1,N+1);
global n;
global fun;
fun="fcharakt";
for (n=0:N)
   f(n+1)=quad("f_razy_xn", 0, 1);
endfor;

I wreszcie obliczamy współczynniki (w odwrotnej kolejności) $\overline{\alpha}=G^{-1}\cdot{}\overline{f}$ :

alfa=G\f'

Uwaga. Ta metoda nie jest dobra dla dużych N, gdyż macierz Hilberta jest źle uwarunkowana.

Stopień wielomianu: 1

Stopień wielomianu: 3

Stopień wielomianu: 5

Stopień wielomianu: 10

Stopień wielomianu: 20

[edytuj] Zobacz też

MATLAB
SCILAB - środowisko podobne do GNU Octave
FreeMat - środowisko podobne do GNU Octave (strona domowa projektu)

[edytuj] Linki zewnętrzne

GNU Octave, strona domowa
Graficzny interfejs Octave Workshop GNU Octave
Repozytorium funkcji dla GNU Octave
Symulacje obliczeń kwantowych w GNU Octave
Cyfrowe przetwarzanie sygnałów w GNU Octave
OpenDX - udostępniony przez IBM wizualizator danych

Źródło: "http://pl.wikipedia.org../../../g/n/u/GNU_Octave_e9e2.html"

Kategorie: Oprogramowanie matematyczne • GNU