Apakah ada aturan sederhana untuk mengatakan apakah itu layak untuk dilakukan SOR daripada Gauss-Seidel? (dan cara yang mungkin bagaimana memperkirakan parameter realxation )$\omega$

Maksud saya hanya dengan melihat pada matriks , atau pengetahuan tentang masalah tertentu yang diwakili oleh matriks?

Saya membaca jawaban atas pertanyaan ini: Apakah ada heuristik untuk mengoptimalkan metode over-relaxation (SOR) yang berurutan? tapi agak terlalu canggih. Saya tidak melihat heuristik sederhana bagaimana memperkirakan jari-jari spektral hanya melihat pada matriks (atau masalah yang diwakilinya).

Saya ingin sesuatu yang jauh lebih sederhana - hanya beberapa contoh matriks (masalah) yang SOR bertemu lebih cepat.

---

Saya sedang bereksperimen dengan SOR untuk matriks raja ini: mana saya adalah matriks identitas, C i j = c ∀ i , j dan R i j s adalah angka acak dari distribusi yang sama sehingga | R i j | < r . Saya berpikir bahwa akan ada beberapa ketergantungan ω optimal pada parameter c , r .$A = I + C + R$$I$$C_{ij}=c$ $\forall i,j$$R_{ij}$$|R_{ij}|<r$$\omega$$c,r$

EDIT: Saya menggunakan sangat kecil , r untuk memastikan bahwa A sangat dominan secara diagonal. ( | c | < 0,1 , r < 2 | c | untuk matriks dimensi 5-10). Saya juga harus mengatakan bahwa A ini nyata dan simetris.$c,r$$A$$|c|<0.1$$r<2|c|$$A$

Namun, saya menemukan bahwa Gauss-Seidel ( ) hampir selalu yang terbaik (?)$\omega=1$ . Apakah ini berarti harus ada korelasi lagi antara untuk mendapatkan keuntungan dari SOR? Atau saya melakukan sesuatu yang salah? $A_{ij}$

---

Saya tahu, bahwa SOR bukan pemecah yang paling efisien (dibandingkan dengan CG, GMRES ...) tetapi sederhana untuk mengimplementasikan dan membuat parsial, dan memodifikasi untuk masalah tertentu. Bagus untuk membuat prototipe.

Konvergensi dari pemecah iteratif klasik untuk sistem linier ditentukan oleh jari-jari spektral dari matriks iterasi, . Untuk sistem linear umum, sulit untuk menentukan parameter SOR yang optimal (atau bahkan bagus) karena kesulitan dalam menentukan jari-jari spektral dari matriks iterasi. Di bawah ini saya telah memasukkan banyak detail tambahan, termasuk contoh masalah nyata di mana berat SOR optimal diketahui.$\rho(\mathbf{G})$

Jari-jari spektral dan konvergensi

Jari-jari spektral didefinisikan sebagai nilai absolut dari nilai eigen magnitudo terbesar. Suatu metode akan konvergen jika dan radius spektral yang lebih kecil berarti konvergensi yang lebih cepat. SOR bekerja dengan mengubah pemisahan matriks yang digunakan untuk menurunkan matriks iterasi berdasarkan pilihan parameter bobot ω , semoga mengurangi jari-jari spektral dari matriks iterasi yang dihasilkan.$\rho<1$$\omega$

Pemecahan matriks

Untuk pembahasan di bawah ini, saya akan menganggap bahwa sistem yang akan dipecahkan diberikan oleh

$$\mathbf{A}x=b,$$

dengan iterasi formulir

$$x^{(k+1)}=v+\mathbf{G}x^{(k)},$$

di mana adalah vektor, dan angka iterasi k dilambangkan x ( k ) .$v$$k$$x^{(k)}$

SOR mengambil rata-rata tertimbang dari iterasi lama dan iterasi Gauss-Seidel. Metode Gauss-Seidel bergantung pada pemisahan matriks formulir

$$\mathbf{A} = \mathbf{D} + \mathbf{L} + \mathbf{U}$$

di mana adalah diagonal A , L adalah matriks segitiga bawah yang berisi semua elemen A tepat di bawah diagonal dan R adalah matriks segitiga atas yang mengandung semua elemen A secara ketat di atas diagonal. Iterasi Gauss-Seidel kemudian diberikan oleh$\mathbf{D}$$\mathbf{A}$$\mathbf{L}$$\mathbf{A}$$\mathbf{R}$$\mathbf{A}$

$$x^{(k+1)}=(\mathbf{D}+\mathbf{L})^{-1}b+\mathbf{G}_{\rm G-S}x^{(k)}$$

dan matriks iterasi adalah

$$\mathbf{G}_{\rm G-S}=-(\mathbf{D}+\mathbf{L})^{-1}\mathbf{U}.$$

SOR kemudian dapat ditulis sebagai

$$x^{(k+1)}=\omega(\mathbf{D}+\omega\mathbf{L})^{-1}b+\mathbf{G}_{\rm SOR}x^{(k)}$$

dimana

$$\mathbf{G}_{\rm SOR} = (\mathbf{D}+\omega\mathbf{L})^{-1}((1-\omega)\mathbf{D}-\omega\mathbf{U}).$$

Menentukan tingkat konvergensi dari skema iteratif benar-benar bermuara pada menentukan jari-jari spektral dari matriks iterasi ini. Secara umum, ini adalah masalah yang sulit kecuali Anda tahu sesuatu yang spesifik tentang struktur matriks. Ada beberapa contoh yang saya tahu di mana koefisien pembobotan yang optimal dapat dihitung. Dalam praktiknya, harus ditentukan dengan cepat berdasarkan pada konvergensi algoritma berjalan yang diamati (diperkirakan). Ini berfungsi dalam beberapa kasus, tetapi gagal dalam kasus lain.$\omega$

SOR optimal

Satu contoh realistis di mana koefisien pembobotan optimal diketahui muncul dalam konteks penyelesaian persamaan Poisson:

$$\nabla^2 u = f~{\rm in}~\Omega\\ u = g~{\rm on}~\partial\Omega$$

Mendiskritkan sistem ini pada domain bujur sangkar dalam 2D ​​menggunakan perbedaan hingga orde kedua dengan spasi grid yang seragam menghasilkan matriks pita simetris dengan 4 pada diagonal, -1 tepat di atas dan di bawah diagonal, dan dua pita lagi -1 agak jauh dari diagonal. Ada beberapa perbedaan karena kondisi batas, tetapi itu adalah struktur dasar. Dengan matriks ini, pilihan optimal terbukti untuk koefisien SOR diberikan oleh

$$\omega = \dfrac{2}{1+\sin(\pi \Delta x/L)}$$

di mana adalah spasi grid dan L adalah ukuran domain. Melakukannya untuk kasus sederhana dengan solusi yang diketahui memberikan kesalahan berikut versus nomor iterasi untuk dua metode ini:$\Delta x$$L$

Seperti yang Anda lihat, SOR mencapai presisi mesin di sekitar 100 iterasi di mana titik Gauss-Seidel sekitar 25 kali lipat lebih buruk. Jika Anda ingin bermain-main dengan contoh ini, saya telah memasukkan kode MATLAB yang saya gunakan di bawah ini.

clear all
close all

%number of iterations:
niter = 150;

%number of grid points in each direction
N = 16;
% [x y] = ndgrid(linspace(0,1,N),linspace(0,1,N));
[x y] = ndgrid(linspace(-pi,pi,N),linspace(-pi,pi,N));
dx = x(2,1)-x(1,1);
L = x(N,1)-x(1,1);

%desired solution:
U = sin(x/2).*cos(y);

% Right hand side for the Poisson equation (computed from U to produce the
% desired known solution)
Ix = 2:N-1;
Iy = 2:N-1;
f = zeros(size(U));
f(Ix,Iy) = (-4*U(Ix,Iy)+U(Ix-1,Iy)+U(Ix+1,Iy)+U(Ix,Iy-1)+U(Ix,Iy+1));

figure(1)
clf
contourf(x,y,U,50,'linestyle','none')
title('True solution')

%initial guess (must match boundary conditions)
U0 = U;
U0(Ix,Iy) = rand(N-2);

%Gauss-Seidel iteration:
UGS = U0; EGS = zeros(1,niter);
for iter=1:niter
    for iy=2:N-1
        for ix=2:N-1
            UGS(ix,iy) = -1/4*(f(ix,iy)-UGS(ix-1,iy)-UGS(ix+1,iy)-UGS(ix,iy-1)-UGS(ix,iy+1));
        end
    end

    %error:
    EGS(iter) = sum(sum((U-UGS).^2))/sum(sum(U.^2));
end

figure(2)
clf
contourf(x,y,UGS,50,'linestyle','none')
title(sprintf('Gauss-Seidel approximate solution, iteration %d', iter))
drawnow

%SOR iteration:
USOR = U0; ESOR = zeros(1,niter);
w = 2/(1+sin(pi*dx/L));
for iter=1:niter
    for iy=2:N-1
        for ix=2:N-1
            USOR(ix,iy) = (1-w)*USOR(ix,iy)-w/4*(f(ix,iy)-USOR(ix-1,iy)-USOR(ix+1,iy)-USOR(ix,iy-1)-USOR(ix,iy+1));
        end
    end

    %error:
    ESOR(iter) = sum(sum((U-USOR).^2))/sum(sum(U.^2));
end

figure(4)
clf
contourf(x,y,USOR,50,'linestyle','none')
title(sprintf('Gauss-Seidel approximate solution, iteration %d', iter))
drawnow


figure(5)
clf
semilogy(EGS,'b')
hold on
semilogy(ESOR,'r')
title('L2 relative error')
xlabel('Iteration number')
legend('Gauss-Seidel','SOR','location','southwest')

Sisi ini sebenarnya bukan keahlian saya, tapi saya rasa ini bukan ujian yang sangat adil untuk banyak aplikasi realistis.

Saya tidak yakin nilai apa yang Anda gunakan untuk c dan r , tapi saya curiga Anda bekerja dengan matriks yang sangat buruk. (Di bawah ini adalah beberapa kode Python yang menunjukkan bahwa ini mungkin bukan matriks yang paling dapat dibalik.)

>>> import numpy
>>> for n in [100, 1000]:
...     for c in [.1, 1, 10]:
...         for r in [.1, 1, 10]:
...             print numpy.linalg.cond(
...                 numpy.eye(n) + 
...                 c * numpy.ones((n, n)) + 
...                 r * numpy.random.random((n, n))
...             )
... 
25.491634739
2034.47889101
2016.33059429
168.220149133
27340.0090644
5532.81258852
1617.33518781
42490.4410689
5326.3865534
6212.01580004
91910.8386417
50543.9269739
24737.6648458
271579.469212
208913.592289
275153.967337
17021788.5576
117365.924601

Jika Anda benar-benar perlu membalikkan matriks yang tidak dikondisikan ini, Anda akan a) menggunakan metode khusus, dan b) mungkin harus mencari bidang baru 😉

Untuk matriks yang dikondisikan dengan baik dari berbagai ukuran, SOR kemungkinan akan lebih cepat. Untuk masalah nyata di mana masalah kecepatan, akan jarang menggunakan SOR - di sisi canggih, ada jauh lebih baik hari ini; di sisi yang lambat namun dapat diandalkan, SOR bukan yang terbaik yang dapat Anda lakukan.

OK, jadi untuk matriks simetris raja ini:

1 t 0 0 0 0 t t 0 0 
t 1 0 0 0 0 0 0 0 0 
0 0 1 0 0 0 0 t 0 t 
0 0 0 1 0 0 0 0 0 t 
0 0 0 0 1 t 0 0 0 0 
0 0 0 0 t 1 0 t 0 0 
t 0 0 0 0 0 1 0 0 0 
t 0 t 0 0 t 0 1 0 0 
0 0 0 0 0 0 0 0 1 t 
0 0 t t 0 0 0 0 t 1 

$t$$t$$t$

$t_i = c + random(-r,r)$

$t$$c=0,r=0.1$$t$

(Ini hanya pengamatan empiris, tidak ada yang ketat)