روش ژاکوبی — به زبان ساده (+ دانلود فیلم آموزش رایگان)

در آموزش‌های قبلی مجله فرادرس، با روش حذفی گاوس برای حل یک دستگاه با $$n$$ معادله خطی و $$n$$ متغیر مجهول به فرم $$ A x = b $$ آشنا شدیم. در این آموزش، با یکی از روش‌های عددی برای حل دستگاه معادلات خطی آشنا خواهیم شد. این الگوریتم تکراری، روش ژاکوبی (Jacobi method) نام دارد.

روش ژاکوبی

فرض کنید $$ x ^ { ( 0 ) } = \begin {bmatrix} x _ 1 ^ { ( 0 ) } \\ x _ 2 ^ { ( 0 ) } \\ \vdots \\ x _ n ^ { ( 0 ) } \end {bmatrix} $$ یک تقریب اولیه برای جواب $$x$$ در دستگاه $$n$$ معادله و $$n $$ مجهول زیر باشد:

$$ \large \begin {align} \quad \begin {matrix} E ( 1 ) : & a _ { 1 1 } x _ 1 & + & a _ { 1 2 } x _ 2 & + & \cdots & + & a _ { 1 n } x _ n & = & b _ 1 \\ E ( 2 ) : & a _ { 2 1 } x _ 1 & + & a _ { 2 2 } x _ 2 & + & \cdots & + & a _ { 2 n } x _ n & = & b _ 2 \\ \vdots & \vdots & & \vdots & & \ddots & & \vdots & & \vdots \\ E ( n ) : & a _ { n 1 } x _ 1 & + & a _ { n 2 } x _ 2 & + & \cdots & + & a _ { n n } x _ n & = & b _ n & \end {matrix} \end {align} $$

در روش تکرار ژاکوبی، از $$i$$اُمین معادله در دستگاه $$ A x = b $$ مقدار هر متغیر $$ x _ i $$ را به دست می‌آوریم. به عبارت دیگر، با فرض $$a_{ii} \neq 0$$، برای $$i = 1, 2, ..., n$$ داریم:

$$ \large \begin {align} x _ 1 = \frac { b _ 1 - \left [ a _ { 1 2 } x _ 2 + a _ { 1 3 } x _ 3 + ... + a _ { 1 n } x _ n \right ] }{ a_ { 1 1 } } \\ x _ 2 = \frac { b _ 2 - \left [ a _ { 2 1 } x _ 1 + a _ { 2 3 } x _ 3 + ... + a _ { 2 n } x _ n \right ] } { a _ { 2 2 } } \\ \quad \quad \quad \quad\quad \quad \vdots \quad \quad \quad \quad \quad \quad \\ x _ n = \frac { b _ n - \left [ a _ { n 1 } x _ 1 + a _ { n 2 } x _ 2 + ... + a _ { n , n - 1 } x _ { n - 1 } \right ] } { a _ { n n } } \end {align} $$

اگر از نماد سیگما استفاده کنیم، برای $$ i = 1, 2, ..., n$$ خواهیم داشت:

$$ \large \begin {align} \quad x _ i = \frac { b _ i - \sum _ { { j = 1 } _{j \neq i} } ^ { n } { a _ { i j } } x _ j} { a _ { i i } } \end {align} $$

برای به دست آوردن تقریب اول $$x^{(1)}$$ با استفاده از روش ژاکوبی معادلات را به صورت جدا نوشته و مقادیر تقریب اولیه $$x^{(0)} $$ را در آن‌ها وارد می‌کنیم:

$$ \large \begin {align} \quad x _ 1 ^ { ( 1) } = \frac { b _ 1 - \left [ a _ { 1 2 } x _ 2 ^ { ( 0 ) } + a _ { 1 3 } x _ 3 ^ { ( 0 ) } + ... + a _ { 1 n } x _ n ^ { ( 0 ) } \right ] } { a _ { 1 1 } } \\ x _ 2 ^ { ( 1 ) } = \frac { b _ 2 - \left [ a _ { 2 1 } x _ 1 ^ { ( 0 ) } + a _ { 2 3 } x _ 3 ^ { ( 0 ) } + ... + a _{ 2 n } x _ n ^ { ( 0 ) } \right ] } { a _ { 2 2 } } \\ \quad \quad \quad \quad \quad \quad \vdots \quad \quad \quad \quad \quad \quad \\ x _ n ^ { ( 1 ) } = \frac { b _ n - \left [ a _ { n 1 } x _ 1 ^ { ( 0 ) } + a _ { n 2 } x _ 2 ^ { ( 0 ) } + ... + a _ { n , n- 1 } x _ { n - 1 } ^ { ( 0 )} \right ] } { a _ { n n } } \end {align} $$

اگر از نماد سیگما برای $$i = 1, 2, ..., n$$ استفاده کنیم، خواهیم داشت:

$$ \large \begin {align} \quad x _ i ^ { ( 1 ) } = \frac { b _ i - \sum _ { j = 1 _ { j \neq i }} ^ { n } a _ { i j } x _ j ^ { ( 0 ) } } { a _ { i i } } \end {align} $$

سپس، به طریق مشابه تقریب $$ x^{(1)} $$، می‌توانیم تقریب $$x^{(2)}$$ را بنویسیم و همین‌گونه ادامه دهیم تا تقریب به جواب واقعی $$ x $$ در $$ A x = b $$ میل کند. برای هر $$k ≥ 1$$ و $$i = 1, 2, ..., n$$، خواهیم داشت:

$$ \large \begin {align} \quad x _ i ^ { ( k ) } = \frac { b _ i - \sum _ { j = 1 _{ j \neq i }} ^ { n } a _ { i j } x _ j ^ { ( k - 1) } } { a _ { i i } } \end {align} $$

به بیان ماتریسی، دستگاه معادلات خطی زیر را در نظر بگیرید:

$$ \large A \mathbf x = \mathbf b $$

که شامل $$n$$ معادله و $$n$$ مجهول است و در آن:

$$ \large A = \begin {bmatrix} a _ { 1 1 } & a _ { 1 2 } & \cdots & a _ { 1 n } \\ a _ { 2 1 } & a _ { 2 2 } & \cdots & a _ { 2 n } \\ \vdots & \vdots & \ddots & \vdots \\a _ { n 1 } & a _ { n 2 } & \cdots & a _ { n n } \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 $$ را می‌توان به ماتریس قطری $$D$$ و ماتریس باقیمانده $$R$$ را به صورت $$A = D + R $$ تجزیه کرد:

$$ \large D = \begin {bmatrix} a _ { 1 1 } & 0 & \cdots & 0 \\ 0 & a _ { 2 2 } & \cdots & 0 \\ \vdots & \vdots & \ddots & \vdots \\0 & 0 & \cdots & a _ { n n } \end {bmatrix} \text{, } \;\;\;\; R = \begin {bmatrix} 0 & a _ { 1 2 } & \cdots & a _ { 1 n } \\ a _ { 2 1 } & 0 & \cdots & a _ { 2 n } \\ \vdots & \vdots & \ddots & \vdots \\ a _ { n 1 } & a _ { n 2 } & \cdots & 0 \end {bmatrix} . $$

در نتیجه، جواب تکراری با معادله زیر به دست می‌آید:

$$ \large \mathbf { x } ^ { ( k + 1 ) } = D ^ { - 1 } ( \mathbf { b } - R \mathbf { x } ^ { ( k ) } ) $$

که در آن، $$ \mathbf{x}^{(k)} $$ همان $$ k$$اُمین تقریب یا تکرار از $$\mathbf{x}$$ و $$\mathbf{x}^{(k+1)}$$ تکرار $$k + 1 $$ است.

همگرایی روش ژاکوبی

شرط همگرایی استاندارد (برای هر روش تکراری) وقتی تحقق می‌یابد که شعاع طیفی ماتریس تکرار کمتر از ۱ باشد:

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

شعاع طیفی: اگر $$\lambda _ 1 $$، .... و $$\lambda _ n $$ مقادیر ویژه ماتریس $$A$$ باشند، شعاع طیفی این ماتریس با $$ \rho (A) \!$$ نشان داده شده و به صورت زیر بیان می‌شود:

$$ \large \rho(A) = \max_i(|\lambda_i|) \! $$

یک شرط لازم (و نه کافی) برای همگرایی این روش، این است که ماتریس $$ A$$ اکیداً غالب قطری باشد. اکیداً غالب قطری سطری، به این معناست که برای هر سطر، قدر مطلق جمله قطری بزرگ‌تر از مجموع قدر مطلق‌های سایر جملات باشد:

$$ \large \left | a _ { i i } \right | > \sum _ { j \ne i } { \left | a _ { i j } \right | } . $$

البته حتی اگر این شرایط نیز برقرار نباشد، گاهی روش ژاکوبی همگرا خواهد شد.

لازم به ذکر است که روش ژاکوبی برای هر ماتریس مثبت معین متقارن همگرا نیست. برای مثال:

$$ \large A =
\begin {pmatrix}
2 9 & 2 & 1 \\
2 & 6 & 1 \\
1 & 1 & \frac { 1 } { 5 }
\end {pmatrix}
\quad \Rightarrow \quad
D ^ { - 1 } R =
\begin {pmatrix}
0 & \frac { 2 } { 2 9 } & \frac { 1 } { 2 9 } \\
\frac { 1 } { 3 } & 0 & \frac { 1 } { 6 } \\
5 & 5 & 0
\end {pmatrix}
\quad \Rightarrow \quad
\rho ( D ^ { - 1 } R ) \approx 1 . 0 6 6 1 \, . $$

مثال‌هایی از روش ژاکوبی

در این بخش، چند مثال را بررسی می‌کنیم.

مثال ۱

جواب دستگاه سه معادله و سه مجهول زیر را با استفاده از روش ژاکوبی بیابید:

$$ \large \begin {align} \quad \begin {matrix} E ( 1 ) : & 3 x _ 1 & + & x _ 2 & + & 4 x _ 3 & = & 1 \\ E ( 2 ) : & x _ 1 & + & 5 x _ 2 & + & 9 x _ 3 & = & 1 \\ E ( 3 ) : & 2 x _ 1 & + & 6 x _ 2 & + & 5 x _ 3 & = & 1 \end {matrix} \end {align} $$

حل: تقریب اولیه $$ x ^ { ( 0 ) } = \begin {bmatrix} 0 . 2 2 \\ 0 . 0 3 \\ 0 . 0 6 \end {bmatrix} $$ را در نظر می‌گیریم. ابتدا متغیرهای $$ x _ i $$ را برای $$ i = 1 , 2 , 3 $$ به صورت زیر می‌نویسیم:

$$ \large \begin {align} \quad x _ 1 = \frac { 1 - x _ 2 - 4 x _ 3 } { 3 } \\ \quad x _ 2 = \frac { 1 - x _ 1 - 9 x _ 3 } { 5 } \\ \quad x _ 3 = \frac { 1 - 2 x _ 1 - 6 x _ 2 } { 5} \end {align} $$

تقریب نخست $$ x^{(1)} $$ با استفاده از روش ژاکوبی برابر است با:

$$ \large \begin {align} \quad x _ 1 ^ { ( 1 ) } = \frac { 1 - ( 0 . 0 3 ) - 4 ( 0 . 0 6 ) } { 3 } = 0 . 2 4 3 \\ \quad x _ 2 ^ { ( 1 ) } = \frac { 1 - ( 0 . 2 2 ) - 9 ( 0 .0 6 ) } { 5 } = 0 . 0 4 8 \\ \quad x _ 3 ^ { ( 1 ) } = \frac { 1 - 2 ( 0 . 2 2 ) - 6 ( 0 . 0 3 ) } { 5 } = 0 . 076 \end {align} $$

اگر همین‌گونه مراحل را ادامه دهیم، به تقریب‌های مناسبی خواهیم رسید. جواب واقعی این دستگاه $$ x = \begin {bmatrix} 0 . 2 \bar { 3 } \\ 0 . 0 \bar { 3 } \\ 0 . 0 \bar { 6 } \end {bmatrix} $$ است.

مثال 2

دستگاه معادلات خطی $$ A x = b $$ با تقریب اولیه $$ x ^ (0)$$ به صورت زیر است:

$$ \large A =
\begin {bmatrix}
2 & 1 \\
5 & 7 \\
\end {bmatrix},
\ b =
\begin {bmatrix}
1 1 \\
1 3 \\
\end {bmatrix}
\ \text {,} \quad x ^ { ( 0 ) } =
\begin {bmatrix}
1 \\
1 \\
\end {bmatrix} . $$

از معادله $$ x ^ { ( k + 1 ) } = D ^ {- 1 } ( b - R x ^ { ( k ) } ) $$ که در بالا آن را معرفی کردیم، مقدار $$ x $$ را تخمین می‌زنیم. ابتدا معادله را کمی ساده‌تر و به فرم $$ D ^ { - 1 } ( b - R x ^ {( k ) } ) = T x ^ { ( k ) } + C $$ می‌نویسیم که در آن، $$T=-D^{-1}R$$ و $$ C = D ^ {-1} b $$. همچنین، $$ R = L + U$$ که $$L$$ و $$U$$، به ترتیب، ماتریس‌های پایین مثلثی و بالامثلثی ماتریس $$A$$ هستند. با توجه به مقادیر مسئله، داریم:

$$ \large D ^ { - 1 } =
\begin {bmatrix}
1 / 2 & 0 \\
0 & 1 / 7 \\
\end {bmatrix} ,
\ L =
\begin {bmatrix}
0 & 0 \\
5 & 0 \\
\end {bmatrix}
, \quad U =
\begin {bmatrix}
0 & 1 \\
0 & 0 \\
\end {bmatrix} . $$

ماتریس $$ T = - D ^ {-1} ( L + U ) $$ را به صورت زیر به دست می‌آوریم:

$$ \large 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} . $$

برای $$C$$ نیز، داریم:

$$ \large C =
\begin {bmatrix}
1 / 2 & 0 \\
0 & 1 / 7 \\
\end {bmatrix}
\begin {bmatrix}
1 1 \\
1 3 \\
\end {bmatrix}
=
\begin {bmatrix}
1 1 / 2 \\
1 3 / 7 \\
\end {bmatrix} . $$

اکنون، با داشتن مقادیرِ محاسبه شده $$T$$ و $$C$$ مقدار $$x$$ را به صورت $$ x ^ {(1)} = T x ^ { ( 0 ) } + C $$ تقریب می‌زنیم:

$$ \large x ^ { ( 1 ) } =
\begin {bmatrix}
0 & - 1 / 2 \\
- 5 / 7 & 0 \\
\end {bmatrix}
\begin {bmatrix}
1 \\
1 \\
\end{bmatrix}
+
\begin {bmatrix}
1 1 / 2 \\
1 3 / 7 \\
\end {bmatrix}
=
\begin {bmatrix}
5 . 0 \\
8 / 7 \\
\end {bmatrix}
\approx
\begin {bmatrix}
5 \\
1 . 1 4 3 \\
\end {bmatrix} . $$

نتیجه تکرار بعدی نیز به صورت زیر است:

$$ \large x ^ { ( 2 ) } =
\begin {bmatrix}
0 & - 1 / 2 \\
- 5 / 7 & 0 \\
\end {bmatrix}
\begin {bmatrix}
5 . 0 \\
8 / 7 \\
\end {bmatrix}
+
\begin {bmatrix}
1 1 / 2 \\
1 3 / 7 \\
\end {bmatrix}
=
\begin {bmatrix}
6 9 / 1 4 \\
- 1 2 / 7 \\
\end {bmatrix}
\approx
\begin {bmatrix}
4 . 9 2 9 \\
- 1 . 7 1 4 \\
\end {bmatrix} . $$

این روال را تا جایی ادامه می‌دهیم که شرط همگرایی $$\|Ax^{(n)} - b\|$$ کوچک شود. جواب بعد از ۲۵ تکرار به صورت زیر است:

$$ \large x = \begin {bmatrix}
7 . 1 1 1 \\
- 3 . 2 2 2
\end {bmatrix} . $$

 مثال ۳

سیستم زیر را در نظر بگیرید:

$$ \large \begin {align}
1 0 x _ 1 - x _ 2 + 2 x _ 3 & = 6 , \\
- x _ 1 + 1 1 x _ 2 - x _ 3 + 3 x _ 4 & = 2 5 , \\
2 x _ 1 - x _ 2 + 1 0 x _ 3 - x _ 4 & = - 1 1 , \\
3 x _ 2 - x _ 3 + 8 x _ 4 & = 1 5 .
\end {align} $$

اگر تقریب اولیه $$ (0,0,0,0)$$ را در نظر بگیریم، جواب تقریبی نخست به صورت زیر است:

$$ \large \begin {align}
x _ 1 & = ( 6 + 0 - ( 2 * 0 ) ) / 1 0 = 0 . 6 , \\
x _ 2 & = ( 2 5 + 0 + 0 - ( 3 * 0 ) ) / 1 1 = 2 5 / 1 1 = 2 . 2 7 2 7 , \\
x _ 3 & = ( - 1 1 - ( 2 * 0 ) + 0 + 0 ) / 1 0 = - 1 . 1 , \\
x _ 4 & = ( 1 5 - ( 3 * 0 ) + 0 ) / 8 = 1 . 8 7 5 .
\end {align} $$

جواب‌های مربوط به ۵ تکرار در جدول زیر آورده شده است:

برنامه پایتون روش ژاکوبی برای این مثال به صورت زیر است:

1import numpy as np
2
3ITERATION_LIMIT = 1000
4
5# initialize the matrix
6A = np.array([[10., -1., 2., 0.],
7              [-1., 11., -1., 3.],
8              [2., -1., 10., -1.],
9              [0.0, 3., -1., 8.]])
10# initialize the RHS vector
11b = np.array([6., 25., -11., 15.])
12
13# prints the system
14print("System:")
15for i in range(A.shape[0]):
16    row = ["{}*x{}".format(A[i, j], j + 1) for j in range(A.shape[1])]
17    print(" + ".join(row), "=", b[i])
18print()
19
20x = np.zeros_like(b)
21for it_count in range(ITERATION_LIMIT):
22    print("Current solution:", x)
23    x_new = np.zeros_like(x)
24
25    for i in range(A.shape[0]):
26        s1 = np.dot(A[i, :i], x[:i])
27        s2 = np.dot(A[i, i + 1:], x[i + 1:])
28        x_new[i] = (b[i] - s1 - s2) / A[i, i]
29
30    if np.allclose(x, x_new, atol=1e-10, rtol=0.):
31        break
32
33    x = x_new
34
35print("Solution:")
36print(x)
37error = np.dot(A, x) - b
38print("Error:")
39print(error)

اگر برنامه بالا را اجرا کنیم، خواهیم داشت:

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]

پیاده‌سازی روش ژاکوبی در متلب

برنامه روش ژاکوبی در متلب به صورت زیر است (توضیحات مربوط به ورودی‌ها و خروجی‌ها در برنامه ذکر شده است):

1function [x,dx] = jacobi(A,b,x,eps,N)
2%
3%  The function jacobi applies Jacobi's method to solve A*x = b.
4%  On entry:
5%    A       coefficient matrix of the linear system;
6%    b       right-hand side vector of the linear system;
7%    x       initial approximation;
8%    eps     accuracy requirement: stop when norm(dx) < eps;
9%    N       maximal number of steps allowed.
10%  On return:
11%    x       approximation for the solution of A*x = b;
12%    dx      vector last used to update x,
13%            if success, then norm(dx) < eps.
14%
15%  Example:
16%    Q = orth(rand(5)); D = diag(rand(1,5));
17%    A = Q*D*Q'; z = rand(5,1); b = A*z;
18%    x = z + rand(5,1)*1e-4;
19%    [x,dx] = jacobi(A,b,x,1e-8,50)
20%
21n = size(A,1);
22fprintf('Running the method of Jacobi...\n');
23for i = 1:N
24   dx = b - A*x;
25   for j = 1:n
26      dx(j) = dx(j)/A(j,j);
27      x(j) = x(j) + dx(j);
28   end;
29   fprintf('  norm(dx) = %.4e\n', norm(dx));
30   if (norm(dx) < eps)
31      fprintf('Succeeded in %d steps\n', i);
32      return;
33   end;
34end;
35fprintf('Failed to reached accuracy requirement in %d steps.\n', N);

اگر این مطلب برایتان مفید بوده است، آموزش‌های زیر نیز به شما پیشنهاد می‌شوند:

• مجموعه آموزش‌های محاسبات عددی
• آموزش محاسبات عددی با MATLAB
• مجموعه آموزش‌های دروس ریاضیات
• آموزش ریاضی پایه دانشگاهی
• معادله خط — به زبان ساده
• دستگاه معادلات دیفرانسیل خطی — به زبان ساده
• جذر یا محاسبه ریشه دوم عدد — به زبان ساده

^^