אנליזה נומרית/פתרון מערכת משוואות לינאריות

מתוך ויקיספר, אוסף ספרי הלימוד והמדריכים החופשי.

בעיות לינאריות בדרך כלל קלות יותר לפתרון מאשר בעיות לא לינאריות. ניתן לגשת לבעיות אלו בצורה אנליטית, ואז הפתרון תלוי בדיוק המחשב, או בצורה נומרית, ואז הפתרון תלוי במידת התכנסות השיטה וגם בדיוק המחשב. השיטות הנומריות בדרך כלל מהירות יותר מהשיטות האנליטיות (אחרת לא היו משתמשים בהן).

מאחר ואנו מתעניינים בבעיות בעלות פתרון (יחיד), נתקל במטריצות ריבועיות בלבד, אשר מקיימות $\ |A|=det(A)\neq 0$ . לשם נוחות, נציג את אופני הכתיבה השונים של מערכת משוואות לינאריות:

$\left\{ \begin{array}{rcrcccrcl} a_{11}x_1 &+& a_{12}x_2 &+& \cdots &+& a_{1n}x_n &=& b_1 \\ a_{21}x_1 &+& a_{22}x_2 &+& \cdots &+& a_{2n}x_n &=& b_2 \\ &&&\vdots&&&&&\vdots \\ a_{n1}x_1 &+& a_{n2}x_2 &+& \cdots &+& a_{nn}x_n &=& b_n \end{array} \right. \quad\Rightarrow\quad \overbrace{\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}}^{\underline{A}} \overbrace{\begin{Bmatrix} x_1 \\ x_2 \\ \vdots \\ x_n \end{Bmatrix}}^{\underline{x}} = \overbrace{\begin{Bmatrix} b_1 \\ b_2 \\ \vdots \\ b_n \end{Bmatrix}}^{\underline{b}}$

$\ (\underline{A})_{ij}=a_{ij} \quad;\quad \underline{A} \cdot \underline{x}= \underline{b} \quad;\quad \sum_{j=1}^n a_{ij}x_j=b_i$

אנליזה נומרית

[עריכה] שיטות אנליטיות

בשיטות האנליטיות מבצעים פעולות שורה אלמנטריות על המטריצה לקבל הפתרון הרצוי, אשר תלוי בדיוק המחשב עליו מתבצע החישוב. נקודת התורפה המרכזית של השיטות האנליטיות היא גודל המערכת המירבי בו ניתן לטפל.

[עריכה] שיטת גאוס

בשיטת גאוס (הנקראת גם "שיטת הדרוג" או "שיטת החילוץ") מדרגים את המטריצה המורחבת [A|b] ובסוף מציבים לאחור. עבור מטריצה מסדר n יידרשו $\ {n^3 \over 3} +n^2 -{n\over 3}$ פעולות.

[עריכה] שיטת גאוס-ג'ורדן

עבור מטריצה מסדר n יידרשו $\ {n^3 \over 2} +n^2 -{n\over 2}$ פעולות.

[עריכה] פירוק LU

מאלגברה לינארית ידוע, כי כל מטרציה שניתן להביא אותה לצורה משולשת-עליונה על ידי פעולות שורה אלמנטריות, ניתן ליצג באמצעות מכפלת מטריצה תחתונה (L) במטריצה עליונה (U). נדגים זאת באמצעות מטריצה מסדר 3:

$\overbrace{\begin{bmatrix} a_{11} & a_{12} & a_{13} \\ a_{21} & a_{22} & a_{23} \\ a_{31} & a_{32} & a_{33} \\ \end{bmatrix}}^{\underline{A}} = \overbrace{\begin{bmatrix} l_{11} & 0 & 0 \\ l_{12} & l_{22} & 0 \\ l_{13} & l_{23} & l_{33} \\ \end{bmatrix}}^{\underline{L}} \overbrace{\begin{bmatrix} 1 & u_{12} & u_{13} \\ 0 & 1 & u_{23} \\ 0 & 0 & 1 \\ \end{bmatrix}}^{\underline{U}}$

כעת ניתן למצוא את המקדמים l_ij, u_ij כתלות ב-a_ij. אחרי שמצויות בידינו שתי המטריצות L,U ניתן לפתור בקלות יחסית מערכת משוואות מהצורה Ax=B:

$\ Ax=LUx \overbrace{=}^{Ux=y} Ly=B$

מאחר ו-L משושלת תחתונה, ניתן למצוא את y ללא קושי, ואז נשאר לפתור Ux=y, מה שגם מתאפשר בקלות מכיוון ש-U משולשת עליונה.

[עריכה] היפוך מטריצה

על מנת לפתור מערכת משוואות מהצורה Ax=B נוכל למצוא את המטריצה ההופכית של A ולבצע מכפלת מטריצות: x=A^-1B. בדרך כלל לא נשתמש בשיטה זו כי הביטוי האנליטי למטריצה הופכית מערב מציאת קופקטורים (זהו למעשה ה-Adjoint). נציג את הביטוי בכל זאת:

$A^{-1}={1 \over \begin{vmatrix}A\end{vmatrix}}\left(C_{ij}\right)^{T}={1 \over \begin{vmatrix}A\end{vmatrix}} \begin{pmatrix} C_{11} & C_{21} & \cdots & C_{j1} \\ C_{12} & \ddots & & C_{j2} \\ \vdots & & \ddots & \vdots \\ C_{1i} & \cdots & \cdots & C_{ji} \\ \end{pmatrix}$

כאשר Cij הם הקופקטורים של המטריצה A, אשר מתקבלים כזכור, על ידי חישוב מינורי המטריצה A.

שיטה אחרת:
בהינתן מערכת מסדר n, נפתור את n הבעיות $\ \underline{A}\cdot\underline{x}_i= \underline{e}_i$ , כאשר $\ \underline{e}_i$ הם וקטורי הבסיס הסטנדרטי. לשם פתרון n הבעיות הללו ניתן להעזר בשיטות המופיעות בדף זה.

[עריכה] פתרון מערכת תלת-אלכסונית

פתרון מערכת תלת-אלכסונית מתבצע באמצעות אלגוריתם הנקרא TDMA (Tridiagonal matrix algorithm), או אלגוריתם תומס.

$\left[ \begin{matrix} {b_1} & {c_1} & { } & { } & { 0 } \\ {a_2} & {b_2} & {c_2} & { } & { } \\ { } & {a_3} & {b_3} & \ddots & { } \\ { } & { } & \ddots & \ddots & {c_{N-1}}\\ { 0 } & { } & { } & {a_N} & {b_N}\\ \end{matrix} \right] \left\{ \begin{matrix} {x_1 } \\ {x_2 } \\ \cdot \\ \cdot \\ {x_N } \\ \end{matrix} \right\} = \left\{ \begin{matrix} {d_1 } \\ {d_2 } \\ \cdot \\ \cdot \\ {d_N } \\ \end{matrix} \right\} \quad\Rightarrow\quad \begin{array}{lcl} b_1x_1 + c_1x_2 & = & d_1 \\ a_2x_1 + b_2x_2 + c_2x_3 & = & d_2 \\ \vdots & { } & { } \\ a_ix_{i-1} + b_ix_i + c_ix_{i+1} & = & d_i \\ \vdots & { } & { } \\ a_Nx_{N-1} + b_Nx_N & = & d_N \\ \end{array}$

יש לנו אם כן, 3N-2 איברים, אשר חלקם יכולים להיות 0.

על מנת להגיע לפתרון, נבצע דירוג של שורה אחר שורה: נכפיל את השורה הראשונה ב- $\ -{a_2 \over b_1}$ ונוסיף אותה לשורה השנייה, כך שנקבל:

$\ \overbrace{\left( b_2 - {a_2c_1 \over b_1} \right)}^{\beta_2} x_2 +c_2x_3= \overbrace{d_2- {a_2d_2 \over b_1}}^{\delta_2} \quad\Rightarrow\quad \beta_2x_2+ c_2x_3= \delta_2$

כעת נכפיל את את המשוואה ב- $\ -{a_3 \over \beta_2}$ ונוסיף אותה לשורה השלישית, כך שנקבל:

$\ \overbrace{\left( b_3 - {a_3c_2 \over \beta_2} \right)}^{\beta_3} x_3 +c_3x_4= \overbrace{d_3- {a_3\delta_2 \over \beta_2}}^{\delta_3} \quad\Rightarrow\quad \beta_3x_3+ c_3x_4= \delta_3$

כעת, אחרי שעברנו לכתיב β,δ ניתן להכליל ולומר שהמשוואות ה-(i-1)-ית וה-i-ית הן מן הצורה:

$\ \beta_{i-1}x_{i-1}+ c_{i-1}x_i= \delta_{i-1}$

$\ a_ix_{i-1}+ b_ix_i+ c_ix_{i+1}= d_i$

בהתאמה, כך שכאשר נכפיל את המשוואה ה-(i-1)-ית ב- $\ -{a_i \over \beta_{i-1}}$ , ונוסיף למשוואה ה-i-ית, נקבל:

$\ \overbrace{\left( b_i - {a_ic_{i-1} \over \beta_{i-1}} \right)}^{\beta_i} x_i +c_ix_{i+1}= \overbrace{d_i- {a_i\delta_{i-1} \over \beta_{i-1}}}^{\delta_i} \quad\Rightarrow\quad \beta_ix_i+ c_ix_{i+1}= \delta_i$

לשם השלמת התמונה, נביט בשתי המשוואות האחרונות:

$\ \beta_{N-1}x_{N-1}+ c_{N-1}x_N= \delta_{N-1}$

$\ a_Nx_{N-1}+ b_Nx_N= d_N$

נכפיל את המשוואה ה-(N-1)-ית ב- $\ -{a_N \over \beta_{N-1}}$ , ונוסיף למשוואה ה-N-ית, כך שנקבל:

$\ \overbrace{\left( b_N - {a_Nc_{N-1} \over \beta_{N-1}} \right)}^{\beta_N} x_N= \overbrace{d_N- {a_N\delta_{N-1} \over \beta_{N-1}}}^{\delta_N} \quad\Rightarrow\quad \beta_Nx_N= \delta_N$

מכאן מחלצים את x_N ומקבלים את שאר הנעלמים על ידי הצבה לאחור.

כאשר נצטרך לכתוב תכנית מחשב, נרצה להכליל את כתיב β,δ גם על השורה הראשונה. לשם כך נוכל להגדיר:

$\ \delta_1=d_1,\ \beta_1=b_1$ , ואז האלגוריתם ימשיך:
$\ \beta_i= \left( b_i - {a_ic_{i-1} \over \beta_{i-1}} \right) \quad,\quad \delta_i= d_i- {a_i\delta_{i-1} \over \beta_{i-1}}$ עבור $\ i=2,3,...,N$ .
$\ x_N= {\delta_N \over \beta_N}$
$\ x_i= -{c_ix_{i+1}+ \delta_i \over \beta_i}$ עבור $\ i=N-1,N-2,...,2,1$ .

בהערכה פשוטה מתקבל כי מספר הפעולות המקסימלי לשיטה זו הינו 5N-4.

שימושים:

בהינתן מד"ר, אם נכתוב את משוואת הפרשים עבורה, נקבל מטריצה תלת-אלכסונית.

[עריכה] קישורים חיצוניים

ערך בוויקיפדיה: צורת ז'ורדן, פירוק LU (אנגלית), מטריצה הפיכה, אלגוריתם תומס

קוד מקור לאלגוריתם תומס.

[עריכה] שיטות איטרטיביות

מאחר והתרגלנו לסמן את מספר האיטרציה ב"n", נסמן את סדר המערכת ב-N.

[עריכה] שיטת Jacobi

שיטה זו משתמשת באיברי המטריצה A על מנת להתכנס לפתרון מתוך ניחוש התחלתי כלשהו. מבצעים איטרציות עד להתכנסות של כל וקטור הנעלמים. כלומר: לא מתבצעות איטרציות עבור כל אחד ואחד מהנעלמים בנפרד. מתוך הסכום : $\ \sum_{j=1}^N a_{ij}x_j=b_i,\ i=1,..,N$ נבודד את הנעלם x_i:

$\ a_{ii}x_i+ \sum_{j=1,j\neq i}^N a_{ij}x_j= b_i \quad\Rightarrow\ x_i= -\sum_{j=1,j\neq i}^N {a_{ij}\over a_{ii}} x_j+ {b_i\over a_{ii}}$

ואז השיטה האיטרטיבית היא:

$\ x_i^{(n+1)}= -\sum_{j=1,j\neq i}^N {a_{ij}\over a_{ii}} x_j^{(n)}+ {b_i\over a_{ii}}\ ,\quad i=1,..,N;\ n=0,1,2,...$

כאשר מספקים ניחוש התחלתי כלשהו:

$\ \underline{x}^{(0)}= \begin{Bmatrix} x_1^{(0)} \\ x_2^{(0)} \\ \vdots \\ x_N^{(0)} \end{Bmatrix}$

מה שמתרחש בפועל הוא 3 לולאות מקוננות אשר משתמשות בוקטור $\ \underline{x}^{(n)}$ על מנת לייצר את הוקטור $\ \underline{x}^{(n+1)}$ (ראו "קישורים חיצוניים" עבור האלגוריתם).

קריטריוני התכנסות
ניתן להשתמש באחד מן הקריטריונים הבאים:

$\ \left| x_i^{(n+1)}- x_i^{(n)} \right| < \epsilon\ ,\quad 1\le i\le N$
$\ \left| 1-\frac{x_i^{(n)}}{x_i^{(n+1)}} \right| < \epsilon\ ,\quad 1\le i\le N$
$\ \sqrt{\sum_{i=1}^n \left( x_i^{(n+1)}-x_i^{(n)} \right)^2} < \epsilon$

בדרך כלל התכנסות השיטה היא איטית ולכן יש לבצע מספר רב של איטרציות. את בדיקת ההתכנסות נהוג לבצע בתום הלולאה עבור כל נעלם בנפרד, ולא עבור כל איטרציה בנפרד.

התנאי להתכנסות
אם $\ \underline{\alpha}=(\alpha_1,\alpha_2...,\alpha_N)$ הוא הפתרון, אז $\ \sum_{j=1}^N a_{ij}\alpha_j \equiv b_i$ . נציב את וקטור השגיאה $\ \epsilon_i^{(n)}= x_i^{(n)}- \alpha_i$ לתוך שיטת Jacobi:

$\ x_i^{(n+1)}= -\sum_{j=1,j\neq i}^N {a_{ij}\over a_{ii}} x_j^{(n)}+ {b_i\over a_{ii}} \quad\Rightarrow\quad \alpha_i+\epsilon_i^{(n+1)}= -\sum_{j=1,j\neq i}^N {a_{ij}\over a_{ii}} \left( \alpha_j+\epsilon_j^{(n)} \right)+ {b_i\over a_{ii}}$

$\ \Rightarrow\quad \epsilon_i^{(n+1)}= -\sum_{j=1,j\neq i}^N {a_{ij}\over a_{ii}} \epsilon_j^{(n)}$

לשם נוחות, נגדיר את השגיאה המקסימלית באיטרציה: $\ \Epsilon^{(n)}= \max_{1\le j\le N, j\neq i} \left\{\left| \epsilon_j^{(n)} \right|\right\}$ ואז:

$\ \left| \epsilon_i^{(n+1)} \right| \le \sum_{j=1,j\neq i}^N \left|{a_{ij}\over a_{ii}}\right| \left|\epsilon_j^{(n)}\right| \le \sum_{j=1,j\neq i}^N \left|{a_{ij}\over a_{ii}}\right| \Epsilon^{(n)} \quad\Rightarrow\quad {\epsilon_i^{(n+1)}\over\Epsilon^{(n)}} \le \sum_{j=1,j\neq i}^N \left|{a_{ij}\over a_{ii}}\right|$

ואז התנאי להתכנסות הוא:

$\ {\epsilon_i^{(n+1)}\over\Epsilon^{(n)}} \le 1 \quad\Rightarrow\quad \sum_{j=1,j\neq i}^N \left|{a_{ij}\over a_{ii}}\right| \le 1 \quad\Rightarrow\quad \sum_{j=1,j\neq i}^N |a_{ij}| \le |a_{ii}|$

כלומר איברי האלכסון בכל שורה במטריצה A צריכים להיות גדולים מסכום כל שאר האיברים באותה השורה, ואז ההתכנסות מובטחת. ניתן להוכיח שנתאי זה מספיק אך לא הכרחי. לתנאי זה קוראים גם בשם "שליטה אלכסונית".

כעת כשאנו מודעים לתנאי ההתכנסות, ננסה לסדר את שורות המטריצה כך שהתנאי יתקיים, לפני הפעלת השיטה.

[עריכה] שיטת Gauss-Seidel

שיטת GS מפצלת את הסכימה לנעלמים לפני הנעלם הנוכחי ולנעלמים אחרי הנעלם הנוכחי:

$\ x_i^{(n+1)}= -\sum_{j=1}^{i-1} {a_{ij}\over a_{ii}} x_j^{(n+1)} -\sum_{j=i+1}^{N} {a_{ij}\over a_{ii}} x_j^{(n)}+ {b_i\over a_{ii}}\ ,\quad i=1,..,N;\ n=0,1,2,...$

בשיטה זו, בכל איטרציה משתמשים בערכים האחרונים שהתקבלו. כלומר כאן i-1 הנעלמים בוקטור הנעלמים מתעדכנים לפי הסדר יחד עם הנעלם ה-i, עם התקדמות הלולאה. מסיבה זו ההתכנסות מהירה פי 2 משיטת Jacobi. בדיקת ההתכנסות תתבצע כמו בשיטה הקודמת.

[עריכה] השוואה בין שיטת יעקובי לשיטת גאוס-זיידל

ננתח את השיטות במקרה של מערכת מסדר 2:

$\ \left[\begin{matrix} a & b \\ c & d \end{matrix}\right] \left\{\begin{matrix} x_1 \\ x_2 \end{matrix}\right\} = \left\{\begin{matrix} e_1 \\ e_2 \end{matrix}\right\}$

במקרה זה, התהליך האיטרטיבי הוא מהצורה:

$\ \begin{matrix} \underline{\mbox{Jacobi}} & \qquad \underline{\mbox{G-S}} \\ ax_1^{(n+1)}= e_1-bx_2^{(n)} & \qquad ax_1^{(n+1)}= e_1-bx_2^{(n)} \\ dx_2^{(n+1)}= e_2-cx_1^{(n)} & \qquad dx_2^{(n+1)}= e_2-cx_1^{(n+1)} \end{matrix}$

כלומר ההבדל היחיד הוא ששיטת גאוס-זיידל משתמשת בערך המעודכן של x₁.

נציב את הביטוי לשגיאה $\ x_i=\epsilon_i+\alpha_i$ ונקבל:

$\ \begin{matrix} \underline{\mbox{Jacobi}} & \qquad \underline{\mbox{G-S}} \\ a \left( \epsilon_1^{(n+1)}+\alpha_1 \right)= e_1-b \left(\epsilon_2^{(n)}+\alpha_2 \right) & \qquad a \left(\epsilon_1^{(n+1)}+\alpha_1 \right)= e_1-b \left(\epsilon_2^{(n)}+\alpha_2 \right) \\ d \left(\epsilon_2^{(n+1)}+\alpha_2 \right)= e_2-c \left(\epsilon_1^{(n)}+\alpha_1 \right) & \qquad d \left(\epsilon_2^{(n+1)}+\alpha_2 \right)= e_2-c \left(\epsilon_1^{(n+1)}+\alpha_1 \right) \end{matrix}$

נזכור כי הוקטור $\ (\alpha_1,\alpha_2)^T$ פותר את המערכת, ואז באמצעות המשוואות הנ"ל נוכל למצוא קשר בין שתי שגיאות עוקבות:

$\ \begin{matrix} \underline{\mbox{Jacobi}} & \qquad \underline{\mbox{G-S}} \\ \epsilon_1^{(n+1)}= -{b\over a}\epsilon_2^{(n)} & \qquad \epsilon_1^{(n+1)}= -{b\over a}\epsilon_2^{(n)} \\ \epsilon_2^{(n+1)}= -{c\over d}\epsilon_1^{(n)} & \qquad \epsilon_2^{(n+1)}= -{c\over d}\epsilon_1^{(n+1)} \\ \end{matrix}$

כך שעבור שיטת יעקובי מתקיים:

$\ \left\{\begin{matrix} \epsilon_1^{(n+1)} \\ \epsilon_2^{(n+1)} \end{matrix}\right\} = \overbrace{{b\over a}{c\over d}}^{\sigma} \left\{\begin{matrix} \epsilon_1^{(n-1)} \\ \epsilon_2^{(n-1)} \end{matrix}\right\} \qquad\Rightarrow \underline{\epsilon}^{(2n)}= \sigma^n \underline{\epsilon}^{(0)}$

ואילו עבור שיטת גאוס-זיידל מתקיים:

$\ \left\{\begin{matrix} \epsilon_1^{(n+1)} \\ \epsilon_2^{(n+1)} \end{matrix}\right\} = \overbrace{{b\over a}{c\over d}}^{\sigma} \left\{\begin{matrix} \epsilon_1^{(n)} \\ \epsilon_2^{(n)} \end{matrix}\right\} \qquad\Rightarrow \underline{\epsilon}^{(n)}= \sigma^n \underline{\epsilon}^{(0)}$

כלומר בעוד שבשיטת יעקובי σ היא השגיאה כעבור כל שתי איטרציות, בשיטת גאוס-זיידל σ היא השגיאה בכל איטרציה בודדת. לכן בשיטת גאוס-זיידל ההתכנסות מהירה פי 2.

[עריכה] שיטת Successive Over-Relaxation

שיטת SOR נועדה על מנת לזרז את התכנסותה של שיטת Gauss-Seidel.

$\ x_i^{(n+1)}= x_i^{(n)}+ \omega \left[ x_{i_{GS}}^{(n+1)}-x_i^{(n)} \right]$

ω הוא פרמטר הרלקסציה אשר מקיים $\ 0<\omega<2$ . על מנת להאיץ תהליך התכנסות איטי לוקחים ω>1 ואילו על מנת להבטיח התכנסות עבור שיטות מתבדרות, לוקחים ω<1. שימו לב כי עבור ω=1 מקבלים חזרה את שיטת GS.

כמו כן, עבור מטריצה A נתונה, קיים ערך אופטימלי של ω אשר מביא להתכנסות המהירה ביותר.

ניתן להוכיח שאם שיטת GS מתכנסת, אז שיטת SOR תתכנס מהר יותר.

דרך אחרת לפיתוח השיטה הוא באמצעות שימוש ייצוג המטריצה A באמצעות המכפלה DLU, כאשר D היא מטריצה אלכסונית (למידע נוסף ראו "קישורים חיצוניים").

[עריכה] קישורים חיצוניים

ערך בוויקיפדיה: שיטת Jacobi (אנגלית), שיטת Gauss-Seidel (אנגלית), שיטת SOR (אנגלית)

הסברים באתר MathWorld: שיטת Jacobi, שיטת Gauss-Seidel, שיטת SOR, אלגוריתם תומס
הסברים באתר אוניברסיטת USCF על שיטת SOR עם קטעי קוד עבור תוכנת Mathematica.

הפרק הקודם:
יציבות נומרית

פתרון מערכת משוואות לינאריות

הפרק הבא:
פתרון משוואות דיפרנציאליות

אנליזה נומרית/פתרון מערכת משוואות לינאריות

מתוך ויקיספר, אוסף ספרי הלימוד והמדריכים החופשי.

תוכן עניינים

[עריכה] שיטות אנליטיות

[עריכה] שיטת גאוס

[עריכה] שיטת גאוס-ג'ורדן

[עריכה] פירוק LU

[עריכה] היפוך מטריצה

[עריכה] פתרון מערכת תלת-אלכסונית

[עריכה] קישורים חיצוניים

[עריכה] שיטות איטרטיביות

[עריכה] שיטת Jacobi

[עריכה] שיטת Gauss-Seidel

[עריכה] השוואה בין שיטת יעקובי לשיטת גאוס-זיידל

[עריכה] שיטת Successive Over-Relaxation

[עריכה] קישורים חיצוניים

צפיות

כלים אישיים

ניווט

קהילה

חיפוש

תיבת כלים

הדפסה/ייצוא