אנליזה נומרית/פתרון משוואות דיפרנציאליות

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

נעבוד עם המקטעים $\ x_i=x_0+ih,\ i\in \mathbb{N}$ . נתרכז במציאת פתרון לבעיה $\ y'(x)=f(x,y)$ , עם תנאי ההתחלה $\ y(x_0)=y_0$ .

מכאן והלאה, הסימול $\ y(x_i)$ יציין הערך האנליטי המדוייק של הפונקציה בנקודה x_i, ואילו $\ y_i$ יציין את הערך המקורב שחושב באמצעות השיטה. שימו לב כי במקרים מיוחדים יתכן שהערכים זהים, כלומר: לפעמים $\ y(x_i)=y_i$ .

[עריכה] פתרון באמצעות טור טיילור

השיטה מתבססת על פיתוח לטור טיילור של ערכי הפונקציה $\ y_1,y_2,...$ סביב הנקודות $\ x_1,x_2,...$ , נקודה אחר נקודה, כך שבכל צעד משתמשים בתוצאה הקודמת. שיטה כזו נקראת שיטה "חד-צעדית" (בהמשך נעסוק גם בשיטות רב-צעדיות).

לשם כך נפתח תחילה את הפונקציה לטור טיילור סביב x₀:

$\ y_1=y_0+hy_0'+ \tfrac{h^2}{2}y_0''+ ...$

השלב הבא הוא לחשב את y₂ באמצעות הערך שהתקבל עבור y₁, ובאופן כללי:

$\ y_{i+1}=y_i+hy_i'+ \tfrac{h^2}{2}y_i''+...$

שימו לב כי רק y₀ הוא ערך מדוייק, כי ערכי y_i עבור i>0 התקבלו באמצעות חישוב מקורב. המטרה היא, אם כן, להקטין ככל הניתן את ההפרש $\ y(x_i)-y_i\ ,\ i\ge1$ .

על מנת למצוא קירוב נומרי, נשתמש, לשם הדגמה, בשלושת האיברים הראשונים של הטור. נמצא את ערכי המקדמים עבור y₁:

$\ y_0=y(x_0)$ (נתון)

$\ y_0'=y'(x_0)=f(x_0,y_0)$

$\ y_0''= \left.\left( \tfrac{d}{dx}y' \right)\right|_{x_0}= \left.\left[ \tfrac{d}{dx}f(x,y(x)) \right]\right|_{x_0}= \left.\left[ \tfrac{\partial f}{\partial x}+ \tfrac{\partial f}{\partial y}\tfrac{dy}{dx} \right]\right|_{x_0}$ (כלל השרשרת)

שימו לב כי ערכו של הביטוי שהתקבל באמצעות כלל השרשרת הינו ידוע ב-x₀. בדרך זו ממשיכים ומוצאים את ערכי y₂, y₃...

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

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

פתרון מד"ר (אוילר):
$\ y_{i+1}= y_i+hf(x_i,y_i)$

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

$\ y_{i+1}=y_i+hy_i'+ R_2= y_i+hf(x_i,y_i)+ O(h^2)$

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

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

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

ערך בוויקיפדיה: שיטת אוילר (אנגלית)

-

[עריכה] פתרון בשיטות רונגה-קוטה

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

כפי שראינו, לפי טור טיילור:

$\ y_{i+1}= y_i+hf(x_i,y_i)+ \frac{h^2}{2}\left[ \frac{\partial f(x_i,y_i)}{\partial x}+ \frac{\partial f(x_i,y_i)}{\partial y}f(x_i,y_i) \right]+O(h^3)$

נשתמש בשיטה הבאה, אשר מזכירה מעט את משפט לגראנז':

$\ y_{i+1}= y_i+\lambda_1hf(x_i,y_i)+ \lambda_2hf\left[x_i+\mu_1h,y_i+\mu_2hf(x_i,y_i)\right]$

ואת הקבועים $\ \lambda_1,\lambda_2,\mu_1,\mu_2$ נמצא על ידי השוואת מקדמים. לשם כך נפתח את השיטה לטור טיילור סביב $\ (x_i,y_i)$ :

$\ y_{i+1}= y_i+\lambda_1hf(x_i,y_i)+ \lambda_2h \left[ f(x_i,y_i)+ \mu_1h \frac{\partial f(x_i,y_i)}{\partial x}+ \mu_2hf(x_i,y_i) \frac{\partial f(x_i,y_i)}{\partial y} +O(h^2) \right]$

$\ \Rightarrow\ \begin{array}{lcl} hf(x_i,y_i) & : & \lambda_1+\lambda_2=1 \\ h^2 \frac{\partial f(x_i,y_i)}{\partial x} & : & \lambda_2\mu_1=\tfrac{1}{2} \\ h^2 f(x_i,y_i) \frac{\partial f(x_i,y_i)}{\partial y} & : & \lambda_2\mu_2=\tfrac{1}{2} \end{array}$

מאחר ויש לנו 4 נעלמים ורק 3 משוואות, נקבע אחד מהם שרירותית וזאת תהיה המשוואה הרביעית. עבור $\ \lambda_1=\tfrac{1}{2}$ נקבל את שיטת Heun, ועבור $\ \lambda_1=0$ נקבל את שיטת improved polygon.

[עריכה] הצורה הכללית של שיטה מסדר 2

אם נסמן $\begin{cases} \lambda_1=1-\kappa \\ \lambda_2=\kappa\neq 0 \\ \mu_1=\mu_2={1\over 2\kappa} \end{cases}$ נקבל את הצורה הכללית של שיטת רונגה קוטה מסדר שני:

$\ y_{i+1}=y_i+ h[(1-\kappa)f(x_i,y_i)+ \kappa f(x_i+ \tfrac{h}{2\kappa}, y_i+\tfrac{h}{2\kappa} f(x_i,y_i) +O(h^3)$

ואז שיטת Heun תתקבל עבור $\ \kappa= \tfrac{1}{2}$ ושיטת improved polygon עבור $\ \kappa=1$ .

[עריכה] שיטת Heun (שיטה מסדר 2)

שיטה זו נקראת גם "שיטת אוילר המשופרת" (improved Euler method). מבחינה גאומטרית, היא משתמשת בממוצע השיפועים בשתי נקודות סמוכות $\ (x_i,y_i),\ (x_i+h,y_i+hy_i')$ המתקבלות משיטת אוילר. שימו לב כי: $\ x_{i+1}=x_i+h,\ y_{i+1}= y_i+hy_i'$ רק בשיטת אוילר, ואילו כאן אלו שתי נקודות שונות. כאן, הנקודה $\ (x_{i+1},y_{i+1})$ מתקבלת מחיתוך הקו הישר היוצא מנקודה $\ (x_i,y_i)$ וששיפועו הוא הממוצע הנ"ל, עם הישר $\ x=x_{i+1}=x+h$ .

שיטת Heun:

$\ y_{i+1}= y_i+ \tfrac{h}{2}(C_i,D_i)$
$\ C_i= f(x_i,y_i)$
$\ D_i=f(x_i+h,y_i+C_i)$

נקח $\ \lambda_1=\tfrac{1}{2}$ ואז $\ \lambda_2=\tfrac{1}{2}\ ,\ \mu_1=\mu_2=1$ , כך שהשיטה המתקבלת היא:

$\ y_{i+1}= y_i+ \tfrac{1}{2}hf(x_i,y_i)+ \tfrac{1}{2}hf[x_i+h,y_i+hf(x_i,y_i)]+ O(h^3)$

[עריכה] שיטת improved polygon (שיטה מסדר 2)

שיטה זו נקראת גם modified Euler method. בניגוד לשיטת Heun, במקום לקחת את ממוצע השיפועים של שתי נקודות, ניקח את השיפוע של נקודת האמצע.

במקרה זה נקח $\ \lambda_1=0$ ואז $\ \lambda_2=1\ ,\ \mu_1=\mu_2= \tfrac{1}{2}$ , כך שהשיטה המתקבלת היא:

$\ y_{i+1}= y_i+ hf[x_i+\tfrac{h}{2},y_i+\tfrac{h}{2}f(x_i,y_i)]+ O(h^3)$

[עריכה] שיקולי דיוק

שני שיקולים ינחו אותנו בעת בחירת מרווח האינטגרציה h:

בהתאם לגודל האופייני של הבעיה הפיזיקלית נחליט בכמה סדרי גודל על h להיות קטן יותר.
ערכו של h יכול להשתנות מנקודה לנקודה בהתאם להתנהגות הפונקציה. אם אין גרדיאנטים חריפים מנקודה לנקדוה, ניתן לקחת מרווחים גדולים יותר בלי לפגוע משמעותית בדיוק, כך ש- $\ h={c \over y'}= {c \over f(x_i,y_i)}$ , כאשר c הוא קבוע פרופורציה כלשהו.

[עריכה] אבחנות

בעוד ששיטת אוילר מחשבת את הערך הבא באמצעות השיפוע בערך הנוכחי, שיטות רונגה-קוטה משתמשות במידע נוסף.
בשיטות רונגה קוטה מסדר שני אנו מגיעים לדיוק מסדר $\ O(h^2)$ באמצעות חישוב ערך הפונקציה f בשתי נקודות בלבד, לעומת שיטת טור-טיילור בה נדרשים שלושה חישובים: $\ f(x_i), \tfrac{\partial f(x_i)}{\partial x}, \tfrac{\partial f(y_i)}{\partial y}$

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

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

הסברים על שיטת Heun באתר אוניברסיטת קוריאה.

[עריכה] שיטות רב-צעדיות

שיטות אלה, בניגוד לנ"ל, משתמשות במספר נקודות קודמות. נפתח את בסיס השיטה:

$\ y'(x)=f(x,y)\quad\Rightarrow\ \int\limits_{x_i}^{x_{i+1}} \frac{dy}{dx}dx= \int\limits_{x_i}^{x_{i+1}} fdx \quad\Rightarrow\ y_{i+1}-y_i= \int\limits_{x_i}^{x_{i+1}} fdx$

נשאר, אם כן, לפתור את האינטגרל $\ \int\limits_{x_i}^{x_{i+1}} fdx$ . לשם כך נבצע החלפת משתנים: $\ x=x_i+ \theta h,\ dx=hd\theta,\ 0\le\theta\le 1$ כך שערכי f המתקבלים בין x_i,x_i+1 ניתנים לתיאור על ידי: $\ f=E^{\theta} f(x_i)$ ולכן האינטגרל מקבל את הצורה:

$\ \int\limits_0^1 E^{\theta} \overbrace{f(x_i,y_i)h}^{const.}d\theta= \int\limits_0^1 (1-\nabla)^{-\theta}f(x_i,y_i)hd\theta=$

$\ =\ \int\limits_0^1 h\left[ 1+\theta\nabla+ \tfrac{\theta(\theta+1)}{2}\nabla^2+... \right]f(x_i,y_i) d\theta=$

$\ = h\left.\left[ \theta+{\theta^2\over 2}\nabla+ \left( {\theta^3\over 6}+{\theta^2\over 4} \right)\nabla^2... \right]f(x_i,y_i)\right|_0^1=$

$\ =\ h\left[ 1+{\nabla\over 2}+ {5 \over 12}\nabla^2+ {9\over 24}\nabla^3+... \right]f(x_i,y_i)$

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

שיטה רב-צעדית:
$\ y_{i+1}= y_i+ h\left[ 1+{\nabla\over 2}+ {5 \over 12}\nabla^2+ {3\over 8}\nabla^3+... \right]f(x_i,y_i)$

שימו לב כי אם ניקח איבר בודד, נקבל את שיטת אוילר.

הערה: אחד המאפיינים של שיטות רב-צעדיות הוא שאינן מאתחלות את עצמן (not self-starting), כלומר יש לספק עבורן מספר תנאי התחלה באמצעות שיטות חד-צעדיות בעלות אותו מרווח h.

[עריכה] שיטת Adams-Bashforth

נפתח את 4 האיברים הראשונים שציינו לעיל:

$\ \begin{align} y_{i+1} & = y_i+{h\over 24} \left[ 55f(x_i,y_i)-59f(x_{i-1},y_{i-1})+ \right. \\ & + \left. 37f(x_{i-2},y_{i-2})-9f(x_{i-3},y_{i-3}) \right]+ O(h^5) \\ \end{align} \ ,\quad i\ge 3$

שימו לב כי לשיטה זו יש לספק 3 תנאי התחלה, אשר בדרך כלל מוצאים אותם בשיטת רונגה-קוטה או בשיטת אוילר. כמו כן, על אותה שיטה להשתמש באותו מרווח h.
נשווה שיטה זו לשיטת רונגה-קוטה מסדר 4, מבחינת מספר פעולות החישוב: $\ \tfrac{R-K}{A-B}= \tfrac{4N}{4\cdot 3+(N-3)}=\tfrac{4N}{N+9}\approx 4$ , כלומר שיטת רונגה-קוטה איטית פי 4.
על מנת להשתמש בשיטה זו עבור מערכת מד"ר, נחליף את הכתיב הסקלרי בכתיב הוקטורי:

$\ \begin{align} \vec y_{i+1} & = \vec y_i+{h\over 24} \left[ 55\vec f(x_i,\vec y_i)-59\vec f(x_{i-1},\vec y_{i-1})+ \right. \\ & + \left. 37\vec f(x_{i-2},\vec y_{i-2})-9\vec f(x_{i-3},\vec y_{i-3}) \right]+ O(h^5)\ \\ \end{align} ,\quad i\ge 3$

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

שיטת אדאמס היא מהצורה:

$\ y_{n+1}=y_n+ {h \over 12} [5y_{n+1}'+8y_n'-y_{n-1}']$

נמצא את הביטוי לשגיאה באמצעות פיתוח האגפים לטור טיילור והשוואת מקדמים:

$\ \mbox{LHS:}\quad y_n+ hy_n'+ \tfrac{h^2}{2}y_n''+ \tfrac{h^3}{3!}y_n'''+ \tfrac{h^4}{4!}y_n^{(4)}+...$

$\ \begin{align} \mbox{RHS:} & \quad y_n+ {h\over 12}\left[\left(5y_n'+ 5hy_n''+ 5\tfrac{h^2}{2}y_n'''+ 5\tfrac{h^3}{3!}y_n^{(4)}+...\right)+ \right. \\ & \left. + (8y_n')- \left(y_n'-hy_n''+ \tfrac{h^2}{2}y_n'''- \tfrac{h^3}{3!}y_n^{(4)}+...\right)\right] +R \\ \end{align}$

$\ \begin{align} y_n: & \quad 1=1 \quad\rightarrow OK\\ y_n': & \quad h=h\left( \tfrac{5}{12}+\tfrac{8}{12}-\tfrac{1}{12} \right)=h \quad\rightarrow OK\\ y_n'': & \quad \tfrac{h^2}{2}=\tfrac{h^2}{12}(5+1)=\tfrac{h^2}{2} \quad\rightarrow OK\\ y_n''': & \quad \tfrac{h^3}{3!}=\tfrac{h^3}{12}\left(\tfrac{5-1}{2}\right)=\tfrac{h^3}{3!} \quad\rightarrow OK\\ y_n^{(4)}: & \quad \tfrac{h^4}{4!} \neq \tfrac{h^4}{12}\left(\tfrac{5+1}{3!}\right)=\tfrac{h^4}{12} \quad\Rightarrow\quad R=y^{(4)}(c)\left( \tfrac{h^4}{24}-\tfrac{h^4}{12} \right)= -\tfrac{h^4}{24}y^{(4)}(c)= O(h^4)\\ \end{align}$

[עריכה] שיטת Cowell-Numerov

שיטה זו היא עבור מד"ר מסדר 2, כאשר f אינה פונקציה של $\ y'$ אלא של $\ y''$ , ותנאי ההתחלה נתונים ב-x=0.

$\ \begin{cases} y_{i+1}=2y_i-y_{i-1}+ {h^2\over 12}[f(x_{i+1},y_{i+1})+10f(x_i,y_i)+f(x_{i-1},y_{i-1})] \\ y''(x)=f(x,y) \end{cases}$

שימו לב כי זוהי שיטה סתומה מאחר ו- $\ y_{i+1}$ תלוי ב-f, אשר תלוי ב- $\ y''$ .

[עריכה] שיטות נוספות

$\ y_{i+1}= y_{i-1}+2hf(x_i,y_i)+ O(h^3)$
$\ y_{i+1}= y_{i-3}+ \tfrac{4h}{3}[2f(x_i,y_i)-f(x_{i-1},y_{i-1})+ 2f(x_{i-2},y_{i-2})] +O(h^5)$

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

ערך בוויקיפדיה: שיטות רב צעדיות (אנגלית)

שיטת אדאמס-באשפורת' באתר MathWorld.
שיטת אדאמס-באשפורת' באתר אוניברסיטת CSUF (כולל אנימציות, הוכחות ודוגמאות שימוש).
אנימציית אדאמס-באשפורת' באתר אוניברסיטת אילנוי.

[עריכה] שיטות מסוג מעריך-מתקן

שיטות מעריך-מתקן (Predictor-Corrector Methods) מתבססות על פתרון האינטגרל על $\ y'(x)= f(x,y)dx$ , או: $\ y_{i+1}= y_i+\int\limits_{x_i}^{x_{i+1}} fdx$ . השלב הבא הוא להשתמש בנוסחה מקורבת לאינטגרל על מנת להתקדם לנקודה הבאה (זכרו כי $\ h=x_{i+1}-x_i$ ).

אם נשתמש בשיטת המלבן להערכת האינטגרל, נקבל:

$\ y_{i+1}= y_i+ \int\limits_{x_i}^{x_{i+1}} fdx= y_i+ hf\left(\tfrac{x_{i+1}+x_i}{2},y_i\right)$

שימו לב כי כל הפרמטרים ידועים, פרט ל-y_i+1, שאותו אנו מחפשים.

אם, לעומת זאת, נשתמש בשיטת הטרפז, נקבל:

$\ y_{i+1}= y_i+ \int\limits_{x_i}^{x_{i+1}} fdx= y_i+ \tfrac{h}{2}[f(x_{i+1},y_{i+1})+f(x_i,y_i)]+ O(h^3)$

שימו לב כי קיבלנו קשר סתום עבור y_i+1, וכביכול אין באפשרותנו להמשיך.

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

כללית, כל שיטת מעריך-מתקן מורכבת משני שלבים:

מעריך: בוחרים בשיטה כלשהי אשר מאפשרת חישוב ישיר של y_i+1.
מתקן איטרטיבי: מציבים את המעריך לתוך השיטה הסתומה.

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

[עריכה] שיקולי התכנסות

ניתן לראות בביטוי הסתום עבור y_i+1 מעין שיטה איטרטיבית (נסמן $\ \hat y=y_{i+1}$ ): $\ \hat y_{n+1}=\Phi(\hat y_n)$ (כל שאר הפרמטרים ידועים). וכזכור, התנאי להתכנסות (מסדר 1) הוא $\ |\Phi'(\hat y)|<1$ . לדוגמה: עבור שיטת הטרפז:

$\ |\Phi'(\hat y)|= \left|{h\over 2}{\partial f\over\partial y}\right|<1\ \Rightarrow\ h < \frac{2}{\tfrac{\partial f}{\partial\hat y}}$

ואז ההתכנסות מובטחת. על מנת להבטיח התכנסות מהירה, יש לדאוג לכך ש- $\ h \ll \frac{2}{\tfrac{\partial f}{\partial\hat y}}$ .

[עריכה] מקרים פרטיים

Adams-Bashforth:

שיטה זו היא למעשה קירוב של פולינום לגראנז' המתקבל על ידי ארבע נקודות. ניתן לאמר שהמעריך הוא:

$\ y_{i+1}^{predictor}= y_i+{h\over 24} \left[ 55f(x_i,y_i)-59f(x_{i-1},y_{i-1})+37f(x_{i-2},y_{i-2})-9f(x_{i-3},y_{i-3}) \right]$

והמתקן יתקבל באמצעות הנקודה החדשה ושלוש הנקודות הקודמות לה:

$\ y_{i+1}^{corrector}= y_i+{h\over 24} \left[ 9f(x_{i+1},y_{i+1})+19f(x_i,y_i)-5f(x_{i-1},y_{i-1})+(x_{i-2},y_{i-2}) \right]$

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

אנימציה של מעריך-מתקן באתר אוניברסיטת אילנוי.
הסברים באתר Dr. Math.
הסברים באתר הפקולטה לפיזיקה של הטכניון.

[עריכה] יציבות של שיטות לפתרון מד"ר

יציבות תלויה במד"ר אותה רוצים לפתור. מניחים כי הפתרון מתנהג כמו משוואת ההפרשים: $\ y_n= \sigma y_{n-1}=...= \sigma^n y_0$ . כאשר נציב ביטוי זה ( $\ y_n=c\cdot\sigma^n$ ) לשיטה, נקבל משוואה עבור σ אשר תלויה ב- $\ f(x,y)$ הנידון. בדרך כלל נקבל שני פתרונות עבור σ, כאשר לאחד נקרא אנליטי ( $\ \sigma_{analytic}$ ) - והוא מבטא את הפתרון למשוואה - ולשני פרזיטי ( $\ \sigma_{parasitic}$ ) - אשר מבטא את השגיאה הנרכשת בכל איטרציה ושאינה קשורה לפתרון האמיתי. התנאים ליציבות הם: