GNU Octave/Dyskretyzacja laplasjanu

Z Wikibooks, biblioteki wolnych podręczników.

[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 LN i porównać czasy wykonania.

Funkcja tworząca macierz LN 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 LN 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: LNu = 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


Spis treści
Źródło „http://pl.wikibooks.org/wiki/GNU_Octave/Dyskretyzacja_laplasjanu”