אנליזה נומרית/אינטגרציה נומרית

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

קפיצה אל: ניווט, חיפוש
כיתוב תמונה

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

\ I_nf(x)=\int\limits_x^{x+nh}f(t)dt \quad ,\quad n\in\mathbb{N}

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

\ DI_nf(x)=\int\limits_x^{x+nh}Df(t)dt= \int\limits_x^{x+nh}{df\over dt}dt= f(x+nh)-f(x)= (E^n-1)f(x)
\ \Rightarrow I_nf(x)=D^{-1}(E^n-1)f(x)


תוכן עניינים

[עריכה] אינטגרציה לפי הפרשים קידמיים

\ I_nf(x)={h\over \ln (1+\Delta)}\left[ (1+\Delta)^n-1 \right]f(x)

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

\ (1+\Delta)^n-1= 1+ n\Delta+{n(n-1)\over 2}\Delta^2+...-1

מאחר ואיננו יודעים את הטור של הכופל השמאלי, נשתמש בקשר \ {1\over \ln (1+\Delta)}\ln (1+\Delta)=1 , כאשר מניחים שהטור של הכופל השמאלי הוא מהצורה: \ {1\over \ln(1+\Delta)}= {a_0\over\Delta}+ a_1+ a_2\Delta+ a_3\Delta^2+.... נבצע השוואת מקדמים של החזקות של Δ ונקבל:

\ \left( {a_0\over\Delta}+ a_1+ a_2\Delta+ a_3\Delta^2+... \right)\left( \Delta- {\Delta^2\over 2}+ {\Delta^3\over 3}-... \right)= 1
\ \Rightarrow\ {1\over \ln (1+\Delta)}={1\over\Delta}+ {1\over 2}- {\Delta\over 12}+ {\Delta^2\over 24}- {19\over 720}\Delta^3+...
\ \Rightarrow I_nf(x)=h\left( {1\over\Delta}+ {1\over 2}- {\Delta\over 12}+... \right)\left[ n\Delta+{n(n-1)\over 2}\Delta^2+... \right]f(x)
\ \Rightarrow \int\limits_x^{x+nh}f(t)dt= nh\left[ 1+{n\over 2}\Delta+ {n(2n-3)\over 12}\Delta^2+ {n(n-2)^2\over 24}\Delta^3+ {n(6n^3-45n^2+110n-90)\over 720}\Delta^4+... \right]f(x)

נשים לב כי באפשרותנו כעת לקבוע שני פרמטרים (קיימות שתי דרגות חופש) - הן מספר המלבנים בחלוקה והן מספר האיברים בטור להתחשב בהם (n ו-h תלויים זה בזה). בחירות שונות שלהם יתנו דיוקים שונים. למשל, אם נבחר n=2 וניקח 3 איברים, נקבל שגיאה לפי \ O(h^5).

[עריכה] n=1

[עריכה] קירוב ראשון: סכום פשוט

LeftRiemann2.png

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

\ I_1f(x)=\int\limits_x^{x+h}f(x)dx= h\left( 1+{\Delta\over 2}- {\Delta^2\over 12}+... \right)f(x)\approx hf(x)+ \left\{ {h^2\over 2}f'(c) \right\}

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

[עריכה] קירוב שני: שיטת הטרפז

Sehnentrapezformel.svg
\ I_1f(x)=hf(x)+{h\Delta f(x)\over 2}+ \left\{ -{h^3\over 12}f''(c) \right\}=
\ =\left[ hf(x)+ {h\over 2}(f(x+h)-f(x)) \right]+ \left\{ -{h^3\over 12}f''(c) \right\}= {h\over 2}\left[ f(x)+f(x+h) \right]+ \left\{ -{h^3\over 12}f''(c) \right\}

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

קירוב זה נקרא קירוב טרפזי לאינטגרל ומסומן, במקרה הכללי, על ידי: \ I_n=\hat{T}(h)=\sum\limits_{i=1}^n{h\over 2}[f(x_i-h)+f(x_i)]. ניתן להראות כי ההתכנסות היא ריבועית, כלומר: \ \hat{T}(h)-I_n=O(h^2)=O\left({1\over n^2}\right).

[עריכה] שיטת הטרפז: הרחבה

נתבונן בשיטת אינטגרציה המבוססת על שיטת הטרפז:

\ \int_0^h f(x)dx= h[Af(0)+Bf(\theta h)]+ R\ ,\quad 0<\theta<1

שימו לב כי θ קובעת גובה הטרפז, וככל ש-θ קרובה ל-1 כך שטח הטרפז גדול יותר.

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

\ \mbox{LHS}: \int_0^h \left( f(0)+ xf'(0)+ \tfrac{x^2}{2}f''(0)+...\right)dx= \left.\left( xf(0)+ \tfrac{x^2}{2}f'(0)+ \tfrac{x^3}{3!}f''(0)+...\right)\right|_0^h= hf(0)+ \tfrac{h^2}{2}f'(0)+ \tfrac{h^3}{3!}f''(0)+...
\ \mbox{RHS}: hAf(0)+ hB \left[ f(0)+\theta h f'(0)+ \tfrac{(\theta h)^2}{2}f''(0)+ \tfrac{(\theta h)^3}{3!}f'''(0)+... \right]+ R

מהשוואת מקדמים של f',f נקבל: \ A=1-\tfrac{1}{2\theta},\ B=\tfrac{1}{2\theta} ואז:

\ \int_0^h f(x)dx= h\left[ \left(1-\tfrac{1}{2p}\right) f(0)+ \tfrac{1}{2p}f(\theta h) \right]+ R

מאחר ומצאנו ערכיהם של שני הקבועים A,B, למעשה קטענו את הטור טיילור של האניטגרל החל מהאיבר השלישי. לכן הביטוי לשארית (שגיאה) של אגף שמאל (באמצעות משפט לגראנז') הוא: \ R_{LHS}= \tfrac{h^3}{3!}f''(c), ואילו השארית (השגיאה) של הטור באגף ימין הוא: \ R_{RHS}= hB\left[\tfrac{(\theta h)^3}{3!}f'''(c)\right]. כעת אם נעביר אגפים על מנת לחלץ את R נקבל:

\ R= R_{LHS}-R_{RHS}= \tfrac{h^3}{3!}f''(c)- hB\left[\tfrac{(\theta h)^3}{3!}f'''(c)\right]= \tfrac{h^3}{12}(2-3\theta)f''(c)

כלומר השיטה היא מסדר שני, ועבור \ \theta=\tfrac{2}{3} נקבל שיטה מסדר 3 (במקרה זה יש לבצע הערכת שגיאה מחדש). כלומר אם נקח בחשבון רק שליש מהמרחק בין שתי נקודות עוקבות, ונשקלל את ערכי הפונקציה בשתי הנקודות העוקבות במשקלים \ \tfrac{1}{4},\tfrac{3}{4} נקבל קירוב טוב יותר.

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

נניח כי חילקנו את תחום האינטגרציה ל-N מלבנים בעלי רוחב h כל אחד. נחשב את האינטגרל המתאים לפי שיטת הטרפז:

\ I=\int\limits_0^{Nh}f(x)dx= \int\limits_0^hf(x)dx+ \int\limits_h^{2h}f(x)dx+ \int\limits_{2h}^{3h}f(x)dx+...+ \int\limits_{(N-1)h}^{Nh}f(x)dx\approx
\ \approx{h\over 2}[f(0)+f(h)]+ {h\over 2}[f(h)+f(2h)]+ {h\over 2}[f(2h)+f(3h)]+...+ {h\over 2}[f((N-1)h)+f(Nh)]

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

\ I=h\left[ {f(0)+f(Nh)\over 2}+ f(h)+ f(2h)+...+ f((N-1)h) \right]

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

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

\ |I_s-I_t|<\epsilon

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

\ I={h\over 2}[f(0)+f(h)]+ {h\over 2}[f(h)+f(2h)]+ {h\over 2}[f(2h)+f(3h)]+...+ {h\over 2}[f((N-1)h)+f(Nh)]+
\ -\left\{ {h^3\over 12}f''(c_1)+ {h^3\over 12}f''(c_2)+...+ {h^3\over 12}f''(c_N) \right\}

כאשר ci היא השגיאה בכל תת-תחום. השגיאה הכוללת היא:

\ R_N= {h^3\over 12}[f''(c_1)+ f''(c_2)+...+ f''(c_N)]

לכן נדרוש:

\ |R_N|\le {h^3\over 12}\left[ \left|f''(c_1)\right|+ \left|f''(c_2)\right|+...+ \left|f''(c_N)\right| \right]< \epsilon

על מנת לפשט את החישוב נגדיר: \ M=\max_{1\le i\le N}\{|f''(c_i)|\}, או לחילופין, כאשר f נתונה והיא אנליטית, ניתן פשוט לקחת את \ M=\max_{a\le x\le b}\{|f''(x)|\}. כך נקבל חסם על Rn ולכן גם על h,N בהתאם לביטוי הבא:

\ |R_N|\le{h^3\over 12}MN<\epsilon,\quad N={b-a\over h}\quad\Rightarrow\quad h< \sqrt{12\epsilon\over M(b-a)}\ ,\quad N>\sqrt{M(b-a)^3\over 12\epsilon}

מקרים פרטיים:

  • אם f מונוטונית עולה אז גם \ f'' מונוטונית עולה ואז \ M=|f''(b)|.
  • אם f מונוטונית יורדת אז גם \ f'' מונוטונית יורדת ואז \ M=|f''(a)|.
  • במקרה אחר (כללי) נמצא את M על ידי חישוב \ |f'''(x)|=0.


דרך אחרת למציאת h המתאים היא חוק ריצ'רדסון לאקסטרפולציה.

כאן מחשבים את ערך האינטרל לפי מלבן ברוחב h ולפי מלבן ברוחב h/2, ולפי השגיאה המתקבלת מחליטים אם לקחת h קטן יותר.

ערך האינטגרל לפי h:

\ \begin{matrix} \ & \bar{I}(h) & \ \\ I= & \overbrace{{h\over 2}[f(x_0)+f(x_1)]} & +\left\{ -{h^3\over 12}f''(c) \right\}\end{matrix}

ערך האינטגרל לפי h/2:

\ \begin{matrix} \ & \bar{I}({h\over 2}) & \ \\ I= & \overbrace{{h/2\over 2}[(f(x_0)+f(x_{1\over 2}))+ (f(x_{1\over 2})+f(x_1))]} & +\left\{ -{(h/2)^3\over 12}[f''(c_1)+ f''(c_2)] \right\}\end{matrix}

לשם פשטות, נניח כי הנגזרת השניה קבועה, ואז: \ k=f''(c)=f''(c_1)=f''(c_2), כך שמתקבלות שתי משוואות בשני נעלמים I, k:

\ \left\{\begin{array}{lcl} I & = & \bar{I}(h)- {h^3\over 12}k \\ I & = & \bar{I}({h\over 2})- {(h/2)^3\over 12}2k \end{array}\right.

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

\ \left\{\begin{array}{lcl} I & = & \bar{I}(h)+ \hat{k}h^R \\ I  & = & \bar{I}({h\over 2})+2\hat{k}\left({h\over 2}\right)^R \end{array}\right. \quad\Rightarrow\quad I=\frac{2^{R-1}\bar{I}\left({h\over 2}\right)- \bar{I}(h)}{2^{R-1} -1}

[עריכה] n=2

\ I_2f(x)=\int\limits_x^{x+2h}f(x)dx= 2h\left[ 1+\Delta+ {1\over 6}\Delta^2+0 -{1\over 180}\Delta^4+... \right]f(x)

[עריכה] קירוב ראשון: סכום פשוט

I_2f(x)= 2hf(x)+ \left\{ 2h^2f'(c) \right\}

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

[עריכה] קירוב שני: סכום פשוט דרך נקודת ביניים

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

Mittelpunktsregel.svg
\ I_2f(x)= [2hf(x)+ 2h\Delta f(x)]+ \left\{ {1\over 3}h^3f''(c) \right\}= 2hf(x+h)+ \left\{ {1\over 3}h^3f''(c) \right\}

[עריכה] קירוב שלישי: שיטת סימפסון

Simpsons method illustration.png
\ I_2f(x)=2h\left[ 1+\Delta+ {1\over 6}\Delta^2 \right]f(x)+ \left\{ -{h^5\over 90}f^{(4)}(c) \right\}=
\ 2h{1\over 6}[f(x)+ 4f(x+h)+ f(x+2h)]+ \left\{ -{1\over 90}h^5f^{(4)}(c) \right\}

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

נשים לב כי 2h הוא אורך תחום האינטגרל, אשר כופל סכום משוקלל של ערכי הפונקציה בתחום, וסכום המשקלים הוא 1. במקרה זה משתמשים ב-3 נקודות במשקלים \ \tfrac{1}{6}, \tfrac{4}{6}, \tfrac{1}{6}. כמו כן, מאחר והקירוב התקבל על ידי פולינום ממעלה 3, אז שיטה זו תתן תוצאה מדוייקת לכל פולינום עד מעלה 3.

שיקולי דיוק: גם כאן ניתן להפעיל את הכללים שראינו עבור n=1.

[עריכה] ראו גם

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


הפרק הקודם:
גזירה נומרית
אינטגרציה נומרית הפרק הבא:
יציבות נומרית