|
Przykładowe zadania z rozwiązaniami w [[w:GNU Octave|GNU Octave]].
===Wprowadzanie danych===
Zakończenie linii średnikiem powoduje, że Octave nie wyświetla wyniku działania komendy :
octave:2> B = rand (3, 2);
Zakończenie linii poprzez naciśnięcie klawisza Enter powoduje wyświetlenie wyniku działania komendy :
octave:3> B = rand (3, 2)
B =
0.515199 0.315178
0.765914 0.018528
0.923735 0.585183
Wklejanie danych z Windows do Octave
*Skopiuj odpowiednie dane z Windows do schowka ( zaznacz, potem Ctrl-C )
*przejdź do okna Octave
*użyj środkowego klawisza myszy
===Generowanie danych===
* Wygenerować wektor długości 10 zawierający losowe liczby rzeczywiste z przedziału [-3, 5]
<pre>
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
</pre>
* Wygenerować macierz 3x5 (3 wiersze, 5 kolumn) zawierającą losowe liczby naturalne z przedziału [1, 99]
<pre>
octave:2> floor(rand(3,5).*99).+1
ans =
80 63 74 19 55
81 15 81 26 3
81 80 70 16 95
</pre>
===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 tica.
Na przykład:
<pre>
octave:85> tic; eig(rand(500,500)); toc
ans = 1.6109
octave:86> tic; toc
ans = 0.00042700
</pre>
===Rozwiązanie układu równań liniowych===
Rozwiązać [[w:układ równań|układ równań]] <math>\begin{cases} 2x-y=1 \\ x+3y=11\end{cases}</math>
A = [2 -1; 1 3];
f = [1 11];
u = A\f'
lub alternatywnie:
u=inv(A)*f'
Wynikiem jest:
u = [2 3]
===Rozwiązanie równania nieliniowego===
Rozwiązać numerycznie równanie <math>(6x^4-x^3-7x^2+x+1)+\sin(x)=0</math>.
Szukamy [[w:miejsce zerowe|miejsc zerowych]] funkcji <math>P(x)=(6x^4-x^3-7x^2+x+1)+sin(x)</math>.
Utwórzmy funkcję <math>P</math>. W tym celu tworzymy osobny plik o nazwie <tt>P.m</tt> 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 [[w:metoda Newtona|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 <tt>P.m</tt>. Najłatwiej to osiągnąć, uruchamiając octave z tego samego katalogu, w którym znajduje się ten plik.
Miejscem zerowym jest <math>x_0\approx{}-0.2709</math>. Sprawdźmy to:
octave:7> P(x0)
ans = 0
====Zadanie tolerancji====
Szukanie miejsc zerowych przy zadanej tolerancji:
fsolve_options("tolerance", 1e-5);
fsolve("P",0.0);
====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
====Szczegóły fsolve====
Funkcja fsolve zwraca również dodatkowe informacje:
#<tt>x</tt> - przybliżone miejsce zerowe
#<tt>info</tt> - informacja o wyniku wykonania, wartość 1 oznacza sukces,
#<tt>msg</tt> - 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
===Rozwiązanie uwikłanego równania nieliniowego===
[[Grafika:Uwiklany uklad rownan.svg|thumb|350px|right|Graficzne rozwiązanie uwikłanego układu równań]]
Rozwiązać numerycznie [[w:funkcja uwikłana|uwikłany]] [[w:układ równań|układ równań]]:
:<math>\begin{cases}
x^2+y^2=1 \\
x+y=0
\end{cases}</math>
Zauważmy, że układ ten opisuje przecięcie okręgu o środku 0 i promieniu 1 z prostą <math>y=-x</math>.
Istnieją zatem dwa rozwiązania i można je obliczyć dokładnie:
:<math>
\begin{cases}
x= \pm \sqrt{2}/2 \approx \pm 0.707107\\
y= \mp \sqrt{2}/2 \approx \mp 0.707107
\end{cases}
</math>
Numerycznie rozwiążemy równanie za pomocą funkcjji '''fsolve'''.
Zdefiniujmy pewną funkcję wektorową <math>h</math>:
:<math>
h\begin{pmatrix}x \\ y\end{pmatrix} = \begin{pmatrix}x^2+y^2-1\\x+y\end{pmatrix}
</math>
Będziemy szukać miejsc zerowych dla h, czyli wektorów <math>\begin{pmatrix}x_0\\y_0\end{pmatrix}</math>
takich, że <math>h\begin{pmatrix}x_0\\y_0\end{pmatrix}=\begin{pmatrix}0\\0\end{pmatrix}</math>
Zdefiniujmy funkcję h w pliku <tt>h.m</tt>:
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 <tt>hprim.m</tt>
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.
===Rysowanie funkcji zadanej w sposób uwikłany===
[[Grafika:Uwiklany okrag.svg|thumb|350px|right|Rysowanie okregu danego za pomocą funkcji uwikłanej <math>x^2+y^2-1=0</math>]]
Funkcja zadana jest w sposób uwikłany, należy narysować jej [[w:wykres funkcji|wykres]] w otoczeniu zadanego punktu. Na przykład: y(x), gdzie
<math>x^2+y^2-1=0</math> na przedziale <math>x\in [-1,1]</math> (górny półokrąg).
Użyć funkcji '''fsolve'''.
Stwórzmy funkcję <tt>okrag</tt> w pliku <tt>okrag.m</tt>.
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 punktów, w których przybliżamy
N=200;
#przedział [A,B], w którym przybliżamy
A=-1;
B=1;
#warunek brzegowy
v=0.0;
#funkcja którą rozwikłujemy
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 przechowująca globalne, aktualnie badane x
global x
[[Grafika:Uwiklana funkcja.svg|thumb|350px|right|Rysowanie wykresu funkcji danej w sposób uwikłany: <math>y+x\sin(y)+\sin(x)=0</math>]]
Obliczamy wektor <tt>wy</tt>. W każdym kroku obliczamy wartość funkcji w punkcie <math>x=A+k \cdot h \in [A,B]</math>.
Jako przybliżoną początkową wartość dla <math>wx(k+1)</math> bierzemy wartość z poprzedniego kroku <math>wy(k)</math>(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 <math>y+x\sin(y)+sin(x)=0</math>
A=-100; B=100;
v=0.005
dla funkcji:
function [z]=u(y)
global x;
z=y+x*sin(y)+sin(x);
endfunction;
====Rysowanie wykresu funkcji odwrotnej====
Szczególnym przypadkiem rysowania funkcji danej w sposób uwikłany jest rysowanie [[w:funkcja odwrotna|funkcji odwrotnej]],
np <math>sin^{-1}</math> (czyli <math>arcsin</math>).
[[Grafika:Arcus sinus.svg|thumb|350px|right|Rysowanie wykresu funkcji odwrotnej do funkcji sinus.]]
Stwórzmy plik <tt>mysin.m</tt>
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;");
===Całkowanie numeryczne===
Narysować wykres [[w:dystrybuanta|dystrybuanty]] 1-wymiarowego [[w:rozkład normalny|rozkładu normalnego]], tj. wykres funkcji:
:<math>
F(x)=\frac{1}{2\pi} \int_{-\infty}^{x}{e^{t^2/2}dt}
</math>
Całki w Octave oblicza się za pomocą '''quad''', na przykład, żeby obliczyć
<math>\int_{0}^{\pi}{\sin(x)dx}</math> możemy wykonać:
octave:2> quad("sin",0,pi)
ans = 2
[[Grafika:Dystrybuanta rozkładu normalnego.svg|thumb|right|350px|Dystrybuanta rozkładu normalnego <math>\mathcal{N}(0,1)</math>.]]
Ż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);");
===Numeryczne rozwiązywanie równań różniczkowych zwyczajnych pierwszego rzędu===
[[Grafika:Funkcja expotencjalna rozwiązanie RRZ.svg|thumb|right|350px|Rozwiązanie równania <math>y'=-ax,\,y(0)=1,\,a=0.1</math>, na przedziale [0,10]]]
Za pomocą funkcji '''lsode''' rozwiązać równania różniczkowe, tj. znaleźć funkcję <tt>y=y(x)</tt>, taką, że:
#<math>\begin{cases}y'=-ay, x\in[0,1], a=\pm{}0.1,\pm{}1,\pm{}100 \\ y(0)=1\end{cases}</math>
#<math>\begin{cases}y'=y^2, x\in[0,3] \\ y(0)=1\end{cases}</math>
#<math>\begin{cases}y'=y\cdot{}x\sin(ax), x\in[0,20] \\ y(0)=1\end{cases}</math>
Rozwiązania powyższych równań można wyznaczyć symbolicznie, są to:
#<math>y=e^{ax}\,</math>
#<math>y=\frac{1}{1-x}\,</math>
#<math>y=e^{ (\sin(ax)-ax\cos(ax))/a^2}</math>, gdzie <math>a=100</math>
[[Grafika:Funkcja hiperboliczna rozwiązanie RRZ.svg|thumb|right|350px|Rozwiązanie równania <math>y'=y^2,\,y(0)=1</math>, 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;
[[Grafika:Funkcja trygonometryczna rozwiązanie RRZ.svg|thumb|right|500px|Rozwiązanie równania <math>y'=y\cdot{}x\sin(x)</math> 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 <tt>a</tt>:
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 <tt>t(1)</tt>, czyli w <tt>t=0</tt>:
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 <math>1/(1-x)</math> wybucha w <math>x=1</math>.
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");
===Numeryczne rozwiązywanie równań różniczkowych zwyczajnych wyższych rzędów===
Rozwiązać [[w:równanie różniczkowe|równanie różniczkowe]] 2. rzędu:
:<math>\begin{cases}
x''=f(x)=-x \\
x(0)=0 \\
x'(0)=1 \\
\end{cases}</math>
Równanie to można scałkować symbolicznie, a jego rozwiązaniem jest <math>x(t)=sin(t)</math>.
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 <math>\overline{x}=[x_1, x_2]=[x(t), x'(t)]</math> ma postać:
:<math>\begin{cases}
x_1' = x_2 \\
x_2' = f(x_1) = -x_1 \\
x_1(0) = 0 \\
x_2(0) = 1
\end{cases}</math>
[[Grafika:Okrag Rozwiazanie 2-wym RRZ.png|thumb|right|350px|Wykres rozwiązania <math>[x,x']</math> równania <math>x''=-x</math> z warunkami początkowymi <math>x(0)=0</math>, <math>x'(0)=1</math>]]
Zdefiniujmy funkcję:
function y=ujfun(x)
y(1)= x(2);
y(2)=-x(1);
endfunction;
Rozwiązujemy równanie <math>\overline{x}'=\texttt{ujfun}(\overline{x})</math>:
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 <math>[1; 0]</math>.
Narysujmy wykres <math>[x_1, x_2]</math>:
plot(Z(:,2),Z(:,1),"-r");
Wynikiem jest okrąg, czyli zbiór <math>\{\,[\sin(t), \cos(t)]\,: t\in [0,2\pi]\}</math>, co zgadza się z rozwiązaniem teoretycznym.
===Numeryczne rozwiązywanie równań różniczkowych zwyczajnych z "odwróconym" warunkiem początkowym===
Rozpatrzmy [[w:zagadnienie Cauchy'ego|zagadnienie Cauchy'ego]] dla [[w:równanie różniczkowe|równania różniczkowego]]:
:<math>\begin{cases}
x' = f(x,t), x \in [0,1] \\
x(0) = s
\end{cases}</math>
Równanie jest sparametryzowane parametrem <math>s\in[0,1]</math> i dla zadanej funkcji <math>f</math>.
Oznaczmy <math>y_s(t)</math> rozwiązanie równania dla ustalonego <math>s</math>.
Dla jakiego <math>s</math> dostaniemy <math>y_s(1)=\frac{1}{2}</math>? Innymi słowy, jeśli oznaczymy
<math>F(s)=y_s(1)</math> należy rozwiązać równanie:
:<math>
F(s)=\frac{1}{2}
</math>
Jak wygląda potok <math>y_s(t)</math> i wykres <math>F(s)</math>?
Rozważmy funkcję <math>f(y,t)=sin(y^2+1)*y*(1-y)</math>.
[[Grafika:Rrz odwrotny war pocz.png|thumb|right|400px|
Wykres potoku dla zagadnienia Cauchy'ego dla równania różniczkowego <math>x'=\sin(x^2+1)x(1-x)</math> z warunkiem początkowym <math>x(0) = s</math> na przedziale <math>t \in [0,1]</math> parametryzowanego parametrem <math>s \in [0, 1]</math>. Linia, która "kończy" się w <math>(1, 0.5)</math> zaczyna się w przybliżeniu w <math>(0, 0.3)</math>.]]
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 <math>s</math> na przedziale <math>[0,1]</math>:
function y=dfunSolveODE(s)
t=[0:0.01:1];
y=lsode("dfun", s, t);
endfunction;
Rysujemy wykres funkcji <math>F(s)=y_s(1)=dfunsolve(s)</math>:
s=[0:0.01:1];
y=dfunSolveODE(s);
plot(s,y);
Rozwiązujemy równanie <math>F(s)-\frac{1}{2}=0</math>. Zdefiniujmy funkcję:
function y=dfunMinusPol(s)
z=dfunSolveODE(s);
y=z(length(z))-0.5;
endfunction;
i znajdźmy jej [[w:miejsce zerowe|miejsce zerowe]]:
fsolve("dfunMinusPol",0.3)
ans = 0.28628
===Numeryczne rozwiązywanie równań różniczkowych zwyczajnych z warunkami brzegowymi ===
Znaleźć funkcję <math>u=u(t)</math> spełniającą [[w:równanie różniczkowe|równanie różniczkowe]] z [[w:warunek brzegowy|warunkami brzegowymi]]:
:<math>
\begin{cases}
u''+u^2=f(t) \\
u(0)=0, u(1)=0
\end{cases}</math>
Dla funkcji <math>f(t)=sin(t)</math>.
Aby rozwiązać to zagadnienie trzeba znaleźć taki parametr <math>s_0</math>, żeby rozwiązanie zagadnienia
:<math>
\begin{cases}
u''+u^2=f(t) \\
u(0)=0, u'(0)=s_0
\end{cases}</math>
spełniało zadany warunek brzegowy w drugim końcu przedziału, tj. <math>u(1)=0</math>.
Można to zrobić metodą z "odwróconym" warunkiem początkowym, opisaną wcześniej.
Przekształćmy zagadnienie na równanie pierwszego rzędu zmiennej <math>[x_1, x_2] = [u(t), u'(t)]</math>:
:<math>
\begin{cases}
x_1' = x_2 \\
x_2' = -x_1^2 + f(t) \\
x_1(0)=0 \\
x_2(0)=s_0
\end{cases}</math>
[[Grafika:RRZ z warunkiem brzegowym, potok fazowy.png|thumb|right|400px|Potok fazowy dla równania <math>u''+u^2=sin(t)</math> z warunkami początkowymi <math>u(0)=0</math> i <math>u'(0)=s</math> dla parametrów <math>s \in [-1,1]</math>.]]
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 [[w:potok fazowy|potok fazowy]] dla tego równania dla wartości <math>s_0 \in [-1, 1]</math>.
(Dokładniej: Funkcja <tt>d1funSolveODE</tt> 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 <math>(0,0)</math> oraz <math>(1,0)</math>
ma w <math>t=0</math> nachylenie bliskie 0.
[[Grafika:RRZ z warunkiem brzegowym, rozwiazanie.png|thumb|right|350px|Rozwiązanie równania <math>u''+u^2=sin(t)</math> z warunkami brzegowymi <math>u(0)=0</math> i <math> u(1)=0</math>.]]
Wyznaczymy teraz tę wartość dokładnie. Zdefiniujmy funkcję <math>u_s(1)</math>
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 <math>(0,0)</math> i <math>(1,0)</math>.
===Rozwiązywanie równań różniczkowych implementując schematy różnicowe===
Rozwiązać numerycznie [[w:zagadnienie Cauchy'ego|zagadnienie Cauchy'ego]]:
:<math>\begin{cases}
x'=f(x)\\
x(a)=x_0
\end{cases}</math>
na przedziale <math>t\in[a,b]</math>. Zaimplementować schemat Eulera i schemat Mid-Point.
====Konstrukcja schematu Eulera====
Załóżmy, że <math>x \in C^1[a,b]</math>. Wówczas ze [[w:wzór Taylora|wzoru Taylora]]:
:<math>x(t+h) = x(t) + h x'(t) + o(h^2) = x(t) + h f(x) + O(h^2) </math>
Oznaczamy <math>x(t+kh)=x_k</math> i pomijamy wyrazy drugiego rzędu:
:<math>x_k=x_{k-1}+h f(x_k)</math>
Schemat ten jest:
*''otwarty'' (kolejny wyraz wylicza się jako nieuwikłana funkcja poprzednich)
*''jednokrokowy'' (do obliczenia <math>x_{k+1}</math> używa tylko <math>x_{k}</math>)
*''samostartujący'' (podajemy <math>x_0</math> i dostajemy wszystkie następne wyrazy)
*''rzędu <math>h^1</math>'' (użyliśmy wzoru Taylora z dokładnością do wyrazów rzędu 1)
====Konstrukcja schematu Mid-point====
Załóżmy, że <math>x \in C^2[a,b]</math>. Wówczas podobnie:
:<math>x(t+h) = x(t) + h x'(t) + h^2 x''(t) / 2 + O(h^3)</math>
:<math>x(t-h) = x(t) - h x'(t) + h^2 x''(t) / 2 + O(h^3)</math>
Odejmujemy stronami i dostajemy:
:<math>x(t+h)=x(t-h)+2 h x'(t) + O(h^3)</math>
Pomijamy wyrazy rzędu 3 i dyskretyzujemy:
:<math>x_{k+1}=x_{k-1}+2 h f(x_k)</math>
Schemat ten jest:
*''otwarty''
*''dwukrokowy'' (do obliczenia <math>x_{k+1}</math> używa <math>x_k</math> i <math>x_{k-1}</math>)
*''nie samostartujący'' (potrzebuje <math>x_0</math> i <math>x_1</math> żeby wystartować)
*''rzędu <math>h^2</math>'' (użyliśmy wzoru Taylora z dokładnością do wyrazów rzędu 2)
====Rozwiązywanie równania za pomocą schematu====
Rozwiązanie zagadnienia Cauchy'ego na przedziale <math>[a, b]</math> za pomocą pewnego schematu:
# Wyznaczamy krok <math>h</math>, który powinien być "wystarczająco" mały.
# Obliczamy ciąg <math>(x_k)_{k=0}^{k=N}</math>, gdzie <math>N=(b-a)/h</math>. Zadajemy <math>x_0=x(a)</math>. Jeśli schemat nie jest samostartujący możemy zadać brakujące wyrazy początkowe wyznaczone na przykład za pomocą schematu Eulera, tj. <math>x_1=x_0+h f(x_0)</math>. Kolejne wyrazy obliczamy według reguły podanej w schemacie.
# Rozwiązaniem jest funkcja, przybliżona w punktach: <math>x(a+k h)=x_k</math>
====Implementacja schematu Eulera====
function y=euler(fun,a,b,x0,N)
h=(b-a)/(N-1);
y=zeros(N,1);
y(1)=x0;
for (k=2:N)
y(k)=y(k-1)+h*feval(fun,y(k-1));
endfor;
endfunction;
====Implementacja schematu Mid-point====
function y=midpoint(fun,a,b,x0,N)
h=(b-a)/(N-1);
y=zeros(N,1);
y(1)=x0;
y(2)=x0+h*feval(fun,x0);
for (k=3:N)
y(k)=y(k-2)+2*h*feval(fun,y(k-1));
endfor;
endfunction;
====Testowanie schematów====
[[Grafika:Schemat Eulera a Midpoint.png|thumb|right|400px|Porównanie schematów Eulera i Midpoint dla funkcji <math>x(t)=exp(-t)</math>.]]
Przetestujmy schematy Euler i Mid-point dla równania:
:<math>\begin{cases}
x'=f(x)=-x, t \in [0, 5]\\
x(0)=1
\end{cases}</math>
Rozwiązaniem tego równania jest <math>x(t)=\exp(-t)</math>.
Zdefiniujmy funkcję:
function y=fun(x)
y=-x;
endfunction;
Zadajmy schematy:
a=0
b=5
x0=1
N=100
Obliczmy schematy:
z=midpoint ("fun",a,b,x0,N);
y=euler ("fun",a,b,x0,N);
Obliczmy rozwiązanie teoretyczne:
t=linspace (a,b,N);
w=1./exp(t);
Rysujmy rozwiązanie teoretyczne i dwa rozwiązania obliczone za pomocą schematów:
plot(t,z,"-r;midpoint;",t,y,"-b;euler;",t,w,"-g;rozwiazanie;");
Widać, że schemat midpoint lepiej przybliża rozwiązanie na początku (jest wyższego rzędu), ale potem "rozjeżdza" się,
w przeciwieństwie do schematu Eulera, który okazuje się być stabilniejszy.
===Rozwiązywanie równań różniczkowych z użyciem gotowych schematów różnicowych===
Rozwiązać za pomocą funkcji '''lsode''' [[w:Sztywny układ równań|sztywny układ równań]] używając [[w:Schemat Adamsa|schematu Adamsa]] i "stiff":
:<math>\begin{cases}
x' = 998 x + 1998 y \\
y' = -999 x - 1999 y \\
x(0) = y(0) = 1
\end{cases}</math>
na odcinku <math>[0,10]</math>. Sprawdzić czy ten układ sztywny?
Sztywność układu sprawdzamy obliczając wartości własne macierzy:
:<math>A=\begin{pmatrix}
998 & 1998 \\
-999 &-1999
\end{pmatrix}</math>
Dostajemy:
A=[998 1998; -999 -1999];
m=eig(A)
m =
-1
-1000
a zatem współczynnik sztywności wynosi wynosi <math>1000</math>, czyli możemy określić układ jako sztywny.
Rozwiążmy go więc napierw metodą "stiff". Zdefiniujmy odpowiednią funkcję:
function y=fun_sztywna(x)
y(1)= 998*x(1)+1998*x(2);
y(2)=-999*x(1)-1999*x(2);
endfunction;
Rozwiązujemy i mierzymy czas:
lsode_options("integration method", "stiff");
tic; z=lsode("fun_sztywna", x0, t); toc
ans = 0.32850
Spróbujmy także schematem Adamsa:
lsode_options("integration method", "adams");
tic; z=lsode("fun_sztywna", x0, t); toc
ans = 26.846
Jak widać mimo, że schemat "stiff" jest schematem zamkniętym (w każdym kroku rozwiązuje równanie), to działa on
znacznie szybciej niż otwarty schemat Adamsa, gdyż rekompensuje sobie długością kroku.
===Chaos w atraktorze Lorentza===
Rozwiązać klasyczny układ równań Lorentza, opisujący uproszczony model pogody:
:<math>
\begin{cases}
\frac{dx}{dt} = \sigma (y-x) \\
\frac{dy}{dt} = rx-y-xz \\
\frac{dz}{dt} = xy-bz \\
x(0)=0, y(0)=1, z(0)=0
\end{cases}
</math>
gdzie <math>\sigma=10</math>, <math>r=8/3</math>, <math>b=28</math>
[[Grafika:Klasyczny uklad lorentza.png|thumb|right|400px|Rozwiązanie klasycznego układu równań różniczkowych Lorentza.]]
Zdefiniujmy funkcję:
function y=lorentz(x)
global sigma;
global r;
global b;
y(1)=sigma*(x(2)-x(1));
y(2)=r*x(1)-x(2)-x(1)*x(3);
y(3)=x(1)*x(2)-b*x(3);
endfunction;
Zadajmy parametry:
global sigma;
global r;
global b;
sigma=10
b=8/3
r=28
Rysujemy ewolucję układu na odcinku <math>[0,300]</math> w 30 krokach, każdy na kolejnym odcinku o długości <math>10</math>.
Punktem początkowym jest punkt końcowy w poprzedniego kroku. Polecenie '''hold on''' wyświetla to, co
jest narysowane do tej pory. W ten sposób możemy obserwować ewolucję układu Lorentza. Polecenie '''hold off''' kończy rysowanie.
x0=[0 1 0]';
for (i=0:30)
t=[i*10:0.01:(i+1)*10];
z=lsode("lorentz", x0, t);
plot(z(:,1),z(:,2));
x0=z(length(z),:)';
hold on
endfor;
hold off
Obserwujemy, że rozwiązanie nie stabilizuje się, zachowuje się "chaotycznie", wciąż
kręci się w atraktorze wpadając to do lewego skrzydła to do prawego.
===Wartości własne===
Obliczyć [[w:wartość własna|wartości własne]], [[w:wektor własny|wektory własne]] i [[w:wielomian charakterystyczny|wielomian charakterystyczny]]
dla zadanej [[w:macierz|macierzy]]
:<math>A=\begin{pmatrix}
3 & 0 & 3 \\
2 & 1 & -1 \\
0 & 0 & 2
\end{pmatrix}</math>
Zadajemy macierz i używamy funkcji '''eig''':
<pre>
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
</pre>
Obliczyliśmy wektory własne (kolumny macierzy v) i wartości własne (wyrazy na diagonali macierzy m).
Sprawdzamy: jeśli <math>X</math>-wektor własny o wartości własnej <math>k</math> to spełnione musi być
równanie <math>AX-kX=0</math>, czyli
<pre>
octave:171> X=v(:,2);
octave:172> k=m(2,2);
octave:173> A*X-k*X
ans =
0
0
0
</pre>
Wielomian charakterystyczny można obliczyć za pomocą funkcji '''poly''':
<pre>
octave:176> poly(A)
ans =
1 -6 11 -6
</pre>
Obliczając ręcznie mamy
<math>\chi(\lambda)=(\lambda-1)(\lambda-2)(\lambda-3)=\lambda^3-6\lambda^2+11\lambda-6</math>.
===Wielomiany===
definicja wielomianu
octave:3> c=[1 1 1]
c = 1 1 1
wyświetlenie wielomianu
octave:4> polyout(c, 'z')
1*z^2 + 1*z^1 + 1
Miejsca zerowe wielomianu czyli rozwiązania równania
<math>1*z^2 + 1*z^1 + 1=0 \,</math>
octave:5> roots(c)
ans =
-0.50000 + 0.86603i
-0.50000 - 0.86603i
Obliczyć wartość wielomianu <math>7x^3-x^2+5 \,</math> w punkcie <math>11\,</math>.
P=[7 -1 0 5];
polyval(P, 11)
Wynikiem jest
9197
===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. [[w:Macierz rzadka|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:
:<math>\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}
</math>
Zadajemy macierz w formacie '''sparse''': wektor współrzędnych x-owych, wektor współrzędnych y-owych, wektor wartości i rozmiary macierzy:
<pre>
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
</pre>
Konwersja do trybu normalnego '''full''' i spowrotem '''sparse''':
<pre>
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
</pre>
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.
===Regresja liniowa===
[[Grafika:Regresja.svg|thumb|right|450 px|Regresja liniowa, stopień dla polyfit=3]]
Zaimplementować [[w:regresja liniowa|regresję liniową]].
Zadajemy punkty <math>(x_1, x_2, .., x_N)</math> i losujemy wartości <math>(y_1, y_2, .., y_N)</math>:
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 <math>(x_k, y_k)</math>
deg=1
O=polyfit(X,Y,deg)
Wykonujemy regresję [[w:metoda najmniejszych kwadratów|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;");
===Interpolacja wielomianowa===
[[Grafika:Interpolacja.svg|thumb|450px|right|Interpolacja wielomianowa dla 21 węzłów]]
Zaimplementować [[w:interpolacja wielomianowa|interpolację wielomianową]]. Znaleźć wielomian interpolacyjny dla węzłów rółnoodległych na odcinku <math>[-5,5]</math> dla funkcji <math>f(x)=x^2+1</math>.
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 [[w:Macierz Vandermonde'a|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
===Dyskretyzacja laplasjanu===
Dyskretyzacją [[w:laplasjan|laplasjanu]] jest macierz:
:<math>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}</math>
(a) Utworzyć różnymi metodami macierz <math>L_N</math> i porównać czasy wykonania.
Funkcja tworząca macierz <math>L_N</math> 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 <math>L_N</math> 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:
<pre>
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
</pre>
Sprawdzamy czas wykonania dla dużych <math>N</math>:
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 [[w:równanie Laplace'a|równanie Laplace'a]] za pomocą odwrotności tej macierzy, policzonej na dwa różne sposoby. Patrz także [[w:zagadnienie własne dla operatora Laplace'a|zagadnienie własne dla operatora Laplace'a]].
Równanie ma postać:
:<math>\begin{cases}
-u''(x) = \sin(\pi{}x) \\
u(0)=u(1)=0
\end{cases}</math>
Jego rozwiązaniem jest funkcja <math>u(x)=\frac{1}{\pi^2}\sin(\pi x)</math>. Dyskretyzacja równania ma postać:
:<math>\begin{cases}
-u_{n-1}+2 u_{n} - u_{n-1} = h^2 \sin(nh) \\
u_0=u_N=0 \\
\end{cases}</math>
gdzie <math>h=\frac{1}{n}</math>.
Zatem musimy rozwiązać równanie: <math>L_N u = f</math>,
gdzie <math>u=[u_1, u_2, \ldots, u_n], f = h^2 sin([h, 2h, \ldots, Nh])</math>.
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
===Metoda kolejnych nadrelaksacji===
[[w:Metoda kolejnych nadrelaksacji|Metoda kolejnych nadrelaksacji]] to iteracyjna metoda rozwiązywania przybliżonego równania Laplace'a na stacjonarny rozkład temperatury w ograniczonym obszarze.
Daną mamy kratę rozmiaru <math>N \times N</math>, 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 <math>u_{i,j}</math>. W metodzie kolejnych nadrelaksacji obliczamy iteracyjnie dla
wewnętrznych punktów kraty:
:<math>
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}
</math>
gdzie <math>\omega \in (0,2)</math> jest pewnym współczynnikiem zależnym od obszaru. Dobre dobranie <math>\omega</math>
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ę <math>(u_{i,j})_{i,j=1}^{N}</math> z wektorem długości <math>N^2</math> 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: <math>(1-\omega)</math> na przekątnej oraz <math>\omega/4</math> w czterech innych
kolumnach.
Na przykład dla <math>N=3</math> ma postać:
:<math>
\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}
</math>
[[Grafika:Metoda Kolejnych Nadrelaksacji.jpg|thumb|right|400px|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 [[w:OpenDx|OpenDx]].]]
Funkcja tworząca tę macierz dla <math>N>2</math>:
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.
===Aproksymacja wielomianowa w przestrzeni Hilberta===
Rozpatrzmy [[w:przestrzeń Hilberta|przestrzeń Hilberta]] [[w:funkcja całkowalna|funkcji całkowalnych]] z drugą potęgą <math>L^2(0,1)</math>, z iloczynem
skalarnym <math><f,g>=\int_{0}^{1}{fg}</math>.
Dana jest <math>f\in L^2(0,1)</math>. Znaleźć wielomian <math>\mathcal{W}=\sum_{i=0}^{N}{\alpha_{i}x^i}</math> stopnia <math>N>1</math>
najlepiej przybliżający <math>f</math>, czyli taki, który minimalizuje normę <math>\|f-\mathcal{W}\|_{L^2}</math>.
W tym celu należy rozwiązać układ równań:
:<math>
G(1,x,x^2,...,x^N) \cdot \overline{\alpha} = \overline{f},
</math>
gdzie
<math>\overline{\alpha}</math> jest szukanym wektorem współczyników <math>\mathcal{W}</math>
macierz G jest [[w:macierz Grama|macierzą Grama]],
a <math>\overline{f}</math> jest wektorem:
:<math>
f_i=\int_{0}^{1}{f(x) x^i dx}
</math>
W naszym przypadku macierz Grama G sprowadza się do [[w:macierz Hilberta|macierzy Hilberta]], gdyż
:<math>
G_{i,j}=<x^i,x^j>=\int_{0}^{1}{x^{i+j}}=\frac{1}{i+j+1}
</math>
Algorytm ten można wykonać w środowisku Octave następująco. Zdefiniujmy funkcję charakterystyczną <math>\chi_{[0,1/2]}\!</math>:
function [y]=fcharakt(x)
y=(x<=0.5);
endfunction;
Do obliczenia wektora <math>\overline{f}</math> potrzebujemy funkcji <math>f\cdot{}x^n</math>:
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 <math>\overline{f}</math>:
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) <math>\overline{\alpha}=G^{-1}\cdot{}\overline{f}</math>:
alfa=G\f'
'''Uwaga'''. Ta metoda nie jest dobra dla dużych N, gdyż macierz Hilberta jest źle [[w:uwarunkowanie macierzy|uwarunkowana]].
{|align=center
|-
|[[Grafika:Aproksymacja wielomianowa funkcji charakterystycznej, stopien 1.svg|thumb|350px|Stopień wielomianu: 1]]
|[[Grafika:Grafika-Aproksymacja wielomianowa funkcji charakterystycznej, stopien 3.svg|thumb|350px|Stopień wielomianu: 3]]
|}
{|align=center
|-
|[[Grafika:Aproksymacja wielomianowa funkcji charakterystycznej, stopien 5.svg|thumb|300px|Stopień wielomianu: 5]]
|[[Grafika:Aproksymacja wielomianowa funkcji charakterystycznej, stopien 10.svg|thumb|300px|Stopień wielomianu: 10]]
|[[Grafika:Aproksymacja wielomianowa funkcji charakterystycznej, stopien 20.svg|thumb|300px|Stopień wielomianu: 20]]
|}
==Tworzenie rysunków==
===Rysunki 2D===
Rysunki dwuwymiarowe można tworzyć za pomocą funkcji '''plot'''. Najczęstsze użycie to:
plot(x, y)
Oba argumenty są wektorami tej samej długości i zawierają współrzędnych odpowiednio x-owe i y-owe.
[[Grafika:Przyblizenie funkcji sinus.png|thumb|right|350px|Wykres funkcji sinus. Rysowanie domyślne.]]
Zatem plot dostaje ciąg punktów: <math>P_i=( x_i, y_i )</math> i rysuje linie:
<math>P_1-P_2, ..., P_{n-1}-P_n</math>. Na przykład:
x=[0:1:10]
y=sin(x)
plot(x,y)
Dostajemy przybliżenie wykresu funkcji sinus.
[[Grafika:Przyblizenie funkcji sinus 2.png|thumb|right|350px|Drugi wykres funkcji sinus. ]]
Punkty domyślnie zostały połączone liniami i narysowane na czerwono. To można zmienić, np.
plot(x,y,"+b")
zaznaczy punkty znakami <math>\times</math> i narysuje je na niebiesko. Żeby lepiej było widać
chcemy też wyświetlić kratę punktów:
grid on
Dostępne kolory i ich oznaczenia:
{| class="wikitable"
!oznaczenie literowe
!oznaczenie numeryczne
!nazwa koloru
|-
| k
| 0
| black
|-
| r
| 1
| red (domyślny)
|-
| g
| 2
| green
|-
| b
| 3
| blue
|-
| m
| 4
| magenta
|-
| c
| 5
| cyan
|-
| w
| 6
| brown
|}
Dostepne style to:
+ * o x - . @ -@
[[Grafika:Funkcje trygonometryczne.png|thumb|right|350px|Wykres funkcji trygonometrycznych. Cztery wykresy na jednym. Każdy innym kolorem i umieszczony w legendzie. Okno ograniczone do <math>[0,10] \times [-\!2,2]</math>.]]
Wykresy można łączyć. Każdy z nich wówczas dobrze jest narysować innym kolorem i opisać w legendzie.
Opis umieszcza się między dwoma średnikami, na przykład:
t = [0:0.1:10];
plot(
t,cos(t),"-r;cos(t);",
t,sin(t),"+m;sin(t);",
t,tan(t),"xc;tan(t);",
t,atan(t),"*g;argtan(t);" );
W tym wykresie funkcja tangens eksplotuje, więc przyda się jeszcze ograniczyć zakres osi:
axis( [0 10 -2 2] )
===Rysunki 3D===
[[Grafika:Wykres domeku.png|thumb|right|350px|Domek narysowany za pomocą '''mesh'''.]]
Narysować domek za pomocą funkcji '''mesh'''.
mesh(x,y,z)
rysuje wykres dla punktów <math>(x_i, y_j, z_{i,j})</math>. Zatem domek o kwadratowej podstawie <math>[0,1]^2</math>,
wysokości <math>1</math> i o dachu wysokości <math>1</math> można narysować za pomocą komendy:
mesh([0 0.5 1], [0 1], [1 2 1; 1 2 1 ])
[[Grafika:Wykres 3d.png|thumb|right|350px|Wykres funkcji <math>sin(x^2+2y^2)</math>.]]
Narysować za pomocą za pomocą funkcji '''mesh''' i '''meshgrid''' wykres funkcji
:<math>f: [-1,1]^2 \to \mathbb{R}</math>
:<math>f(x,y)=sin(x^2+2y^2)\,</math>
Przygotujmy brzegi:
x=[-1:0.1:1];
y=[-1:0.1:1];
Obliczamy wartości:
[X,Y]=meshgrid(x,y);
Z=sin(X.*X.+2.*Y.*Y)
Rysujemy:
mesh(x,y,Z)
===Zapis do pliku===
Sposób zapisu do pliku zależy od wersji GNU Octave.
Jeśli rysujemy za pomocą [[w:Gnuplot|Gnuplota]], można to zrobić "niskopoziomowo" za pomocą '''__gnuplot_raw__''':
__gnuplot_raw__("set term png \n");
__gnuplot_raw__("set output \"plik.png\" \n");
plot(...)
Wynikiem będzie "plik.png" w formacie [[w:PNG|PNG]].
Inna metoda, już odradzana, to użycie funkcji '''gset'''
gset term postscript eps color
gset output "plik.eps"
plot(...)
===Zaawansowany gnuplot===
Z poziomu Octave'a można użyć zaawansowanych opcji [[w:Gnuplot|Gnuplota]]. Wykonując
__gnuplot_raw__("polecenie\n");
to tak, jakbyśmy wykonywali <tt>polecenie</tt> w samym gnuplocie.
Na przykład:
__gnuplot_raw__("set key off\n")
==Źródła==
* [http://www.octave.org/ Strona domowa GNU Octave]
* [http://www.gnu.org/software/octave/doc/interpreter/ Dokumentacja GNU Octave]
* [[w:GNU Octave|GNU Octave na Wikipedii]]
* [http://math.iu-bremen.de/oliver/teaching/iub/resources/octave/octave-intro/octave-intro.html Wprowadzenie do Octave]
[[Kategoria:Programy matematyczne]]
[[en:Octave Programming Tutorial]]
[[fr:Programmation Octave]]
| 