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

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

מאחר ואנו מתעניינים בבעיות בעלות פתרון (יחיד), נתקל במטריצות ריבועיות בלבד, אשר מקיימות $|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_{1}x_{1}+c_{1}x_{2}&=&d_{1}\\a_{2}x_{1}+b_{2}x_{2}+c_{2}x_{3}&=&d_{2}\\\vdots &{}&{}\\a_{i}x_{i-1}+b_{i}x_{i}+c_{i}x_{i+1}&=&d_{i}\\\vdots &{}&{}\\a_{N}x_{N-1}+b_{N}x_{N}&=&d_{N}\\\end{array}}$

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

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

$\ \overbrace {\left(b_{2}-{a_{2}c_{1} \over b_{1}}\right)} ^{\beta _{2}}x_{2}+c_{2}x_{3}=\overbrace {d_{2}-{a_{2}d_{2} \over b_{1}}} ^{\delta _{2}}\quad \Rightarrow \quad \beta _{2}x_{2}+c_{2}x_{3}=\delta _{2}$

כעת נכפיל את את המשוואה ב- $\ -{a_{3} \over \beta _{2}}$ ונוסיף אותה לשורה השלישית, כך שנקבל:

$\ \overbrace {\left(b_{3}-{a_{3}c_{2} \over \beta _{2}}\right)} ^{\beta _{3}}x_{3}+c_{3}x_{4}=\overbrace {d_{3}-{a_{3}\delta _{2} \over \beta _{2}}} ^{\delta _{3}}\quad \Rightarrow \quad \beta _{3}x_{3}+c_{3}x_{4}=\delta _{3}$

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

$\ \beta _{i-1}x_{i-1}+c_{i-1}x_{i}=\delta _{i-1}$

$\ a_{i}x_{i-1}+b_{i}x_{i}+c_{i}x_{i+1}=d_{i}$

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

$\ \overbrace {\left(b_{i}-{a_{i}c_{i-1} \over \beta _{i-1}}\right)} ^{\beta _{i}}x_{i}+c_{i}x_{i+1}=\overbrace {d_{i}-{a_{i}\delta _{i-1} \over \beta _{i-1}}} ^{\delta _{i}}\quad \Rightarrow \quad \beta _{i}x_{i}+c_{i}x_{i+1}=\delta _{i}$

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

$\ \beta _{N-1}x_{N-1}+c_{N-1}x_{N}=\delta _{N-1}$

$\ a_{N}x_{N-1}+b_{N}x_{N}=d_{N}$

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

$\ \overbrace {\left(b_{N}-{a_{N}c_{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 _{N}x_{N}=\delta _{N}$

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

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

1. $\ \delta _{1}=d_{1},\ \beta _{1}=b_{1}$, ואז האלגוריתם ימשיך:
2. $\ \beta _{i}=\left(b_{i}-{a_{i}c_{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$.
3. $\ x_{N}={\delta _{N} \over \beta _{N}}$
4. $\ x_{i}={\delta _{i}-c_{i}x_{i+1} \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,...$

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

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

• $\ \left|x_{i}^{(n+1)}-x_{i}^{(n)}\right|<\epsilon \ ,\quad 1\leq i\leq N$
• $\ \left|1-{\frac {x_{i}^{(n)}}{x_{i}^{(n+1)}}}\right|<\epsilon \ ,\quad 1\leq i\leq 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)}$

לשם נוחות, נגדיר את השגיאה המקסימלית באיטרציה: $\ \mathrm {E} ^{(n)}=\max _{1\leq j\leq N,j\neq i}\left\{\left|\epsilon _{j}^{(n)}\right|\right\}$ ואז:

$\ \left|\epsilon _{i}^{(n+1)}\right|\leq \sum _{j=1,j\neq i}^{N}\left|{a_{ij} \over a_{ii}}\right|\left|\epsilon _{j}^{(n)}\right|\leq \sum _{j=1,j\neq i}^{N}\left|{a_{ij} \over a_{ii}}\right|\mathrm {E} ^{(n)}\quad \Rightarrow \quad {\epsilon _{i}^{(n+1)} \over \mathrm {E} ^{(n)}}\leq \sum _{j=1,j\neq i}^{N}\left|{a_{ij} \over a_{ii}}\right|$

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

$\ {\epsilon _{i}^{(n+1)} \over \mathrm {E} ^{(n)}}\leq 1\quad \Rightarrow \quad \sum _{j=1,j\neq i}^{N}\left|{a_{ij} \over a_{ii}}\right|\leq 1\quad \Rightarrow \quad \sum _{j=1,j\neq i}^{N}|a_{ij}|\leq |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.

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

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

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