שיטות גרפיות ונומריות

דיוק גבוה עם Runge–Kutta

Higher Accuracy with Runge–Kutta

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

יחידת לימוד עצמאיתשיעור 5 מתוך 5כ־70 דקות כולל תרגול4 מטרות למידה

1
נקודת פתיחהמה כבר ידוע ומה השאלה?
2
ניבוימה יקרה במקרה הפשוט?
3
המחשהמה באמת משתנה בתמונה?
4
ניסוח והצדקהשם, נוסחה ולמה היא נכונה
5
בדיקהדוגמה, תרגול ומשוב
6
העברהמקרה חדש והמשך המסלול

שאלה לפני נוסחה

מה אנחנו מנסים לראות?

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

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

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

לראות לפני שמחשבים · המקרה הפשוט ביותר

חציית הצעד חושפת את סדר השיטה

עקבו מהצורה הראשונית אל התוצאה

ההמחשה משווה בסולם לוגריתמי את השגיאה של Euler,‏ midpoint/RK2 ו־RK4 כאשר מקטינים את h בפתרון y′=y עד t=1.

מה כדאי לגלות: עקבו אחר קצב הירידה: בכל חצייה של h השגיאה קטנה בקירוב פי 2, פי 4 או פי 16. השיפוע התלול של RK4 הוא המשמעות המעשית של סדר רביעי.

אחרי שראינו

בסוף השיעור תוכלו

לחשב את ארבעת שלבי RK4 לפי הסדר ולבצע צעד מלא עבור משוואה סקלרית או מערכת.
להסביר כיצד משקלי 1,2,2,1 משלבים שיפועים מתחילת הצעד, מאמצעו ומסופו.
להבחין בין שגיאה מקומית O(h⁵) לשגיאה גלובלית O(h⁴) ולבדוק סדר באמצעות חציית הצעד.
להשתמש בפולינום היציבות כדי להבין מדוע סדר גבוה אינו מבטל מגבלות יציבות.

נותנים שמות למה שראינו

מושגים שנשתמש בהם

שלב: הערכת שיפוע פנימית kᵢ שנבנית מנקודת זמן ומצב זמני בתוך הצעד.
midpoint / RK2: שיטה מפורשת בת שני שלבים: מודדים k₁ בתחילה, משתמשים בו כדי להגיע לאמצע, ואז מעדכנים בעזרת שיפוע האמצע k₂; הסדר הגלובלי שלה הוא 2.
RK4: שיטת Runge–Kutta מפורשת בת ארבעה שלבים ובעלת סדר גלובלי רביעי.
סדר שיטה: החזקה p בקשר E(h)=O(hᵖ) לשגיאה הגלובלית עד זמן קבוע.
סימון O(hᵖ): חסם אסימפטוטי: עבור h קטן מספיק, גודל השגיאה אינו עולה על C·hᵖ עבור קבוע C שאינו תלוי ב־h.
שגיאת קטיעה מקומית: הפער בצעד אחד כאשר מתחילים מן הפתרון המדויק; ב־RK4 היא O(h⁵).
פולינום יציבות: הפונקציה R(z) שמכפילה את הערך בצעד אחד של y′=λy.
כלל סימפסון: קירוב אינטגרל במשקלים 1,4,1 בנקודת ההתחלה, האמצע והסוף; הוא מופיע מ־RK4 כאשר f תלויה רק בזמן.

מהמילים אל הסימנים

כותבים את הפעולה שכבר ראינו

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

להצגת נוסחת המפתח

קראו כל סימן כפעולה, לא כמשהו שצריך לשנןyₙ₊₁ = yₙ + h(k₁+2k₂+2k₃+k₄)/6

מרכיבים את הפעולות

מהתמונה אל הנוסחה — הגזירה, בלי לדלג על הצעדים

גשר קצר: שיטת האמצע היא RK2
לפני ארבעה שלבים, אפשר לשפר את אוילר בעזרת שניים: k₁ מנבא את מצב האמצע ו־k₂ הוא השיפוע שבו מעדכנים את כל הצעד. זו שיטה מפורשת מסדר שני, והיא קו ההשוואה שבתרשים.
k₁=f(tₙ,yₙ), k₂=f(tₙ+h/2,yₙ+hk₁/2), yₙ₊₁=yₙ+hk₂
השיפוע הראשון יוצא מן המצב הידוע
k₁ הוא אותו שיפוע שבו משתמש אוילר קדמי. הוא יוצר תחזית ראשונה לנקודת האמצע.
k₁=f(tₙ,yₙ)
השיפוע השני בודק את אמצע הצעד
מתקדמים חצי צעד בזמן וחצי תיקון במצב על פי k₁. כך k₂ כבר מרגיש שינוי בשיפוע לאורך המסלול.
k₂=f(tₙ+h/2, yₙ+hk₁/2)
השיפוע השלישי מתקן את האמצע
במקום לסמוך על תחזית האמצע הראשונה, בונים תחזית חדשה בעזרת k₂. זהו שלב נפרד גם אם זמן ההערכה זהה.
k₃=f(tₙ+h/2, yₙ+hk₂/2)
השיפוע הרביעי מביט בקצה
k₃ משמש לתחזית של צעד מלא, ושם מודדים את k₄. ארבעת השיפועים כעת דוגמים את כל הצעד.
k₄=f(tₙ+h, yₙ+hk₃)
ממוצע משוקלל מעדכן את הפתרון
שני שיפועי האמצע מקבלים משקל כפול. המשקלים נבחרו כך שפיתוח טיילור של העדכון יתאים לפתרון המדויק עד איברי h⁴.
yₙ₊₁=yₙ+(h/6)(k₁+2k₂+2k₃+k₄)
משוואת הבדיקה חושפת את פולינום היציבות
עבור y′=λy כל שלב הוא כפל אלגברי. מתקבל פולינום שהוא חיתוך טור eᶻ לאחר החזקה הרביעית, כאשר z=hλ.
R₄(z)=1+z+z²/2+z³/6+z⁴/24
מקומי מסדר חמש, גלובלי מסדר ארבע
בצעד יחיד חסר האיבר הבא בטור ולכן השגיאה O(h⁵). עד זמן T מבצעים בערך T/h צעדים; הצטברות יציבה מורידה חזקה אחת ונותנת O(h⁴).
τₙ=O(h⁵), E(T)=O(h⁴)
חציית הצעד היא ניסוי סדר
אם האיבר המוביל הוא Ch⁴, מעבר מ־h ל־h/2 מקטין את השגיאה בערך פי 2⁴=16. עבור צעדים לא קטנים היחס רק מתקרב ל־16.
E(h)/E(h/2)≈16
הקשר לכלל סימפסון ולמערכות
כאשר f תלויה רק בזמן, שלבי האמצע מתאחדים והמשקלים נעשים 1,4,1 של סימפסון. עבור מערכת וקטורית אותן נוסחאות פועלות רכיב־רכיב עם kᵢ וקטוריים.
∫[tₙ→tₙ+h]f(t)dt ≈ (h/6)[f(tₙ)+4f(tₙ+h/2)+f(tₙ+h)]

עכשיו רואים את השיטה בפעולה

דוגמה פתורה, צעד אחר צעד

דוגמה 1

צעד RK4 מלא ביד

קרבו את y(1) עבור y′=−y,‏ y(0)=1 בצעד יחיד h=1.

k₁=f(0,1)=−1.
k₂=f(1/2,1+k₁/2)=−(1−1/2)=−1/2.
k₃=f(1/2,1+k₂/2)=−(1−1/4)=−3/4.
k₄=f(1,1+k₃)=−(1−3/4)=−1/4.
העדכון הוא y₁=1+(1/6)[−1+2(−1/2)+2(−3/4)−1/4]=3/8=0.375.
הערך המדויק e^(−1)≈0.367879, ולכן השגיאה המוחלטת היא כ־0.00712.

התוצאה ומה היא אומרת

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

y₁^RK4=0.375, e^(−1)≈0.367879

איפה זה פוגש את העולם

יישום

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

מלכודת נפוצה

מה כדאי לבדוק

לשכוח שכבר מכפילים את kᵢ ב־h רק בעדכון ובנקודות הביניים, או להשתמש באותו y לכל ארבעת השיפועים. סדר רביעי גם אינו מבטיח יציבות לצעד גדול במערכת קשיחה.

בודקים שהרעיון יושב

תרגול ומשוב

בדיקת הבנה

בחרו תשובה ורק אז בדקו את עצמכם

איזו נוסחה מגדירה נכון את k₃ ב־RK4 הקלאסית?

f(tₙ+h/2,yₙ+hk₂/2)f(tₙ+h,yₙ+hk₂)f(tₙ+h/2,yₙ+hk₁/2)

בחרו אחת מהאפשרויות כדי להמשיך.

עכשיו אתם

אתגר קצר

פתרו את y′=−y עד T=1 עם RK4 בצעדים h=0.5,‏ 0.25,‏ 0.125 ו־0.0625. חשבו את השגיאה מול e^(−1) ואת היחס בין שגיאות עוקבות. כיצד התוצאות מאמתות סדר רביעי?

רמז ראשון

למשוואה הזאת כל צעד מכפיל ב־R₄(−h). לכן y(1)=[R₄(−h)]^(1/h), ואין צורך לחשב בנפרד את כל ארבעת השלבים בכל פעם.

הצגת דרך הפתרון

ב־h=0.5 מתקבל y≈0.368170844 ושגיאה 2.914×10⁻⁴.
ב־h=0.25 מתקבל y≈0.367894199 ושגיאה 1.476×10⁻⁵; היחס הקודם הוא כ־19.7.
ב־h=0.125 מתקבל y≈0.367880272 ושגיאה 8.308×10⁻⁷; היחס הוא כ־17.8.
ב־h=0.0625 מתקבל y≈0.367879490 ושגיאה 4.928×10⁻⁸; היחס הוא כ־16.9.
ככל שהצעד קטן, איברי השגיאה הגבוהים נעשים זניחים והיחס מתקרב ל־16=2⁴, כצפוי משגיאה גלובלית מסדר רביעי.

מימוש הנושא

קוד שאפשר לקחת לניסוי

הקוד מיישם את אותו רעיון מתמטי. ב־JavaScript אפשר להריץ אותו כאן; גרסת Python זמינה להשוואה ולהעתקה.

הקשר למתמטיקה

rk4Step מממש ישירות את ארבעת השלבים, ו־integrateRK4 חוזרת על הצעד עד זמן היעד. השוואת כמה ערכי h הופכת את טענת הסדר לניסוי מספרי.

פלט צפוי

עבור h=0.5 יתקבל y(1)≈0.368170844. כאשר חוצים שוב ושוב את h, יחסי השגיאה יתקרבו ל־16.

ניסוי המשך

החליפו את האגף הימני ב־−20y והריצו בצעדים h=0.2 ו־0.05. האם שיטה מסדר רביעי יכולה עדיין להיות לא־יציבה בצעד גדול?

מהרעיון אל הקוד

מימוש חישובי

function rk4Step(f, t, y, h) {
  const k1 = f(t, y);
  const k2 = f(t + h/2, y + h*k1/2);
  const k3 = f(t + h/2, y + h*k2/2);
  const k4 = f(t + h, y + h*k3);
  return y + h*(k1 + 2*k2 + 2*k3 + k4)/6;
}

function integrateRK4(f, t0, y0, h, steps) {
  const points = [{ t: t0, y: y0 }];
  let t = t0;
  let y = y0;
  for (let step = 0; step < steps; step += 1) {
    y = rk4Step(f, t, y, h);
    t += h;
    points.push({ t, y });
  }
  return points;
}

const exact = Math.exp(-1);
const results = [.5, .25, .125, .0625].map((h) => {
  const points = integrateRK4((_t, y) => -y, 0, 1, h, Math.round(1/h));
  const y = points.at(-1).y;
  return { h, y, error: Math.abs(y - exact) };
});

console.table(results.map((row, index) => ({
  ...row,
  errorRatio: index === 0 ? null : results[index - 1].error / row.error,
})));

from math import exp

def rk4_step(f, t, y, h):
    k1 = f(t, y)
    k2 = f(t+h/2, y+h*k1/2)
    k3 = f(t+h/2, y+h*k2/2)
    k4 = f(t+h, y+h*k3)
    return y + h*(k1+2*k2+2*k3+k4)/6

def integrate_rk4(f, t0, y0, h, steps):
    points = [(t0, y0)]
    t, y = t0, y0
    for _ in range(steps):
        y = rk4_step(f, t, y, h)
        t += h
        points.append((t, y))
    return points

exact = exp(-1)
results = []
for h in (0.5, 0.25, 0.125, 0.0625):
    points = integrate_rk4(lambda _t, y: -y, 0, 1, h, round(1/h))
    y = points[-1][1]
    results.append({"h": h, "y": y, "error": abs(y-exact)})

for index, row in enumerate(results):
    ratio = None if index == 0 else results[index-1]["error"]/row["error"]
    print({**row, "error_ratio": ratio})

לפני שממשיכים

מה לוקחים מהשיעור

RK4 מעריכה ארבעה שיפועים תלויים זה בזה בתוך כל צעד.
הממוצע במשקלים 1,2,2,1 מתאים לטיילור עד סדר רביעי בלי לגזור את f.
השגיאה המקומית היא O(h⁵), אך השגיאה המצטברת עד זמן קבוע היא O(h⁴).
חציית h אמורה להקטין את השגיאה הגלובלית בערך פי 16 באזור האסימפטוטי.
סדר גבוה אינו מבטל מגבלת יציבות; RK4 עדיין שיטה מפורשת ואינה הבחירה הטבעית לכל בעיה קשיחה.

משימת מעבדה לשיעור הזה

תחנה 4 — ארבעה שלבים וסדר התכנסות

עברו שלב־שלב בתוך צעד RK4, ואז חצו את h ובדקו אם יחס השגיאות מתקרב לערך 16 שמנבא הסדר הרביעי.

לפתיחת התחנה

מקורות להעמקה בלבד — השיעור עומד בפני עצמו

מדריך המרצים: עמ׳ 95–98

ספר הלימוד: עמ׳ 191–196

מה אנחנו מנסים לראות?

חציית הצעד חושפת את סדר השיטה

בסוף השיעור תוכלו

מושגים שנשתמש בהם

כותבים את הפעולה שכבר ראינו

מהתמונה אל הנוסחה — הגזירה, בלי לדלג על הצעדים

גשר קצר: שיטת האמצע היא RK2

השיפוע הראשון יוצא מן המצב הידוע

השיפוע השני בודק את אמצע הצעד

השיפוע השלישי מתקן את האמצע

השיפוע הרביעי מביט בקצה

ממוצע משוקלל מעדכן את הפתרון

משוואת הבדיקה חושפת את פולינום היציבות

מקומי מסדר חמש, גלובלי מסדר ארבע

חציית הצעד היא ניסוי סדר

הקשר לכלל סימפסון ולמערכות