שיטות גרפיות ונומריות
שיטות אוילר
The Basic Euler Methods
צעד קדימה, צעד אחורה ו־trapezoid — דיוק מול יציבות.
- 1נקודת פתיחהמה כבר ידוע ומה השאלה?
- 2ניבוימה יקרה במקרה הפשוט?
- 3המחשהמה באמת משתנה בתמונה?
- 4ניסוח והצדקהשם, נוסחה ולמה היא נכונה
- 5בדיקהדוגמה, תרגול ומשוב
- 6העברהמקרה חדש והמשך המסלול
שאלה לפני נוסחה
מה אנחנו מנסים לראות?
להשוות בין אוילר קדמי, אוילר אחורי ושיטת הטרפז, ולבחון בכל אחת דיוק, עלות ויציבות לפני שמפעילים אותה על מערכת דינמית.
כל צעד מחליף קטע קצר של עקומה בקו ישר. אוילר קדמי משתמש בשיפוע שכבר ידוע, אוילר אחורי דורש את השיפוע בקצה החדש, ושיטת הטרפז מאזנת בין שניהם.
לפני ההמחשה: קראו רק את הכותרת “אותו צעד, שלוש בחירות לשיפוע”. מה לדעתכם ישתנה, מה יישאר קבוע, ואיזה סימן בתמונה יוכל להפריך את הניבוי?
לראות לפני שמחשבים · המקרה הפשוט ביותר
אותו צעד, שלוש בחירות לשיפוע
ההמחשה משווה על y′=−10y את אוילר הקדמי, אוילר האחורי והטרפז לצד העקומה המדויקת ומציגה את גורם ההגברה של כל אחת.
מה כדאי לגלות: הפרידו בין שתי שאלות: האם הנקודות נשארות חסומות, והאם הן קרובות לעקומה המדויקת? אלה אינן אותה שאלה.
אחרי שראינו
בסוף השיעור תוכלו
- לגזור את נוסחאות אוילר הקדמי, אוילר האחורי והטרפז מן השינוי המצטבר בצעד זמן אחד.
- לבצע ידנית כמה צעדים ולפתור את המשוואה הסמויה של אוילר אחורי במקרה לינארי.
- להבחין בין סדר דיוק לבין יציבות ולחשב את גורם ההגברה של כל שיטה.
- לבחור שיטה וגודל צעד מתאימים לבעיה הכוללת קצב דעיכה מהיר.
נותנים שמות למה שראינו
מושגים שנשתמש בהם
- רשת זמן
- הנקודות tₙ=t₀+nh שבהן נשמרים הקירובים yₙ≈y(tₙ).
- גודל צעד
- המרחק h בין שתי נקודות זמן עוקבות; צעד קטן לרוב משפר דיוק אך מגדיל את העלות.
- שיטה מפורשת
- עדכון שמחשב את yₙ₊₁ ישירות מנתונים שכבר ידועים, כמו אוילר קדמי.
- שיטה סמויה
- עדכון שבו yₙ₊₁ מופיע גם בתוך אגף ימין ולכן יש לפתור משוואה בכל צעד.
- שגיאה מקומית
- השגיאה שנוצרת בצעד יחיד בהנחה שנכנסנו אליו עם הערך המדויק.
- שגיאה גלובלית
- הפער המצטבר לאחר צעדים רבים עד זמן קבוע.
- סימון O(hᵖ)
- טענה שקיים קבוע C שאינו תלוי ב־h כך שהשגיאה קטנה לכל היותר מ־C·hᵖ כאשר h קטן מספיק.
- גורם הגברה
- המספר R(z) שמכפיל את yₙ בצעד אחד של משוואת הבדיקה y′=λy, כאשר z=hλ.
- A-יציבות
- תכונה שלפיה כל z בעל Re(z)≤0 נמצא בתחום היציבות |R(z)|≤1.
- בעיה קשיחה
- בעיה שבה קצבים מהירים ואיטיים יחד כופים על שיטה מפורשת צעד קטן בגלל יציבות, גם כשהפתרון הרצוי משתנה לאט.
מהמילים אל הסימנים
כותבים את הפעולה שכבר ראינו
פתחו את הנוסחה רק אחרי שאתם יכולים לומר במילים מה כל פעולה אמורה לעשות. בשלב הבא נצדיק אותה צעד אחר צעד.
להצגת נוסחת המפתח
קראו כל סימן כפעולה, לא כמשהו שצריך לשנןyₙ₊₁ = yₙ + h f(tₙ,yₙ)מרכיבים את הפעולות
מהתמונה אל הנוסחה — הגזירה, בלי לדלג על הצעדים
הצעד המדויק הוא אינטגרל
כל שיטת צעד אחד מקרבת את השטח מתחת לשיפוע בין tₙ ל־tₙ₊₁. ההבדל בין השיטות הוא בנקודות שבהן הן דוגמות את f.
y(tₙ₊₁)=y(tₙ)+∫[tₙ→tₙ₊₁] f(t,y(t))dtאוילר קדמי: שיפוע מתחילת הצעד
מחליפים את האינטגרל במלבן שגובהו השיפוע הידוע f(tₙ,yₙ). זהו עדכון מפורש וזול.
yₙ₊₁=yₙ+h f(tₙ,yₙ)אוילר אחורי: שיפוע מסוף הצעד
גובה המלבן נמדד במצב החדש, ולכן yₙ₊₁ מופיע בשני האגפים. במשוואה לינארית פותרים אלגברית; במשוואה לא־לינארית נדרש בדרך כלל פותר שורשים.
yₙ₊₁=yₙ+h f(tₙ₊₁,yₙ₊₁)שיטת הטרפז: ממוצע שני הקצוות
ממוצע השיפוע בתחילה ובסוף מקרב את השטח בטרפז. גם זו שיטה סמויה, אך הסימטריה מבטלת את איבר השגיאה מסדר שני.
yₙ₊₁=yₙ+(h/2)[f(tₙ,yₙ)+f(tₙ₊₁,yₙ₊₁)]סדר דיוק הוא קצב ירידת השגיאה
לאוילר קדמי ואחורי שגיאה מקומית O(h²) וגלובלית O(h). לטרפז שגיאה מקומית O(h³) וגלובלית O(h²), בהנחת פתרון חלק וחישוב יציב.
Euler: global O(h); trapezoid: global O(h²)יציבות נבדקת על אקספוננט
במשוואת הבדיקה y′=λy כל שיטה נעשית כפל ב־R(z), כאשר z=hλ. התנאי |R(z)|≤1 מבטיח חסימות; לדעיכה אסימפטוטית ממש נדרש |R(z)|<1.
R_F=1+z, R_B=1/(1−z), R_T=(1+z/2)/(1−z/2)גבול אוילר הקדמי
עבור λ ממשי שלילי, התנאי |1+hλ|≤1 שקול ל־0≤h|λ|≤2. קצב מהיר בעל |λ| גדול עלול לכפות צעד זעיר אף שכבר דעך כמעט לאפס.
λ<0: h≤2/|λ|A-stable אינו זהה לריסון חזק
אוילר אחורי והטרפז יציבים לכל Re(z)≤0. אך כאשר z→−∞, גורם אוילר אחורי שואף ל־0 ואילו גורם הטרפז שואף ל־−1; הטרפז עלול להשאיר תנודה מחליפת סימן של מצב מהיר במקום למחוק אותו.
z→−∞: R_B(z)→0, R_T(z)→−1עוצרים לחשוב: יציבות מונעת פיצוץ; היא אינה מבטיחה דיוק ואינה מבטיחה ריסון מהיר.
עכשיו רואים את השיטה בפעולה
דוגמה פתורה, צעד אחר צעד
פתרון דועך ושיטה שמתפוצצת
פתרו נומרית y′=−10y, y(0)=1 בארבעה צעדים של h=0.25 בשלוש השיטות.
כאן z=hλ=−2.5. אוילר קדמי נותן R_F=1−2.5=−1.5, ולכן |R_F|>1.
אחרי ארבעה צעדים y₄^F=(−1.5)⁴=5.0625: הערכים מחליפים סימן וגדלים, אף שהפתרון האמיתי דועך.
אוילר אחורי נותן R_B=1/(1+2.5)=2/7, ולכן y₄^B=(2/7)⁴≈0.00666.
הטרפז נותן R_T=(1−1.25)/(1+1.25)=−1/9, ולכן y₄^T=(−1/9)⁴≈0.000152.
בזמן t=1 הפתרון המדויק הוא e^(−10)≈0.0000454. שתי השיטות הסמויות יציבות, אך עדיין אינן מדויקות באותה מידה.
אוילר קדמי נכשל בגלל יציבות; אוילר אחורי והטרפז נשארים חסומים, אך יציבות לבדה אינה מבטיחה התאמה מלאה לערך המדויק.
R_F=−1.5, R_B=2/7, R_T=−1/9איפה זה פוגש את העולם
יישום
במודלים קשיחים מופיעים יחד קצבי דעיכה איטיים ומהירים. צעד שמתאים לתהליך האיטי עלול לפוצץ את אוילר הקדמי, בעוד שיטה סמויה שומרת את המצב המהיר תחת שליטה.
מלכודת נפוצה
מה כדאי לבדוק
לבלבל בין יציבות לדיוק: שיטה יכולה להיות יציבה אך לא מדויקת. באוילר אחורי אסור להציב f(tₙ,yₙ) במקום f(tₙ₊₁,yₙ₊₁), ובבעיה לא־לינארית ייתכן שצריך לפתור משוואה ואף לבחור בין כמה שורשים.
בודקים שהרעיון יושב
תרגול ומשוב
בדיקת הבנה
בחרו תשובה ורק אז בדקו את עצמכם
בחרו אחת מהאפשרויות כדי להמשיך.
עכשיו אתם
אתגר קצר
המערכת הקשיחה u′=−u, v′=−50v מתחילה ב־(1,1). עבור h=0.1 חשבו את גורמי ההגברה של שלוש השיטות לכל רכיב, וקבעו איזו שיטה מרסנת היטב את המצב המהיר.
רמז ראשון
הציבו פעם λ=−1 ופעם λ=−50 בנוסחאות R_F, R_B ו־R_T. השוו לא רק |R|≤1 אלא גם את גודל |R| עבור המצב המהיר.
הצגת דרך הפתרון
- לרכיב האיטי מתקבל hλ=−0.1: לכן R_F=0.9, R_B=1/1.1≈0.909 ו־R_T=0.95/1.05≈0.905; שלושתן יציבות.
- לרכיב המהיר מתקבל hλ=−5: אוילר קדמי נותן R_F=−4 ולכן מתפוצץ.
- אוילר אחורי נותן R_B=1/6≈0.167, ולכן מוחק במהירות את המצב שכבר אמור לדעוך.
- הטרפז נותן R_T=(1−2.5)/(1+2.5)=−3/7≈−0.429: הוא יציב אך משאיר תנודה מחליפת סימן.
- לבעיה הזאת אוילר אחורי הוא הבחירה המרוסנת ביותר בצעד הנתון; אם הדיוק במצב האיטי חשוב, עדיין יש לבדוק את השגיאה או להקטין h.
מימוש הנושא
קוד שאפשר לקחת לניסוי
הקוד מיישם את אותו רעיון מתמטי. ב־JavaScript אפשר להריץ אותו כאן; גרסת Python זמינה להשוואה ולהעתקה.
linearEuler מחשבת פעם אחת את גורמי ההגברה ואז מחילה אותם שוב ושוב, בדיוק כפי שהנוסחאות האנליטיות מנבאות.
עבור λ=−10, h=0.25 יודפסו הגורמים −1.5, 0.285714 ו־−0.111111; אחרי ארבעה צעדים הערכים יהיו בקירוב 5.0625, 0.00666 ו־0.000152.
הריצו שוב עם h=0.2 ואז h=0.1. בדקו מתי אוילר קדמי רק מחליף סימן, ומתי הוא גם דועך ומתקרב לפתרון המדויק.
מהרעיון אל הקוד
מימוש חישובי
function linearEuler(lambda, h, steps, y0=1) {
let forward = y0;
let backward = y0;
let trapezoid = y0;
const z = h * lambda;
const factors = {
forward: 1 + z,
backward: 1 / (1 - z),
trapezoid: (1 + z/2) / (1 - z/2),
};
const values = [{ step: 0, forward, backward, trapezoid }];
for (let step = 1; step <= steps; step += 1) {
forward *= factors.forward;
backward *= factors.backward;
trapezoid *= factors.trapezoid;
values.push({ step, forward, backward, trapezoid });
}
return { factors, values };
}
console.log(linearEuler(-10, .25, 4));def linear_euler(rate, h, steps, y0=1.0):
z = h*rate
factors = {
"forward": 1+z,
"backward": 1/(1-z),
"trapezoid": (1+z/2)/(1-z/2),
}
current = {name: y0 for name in factors}
values = [{"step": 0, **current}]
for step in range(1, steps+1):
current = {name: current[name]*factor
for name, factor in factors.items()}
values.append({"step": step, **current})
return factors, values
print(linear_euler(-10, .25, 4))לפני שממשיכים
מה לוקחים מהשיעור
- שלוש השיטות מקרבות את אותו אינטגרל באמצעות שיפוע בתחילה, בסוף או בממוצע.
- אוילר קדמי ואחורי הם מסדר ראשון; שיטת הטרפז היא מסדר שני.
- דיוק ויציבות הם שני מבחנים נפרדים: שיטה יציבה יכולה להיות לא מדויקת.
- אוילר קדמי מוגבל על דעיכה מהירה; אוילר אחורי והטרפז A-stable.
- בגבול קשיח אוילר אחורי מרסן מצבים מהירים, בעוד גורם הטרפז מתקרב ל־−1.
משימת מעבדה לשיעור הזה
תחנה 3 — דיוק מול יציבות
בחרו קצב דעיכה וגודל צעד, נבאו את גורם ההגברה של כל שיטה והשוו את התחזית למסלולים הנומריים.
מקורות להעמקה בלבד — השיעור עומד בפני עצמו
מדריך המרצים: עמ׳ 92–95
ספר הלימוד: עמ׳ 184–190