GNU Octave/Metoda kolejnych nadrelaksacji

Metoda kolejnych nadrelaksacjiEdytuj

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  , 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  . W metodzie kolejnych nadrelaksacji obliczamy iteracyjnie dla wewnętrznych punktów kraty:

 

gdzie   jest pewnym współczynnikiem zależnym od obszaru. Dobre dobranie   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ę   z wektorem długości   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:   na przekątnej oraz   w czterech innych

kolumnach. Na przykład dla   ma postać:

 
 
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  :

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.