הצגת המקור של הדף "אנליזה נומרית/אינטגרציה נומרית"

[[תמונה:Finite Differences.png|שמאל|ממוזער|כיתוב תמונה]]
כרגיל, תחום האינטגרציה מחולק למלבנים אשר יש לחשב את שטחם. אם מסכמים n מלבנים, נגדיר את הפעולה הבאה:
:<math display=block>I_nf(x)=\int\limits_x^{x+nh}f(t)dt\quad ,\quad n\in\mathbb{N}</math>
עלינו למצוא קשר אופרטורי בין <math>\ I_n</math> לבין האופרטורים שפיתחנו, על מנת לייצג סכום כנ"ל. נבצע לכן את הפעולות הבאות:
:<math display=block>\begin{align}DI_nf(x)&=\int\limits_x^{x+nh}Df(t)dt=\int\limits_x^{x+nh}\frac{df}{dt}dt=f(x+nh)-f(x)=(E^n-1)f(x)\\&\Rarr I_nf(x)=D^{-1}(E^n-1)f(x)\end{align}</math>
{{אנליזה נומרית}}

==אינטגרציה לפי הפרשים קדמיים==
:<math display=block>I_nf(x)=\frac{h}{\ln(1+\Delta)}\bigl[(1+\Delta)^n-1\bigr]f(x)</math>
ההמשך הוא לפתח לטור טיילור ולקחת מקדמים ראשונים כקירוב. הטור של הכופל הימני הוא:
:<math display=block>(1+\Delta)^n-1=1+n\Delta+\frac{n(n-1)}{2}\Delta^2+\cdots-1</math>
מאחר ואיננו יודעים את הטור של הכופל השמאלי, נשתמש בקשר <math>\frac{1}{\ln(1+\Delta)}\ln(1+\Delta)=1</math> , כאשר מניחים שהטור של הכופל השמאלי הוא מהצורה:
:<math display=block>\frac{1}{\ln(1+\Delta)}=\frac{a_0}{\Delta}+a_1+a_2\Delta+a_3\Delta^2+\cdots</math>
נבצע השוואת מקדמים של החזקות של <math>\Delta</math> ונקבל:
:<math display=block>\begin{align}
&\left(\frac{a_0}{\Delta}+a_1+a_2\Delta+a_3\Delta^2+\cdots\right)\left(\Delta-\frac{\Delta^2}{2}+\frac{\Delta^3}{3}-\cdots\right)=1\\
&\Rarr\quad\frac{1}{\ln(1+\Delta)}=\frac{1}{\Delta}+\frac12-\frac{\Delta}{12}+\frac{\Delta^2}{24}-\frac{19}{720}\Delta^3+\cdots\\
&\Rarr\quad I_nf(x)=h\left(\frac{1}{\Delta}+\dfrac12-\frac{\Delta}{12}+\cdots\right)\left[n\Delta+\frac{n(n-1)}{2}\Delta^2+\cdots\right]f(x)\\
&\Rarr\quad\int\limits_x^{x+nh}f(t)dt=nh\left[1+\frac{n}{2}\Delta+\frac{n(2n-3)}{12}\Delta^2+\frac{n(n-2)^2}{24}\Delta^3+\frac{n(6n^3-45n^2+110n-90)}{720}\Delta^4+\cdots\right]f(x)
\end{align}</math>
נשים לב כי באפשרותנו כעת לקבוע שני פרמטרים (קיימות שתי דרגות חופש) – הן מספר המלבנים בחלוקה והן מספר האיברים בטור להתחשב בהם (<math>h,n</math> תלויים זה בזה). בחירות שונות שלהם יתנו דיוקים שונים. למשל, אם נבחר <math>n=2</math> וניקח 3 איברים, נקבל שגיאה לפי <math>O(h^5)</math> .

==n=1==
===קירוב ראשון: סכום פשוט===
[[תמונה:LeftRiemann2.png|שמאל|ממוזער|250 פיקסלים]]
כרגיל, נבטא את השגיאה באמצעות משפט לגראנז' ונסמנה בסוגריים מסולסלים:
:<math>I_1f(x)=\int\limits_x^{x+h}f(x)dx=h\left(1+\frac{\Delta}{2}-\frac{\Delta^2}{12}+\cdots\right)f(x)\approx hf(x)+\left\{\frac{h^2}{2}f'(c)\right\}</math>
המשמעות הגאומטרית היא מלבן אחד הנלקח לפי הנקודה השמאלית. אם הפונקציה מונוטונית עולה – זהו סכום תחתון, ואם היא מונוטונית יורדת – סכום עליון.

===קירוב שני: שיטת הטרפז===
[[תמונה:Sehnentrapezformel.svg|שמאל|ממוזער|250 פיקסלים]]
:<math display=block>\begin{align}I_1f(x)&=hf(x)+\frac{h\Delta f(x)}{2}+\left\{-\frac{h^3}{12}f''(c)\right\}\\&=\left[hf(x)+\frac{h}{2}\bigl[f(x+h)-f(x)\bigr]\right]+\left\{-\frac{h^3}{12}f''(c)\right\}\\&=\frac{h}{2}\bigl[f(x)+f(x+h)\bigr]+\left\{-\frac{h^3}{12}f''(c)\right\}\end{align}</math>
נשים לב כי <math>h</math> הוא אורך תחום האינטגרל, אשר כופל סכום משוקלל (במקרה זה יש שתי נקודות ולכן המשקלים שווים ל־0.5) של ערכי הפונקציה בתחום, וסכום המשקלים הוא 1.

קירוב זה נקרא '''קירוב טרפזי לאינטגרל''' ומסומן, במקרה הכללי, על ידי:
:<math>I_n=\hat T(h)=\sum_{i=1}^n\frac{h}{2}\bigl[f(x_i-h)+f(x_i)\bigr]</math>
ניתן להראות כי ההתכנסות היא ריבועית, כלומר <math>\hat T(h)-I_n=O(h^2)=O\left(\frac{1}{n^2}\right)</math> .

===שיטת הטרפז: הרחבה===
נתבונן בשיטת אינטגרציה המבוססת על שיטת הטרפז:
:<math>\int\limits_0^hf(x)dx=h\bigl[Af(0)+Bf(\theta h)\bigr]+R\ ,\ 0<\theta<1</math>
שימו לב כי <math>\theta</math> קובעת גובה הטרפז, וככל שהיא מתקרבת ל־1 כך שטח הטרפז גדול יותר.

ננסה למצוא את הקבועים כפונקציה של <math>\theta</math> . לשם כך נפתח את את שני האגפים לטור טיילור ונשווה מקדמים:
:<math display=block>\begin{align}\text{LHS}:&\quad\int\limits_0^h\left(f(0)+xf'(0)+\frac{x^2}{2!}f''(0)+\cdots\right)dx=\left(xf(0)+\frac{x^2}{2!}f'(0)+\frac{x^3}{3!}f''(0)+\cdots\right)\Bigg|_0^h=hf(0)+\frac{h^2}{2!}f'(0)+\frac{h^3}{3!}f''(0)+\cdots\\\text{RHS}:&\quad hAf(0)+hB\left[f(0)+\theta hf'(0)+\frac{(\theta h)^2}{2!}f''(0)+\frac{(\theta h)^3}{3!}f'''(0)+\cdots\right]+R\end{align}</math>
מהשוואת מקדמים של f',f נקבל: <math>A=1-\tfrac{1}{2\theta},\ B=\tfrac{1}{2\theta}</math> ואז:
:<math>\int\limits_0^hf(x)dx=h\left[\left(1-\tfrac{1}{2p}\right)f(0)+\tfrac{1}{2p}f(\theta h)\right]+R</math>
מאחר ומצאנו ערכיהם של שני הקבועים A,B, למעשה קטענו את הטור טיילור של האניטגרל החל מהאיבר השלישי. לכן הביטוי לשארית (שגיאה) של אגף שמאל (באמצעות משפט לגראנז') הוא: <math>\ R_{LHS}= \tfrac{h^3}{3!}f''(c)</math>, ואילו השארית (השגיאה) של הטור באגף ימין הוא: <math>\ R_{RHS}= hB\left[\tfrac{(\theta h)^3}{3!}f'''(c)\right]</math>.
כעת אם נעביר אגפים על מנת לחלץ את R נקבל:
:<math>\ 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)</math>

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

===שיקולי יעילות===
נניח כי חילקנו את תחום האינטגרציה ל-N מלבנים בעלי רוחב h כל אחד. נחשב את האינטגרל המתאים לפי שיטת הטרפז:
<div style="direction: ltr;">
:<math>\ 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</math>
:<math>\ \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)]</math>
</div>
נשים לב כי בביטוי הנ"ל, פרט לשתי הנקודות הקיצוניות, מחשבים כל נקודה פעמיים. לכן על מנת לחסוך כמחצית מזמן החישוב, נכתוב ביטוי יעיל יותר:
<div style="direction: ltr;">
<math>\ I=h\left[ {f(0)+f(Nh)\over 2}+ f(h)+ f(2h)+...+ f((N-1)h) \right]</math>
</div>

===שיקולי דיוק===
כאשר נדרש דיוק &epsilon; בחישוב האינטגרל, עלינו לבחור ב-h (או לחילופין - N) מתאים.
הצורה הכללית ביותר לבדיקת התכנסות היא אנאלוגית לבדיקת התכנסות של סדרת סכומים חלקיים:
:<math>\ |I_s-I_t|<\epsilon</math>
שיטה זו היא מסורבלת ולכן לשם הדיון, נפתח שיטה ישירה אחרת אשר תתייחס לשיטת הטרפז.
:<math>\ 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)]+</math>
:<math>\ -\left\{ {h^3\over 12}f''(c_1)+ {h^3\over 12}f''(c_2)+...+ {h^3\over 12}f''(c_N) \right\}</math>
כאשר c<sub>i</sub> היא השגיאה בכל תת-תחום. השגיאה הכוללת היא:
<div style="direction: ltr;">
:<math>\ R_N= {h^3\over 12}[f''(c_1)+ f''(c_2)+...+ f''(c_N)]</math>
</div>
לכן נדרוש:
<div style="direction: ltr;">
:<math>\ |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</math>
</div>
על מנת לפשט את החישוב נגדיר: <math>\ M=\max_{1\le i\le N}\{|f''(c_i)|\}</math>, או לחילופין, כאשר f נתונה והיא אנליטית, ניתן פשוט לקחת את  <math>\ M=\max_{a\le x\le b}\{|f''(x)|\}</math>. כך נקבל חסם על R<sub>n</sub> ולכן גם על h,N בהתאם לביטוי הבא:
<div style="direction: ltr;">
:<math>\ |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}</math>
</div>

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


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

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

ערך האינטגרל לפי h:
<div style="direction: ltr;">
:<math>\ \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}</math>
</div>
ערך האינטגרל לפי h/2:
<div style="direction: ltr;">
:<math>\ \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}</math>
</div>
לשם פשטות, נניח כי הנגזרת השניה קבועה, ואז: <math>\ k=f''(c)=f''(c_1)=f''(c_2)</math>, כך שמתקבלות שתי משוואות בשני נעלמים I, k:
<div style="direction: ltr;">
:<math>\ \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.</math>
</div>
כאשר נפתור את המשוואות, נקבל את ערכם של k,I, ולכן נידע גם מהי השגיאה שהתקבלה בחישוב. מכאן נחליט אם הדיוק משביע רצון, אחרת ניתן לבחור h אחר. נפתור עבור מקרה כללי:
<div style="direction: ltr;">
:<math>\ \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}</math>
</div>

==n=2==
:<math>\ 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)</math>

===קירוב ראשון: סכום פשוט===
:<math>I_2f(x)= 2hf(x)+ \left\{ 2h^2f'(c) \right\}</math>
קירוב זה משתמש בשתי נקודות. זהו קירוב גס, מכיוון שלוקחים מלבן אחד לאורך של שתי נקודות, אשר גובהו נקבע לפי הנקודה הראשונה.

===קירוב שני: סכום פשוט דרך נקודת ביניים===
קירוב זה משתמש בשלוש נקודות, כאשר הנקודה האמצעית קובעת את גובה המלבן, אשר נפרש לאורך שלושת הנקודות.
[[תמונה:Mittelpunktsregel.svg|left|100px]]
:<math>\ 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\}</math>

===קירוב שלישי: שיטת סימפסון===
[[תמונה:Simpsons method illustration.svg|left|100px]]
<div style="direction: ltr;">
:<math>\ I_2f(x)=2h\left[ 1+\Delta+ {1\over 6}\Delta^2 \right]f(x)+ \left\{ -{h^5\over 90}f^{(4)}(c) \right\}=</math>
:<math>\ 2h{1\over 6}[f(x)+ 4f(x+h)+ f(x+2h)]+ \left\{ -{1\over 90}h^5f^{(4)}(c) \right\}</math>
</div>
לפי הביטויים שקיבלנו אפשר לראות כי ניתן להשתמש בשיטה רק כאשר יש לנו מספר זוגי של מלבנים.

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

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

==ראו גם==
* [[מתמטיקה תיכונית/אלגברה תיכונית/משוואות/משוואות בשני נעלמים או יותר]]
* [[חשבון דיפרנציאלי ואינטגרלי/אינטגרלים]]
* [[חשבון אינפיניטסימלי]]

==קישורים חיצוניים==
{{מיזמים|ויקיפדיה=אינטגרל|ויקיפדיה 2=שיטות נומריות לחישוב אינטגרלים מסוימים}}

{{אנליזה נומרית|מוגבל=כן}}

[[קטגוריה:אנליזה נומרית]]