Jacobi metodu

Jacobi metodu (diğer adıyla Jacobi yinelemeli metodu^[1]), sayısal lineer cebirde lineer denklemlerin diyagonal olarak baskın sistemlerin çözümlerinin belirlenmesi için oluşturulmuş bir algoritmadır. Her diyagonal eleman tek tek çözülür ve yaklaşık bir değer olarak alınır. Bu aşama onlar yakınsayana kadar tekrarlanır. Bu algoritma matris köşegenleştirilmesi Jacobi dönüşüm metodunun (diğer adıyla Jacobi özdeğer algoritmasının) sadeleştirilmiş şeklidir. Bu metot daha sonra Carl Gustav Jacob Jacobi olarak isimlendirilmiştir.

Açıklama

n lineer denklemli bir kare sistemi verilmiş olsun.

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

$A={\begin{bmatrix}a_{11}&a_{12}&\cdots &a_{1n}\\a_{21}&a_{22}&\cdots &a_{2n}\\\vdots &\vdots &\ddots &\vdots \\a_{n1}&a_{n2}&\cdots &a_{nn}\end{bmatrix}},\qquad \mathbf {x} ={\begin{bmatrix}x_{1}\\x_{2}\\\vdots \\x_{n}\end{bmatrix}},\qquad \mathbf {b} ={\begin{bmatrix}b_{1}\\b_{2}\\\vdots \\b_{n}\end{bmatrix}}.$ A matrisi diyagonal (D) ve geriye kalan (R) bileşenlerine ayrılır.

A=D+R\qquad {\text{where}}\qquad D={\begin{bmatrix}a_{11}&0&\cdots &0\\0&a_{22}&\cdots &0\\\vdots &\vdots &\ddots &\vdots \\0&0&\cdots &a_{nn}\end{bmatrix}}{\text{ and }}R={\begin{bmatrix}0&a_{12}&\cdots &a_{1n}\\a_{21}&0&\cdots &a_{2n}\\\vdots &\vdots &\ddots &\vdots \\a_{n1}&a_{n2}&\cdots &0\end{bmatrix}}.

Bu çözüm daha sonra : $\mathbf {x} ^{(k+1)}=D^{-1}(\mathbf {b} -R\mathbf {x} ^{(k)}).$ aracılığıyla yinelemeli olarak elde edilir. Eleman tabanlı formül böylece : $x_{i}^{(k+1)}={\frac {1}{a_{ii}}}\left(b_{i}-\sum _{j\neq i}a_{ij}x_{j}^{(k)}\right),\quad i=1,2,\ldots ,n.$ olur. xi(k+1) hesaplanması kendisi dışında x(k)'daki her elemanı gerektirir. Gauss-Seidel yönteminin aksine, xi(k) ile xi(k+1)'yi, hesaplamanın geri kalanı tarafından ihtiyaç duyulacak değer olarak üzerine yazamayız. Minimum saklama miktarı, büyüklüğü n olan iki vektördür.

Algoritma

x^{(0)}

için bir başlangıç tahmini seçelim.

k=0

while (yakınsama ulaşamamış) do

for (i:=1; i<n) do

\sigma =0

for (j:=1; j<n) do

if ( j≠i) then

\sigma =\sigma +a_{ij}x_{j}^{(k)}

end if

end (j döngüsü)

check (yakınsama ulaşmış?)

k=k+1

end (while döngüsü)

Yakınsama

Standart yakınsama durumu (her tekrarlamalı metot için) yineleme matrisinin spektral yarıçapı en az bir olmasıdır.

\rho (D^{-1}R)<1.

A matrisi kesin ve indirgenemez bir şekilde diyagonal olarak baskın ise bu yöntem yakınsama garantilidir. Kesin olan satırın diyagonal baskınlığı, her satır için diyagonal terimin mutlak değerinin, diğer terimlerin mutlak değerlerinin toplamından büyük olmasıdır.

\left|a_{ii}\right|>\sum _{j\neq i}{\left|a_{ij}\right|}.

Jacobi metodu bazen bu şartlar sağlanmasa bile yakınsar.

Örnek

$x^{(0)}$ başlangıç değerli $Ax=b$ formundaki lineer sisteminde

A={\begin{bmatrix}2&1\\5&7\\\end{bmatrix}},\ b={\begin{bmatrix}11\\13\\\end{bmatrix}}\quad {\text{and}}\quad x^{(0)}={\begin{bmatrix}1\\1\\\end{bmatrix}}.

olarak veriliyor. X' i tahmin etmek için yukarıda tanımlanmış $x^{(k+1)}=D^{-1}(b-Rx^{(k)})$ , denklemini kullanırız. İlk olarak; denklemi $T=-D^{-1}R$ ve $C=D^{-1}b$ olduğu yerde daha uygun formda yazarız: $D^{-1}(b-Rx^{(k)})=Tx^{(k)}+C$ L ve U, A matrisinin kesin alt ve üst kısımları olduğu yerde $R=L+U$ 'dur. Bilinen değerlerden; : $D^{-1}={\begin{bmatrix}1/2&0\\0&1/7\\\end{bmatrix}},\ L={\begin{bmatrix}0&0\\5&0\\\end{bmatrix}}\quad {\text{and}}\quad U={\begin{bmatrix}0&1\\0&0\\\end{bmatrix}}.$ $T=-D^{-1}(L+U)$ olarak tanımlayabiliriz.

T={\begin{bmatrix}1/2&0\\0&1/7\\\end{bmatrix}}\left\{{\begin{bmatrix}0&0\\-5&0\\\end{bmatrix}}+{\begin{bmatrix}0&-1\\0&0\\\end{bmatrix}}\right\}={\begin{bmatrix}0&-1/2\\-5/7&0\\\end{bmatrix}}.

Ayrıca, : $C={\begin{bmatrix}1/2&0\\0&1/7\\\end{bmatrix}}{\begin{bmatrix}11\\13\\\end{bmatrix}}={\begin{bmatrix}11/2\\13/7\\\end{bmatrix}}.$ olarak buluruz. Hesaplanan T ve C ile x'i $x$ as $x^{(1)}=Tx^{(0)}+C$ : olarak tahmin ederiz.

x^{(1)}={\begin{bmatrix}0&-1/2\\-5/7&0\\\end{bmatrix}}{\begin{bmatrix}1\\1\\\end{bmatrix}}+{\begin{bmatrix}11/2\\13/7\\\end{bmatrix}}={\begin{bmatrix}5.0\\8/7\\\end{bmatrix}}\approx {\begin{bmatrix}5\\1.143\\\end{bmatrix}}.

Bir sonraki yineleme bize

x^{(2)}={\begin{bmatrix}0&-1/2\\-5/7&0\\\end{bmatrix}}{\begin{bmatrix}5.0\\8/7\\\end{bmatrix}}+{\begin{bmatrix}11/2\\13/7\\\end{bmatrix}}={\begin{bmatrix}69/14\\-12/7\\\end{bmatrix}}\approx {\begin{bmatrix}4.929\\-1.714\\\end{bmatrix}}.

verir. Bu aşama yakınsayana kadar tekrarlanır ( $\|Ax^{(n)}-b\|$ küçük değer alana kadar). 25 yineleme sonra çözüm:

x={\begin{bmatrix}7.111\\-3.222\end{bmatrix}}.

Başka bir örnek

Aşağıdaki lineer sistemin verildiğini varsayalım:

{\begin{aligned}10x_{1}-x_{2}+2x_{3}&=6,\\-x_{1}+11x_{2}-x_{3}+3x_{4}&=25,\\2x_{1}-x_{2}+10x_{3}-x_{4}&=-11,\\3x_{2}-x_{3}+8x_{4}&=15.\end{aligned}}

Başlangıç tahmini olarak (0, 0, 0, 0) alalım, ardından ilk tahminimiz:

{\begin{aligned}x_{1}&=(6-0-0)/10=0.6,\\x_{2}&=(25-0-0)/11=25/11=2.2727\\x_{3}&=(-11-0-0)/10=-1.1,\\x_{4}&=(15-0-0)/8=1.875.\end{aligned}}

olacaktır. Elde edilen tahminleri kullanarak, istenilen doğruluğa ulaşıncaya kadar yineleme işlemi tekrarlanır. Beş iterasyon sonra yaklaşık çözümler şunlardır:

$x_{1}$	$x_{2}$	$x_{3}$	$x_{4}$
$0.6$	$2.27272$	$-1.1$	$1.875$
$1.04727$	$1.7159$	$-0.80522$	$0.88522$
$0.93263$	$2.05330$	$-1.0493$	$1.13088$
$1.01519$	$1.95369$	$-0.9681$	$0.97384$
$0.98899$	$2.0114$	$-1.0102$	$1.02135$

Sistemin tam çözümü ise (1, 2, −1, 1) olur.

Python 3 ve Numpy kullanılarak yapılmış bir örnek

Aşağıdaki sayısal işlem çözüm vektörü oluşturmak için yinelenir.

import numpy as np

ITERATION_LIMIT = 1000

# initialize the matrix
A = np.array([[10., -1., 2., 0.],
              [-1., 11., -1., 3.],
              [2., -1., 10., -1.],
              [0.0, 3., -1., 8.]])
# initialize the RHS vector
b = np.array([6., 25., -11., 15.])

# prints the system
print("System:")
for i in range(A.shape[0]):
    row = ["{}*x{}".format(A[i, j], j + 1) for j in range(A.shape[1])]
    print(" + ".join(row), "=", b[i])
print()

x = np.zeros_like(b)
for it_count in range(ITERATION_LIMIT):
    print("Current solution:", x)
    x_new = np.zeros_like(x)

    for i in range(A.shape[0]):
        s1 = np.dot(A[i, :i], x[:i])
        s2 = np.dot(A[i, i + 1:], x[i + 1:])
        x_new[i] = (b[i] - s1 - s2) / A[i, i]

    if np.allclose(x, x_new, rtol=1e-10):
        break

    x = x_new

print("Solution:")
print(x)
error = np.dot(A, x) - b
print("Error:")
print(error)

Aşağıdaki çıktıyı üretir:

System:
10.0*x1 + -1.0*x2 + 2.0*x3 + 0.0*x4 = 6.0
-1.0*x1 + 11.0*x2 + -1.0*x3 + 3.0*x4 = 25.0
2.0*x1 + -1.0*x2 + 10.0*x3 + -1.0*x4 = -11.0
0.0*x1 + 3.0*x2 + -1.0*x3 + 8.0*x4 = 15.0

Current solution: [ 0.  0.  0.  0.]
Current solution: [ 0.6         2.27272727 -1.1         1.875     ]
Current solution: [ 1.04727273  1.71590909 -0.80522727  0.88522727]
Current solution: [ 0.93263636  2.05330579 -1.04934091  1.13088068]
Current solution: [ 1.01519876  1.95369576 -0.96810863  0.97384272]
Current solution: [ 0.9889913   2.01141473 -1.0102859   1.02135051]
Current solution: [ 1.00319865  1.99224126 -0.99452174  0.99443374]
Current solution: [ 0.99812847  2.00230688 -1.00197223  1.00359431]
Current solution: [ 1.00062513  1.9986703  -0.99903558  0.99888839]
Current solution: [ 0.99967415  2.00044767 -1.00036916  1.00061919]
Current solution: [ 1.0001186   1.99976795 -0.99982814  0.99978598]
Current solution: [ 0.99994242  2.00008477 -1.00006833  1.0001085 ]
Current solution: [ 1.00002214  1.99995896 -0.99996916  0.99995967]
Current solution: [ 0.99998973  2.00001582 -1.00001257  1.00001924]
Current solution: [ 1.00000409  1.99999268 -0.99999444  0.9999925 ]
Current solution: [ 0.99999816  2.00000292 -1.0000023   1.00000344]
Current solution: [ 1.00000075  1.99999868 -0.99999899  0.99999862]
Current solution: [ 0.99999967  2.00000054 -1.00000042  1.00000062]
Current solution: [ 1.00000014  1.99999976 -0.99999982  0.99999975]
Current solution: [ 0.99999994  2.0000001  -1.00000008  1.00000011]
Current solution: [ 1.00000003  1.99999996 -0.99999997  0.99999995]
Current solution: [ 0.99999999  2.00000002 -1.00000001  1.00000002]
Current solution: [ 1.          1.99999999 -0.99999999  0.99999999]
Current solution: [ 1.  2. -1.  1.]
Solution:
[ 1.  2. -1.  1.]
Error:
[ -2.81440107e-08   5.15706873e-08  -3.63466359e-08   4.17092547e-08]

Ağırlıklı Jacobi yöntemi

Ağırlıklı Jacobi yinelemesi : $\mathbf {x} ^{(k+1)}=\omega D^{-1}(\mathbf {b} -R\mathbf {x} ^{(k)})+\left(1-\omega \right)\mathbf {x} ^{(k)}$ yinelemesini olağan seçenek olan $\omega =2/3$ ile hesaplarkenw parametresini kullanır.^[2]

Son gelişmeler

2014 yılında planlanmış relaksasyon Jacobi yöntemi adıyla anılan bir algoritma düzeltmesi yayımlandı.^[1]^[3] Bu yeni yöntem bir üst ve alt relaksasyon programı kullanır ve büyük iki ve üç boyutlu Kartezyen örgülerle ayrıştırılmış eliptik denklemlerin çözümünde iki yüz kat performans artışı sağlar.

Kaynakça

^ ^a ^b "19th century math tactic gets a makeover—and yields answers up to 200 times faster". Douglas, Man Adası: Phys.org. 30 Haziran 2014. 6 Temmuz 2014 tarihinde kaynağından arşivlendi. Erişim tarihi: 1 Temmuz 2014.
^ Saad, Yousef (2003). Iterative Methods for Sparse Linear Systems (2 bas.). SIAM. s. 414. ISBN 0898715342.
^ Yang, Xiang; Mittal, Rajat (27 Haziran 2014). "Acceleration of the Jacobi iterative method by factors exceeding 100 using scheduled relaxation". Journal of Computational Physics. doi:10.1016/j.jcp.2014.06.010.

[JHU_Jacobi-1] "19th century math tactic gets a makeover—and yields answers up to 200 times faster". Douglas, Man Adası: Phys.org. 30 Haziran 2014. 6 Temmuz 2014 tarihinde kaynağından arşivlendi. Erişim tarihi: 1 Temmuz 2014.

[2] Saad, Yousef (2003). Iterative Methods for Sparse Linear Systems (2 bas.). SIAM. s. 414. ISBN 0898715342.

[3] Yang, Xiang; Mittal, Rajat (27 Haziran 2014). "Acceleration of the Jacobi iterative method by factors exceeding 100 using scheduled relaxation". Journal of Computational Physics. doi:10.1016/j.jcp.2014.06.010.

[1]

[2]

[3]