מתמטיקה שימושית ו־AᵀA

תנאי שפה והבדלים סופיים

Boundary Conditions Replace Initial Conditions

בבעיית שפה שני ערכים בקצות התחום מחליפים את נתוני ההתחלה. מטריצת הפרשים ראשונים A יוצרת את K=AᵀA, קירוב סימטרי ל־−u″.

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

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

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

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

להבחין בין תנאי התחלה לתנאי שפה, לפתור בעיית שפה רציפה, ולבנות ממנה מערכת הבדלים סופיים K=AᵀA תוך מעקב אחר h² וערכי הקצה.

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

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

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

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

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

האיור משווה את הפתרון המדויק של −u″=12x² לערכי ההבדלים הסופיים בין שני קצוות קבועים.

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

אחרי שראינו

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

לפתור −u″=f עם שני תנאי קצה באמצעות פתרון פרטי והומוגני.
לגזור את נוסחת ההפרש השני ולבנות את המטריצה התלת־אלכסונית.
להכניס ערכי שפה לצד ימין, לבדוק שגיאת רשת ולהסביר מתי תנאי Neumann יוצרים סינגולריות.

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

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

בעיית שפה: משוואה דיפרנציאלית עם תנאים בקצות התחום, למשל u(0),u(1).
תנאי Dirichlet: קביעת ערך הפונקציה בקצה, למשל u(0)=0.
תנאי Neumann: קביעת השיפוע או השטף בקצה, למשל u′(1)=0.
הבדל סופי: צירוף ערכי רשת המקרב נגזרת; ההפרש השני מקרב u″.
גודל רשת h: המרחק בין נקודות סמוכות; כאשר h קטן, הקירוב משתפר בדרך כלל.
מטריצת קשיחות K: המטריצה התלת־אלכסונית של −ההפרש השני; בתנאים מתאימים K=AᵀA.

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

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

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

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

קראו כל סימן כפעולה, לא כמשהו שצריך לשנן−u″=f; −uᵢ₋₁+2uᵢ−uᵢ₊₁=h²fᵢ

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

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

פותרים תחילה ברצף
אינטגרציה פעמיים של −u″=f נותנת פתרון פרטי ועוד Ax+B. שני תנאי השפה קובעים את שני הקבועים.
u=uₚ+A+Bx
בונים קירוב לנגזרת שנייה
פיתוחי טיילור של u(xᵢ+h) ושל u(xᵢ−h) מבטלים את הנגזרת הראשונה ומשאירים את הנגזרת השנייה.
u″(xᵢ)≈(uᵢ₋₁−2uᵢ+uᵢ₊₁)/h²
מתאימים את סימן המשוואה
עבור −u″=f הופכים את הסימנים ומכפילים ב־h². לכן האלכסון הוא 2 והשכנים מקבלים −1.
−uᵢ₋₁+2uᵢ−uᵢ₊₁=h²f(xᵢ)
מפרידים נעלמים מערכי שפה
רק ערכי הפנים הם נעלמים. כאשר שכן נמצא על השפה, ערכו הידוע עובר לצד ימין בסימן חיובי.
KU=h²F+boundary contributions
רואים את AᵀA
A מודדת הפרשים ראשונים בין נקודות סמוכות. האנרגיה ‖AU‖² היא סכום ריבועי השיפועים, ו־AᵀA מחזירה את ההפרש השני.
K=AᵀA; UᵀKU=‖AU‖²
בודקים סדר דיוק
שגיאת הקיטוע של ההפרש המרכזי היא מסדר h² עבור פונקציה חלקה. חציית h אמורה להקטין את השגיאה בקירוב פי 4.
max error=O(h²)
מבינים סינגולריות בתנאי שיפוע
אם קובעים רק שיפועים בשני הקצוות, הוספת קבוע ל־u אינה משנה את המשוואה או את התנאים. לכן K מחזיקה את וקטור הקבועים בגרעין, ונדרש תנאי התאמה למקור.
K𝟙=0; ∫f=boundary flux difference

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

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

דוגמה 1

פתרון רציף מדוגמת המקור

פתרו −u″=12x² עם u(0)=u(1)=0.

פתרון פרטי הוא uₚ=−x⁴, כי −uₚ″=12x².
הפתרון ההומוגני של −u″=0 הוא A+Bx.
מתנאי u(0)=0 מקבלים A=0.
מתנאי u(1)=0 מקבלים −1+B=0, ולכן B=1.
בדיקה: u=−x⁴+x,‏ u(0)=u(1)=0 ו־−u″=12x².

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

u(x)=−x⁴+x.

דוגמה 2

שלוש נקודות פנים

קרבו את אותה בעיה בנקודות x=1/4,1/2,3/4.

כאן h=1/4 ו־U=(u₁,u₂,u₃).
K=[[2,−1,0],[−1,2,−1],[0,−1,2]].
צד ימין הוא h²·12(x₁²,x₂²,x₃²), מפני ששני ערכי השפה אפס.
פותרים KU=h²F ומקבלים קירוב לערכי −x⁴+x בשלוש הנקודות.
הצבה חזרה ב־KU−h²F בודקת שהשארית האלגברית כמעט אפס; ההפרש מן העקומה המדויקת הוא שגיאת הדיסקרטיזציה.

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

מערכת תלת־אלכסונית מחליפה בעיית שפה רציפה במערכת לינארית סימטרית וחיובית.

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

יישום

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

מלכודת נפוצה

מה כדאי לבדוק

לשכוח את h² בצד ימין. המטריצה [−1,2,−1] מייצגת h²(−u″), לא את −u″ עצמו. בלי הגורם, שינוי מספר נקודות הרשת משנה את הבעיה הפיזיקלית.

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

תרגול ומשוב

בדיקת הבנה

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

במשוואה −u″=f, מדוע האלכסון של מטריצת ההפרשים הוא 2?

כי −u″≈(−uᵢ₋₁+2uᵢ−uᵢ₊₁)/h²כי יש שני תנאי שפהכי f מוכפל ב־2

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

עכשיו אתם

אתגר קצר

פתרו −u″=2 עם u(0)=1,u(1)=0 ברצף, ואז כתבו את מערכת ההפרשים לשתי נקודות פנים.

רמז ראשון

אינטגרציה נותנת u=−x²+A+Bx. ברשת h=1/3, והערך השמאלי 1 נוסף לרכיב הראשון של צד ימין.

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

מן התנאים מתקבל A=1 ו־B=0, ולכן u(x)=1−x².
K=[[2,−1],[−1,2]] ו־h²F=(2/9,2/9).
תרומת השפה היא (1,0), ולכן KU=(11/9,2/9).
פתרון המערכת נותן U=(8/9,5/9), בדיוק ערכי 1−x² בנקודות 1/3 ו־2/3.

מימוש הנושא

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

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

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

הקוד בונה את K לכל מספר נקודות פנים, פותר באלימינציה ומשווה לפתרון המדויק −x⁴+x.

פלט צפוי

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

ניסוי המשך

הריצו עם 3,7 ו־15 נקודות פנים. האם חציית h מקטינה את השגיאה בקירוב פי 4?

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

מימוש חישובי

function solve(A,b) {
  const M=A.map((row,i)=>[...row,b[i]]), n=A.length;
  for(let c=0;c<n;c++) {
    let p=c; for(let r=c+1;r<n;r++) if(Math.abs(M[r][c])>Math.abs(M[p][c])) p=r;
    [M[c],M[p]]=[M[p],M[c]];
    const d=M[c][c]; for(let j=c;j<=n;j++) M[c][j]/=d;
    for(let r=0;r<n;r++) if(r!==c) { const f=M[r][c]; for(let j=c;j<=n;j++) M[r][j]-=f*M[c][j]; }
  }
  return M.map(row=>row[n]);
}

function boundarySolution(N) {
  const h=1/(N+1), x=Array.from({length:N},(_,i)=>(i+1)*h);
  const K=Array.from({length:N},(_,i)=>Array.from({length:N},(_,j)=>i===j?2:Math.abs(i-j)===1?-1:0));
  const rhs=x.map(value=>h*h*12*value*value);
  const u=solve(K,rhs), exact=x.map(value=>-(value**4)+value);
  return {h,x,u,exact,maxError:Math.max(...u.map((value,i)=>Math.abs(value-exact[i])))};
}

console.log(boundarySolution(7));

def solve(A,b):
    M=[row[:] + [b[i]] for i,row in enumerate(A)]
    n=len(A)
    for c in range(n):
        p=max(range(c,n),key=lambda r:abs(M[r][c]))
        M[c],M[p]=M[p],M[c]
        d=M[c][c]
        M[c]=[value/d for value in M[c]]
        for r in range(n):
            if r!=c:
                f=M[r][c]
                M[r]=[M[r][j]-f*M[c][j] for j in range(n+1)]
    return [row[n] for row in M]

def boundary_solution(N):
    h=1/(N+1)
    x=[(i+1)*h for i in range(N)]
    K=[[2 if i==j else -1 if abs(i-j)==1 else 0 for j in range(N)] for i in range(N)]
    rhs=[h*h*12*value*value for value in x]
    u=solve(K,rhs)
    exact=[-value**4+value for value in x]
    return {"h":h,"x":x,"u":u,"exact":exact,
            "max_error":max(abs(value-exact[i]) for i,value in enumerate(u))}

print(boundary_solution(7))

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

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

בעיית שפה משתמשת בשני קצוות; סוג התנאים קובע קיום ויחידות.
הפרש שני מרכזי יוצר מטריצה סימטרית מקומית עם 2 באלכסון ו־−1 לשכנים.
K=AᵀA מחברת בין שיפועים, אנרגיה ועקמומיות.
ערכי שפה עוברים לצד ימין, ו־h² חייב להישמר כדי לקרב את אותה משוואה.

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

תחנה 3 — מצופפים את הרשת

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

לפתיחת התחנה

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

מדריך המרצים: עמ׳ 204–207

ספר הלימוד: עמ׳ 406–415