Tekst źródłowy strony GNU Octave/Dyskretyzacja laplasjanu


===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

<noinclude>
{{Nawigacja|GNU Octave|
[[../Interpolacja wielomianowa/]]|
[[../Metoda kolejnych nadrelaksacji/]]|
}}</noinclude>