GNU Octave/Metoda kolejnych nadrelaksacji

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 $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 $\omega$ 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:

1. W wierszu odpowiadającym punktom punktom brzegowym kraty: 1 na przekątnej, pozostałe: 0
2. W wierszach odpowiadających punktom wewnętrzmym: $(1-\omega )$ na przekątnej oraz $\omega /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.