GNU Octave/Dyskretyzacja laplasjanu

Dyskretyzacja laplasjanu

Dyskretyzacją laplasjanu jest macierz:

$L_N=\left[\begin{array}{cccccc}2& -1& 0& \cdots & 0& 0\\ -1& 2& -1& \cdots & 0& 0\\ 0& -1& 2& \cdots & 0& 0\\ ⋮& ⋮& \ddots & \ddots & ⋮& ⋮\\ 0& \cdots & \cdots & 2& -1& 0\\ 0& \cdots & \cdots & -1& 2& -1\\ 0& \cdots & \cdots & 0& -1& 2\end{array}\right]$

(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{array}{c}-u^{″}(x)=\sin (\pi x)\\ u(0)=u(1)=0\end{array}$

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

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

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

Zatem musimy rozwiązać równanie: $L_{N}u=f$, gdzie $u=[u_1,u_2,\dots ,u_n],f=h^{2}sin([h,2h,\dots ,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

Szablon:Nawigacja