עריכה
השיטה מתבססת על פיתוח לטור טיילור של ערכי הפונקציה
y
1
,
y
2
,
.
.
.
{\displaystyle \ y_{1},y_{2},...}
x
1
,
x
2
,
.
.
.
{\displaystyle \ x_{1},x_{2},...}
לשם כך נפתח תחילה את הפונקציה לטור טיילור סביב x0 :
y
1
=
y
0
+
h
y
0
′
+
h
2
2
y
0
″
+
.
.
.
{\displaystyle \ y_{1}=y_{0}+hy_{0}'+{\tfrac {h^{2}}{2}}y_{0}''+...}
השלב הבא הוא לחשב את y2 באמצעות הערך שהתקבל עבור y1 , ובאופן כללי:
y
i
+
1
=
y
i
+
h
y
i
′
+
h
2
2
y
i
″
+
.
.
.
{\displaystyle \ y_{i+1}=y_{i}+hy_{i}'+{\tfrac {h^{2}}{2}}y_{i}''+...}
שימו לב כי רק y0 הוא ערך מדוייק, כי ערכי yi עבור i>0 התקבלו באמצעות חישוב מקורב. המטרה היא, אם כן, להקטין ככל הניתן את ההפרש
y
(
x
i
)
−
y
i
,
i
≥
1
{\displaystyle \ y(x_{i})-y_{i}\ ,\ i\geq 1}
על מנת למצוא קירוב נומרי, נשתמש, לשם הדגמה, בשלושת האיברים הראשונים של הטור. נמצא את ערכי המקדמים עבור y1 :
y
0
=
y
(
x
0
)
{\displaystyle \ y_{0}=y(x_{0})}
y
0
′
=
y
′
(
x
0
)
=
f
(
x
0
,
y
0
)
{\displaystyle \ y_{0}'=y'(x_{0})=f(x_{0},y_{0})}
y
0
″
=
(
d
d
x
y
′
)
|
x
0
=
[
d
d
x
f
(
x
,
y
(
x
)
)
]
|
x
0
=
[
∂
f
∂
x
+
∂
f
∂
y
d
y
d
x
]
|
x
0
{\displaystyle \ y_{0}''=\left.\left({\tfrac {d}{dx}}y'\right)\right|_{x_{0}}=\left.\left[{\tfrac {d}{dx}}f(x,y(x))\right]\right|_{x_{0}}=\left.\left[{\tfrac {\partial f}{\partial x}}+{\tfrac {\partial f}{\partial y}}{\tfrac {dy}{dx}}\right]\right|_{x_{0}}}
שימו לב כי ערכו של הביטוי שהתקבל באמצעות כלל השרשרת הינו ידוע ב-x0 . בדרך זו ממשיכים ומוצאים את ערכי y2 , y3 ...
ככל שרוצים לקבל דיוק רב יותר, יש צורך לחשב נגזרות מסדר גבוה יותר, ולפעמים הדבר בלתי אפשרי.
פתרון מד"ר (אוילר):
y
i
+
1
=
y
i
+
h
f
(
x
i
,
y
i
)
{\displaystyle \ y_{i+1}=y_{i}+hf(x_{i},y_{i})}
שיטת אוילר היא פתרון באמצעות טור טיילור, כאשר לוקחים שני איברים ראשונים בלבד, דבר אשר חוסך חישוב נגזרות על חשבון פגיעה בדיוק הפתרון:
y
i
+
1
=
y
i
+
h
y
i
′
+
R
2
=
y
i
+
h
f
(
x
i
,
y
i
)
+
O
(
h
2
)
{\displaystyle \ y_{i+1}=y_{i}+hy_{i}'+R_{2}=y_{i}+hf(x_{i},y_{i})+O(h^{2})}
כאמור, בשיטה זו הדיוק אינו רב אך יתרה מזאת, השיטה לעתים קרובות לא יציבה מכיוון שעבור פתרונות שאינם די מתונים, שגיאות קטנות מוגברות ככל שמתקדמים ב-x.
כפי שנראה בהמשך, שיטת אוילר היא שיטת רונגה-קוטה מסדר ראשון.
-
עריכה
אלו הן קבוצה של שיטות רבות אשר נבדלות זו מזו בכמות החישובים שיש לבצע, ולכן גם בדיוק. שיטות אלה שימושיות במיוחד, עקב נוחיותן בכך שאינן מצריכות חישוב כל נגזרת של הפונקציה המדוברת, אלא רק חישוב ערכי הפונקציה עצמה.
כפי שראינו, לפי טור טיילור:
y
i
+
1
=
y
i
+
h
f
(
x
i
,
y
i
)
+
h
2
2
[
∂
f
(
x
i
,
y
i
)
∂
x
+
∂
f
(
x
i
,
y
i
)
∂
y
f
(
x
i
,
y
i
)
]
+
O
(
h
3
)
{\displaystyle \ y_{i+1}=y_{i}+hf(x_{i},y_{i})+{\frac {h^{2}}{2}}\left[{\frac {\partial f(x_{i},y_{i})}{\partial x}}+{\frac {\partial f(x_{i},y_{i})}{\partial y}}f(x_{i},y_{i})\right]+O(h^{3})}
נשתמש בשיטה הבאה, אשר מזכירה מעט את משפט לגראנז':
y
i
+
1
=
y
i
+
λ
1
h
f
(
x
i
,
y
i
)
+
λ
2
h
f
[
x
i
+
μ
1
h
,
y
i
+
μ
2
h
f
(
x
i
,
y
i
)
]
{\displaystyle \ y_{i+1}=y_{i}+\lambda _{1}hf(x_{i},y_{i})+\lambda _{2}hf\left[x_{i}+\mu _{1}h,y_{i}+\mu _{2}hf(x_{i},y_{i})\right]}
ואת הקבועים
λ
1
,
λ
2
,
μ
1
,
μ
2
{\displaystyle \ \lambda _{1},\lambda _{2},\mu _{1},\mu _{2}}
(
x
i
,
y
i
)
{\displaystyle \ (x_{i},y_{i})}
y
i
+
1
=
y
i
+
λ
1
h
f
(
x
i
,
y
i
)
+
λ
2
h
[
f
(
x
i
,
y
i
)
+
μ
1
h
∂
f
(
x
i
,
y
i
)
∂
x
+
μ
2
h
f
(
x
i
,
y
i
)
∂
f
(
x
i
,
y
i
)
∂
y
+
O
(
h
2
)
]
{\displaystyle \ y_{i+1}=y_{i}+\lambda _{1}hf(x_{i},y_{i})+\lambda _{2}h\left[f(x_{i},y_{i})+\mu _{1}h{\frac {\partial f(x_{i},y_{i})}{\partial x}}+\mu _{2}hf(x_{i},y_{i}){\frac {\partial f(x_{i},y_{i})}{\partial y}}+O(h^{2})\right]}
⇒
h
f
(
x
i
,
y
i
)
:
λ
1
+
λ
2
=
1
h
2
∂
f
(
x
i
,
y
i
)
∂
x
:
λ
2
μ
1
=
1
2
h
2
f
(
x
i
,
y
i
)
∂
f
(
x
i
,
y
i
)
∂
y
:
λ
2
μ
2
=
1
2
{\displaystyle \ \Rightarrow \ {\begin{array}{lcl}hf(x_{i},y_{i})&:&\lambda _{1}+\lambda _{2}=1\\h^{2}{\frac {\partial f(x_{i},y_{i})}{\partial x}}&:&\lambda _{2}\mu _{1}={\tfrac {1}{2}}\\h^{2}f(x_{i},y_{i}){\frac {\partial f(x_{i},y_{i})}{\partial y}}&:&\lambda _{2}\mu _{2}={\tfrac {1}{2}}\end{array}}}
מאחר ויש לנו 4 נעלמים ורק 3 משוואות, נקבע אחד מהם שרירותית וזאת תהיה המשוואה הרביעית. עבור
λ
1
=
1
2
{\displaystyle \ \lambda _{1}={\tfrac {1}{2}}}
λ
1
=
0
{\displaystyle \ \lambda _{1}=0}
עריכה
אם נסמן
{
λ
1
=
1
−
κ
λ
2
=
κ
≠
0
μ
1
=
μ
2
=
1
2
κ
{\displaystyle {\begin{cases}\lambda _{1}=1-\kappa \\\lambda _{2}=\kappa \neq 0\\\mu _{1}=\mu _{2}={1 \over 2\kappa }\end{cases}}}
y
i
+
1
=
y
i
+
h
[
(
1
−
κ
)
f
(
x
i
,
y
i
)
+
κ
f
(
x
i
+
h
2
κ
,
y
i
+
h
2
κ
f
(
x
i
,
y
i
)
)
]
+
O
(
h
3
)
{\displaystyle \ y_{i+1}=y_{i}+h[(1-\kappa )f(x_{i},y_{i})+\kappa f(x_{i}+{\tfrac {h}{2\kappa }},y_{i}+{\tfrac {h}{2\kappa }}f(x_{i},y_{i}))]+O(h^{3})}
ואז שיטת Heun תתקבל עבור
κ
=
1
2
{\displaystyle \ \kappa ={\tfrac {1}{2}}}
κ
=
1
{\displaystyle \ \kappa =1}
עריכה
שיטה זו נקראת גם "שיטת אוילר המשופרת" (improved Euler method). מבחינה גאומטרית, היא משתמשת בממוצע השיפועים בשתי נקודות סמוכות
(
x
i
,
y
i
)
,
(
x
i
+
h
,
y
i
+
h
y
i
′
)
{\displaystyle \ (x_{i},y_{i}),\ (x_{i}+h,y_{i}+hy_{i}')}
x
i
+
1
=
x
i
+
h
,
y
i
+
1
=
y
i
+
h
y
i
′
{\displaystyle \ x_{i+1}=x_{i}+h,\ y_{i+1}=y_{i}+hy_{i}'}
(
x
i
+
1
,
y
i
+
1
)
{\displaystyle \ (x_{i+1},y_{i+1})}
(
x
i
,
y
i
)
{\displaystyle \ (x_{i},y_{i})}
x
=
x
i
+
1
=
x
+
h
{\displaystyle \ x=x_{i+1}=x+h}
נקח
λ
1
=
1
2
{\displaystyle \ \lambda _{1}={\tfrac {1}{2}}}
λ
2
=
1
2
,
μ
1
=
μ
2
=
1
{\displaystyle \ \lambda _{2}={\tfrac {1}{2}}\ ,\ \mu _{1}=\mu _{2}=1}
y
i
+
1
=
y
i
+
1
2
h
f
(
x
i
,
y
i
)
+
1
2
h
f
[
x
i
+
h
,
y
i
+
h
f
(
x
i
,
y
i
)
]
+
O
(
h
3
)
{\displaystyle \ y_{i+1}=y_{i}+{\tfrac {1}{2}}hf(x_{i},y_{i})+{\tfrac {1}{2}}hf[x_{i}+h,y_{i}+hf(x_{i},y_{i})]+O(h^{3})}
את הביטוי האחרון נהוג לבטא במקדמים
K
0
,
K
1
{\displaystyle K_{0},K_{1}}
y
i
+
1
=
y
i
+
h
2
[
f
(
x
i
,
y
i
)
⏞
K
0
+
f
(
x
i
+
h
,
y
i
+
h
f
(
x
i
,
y
i
)
⏞
K
0
)
⏞
K
1
]
+
O
(
h
3
)
=
y
i
+
h
2
(
K
0
+
K
1
)
+
O
(
h
3
)
{\displaystyle \ y_{i+1}=y_{i}+{\frac {h}{2}}\left[\overbrace {f(x_{i},y_{i})} ^{K_{0}}+\overbrace {f\left(x_{i}+h,y_{i}+h\overbrace {f(x_{i},y_{i})} ^{K_{0}}\right)} ^{K_{1}}\right]+O(h^{3})=y_{i}+{\tfrac {h}{2}}(K_{0}+K_{1})+O(h^{3})}
(כאשר בשיטות מסדר גבוה יותר המקדמים
K
i
{\displaystyle K_{i}}
עריכה
שיטה זו נקראת גם modified Euler method. בניגוד לשיטת Heun, במקום לקחת את ממוצע השיפועים של שתי נקודות, ניקח את השיפוע של נקודת האמצע.
במקרה זה נקח
λ
1
=
0
{\displaystyle \ \lambda _{1}=0}
λ
2
=
1
,
μ
1
=
μ
2
=
1
2
{\displaystyle \ \lambda _{2}=1\ ,\ \mu _{1}=\mu _{2}={\tfrac {1}{2}}}
y
i
+
1
=
y
i
+
h
f
[
x
i
+
h
2
,
y
i
+
h
2
f
(
x
i
,
y
i
)
]
+
O
(
h
3
)
{\displaystyle \ y_{i+1}=y_{i}+hf[x_{i}+{\tfrac {h}{2}},y_{i}+{\tfrac {h}{2}}f(x_{i},y_{i})]+O(h^{3})}
שני שיקולים ינחו אותנו בעת בחירת מרווח האינטגרציה h:
בהתאם לגודל האופייני של הבעיה הפיזיקלית נחליט בכמה סדרי גודל על h להיות קטן יותר.
ערכו של h יכול להשתנות מנקודה לנקודה בהתאם להתנהגות הפונקציה. אם אין גרדיאנטים חריפים מנקודה לנקדוה, ניתן לקחת מרווחים גדולים יותר בלי לפגוע משמעותית בדיוק, כך ש-
h
=
c
y
′
=
c
f
(
x
i
,
y
i
)
{\displaystyle \ h={c \over y'}={c \over f(x_{i},y_{i})}}
בעוד ששיטת אוילר מחשבת את הערך הבא באמצעות השיפוע בערך הנוכחי, שיטות רונגה-קוטה משתמשות במידע נוסף.
בשיטות רונגה קוטה מסדר שני אנו מגיעים לדיוק מסדר
O
(
h
2
)
{\displaystyle \ O(h^{2})}
f
(
x
i
)
,
∂
f
(
x
i
)
∂
x
,
∂
f
(
y
i
)
∂
y
{\displaystyle \ f(x_{i}),{\tfrac {\partial f(x_{i})}{\partial x}},{\tfrac {\partial f(y_{i})}{\partial y}}}
עריכה
בשיטה זו לוקחים
K
0
=
f
(
x
i
,
y
i
)
K
1
=
f
(
x
i
+
h
2
,
y
i
+
h
2
K
0
)
K
2
=
f
(
x
i
+
h
2
,
y
i
+
h
2
K
1
)
K
3
=
f
(
x
i
+
h
,
y
i
+
h
K
2
)
y
i
+
1
=
y
i
+
h
6
(
K
0
+
2
K
1
+
2
K
2
+
K
3
)
{\displaystyle {\begin{array}{rcl}K_{0}&=&f(x_{i},y_{i})\\K_{1}&=&f(x_{i}+{\tfrac {h}{2}},y_{i}+{\tfrac {h}{2}}K_{0})\\K_{2}&=&f(x_{i}+{\tfrac {h}{2}},y_{i}+{\tfrac {h}{2}}K_{1})\\K_{3}&=&f(x_{i}+h,y_{i}+hK_{2})\\y_{i+1}&=&y_{i}+{\tfrac {h}{6}}\left(K_{0}+2K_{1}+2K_{2}+K_{3}\right)\end{array}}}
עריכה
בשיטה זו לוקחים
K
0
=
f
(
x
i
,
y
i
)
K
1
=
f
(
x
i
+
1
5
h
,
y
i
+
1
5
h
K
0
)
K
2
=
f
(
x
i
+
3
10
h
,
y
i
+
3
40
h
K
0
+
9
40
h
K
1
)
K
3
=
f
(
x
i
+
4
5
h
,
y
i
+
44
45
h
K
0
−
56
15
h
K
1
+
32
9
h
K
2
)
K
4
=
f
(
x
i
+
8
9
h
,
y
i
+
19372
6561
h
K
0
−
25360
2187
h
K
1
+
64448
6561
h
K
2
−
212
729
h
K
3
)
K
5
=
f
(
x
i
+
h
,
y
i
+
9017
3168
h
K
0
−
355
33
h
K
1
+
46732
5247
h
K
2
+
49
176
h
K
3
−
5103
18656
h
K
4
)
K
6
=
f
(
x
i
+
h
,
y
i
+
35
384
h
K
0
+
500
1113
h
K
2
+
125
192
h
K
3
−
2187
6784
h
K
4
+
11
84
h
K
5
)
y
i
+
1
=
y
i
+
h
(
5179
57600
K
0
+
7571
16695
K
2
+
393
640
K
3
−
92097
339200
K
4
+
187
2100
K
5
+
1
40
K
6
)
+
O
(
h
5
)
y
i
+
1
=
y
i
+
h
(
35
384
K
0
+
500
1113
K
2
+
125
192
K
3
−
2187
6784
K
4
+
11
84
K
5
)
+
O
(
h
6
)
{\displaystyle {\begin{array}{rcl}K_{0}&=&f(x_{i},y_{i})\\K_{1}&=&f(x_{i}+{\tfrac {1}{5}}h,y_{i}+{\tfrac {1}{5}}hK_{0})\\K_{2}&=&f(x_{i}+{\tfrac {3}{10}}h,y_{i}+{\tfrac {3}{40}}hK_{0}+{\tfrac {9}{40}}hK_{1})\\K_{3}&=&f(x_{i}+{\tfrac {4}{5}}h,y_{i}+{\tfrac {44}{45}}hK_{0}-{\tfrac {56}{15}}hK_{1}+{\tfrac {32}{9}}hK_{2})\\K_{4}&=&f(x_{i}+{\tfrac {8}{9}}h,y_{i}+{\tfrac {19372}{6561}}hK_{0}-{\tfrac {25360}{2187}}hK_{1}+{\tfrac {64448}{6561}}hK_{2}-{\tfrac {212}{729}}hK_{3})\\K_{5}&=&f(x_{i}+h,y_{i}+{\tfrac {9017}{3168}}hK_{0}-{\tfrac {355}{33}}hK_{1}+{\tfrac {46732}{5247}}hK_{2}+{\tfrac {49}{176}}hK_{3}-{\tfrac {5103}{18656}}hK_{4})\\K_{6}&=&f(x_{i}+h,y_{i}+{\tfrac {35}{384}}hK_{0}+{\tfrac {500}{1113}}hK_{2}+{\tfrac {125}{192}}hK_{3}-{\tfrac {2187}{6784}}hK_{4}+{\tfrac {11}{84}}hK_{5})\\&&\\y_{i+1}&=&y_{i}+h\left({\tfrac {5179}{57600}}K_{0}+{\tfrac {7571}{16695}}K_{2}+{\tfrac {393}{640}}K_{3}-{\tfrac {92097}{339200}}K_{4}+{\tfrac {187}{2100}}K_{5}+{\tfrac {1}{40}}K_{6}\right)+O(h^{5})\\y_{i+1}&=&y_{i}+h\left({\tfrac {35}{384}}K_{0}+{\tfrac {500}{1113}}K_{2}+{\tfrac {125}{192}}K_{3}-{\tfrac {2187}{6784}}K_{4}+{\tfrac {11}{84}}K_{5}\right)+O(h^{6})\end{array}}}
שיטות אלה, בניגוד לנ"ל, משתמשות במספר נקודות קודמות. נפתח את בסיס השיטה:
y
′
(
x
)
=
f
(
x
,
y
)
⇒
∫
x
i
x
i
+
1
d
y
d
x
d
x
=
∫
x
i
x
i
+
1
f
d
x
⇒
y
i
+
1
−
y
i
=
∫
x
i
x
i
+
1
f
d
x
{\displaystyle \ y'(x)=f(x,y)\quad \Rightarrow \ \int \limits _{x_{i}}^{x_{i+1}}{\frac {dy}{dx}}dx=\int \limits _{x_{i}}^{x_{i+1}}fdx\quad \Rightarrow \ y_{i+1}-y_{i}=\int \limits _{x_{i}}^{x_{i+1}}fdx}
נשאר, אם כן, לפתור את האינטגרל
∫
x
i
x
i
+
1
f
d
x
{\displaystyle \ \int \limits _{x_{i}}^{x_{i+1}}fdx}
x
=
x
i
+
θ
h
,
d
x
=
h
d
θ
,
0
≤
θ
≤
1
{\displaystyle \ x=x_{i}+\theta h,\ dx=hd\theta ,\ 0\leq \theta \leq 1}
i ,xi+1 ניתנים לתיאור על ידי:
f
=
E
θ
f
(
x
i
)
{\displaystyle \ f=E^{\theta }f(x_{i})}
∫
0
1
E
θ
f
(
x
i
,
y
i
)
h
⏞
c
o
n
s
t
.
d
θ
=
∫
0
1
(
1
−
∇
)
−
θ
f
(
x
i
,
y
i
)
h
d
θ
=
{\displaystyle \ \int \limits _{0}^{1}E^{\theta }\overbrace {f(x_{i},y_{i})h} ^{const.}d\theta =\int \limits _{0}^{1}(1-\nabla )^{-\theta }f(x_{i},y_{i})hd\theta =}
=
∫
0
1
h
[
1
+
θ
∇
+
θ
(
θ
+
1
)
2
∇
2
+
.
.
.
]
f
(
x
i
,
y
i
)
d
θ
=
{\displaystyle \ =\ \int \limits _{0}^{1}h\left[1+\theta \nabla +{\tfrac {\theta (\theta +1)}{2}}\nabla ^{2}+...\right]f(x_{i},y_{i})d\theta =}
=
h
[
θ
+
θ
2
2
∇
+
(
θ
3
6
+
θ
2
4
)
∇
2
.
.
.
]
f
(
x
i
,
y
i
)
|
0
1
=
{\displaystyle \ =h\left.\left[\theta +{\theta ^{2} \over 2}\nabla +\left({\theta ^{3} \over 6}+{\theta ^{2} \over 4}\right)\nabla ^{2}...\right]f(x_{i},y_{i})\right|_{0}^{1}=}
=
h
[
1
+
∇
2
+
5
12
∇
2
+
9
24
∇
3
+
.
.
.
]
f
(
x
i
,
y
i
)
{\displaystyle \ =\ h\left[1+{\nabla \over 2}+{5 \over 12}\nabla ^{2}+{9 \over 24}\nabla ^{3}+...\right]f(x_{i},y_{i})}
כך שככל שנקח יותר איברים, נקבל ביטוי התלוי ביותר נקודות קודמות. שימו לב כי השתמשנו באופרטור הפרשים אחוריים על מנת לקבל תלות בנקודות הקודמות.
שיטה רב-צעדית:
y
i
+
1
=
y
i
+
h
[
1
+
∇
2
+
5
12
∇
2
+
3
8
∇
3
+
.
.
.
]
f
(
x
i
,
y
i
)
{\displaystyle \ y_{i+1}=y_{i}+h\left[1+{\nabla \over 2}+{5 \over 12}\nabla ^{2}+{3 \over 8}\nabla ^{3}+...\right]f(x_{i},y_{i})}
שימו לב כי אם ניקח איבר בודד, נקבל את שיטת אוילר.
הערה : אחד המאפיינים של שיטות רב-צעדיות הוא שאינן מאתחלות את עצמן (not self-starting), כלומר יש לספק עבורן מספר תנאי התחלה באמצעות שיטות חד-צעדיות בעלות אותו מרווח h.
עריכה
נפתח את 4 האיברים הראשונים שציינו לעיל:
y
i
+
1
=
y
i
+
h
24
[
55
f
(
x
i
,
y
i
)
−
59
f
(
x
i
−
1
,
y
i
−
1
)
+
+
37
f
(
x
i
−
2
,
y
i
−
2
)
−
9
f
(
x
i
−
3
,
y
i
−
3
)
]
+
O
(
h
5
)
,
i
≥
3
{\displaystyle \ {\begin{aligned}y_{i+1}&=y_{i}+{h \over 24}\left[55f(x_{i},y_{i})-59f(x_{i-1},y_{i-1})+\right.\\&+\left.37f(x_{i-2},y_{i-2})-9f(x_{i-3},y_{i-3})\right]+O(h^{5})\\\end{aligned}}\ ,\quad i\geq 3}
שימו לב כי לשיטה זו יש לספק 3 תנאי התחלה, אשר בדרך כלל מוצאים אותם בשיטת רונגה-קוטה או בשיטת אוילר. כמו כן, על אותה שיטה להשתמש באותו מרווח h.
נשווה שיטה זו לשיטת רונגה-קוטה מסדר 4, מבחינת מספר פעולות החישוב:
R
−
K
A
−
B
=
4
N
4
⋅
3
+
(
N
−
3
)
=
4
N
N
+
9
≈
4
{\displaystyle \ {\tfrac {R-K}{A-B}}={\tfrac {4N}{4\cdot 3+(N-3)}}={\tfrac {4N}{N+9}}\approx 4}
על מנת להשתמש בשיטה זו עבור מערכת מד"ר , נחליף את הכתיב הסקלרי בכתיב הוקטורי:
y
→
i
+
1
=
y
→
i
+
h
24
[
55
f
→
(
x
i
,
y
→
i
)
−
59
f
→
(
x
i
−
1
,
y
→
i
−
1
)
+
+
37
f
→
(
x
i
−
2
,
y
→
i
−
2
)
−
9
f
→
(
x
i
−
3
,
y
→
i
−
3
)
]
+
O
(
h
5
)
,
i
≥
3
{\displaystyle \ {\begin{aligned}{\vec {y}}_{i+1}&={\vec {y}}_{i}+{h \over 24}\left[55{\vec {f}}(x_{i},{\vec {y}}_{i})-59{\vec {f}}(x_{i-1},{\vec {y}}_{i-1})+\right.\\&+\left.37{\vec {f}}(x_{i-2},{\vec {y}}_{i-2})-9{\vec {f}}(x_{i-3},{\vec {y}}_{i-3})\right]+O(h^{5})\ \\\end{aligned}},\quad i\geq 3}
שיטת אדאמס היא מהצורה:
y
n
+
1
=
y
n
+
h
12
[
5
y
n
+
1
′
+
8
y
n
′
−
y
n
−
1
′
]
{\displaystyle \ y_{n+1}=y_{n}+{h \over 12}[5y_{n+1}'+8y_{n}'-y_{n-1}']}
נמצא את הביטוי לשגיאה באמצעות פיתוח האגפים לטור טיילור והשוואת מקדמים:
עריכה
שיטה זו היא עבור מד"ר מסדר 2, כאשר f אינה פונקציה של
y
′
{\displaystyle \ y'}
y
″
{\displaystyle \ y''}
{
y
i
+
1
=
2
y
i
−
y
i
−
1
+
h
2
12
[
f
(
x
i
+
1
,
y
i
+
1
)
+
10
f
(
x
i
,
y
i
)
+
f
(
x
i
−
1
,
y
i
−
1
)
]
y
″
(
x
)
=
f
(
x
,
y
)
{\displaystyle \ {\begin{cases}y_{i+1}=2y_{i}-y_{i-1}+{h^{2} \over 12}[f(x_{i+1},y_{i+1})+10f(x_{i},y_{i})+f(x_{i-1},y_{i-1})]\\y''(x)=f(x,y)\end{cases}}}
שימו לב כי זוהי שיטה סתומה מאחר ו-
y
i
+
1
{\displaystyle \ y_{i+1}}
y
″
{\displaystyle \ y''}
y
i
+
1
=
y
i
−
1
+
2
h
f
(
x
i
,
y
i
)
+
O
(
h
3
)
{\displaystyle \ y_{i+1}=y_{i-1}+2hf(x_{i},y_{i})+O(h^{3})}
y
i
+
1
=
y
i
−
3
+
4
h
3
[
2
f
(
x
i
,
y
i
)
−
f
(
x
i
−
1
,
y
i
−
1
)
+
2
f
(
x
i
−
2
,
y
i
−
2
)
]
+
O
(
h
5
)
{\displaystyle \ y_{i+1}=y_{i-3}+{\tfrac {4h}{3}}[2f(x_{i},y_{i})-f(x_{i-1},y_{i-1})+2f(x_{i-2},y_{i-2})]+O(h^{5})}
עריכה
שיטות מעריך-מתקן (Predictor-Corrector Methods) מתבססות על פתרון האינטגרל על
y
′
(
x
)
=
f
(
x
,
y
)
d
x
{\displaystyle \ y'(x)=f(x,y)dx}
y
i
+
1
=
y
i
+
∫
x
i
x
i
+
1
f
d
x
{\displaystyle \ y_{i+1}=y_{i}+\int \limits _{x_{i}}^{x_{i+1}}fdx}
h
=
x
i
+
1
−
x
i
{\displaystyle \ h=x_{i+1}-x_{i}}
אם נשתמש בשיטת המלבן להערכת האינטגרל, נקבל:
y
i
+
1
=
y
i
+
∫
x
i
x
i
+
1
f
d
x
=
y
i
+
h
f
(
x
i
+
1
+
x
i
2
,
y
i
)
{\displaystyle \ y_{i+1}=y_{i}+\int \limits _{x_{i}}^{x_{i+1}}fdx=y_{i}+hf\left({\tfrac {x_{i+1}+x_{i}}{2}},y_{i}\right)}
שימו לב כי כל הפרמטרים ידועים, פרט ל-yi+1 , שאותו אנו מחפשים. אם, לעומת זאת, נשתמש בשיטת הטרפז, נקבל:
y
i
+
1
=
y
i
+
∫
x
i
x
i
+
1
f
d
x
=
y
i
+
h
2
[
f
(
x
i
+
1
,
y
i
+
1
)
+
f
(
x
i
,
y
i
)
]
+
O
(
h
3
)
{\displaystyle \ y_{i+1}=y_{i}+\int \limits _{x_{i}}^{x_{i+1}}fdx=y_{i}+{\tfrac {h}{2}}[f(x_{i+1},y_{i+1})+f(x_{i},y_{i})]+O(h^{3})}
שימו לב כי קיבלנו קשר סתום עבור yi+1 , וכביכול אין באפשרותנו להמשיך. כאן, ובכלל, במקרים בהם יש צורך להשתמש בערכים שאינם קיימים עדיין לרשותנו, נכנסת לתמונה שיטת המעריך-מתקן.
כללית, כל שיטת מעריך-מתקן מורכבת משני שלבים:
מעריך : בוחרים בשיטה כלשהי אשר מאפשרת חישוב ישיר של yi+1 .מתקן איטרטיבי : מציבים את המעריך לתוך השיטה הסתומה.אם נבחר בשיטה כלשהי על מנת להעריך את ערכו של yi+1 , כמו למשל בשיטת המלבן, נוכל להציב איטרטיבית את התוצאה בשיטת הטרפז. שיטת הטרפז תניב ערך מתוקן אשר קרוב יותר לאמיתי. ניתן להפעיל שוב את השלב האיטרטיבי, אולם לא נהוג להפעילה יותר מפעמיים. אם נדרש דיוק רב יותר, ניתן להקטין את המרווח h, ואף ניתן לעשות זאת בכל איטרציה, אך במקרה זה יש לבצע מחדש את שלב המעריך.
ניתן לראות בביטוי הסתום עבור yi+1 מעין שיטה איטרטיבית (נסמן
y
^
=
y
i
+
1
{\displaystyle \ {\hat {y}}=y_{i+1}}
y
^
n
+
1
=
Φ
(
y
^
n
)
{\displaystyle \ {\hat {y}}_{n+1}=\Phi ({\hat {y}}_{n})}
|
Φ
′
(
y
^
)
|
<
1
{\displaystyle \ |\Phi '({\hat {y}})|<1}
|
Φ
′
(
y
^
)
|
=
|
h
2
∂
f
∂
y
|
<
1
⇒
h
<
2
∂
f
∂
y
^
{\displaystyle \ |\Phi '({\hat {y}})|=\left|{h \over 2}{\partial f \over \partial y}\right|<1\ \Rightarrow \ h<{\frac {2}{\tfrac {\partial f}{\partial {\hat {y}}}}}}
ואז ההתכנסות מובטחת. על מנת להבטיח התכנסות מהירה, יש לדאוג לכך ש-
h
≪
2
∂
f
∂
y
^
{\displaystyle \ h\ll {\frac {2}{\tfrac {\partial f}{\partial {\hat {y}}}}}}
Adams-Bashforth :
שיטה זו היא למעשה קירוב של פולינום לגראנז' המתקבל על ידי ארבע נקודות. ניתן לאמר שהמעריך הוא:
y
i
+
1
p
r
e
d
i
c
t
o
r
=
y
i
+
h
24
[
55
f
(
x
i
,
y
i
)
−
59
f
(
x
i
−
1
,
y
i
−
1
)
+
37
f
(
x
i
−
2
,
y
i
−
2
)
−
9
f
(
x
i
−
3
,
y
i
−
3
)
]
{\displaystyle \ y_{i+1}^{predictor}=y_{i}+{h \over 24}\left[55f(x_{i},y_{i})-59f(x_{i-1},y_{i-1})+37f(x_{i-2},y_{i-2})-9f(x_{i-3},y_{i-3})\right]}
והמתקן יתקבל באמצעות הנקודה החדשה ושלוש הנקודות הקודמות לה:
y
i
+
1
c
o
r
r
e
c
t
o
r
=
y
i
+
h
24
[
9
f
(
x
i
+
1
,
y
i
+
1
)
+
19
f
(
x
i
,
y
i
)
−
5
f
(
x
i
−
1
,
y
i
−
1
)
+
(
x
i
−
2
,
y
i
−
2
)
]
{\displaystyle \ y_{i+1}^{corrector}=y_{i}+{h \over 24}\left[9f(x_{i+1},y_{i+1})+19f(x_{i},y_{i})-5f(x_{i-1},y_{i-1})+(x_{i-2},y_{i-2})\right]}
עריכה
![{\displaystyle \ y_{0}''=\left.\left({\tfrac {d}{dx}}y'\right)\right|_{x_{0}}=\left.\left[{\tfrac {d}{dx}}f(x,y(x))\right]\right|_{x_{0}}=\left.\left[{\tfrac {\partial f}{\partial x}}+{\tfrac {\partial f}{\partial y}}{\tfrac {dy}{dx}}\right]\right|_{x_{0}}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/955251e57c3e4d6d90c38f6b53dcc07c3f305d3a)
![{\displaystyle \ y_{i+1}=y_{i}+hf(x_{i},y_{i})+{\frac {h^{2}}{2}}\left[{\frac {\partial f(x_{i},y_{i})}{\partial x}}+{\frac {\partial f(x_{i},y_{i})}{\partial y}}f(x_{i},y_{i})\right]+O(h^{3})}](https://wikimedia.org/api/rest_v1/media/math/render/svg/420c206da6fd7568f44240534c156b9a7053f076)
![{\displaystyle \ y_{i+1}=y_{i}+\lambda _{1}hf(x_{i},y_{i})+\lambda _{2}hf\left[x_{i}+\mu _{1}h,y_{i}+\mu _{2}hf(x_{i},y_{i})\right]}](https://wikimedia.org/api/rest_v1/media/math/render/svg/b1599faa671928d1c3ad6bcedfe647e6de7d3893)
![{\displaystyle \ y_{i+1}=y_{i}+\lambda _{1}hf(x_{i},y_{i})+\lambda _{2}h\left[f(x_{i},y_{i})+\mu _{1}h{\frac {\partial f(x_{i},y_{i})}{\partial x}}+\mu _{2}hf(x_{i},y_{i}){\frac {\partial f(x_{i},y_{i})}{\partial y}}+O(h^{2})\right]}](https://wikimedia.org/api/rest_v1/media/math/render/svg/71525f8aba30603c3b4734c3df27bf0a4d7b8743)
![{\displaystyle \ y_{i+1}=y_{i}+h[(1-\kappa )f(x_{i},y_{i})+\kappa f(x_{i}+{\tfrac {h}{2\kappa }},y_{i}+{\tfrac {h}{2\kappa }}f(x_{i},y_{i}))]+O(h^{3})}](https://wikimedia.org/api/rest_v1/media/math/render/svg/fef754cfb8be525fe9a75bbdbeedfa8198bd6705)
![{\displaystyle \ y_{i+1}=y_{i}+{\tfrac {1}{2}}hf(x_{i},y_{i})+{\tfrac {1}{2}}hf[x_{i}+h,y_{i}+hf(x_{i},y_{i})]+O(h^{3})}](https://wikimedia.org/api/rest_v1/media/math/render/svg/2e7220bc01d3d492bea62ee8b92dfdfd7f0b0e37)
![{\displaystyle \ y_{i+1}=y_{i}+{\frac {h}{2}}\left[\overbrace {f(x_{i},y_{i})} ^{K_{0}}+\overbrace {f\left(x_{i}+h,y_{i}+h\overbrace {f(x_{i},y_{i})} ^{K_{0}}\right)} ^{K_{1}}\right]+O(h^{3})=y_{i}+{\tfrac {h}{2}}(K_{0}+K_{1})+O(h^{3})}](https://wikimedia.org/api/rest_v1/media/math/render/svg/6d29710a8e2b8926453b604f46d9c04eaaf0a2cf)
![{\displaystyle \ y_{i+1}=y_{i}+hf[x_{i}+{\tfrac {h}{2}},y_{i}+{\tfrac {h}{2}}f(x_{i},y_{i})]+O(h^{3})}](https://wikimedia.org/api/rest_v1/media/math/render/svg/8354b45aa9cb4e753e9baab52cf5fd7ad831995a)
![{\displaystyle \ =\ \int \limits _{0}^{1}h\left[1+\theta \nabla +{\tfrac {\theta (\theta +1)}{2}}\nabla ^{2}+...\right]f(x_{i},y_{i})d\theta =}](https://wikimedia.org/api/rest_v1/media/math/render/svg/c667d822e0f89472e465539f295552f942bbfe18)
![{\displaystyle \ =h\left.\left[\theta +{\theta ^{2} \over 2}\nabla +\left({\theta ^{3} \over 6}+{\theta ^{2} \over 4}\right)\nabla ^{2}...\right]f(x_{i},y_{i})\right|_{0}^{1}=}](https://wikimedia.org/api/rest_v1/media/math/render/svg/6d94561b066c5826521b393e0e963a2e32721cb0)
![{\displaystyle \ =\ h\left[1+{\nabla \over 2}+{5 \over 12}\nabla ^{2}+{9 \over 24}\nabla ^{3}+...\right]f(x_{i},y_{i})}](https://wikimedia.org/api/rest_v1/media/math/render/svg/f6256f5fe8c90237175a51391648e0b7ec896390)
![{\displaystyle \ y_{i+1}=y_{i}+h\left[1+{\nabla \over 2}+{5 \over 12}\nabla ^{2}+{3 \over 8}\nabla ^{3}+...\right]f(x_{i},y_{i})}](https://wikimedia.org/api/rest_v1/media/math/render/svg/ee8a5efe9da87cf447172d0c4cdce3aa7ffb9d85)
![{\displaystyle \ {\begin{aligned}y_{i+1}&=y_{i}+{h \over 24}\left[55f(x_{i},y_{i})-59f(x_{i-1},y_{i-1})+\right.\\&+\left.37f(x_{i-2},y_{i-2})-9f(x_{i-3},y_{i-3})\right]+O(h^{5})\\\end{aligned}}\ ,\quad i\geq 3}](https://wikimedia.org/api/rest_v1/media/math/render/svg/c63115552f40d0ab8def2dd31da319b430ea8a70)
![{\displaystyle \ {\begin{aligned}{\vec {y}}_{i+1}&={\vec {y}}_{i}+{h \over 24}\left[55{\vec {f}}(x_{i},{\vec {y}}_{i})-59{\vec {f}}(x_{i-1},{\vec {y}}_{i-1})+\right.\\&+\left.37{\vec {f}}(x_{i-2},{\vec {y}}_{i-2})-9{\vec {f}}(x_{i-3},{\vec {y}}_{i-3})\right]+O(h^{5})\ \\\end{aligned}},\quad i\geq 3}](https://wikimedia.org/api/rest_v1/media/math/render/svg/b3a4d862e8330d1df53857723235a6a5f9ba6b25)
![{\displaystyle \ y_{n+1}=y_{n}+{h \over 12}[5y_{n+1}'+8y_{n}'-y_{n-1}']}](https://wikimedia.org/api/rest_v1/media/math/render/svg/aaf18afa1bf78d6a014f70c17433d84e1348253c)
![{\displaystyle \ {\begin{aligned}{\mbox{RHS:}}&\quad y_{n}+{h \over 12}\left[\left(5y_{n}'+5hy_{n}''+5{\tfrac {h^{2}}{2}}y_{n}'''+5{\tfrac {h^{3}}{3!}}y_{n}^{(4)}+...\right)+\right.\\&\left.+(8y_{n}')-\left(y_{n}'-hy_{n}''+{\tfrac {h^{2}}{2}}y_{n}'''-{\tfrac {h^{3}}{3!}}y_{n}^{(4)}+...\right)\right]+R\\\end{aligned}}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/9739c32f1224a926a88e4d0f89e8aa06b48d798d)
![{\displaystyle \ {\begin{cases}y_{i+1}=2y_{i}-y_{i-1}+{h^{2} \over 12}[f(x_{i+1},y_{i+1})+10f(x_{i},y_{i})+f(x_{i-1},y_{i-1})]\\y''(x)=f(x,y)\end{cases}}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/98b954fea8d4f2047e22b3208b179380d8dba095)
![{\displaystyle \ y_{i+1}=y_{i-3}+{\tfrac {4h}{3}}[2f(x_{i},y_{i})-f(x_{i-1},y_{i-1})+2f(x_{i-2},y_{i-2})]+O(h^{5})}](https://wikimedia.org/api/rest_v1/media/math/render/svg/f24616b6370c220f95dc279bd78fd2a8395771cb)
![{\displaystyle \ y_{i+1}=y_{i}+\int \limits _{x_{i}}^{x_{i+1}}fdx=y_{i}+{\tfrac {h}{2}}[f(x_{i+1},y_{i+1})+f(x_{i},y_{i})]+O(h^{3})}](https://wikimedia.org/api/rest_v1/media/math/render/svg/4cc0b8c39bba817f7f83de2cf2d2e7169f2e5003)
![{\displaystyle \ y_{i+1}^{predictor}=y_{i}+{h \over 24}\left[55f(x_{i},y_{i})-59f(x_{i-1},y_{i-1})+37f(x_{i-2},y_{i-2})-9f(x_{i-3},y_{i-3})\right]}](https://wikimedia.org/api/rest_v1/media/math/render/svg/791f2a8d0a5f05ed07b0e7e1b46af8b9c3b16a0d)
![{\displaystyle \ y_{i+1}^{corrector}=y_{i}+{h \over 24}\left[9f(x_{i+1},y_{i+1})+19f(x_{i},y_{i})-5f(x_{i-1},y_{i-1})+(x_{i-2},y_{i-2})\right]}](https://wikimedia.org/api/rest_v1/media/math/render/svg/dbd80b003d7a40a63758d0fa3e3c58429997074e)
