Применяя формулу численного дифференцирования
править
y
(
t
i
+
1
)
−
y
(
t
i
)
h
=
f
(
t
i
,
y
(
t
i
)
)
⇒
y
i
+
1
=
y
i
+
h
⋅
f
(
t
i
,
y
i
)
{\displaystyle {\frac {y(t_{i+1})-y(t_{i})}{h}}=f(t_{i},y(t_{i}))\quad \Rightarrow \quad y_{i+1}=y_{i}+h\cdot f(t_{i},y_{i})}
{
y
′
=
f
(
t
,
y
(
t
)
)
,
t
i
=
t
0
+
i
h
,
i
=
1
,
n
¯
y
(
t
0
)
=
y
0
{\displaystyle {\begin{cases}y'=f(t,y(t)),&\quad t_{i}=t_{0}+ih,i={\overline {1,n}}\\y(t_{0})=y_{0}\end{cases}}}
Дискретизация задачи (переход от непрерывной функции к дискретной)
y
(
t
i
)
{\displaystyle y(t_{i})}
- точное решение.
y
i
{\displaystyle y_{i}}
- приближенное.
y
′
(
t
i
)
≈
h
−
1
(
y
i
+
1
−
y
i
)
{\displaystyle y'(t_{i})\approx h^{-1}(y_{i+1}-y_{i})}
,
y
i
+
1
−
y
i
h
=
f
(
t
i
,
y
i
)
{\displaystyle {\frac {y_{i+1}-y_{i}}{h}}=f(t_{i},y_{i})}
y
n
=
(
y
0
⋮
y
N
)
{\displaystyle y^{n}={\begin{pmatrix}y_{0}\\\vdots \\y_{N}\end{pmatrix}}}
- в ответе получаем вектор.
{
y
i
+
1
=
y
i
+
h
f
(
t
i
,
y
i
)
y
0
- задано
{\displaystyle {\begin{cases}y_{i+1}=y_{i}+hf(t_{i},y_{i})&\\y_{0}{\text{ - задано}}\end{cases}}\qquad }
шаг обычно берут h = 0.1}
Геометрическая интерпретация.
Геометрическая интерпретация метода Эйлера
l
(
t
)
=
y
(
t
i
)
+
y
′
(
t
i
)
(
t
−
t
i
)
{\displaystyle l(t)=y(t_{i})+y'(t_{i})(t-t_{i})}
- уравнение касательной.
Погрешности:
Определение. Погрешность на шаге определяется по формуле
r
i
=
y
(
t
i
+
1
)
−
y
i
+
1
{\displaystyle r_{i}=y(t_{i+1})-y_{i+1}}
, где
y
(
t
i
)
{\displaystyle y(t_{i})}
- точное решение, а
y
i
{\displaystyle y_{i}}
- приближенное.
Определение. Локальная погрешность:
l
i
=
y
(
t
i
)
−
y
i
{\displaystyle l_{i}=y(t_{i})-y_{i}}
, если
y
i
−
1
≡
y
(
t
i
−
1
)
{\displaystyle y_{i-1}\equiv y(t_{i-1})}
- точное.
Определение. Глобальная погрешность:
E
(
h
)
=
max
0
≤
i
≤
N
−
1
‖
r
i
‖
{\displaystyle E(h)=\max {0\leq i\leq N-1}\|r_{i}\|}
Определение. Будем называть метод
p
{\displaystyle p}
-го порядка точности по
h
{\displaystyle h}
, если
E
(
h
)
≤
c
h
p
{\displaystyle E(h)\leq ch^{p}}
, где
c
{\displaystyle c}
- константа.
y
i
+
1
=
y
i
+
h
k
i
{\displaystyle y_{i+1}=y_{i}+hk_{i}}
. изменяя
k
i
{\displaystyle k_{i}}
можно легко получить меньшую погрешность.
k
i
=
f
(
t
i
+
1
2
,
y
i
+
1
2
{\displaystyle k_{i}=f(t_{i+{\frac {1}{2}}},y_{i+{\frac {1}{2}}}}
) - усовершенствованный метод Эйлера.
{
y
i
+
1
2
=
y
i
+
h
2
f
(
t
i
,
y
i
)
y
i
+
1
=
y
i
+
h
f
(
t
i
+
1
2
,
y
i
+
1
2
)
⇔
y
i
+
1
=
y
i
+
h
k
i
,
k
i
=
f
(
t
i
+
1
2
,
y
i
+
h
2
f
(
t
i
,
y
i
)
)
{\displaystyle {\begin{cases}y_{i+{\frac {1}{2}}}=y_{i}+{\frac {h}{2}}f(t_{i},y_{i})\\y_{i+1}=y_{i}+hf(t_{i+{\frac {1}{2}}},y_{i+{\frac {1}{2}}})\end{cases}}\Leftrightarrow y_{i+1}=y_{i}+hk_{i},k_{i}=f(t_{i+{\frac {1}{2}}},y_{i}+{\frac {h}{2}}f(t_{i},y_{i}))}
(за угол наклона берется наклон касательной в средней точке).
E
(
h
)
=
c
h
2
{\displaystyle E(h)=ch^{2}}
y
i
+
1
+
h
k
i
,
k
i
=
k
i
1
+
k
i
2
2
,
k
i
1
=
f
(
t
i
,
y
i
)
,
k
i
2
=
f
(
t
i
+
h
,
y
i
+
h
k
i
2
)
{\displaystyle y_{i+1}+hk_{i},\quad k_{i}={\frac {k_{i}^{1}+k_{i}^{2}}{2}},\quad k_{i}^{1}=f(t_{i},y_{i}),k_{i}^{2}=f(t_{i}+h,y_{i}+hk_{i}^{2})}
.
E
(
h
)
=
c
h
2
{\displaystyle E(h)=ch^{2}}
.
Метод Эйлера - это метод Рунге-Кутты I-го порядка точности.
y
i
+
1
=
y
i
+
h
k
i
,
k
i
{\displaystyle y_{i+1}=y_{i}+hk_{i},\quad k_{i}}
- угловой коэффициент секущей.
m
{\displaystyle m}
-этапный метод Рунге-Кутты
y
i
=
y
(
t
i
)
{\displaystyle y_{i}=y(t_{i})}
Метод Рунге-Кутты
t
i
(
1
)
=
t
i
{\displaystyle t_{i}^{(1)}=t_{i}}
t
i
(
2
)
=
t
i
(
1
)
+
α
2
h
0
<
α
2
<
1
t
i
(
m
)
=
t
i
(
1
)
+
α
m
h
0
<
α
m
<
1
{\displaystyle t_{i}^{(2)}=t_{i}^{(1)}+\alpha _{2}h\quad 0<\alpha _{2}<1\quad t_{i}^{(m)}=t_{i}^{(1)}+\alpha _{m}h\quad 0<\alpha _{m}<1}
k
i
(
1
)
=
f
(
t
i
(
1
)
,
y
i
(
1
)
)
=
f
(
t
i
,
y
i
)
{\displaystyle k_{i}^{(1)}=f(t_{i}^{(1)},y_{i}^{(1)})=f(t_{i},y_{i})}
k
i
(
2
)
=
f
(
t
i
(
1
)
+
α
2
h
,
y
i
+
h
β
21
k
i
(
1
)
)
{\displaystyle k_{i}^{(2)}=f(t_{i}^{(1)}+\alpha _{2}h,y_{i}+h\beta _{21}k_{i}^{(1)})}
k
i
(
3
)
=
f
(
t
i
(
1
)
+
α
3
h
,
y
i
+
h
(
β
31
k
i
(
1
)
+
β
32
k
i
(
2
)
)
)
{\displaystyle k_{i}^{(3)}=f(t_{i}^{(1)}+\alpha _{3}h,y_{i}+h(\beta _{31}k_{i}^{(1)}+\beta _{32}k_{i}^{(2)}))}
…
{\displaystyle \dots }
k
i
(
m
)
=
f
(
t
i
(
1
)
+
α
m
h
,
y
i
+
h
(
β
m
1
k
i
(
1
)
+
.
.
.
+
β
m
,
m
−
1
k
i
(
m
−
1
)
)
)
{\displaystyle k_{i}^{(m)}=f(t_{i}^{(1)}+\alpha _{m}h,y_{i}+h(\beta _{m1}k_{i}^{(1)}+...+\beta _{m,m-1}k_{i}^{(m-1)}))}
Все коэффициенты (
α
{\displaystyle \alpha }
и
β
{\displaystyle \beta }
) подлежат определению:
k
i
=
c
1
k
i
(
1
)
+
c
2
k
i
(
2
)
+
.
.
.
+
c
m
k
i
(
m
)
{\displaystyle k_{i}=c_{1}k_{i}^{(1)}+c_{2}k_{i}^{(2)}+...+c_{m}k_{i}^{(m)}}
коэффициенты
c
i
{\displaystyle c_{i}}
также неизвестны.
y
(
t
+
h
)
=
y
(
t
)
+
h
y
′
(
t
)
+
h
2
2
y
″
(
t
)
+
.
.
.
+
h
p
p
!
y
(
p
)
(
t
)
+
O
¯
¯
(
h
p
+
1
)
{\displaystyle y(t+h)=y(t)+hy'(t)+{\frac {h^{2}}{2}}y''(t)+...+{\frac {h^{p}}{p!}}y^{(p)}(t)+{\overline {\overline {O}}}(h^{p+1})}
Если хотим построить метод
p
{\displaystyle p}
-го порядка точности по
h
{\displaystyle h}
,
y
i
+
1
=
y
i
+
h
k
i
{\displaystyle y_{i+1}=y_{i}+hk_{i}}
. Возьмем разложение для точного решения. Правую часть раскладываем в
t
i
{\displaystyle t_{i}}
. Коэффициенты
α
,
β
{\displaystyle \alpha ,\beta }
и
c
{\displaystyle c}
подбираем так, чтобы разложение совпало до
p
{\displaystyle p}
-го порядка включительно. Построим параметрическое семейство двухэтапных методов Рунге-Кутты (2-го порядка точности по
h
{\displaystyle h}
)
t
i
(
1
)
=
t
i
,
t
i
(
2
)
=
i
i
+
α
h
{\displaystyle t_{i}^{(1)}=t_{i},\quad t_{i}^{(2)}=i_{i}+\alpha h}
k
i
(
1
)
=
f
(
t
i
,
y
i
)
,
k
i
(
2
)
=
f
(
t
i
+
α
h
,
y
i
+
h
β
k
i
(
1
)
,
k
i
=
c
1
k
i
(
1
)
+
c
2
k
i
(
2
)
{\displaystyle k_{i}^{(1)}=f(t_{i},y_{i}),\quad k_{i}^{(2)}=f(t_{i}+\alpha h,y_{i}+h\beta k_{i}^{(1)},\quad k_{i}=c_{1}k_{i}^{(1)}+c_{2}k_{i}^{(2)}}
y
i
+
1
=
y
i
+
h
k
i
=
y
i
+
h
[
c
1
f
(
t
i
,
y
i
)
+
c
2
f
(
t
i
+
α
h
,
y
i
+
h
β
k
i
(
1
)
)
]
{\displaystyle y_{i+1}=y_{i}+hk_{i}=y_{i}+h\left[c_{1}f(t_{i},y_{i})+c_{2}f(t_{i}+\alpha h,y_{i}+h\beta k_{i}^{(1)})\right]}
f
(
t
i
+
α
h
,
y
i
+
h
β
k
i
(
1
)
)
=
f
(
t
i
,
y
i
)
+
f
t
′
(
t
i
,
y
i
)
α
h
f
y
′
(
t
i
,
y
i
)
h
β
k
i
(
1
)
+
O
¯
¯
(
h
2
)
{\displaystyle f(t_{i}+\alpha h,y_{i}+h\beta k_{i}^{(1)})=f(t_{i},y_{i})+f'_{t}(t_{i},y_{i})\alpha hf'_{y}(t_{i},y_{i})h\beta k_{i}^{(1)}+{\overline {\overline {O}}}(h^{2})}
y
i
+
1
=
y
i
+
h
[
c
1
f
(
t
i
,
y
i
)
+
c
2
(
f
(
t
i
,
y
i
)
+
f
t
′
(
t
i
,
y
i
)
α
h
+
h
β
f
y
′
(
t
i
,
y
i
)
k
i
(
1
)
)
]
+
O
¯
¯
(
h
3
)
{\displaystyle y_{i+1}=y_{i}+h\left[c_{1}f(t_{i},y_{i})+c_{2}(f(t_{i},y_{i})+f'_{t}(t_{i},y_{i})\alpha h+h\beta f'_{y}(t_{i},y_{i})k_{i}^{(1)})\right]+{\overline {\overline {O}}}(h^{3})}
- разложение приблзительного решения.
В полученное разложение будем подставлять точное решение, будем считать
α
=
β
{\displaystyle \alpha =\beta }
y
″
(
t
i
)
=
f
t
′
(
t
i
,
y
(
t
i
)
)
+
f
y
′
(
t
−
i
,
y
(
t
i
)
)
f
(
t
i
,
y
i
(
t
i
)
)
{\displaystyle y''(t_{i})=f'_{t}(t_{i},y(t_{i}))+f'_{y}(t-i,y(t_{i}))f(t_{i},y_{i}(t_{i}))}
y
(
t
i
+
1
)
=
y
(
t
i
)
+
h
(
c
1
+
c
2
)
f
(
t
i
,
y
(
t
i
)
)
+
c
2
α
h
2
[
f
t
′
(
t
i
,
y
(
t
i
)
)
+
f
y
′
(
t
i
,
y
(
t
i
)
)
f
(
t
i
,
y
(
t
i
)
)
]
⇒
y
(
t
i
+
1
)
=
y
(
t
i
)
+
h
(
c
1
+
c
2
)
y
′
+
α
c
2
h
2
y
″
(
t
i
)
+
O
¯
¯
(
h
3
)
{\displaystyle y(t_{i+1})=y(t_{i})+h(c_{1}+c_{2})f(t_{i},y(t_{i}))+c_{2}\alpha h^{2}\left[f'_{t}(t_{i},y(t_{i}))+f'_{y}(t_{i},y(t_{i}))f(t_{i},y(t_{i}))\right]\Rightarrow y(t_{i+1})=y(t_{i})+h(c_{1}+c_{2})y'+\alpha c_{2}h^{2}y''(t_{i})+{\overline {\overline {O}}}(h^{3})}
{
c
1
+
c
2
=
1
c
2
α
=
1
2
⇒
{
c
2
=
1
2
α
c
1
=
1
−
1
2
α
,
α
{\displaystyle {\begin{cases}c_{1}+c_{2}=1\\c_{2}\alpha ={\frac {1}{2}}\end{cases}}\Rightarrow {\begin{cases}c_{2}={\frac {1}{2\alpha }}\\c_{1}=1-{\frac {1}{2\alpha }}\end{cases}},\quad \alpha }
- параметр
Имеем:
{
y
i
+
1
=
y
+
h
k
i
k
i
(
1
)
=
f
(
t
i
,
y
i
)
k
i
(
2
)
=
f
(
t
i
+
α
h
,
y
i
+
α
h
k
i
(
1
)
)
k
i
=
(
1
−
1
2
α
)
k
i
(
1
)
+
1
2
α
k
i
(
2
)
{\displaystyle {\begin{cases}y_{i+1}=y_{+}hk_{i}\\k_{i}^{(1)}=f(t_{i},y_{i})\\k_{i}^{(2)}=f(t_{i}+\alpha h,y_{i}+\alpha hk_{i}^{(1)})\\k_{i}=(1-{\frac {1}{2\alpha }})k_{i}^{(1)}+{\frac {1}{2\alpha }}k_{i}^{(2)}\end{cases}}}
При
α
=
1
{\displaystyle \alpha =1}
получаем метод Эйлера-Коши,
α
=
1
2
{\displaystyle \alpha ={\frac {1}{2}}}
- усовершенствованный метод Эйлера.
По построению: локальная погрешность этих методов для
p
{\displaystyle p}
-го порядка точности
l
=
c
h
p
+
1
{\displaystyle l=ch^{p+1}}
Каноническая форма записи методов численного интегрирования
править
y
i
+
1
=
y
i
+
h
f
(
t
i
,
y
i
)
→
y
i
+
1
−
y
i
h
=
f
(
t
i
,
y
i
)
{\displaystyle y_{i+1}=y_{i}+hf(t_{i},y_{i})\to {\frac {y_{i+1}-y_{i}}{h}}=f(t_{i},y_{i})}
- явный одношаговый метод.
1
h
∑
j
=
0
k
α
j
y
i
+
1
−
j
=
Φ
(
t
i
,
y
i
+
1
−
k
,
.
.
.
,
y
i
,
y
i
+
1
,
h
)
,
α
0
≠
0
{\displaystyle {\frac {1}{h}}\sum _{j=0}^{k}\alpha _{j}y_{i+1-j}=\Phi (t_{i},y_{i+1-k},...,y_{i},y_{i+1},h),\alpha _{0}\neq 0}
- канонический вид формулы численного интегрирования. (неявный
k
{\displaystyle k}
-шаговый метод)
Если правая часть этой формулы не содержит
y
i
+
1
{\displaystyle y_{i+1}}
, тогда это явный
k
{\displaystyle k}
-шаговый метод. Все методы Рунге-Кутты - явные, одношаговые.
y
i
+
1
−
y
i
f
=
Φ
(
t
i
,
y
i
,
h
)
{\displaystyle {\frac {y_{i+1}-y_{i}}{f}}=\Phi (t_{i},y_{i},h)}
- общая формула методов Р-К (каноническая формула явных одношаговых методов).
Определение. Аппроксимацией будем называть сеточную функцию
Ψ
i
h
{\displaystyle \Psi _{i}^{h}}
, определенную так: \\
Ψ
i
h
=
1
h
∑
j
=
0
k
α
i
y
(
t
i
+
1
−
j
)
−
Φ
(
t
i
,
y
(
t
i
+
1
−
k
)
,
.
.
.
,
y
(
t
i
)
,
y
(
t
i
+
1
)
,
h
)
{\displaystyle \Psi _{i}^{h}={\frac {1}{h}}\sum _{j=0}^{k}\alpha _{i}y(t_{i+1-j})-\Phi (t_{i},y(t_{i+1-k}),...,y(t_{i}),y(t_{i+1}),h)}
Определение. Будем говорить, что метод сходится с порядком
p
{\displaystyle p}
, если величина
E
(
h
)
=
max
0
≤
1
≤
N
|
y
i
−
y
(
t
i
)
|
≤
c
h
p
{\displaystyle E(h)=\max {0\leq 1\leq N}|y_{i}-y(t_{i})|\leq ch^{p}}
Погрешность аппроксимации метода Эйлера:
Ψ
i
h
=
h
2
2
y
″
(
t
)
{\displaystyle \Psi _{i}^{h}={\frac {h^{2}}{2}}y''(t)}
{
y
′
=
f
(
t
,
y
)
y
(
t
0
)
=
y
0
{
y
~
′
=
f
(
t
,
y
~
)
−
ϕ
y
~
(
t
0
)
=
y
~
0
{\displaystyle {\begin{cases}y'=f(t,y)\\y(t_{0})=y_{0}\end{cases}}\quad {\begin{cases}{\tilde {y}}'=f(t,{\tilde {y}})-\phi \\{\tilde {y}}(t_{0})={\tilde {y}}_{0}\end{cases}}\quad }
- возмущенная задача
ε
(
t
)
=
y
(
t
)
−
y
~
(
t
)
⇒
{
ε
′
=
f
(
t
,
y
)
−
f
(
t
,
y
~
)
+
ϕ
ε
(
t
0
)
=
y
0
−
y
~
0
=
ε
0
{\displaystyle \varepsilon (t)=y(t)-{\tilde {y}}(t)\Rightarrow {\begin{cases}\varepsilon '=f(t,y)-f(t,{\tilde {y}})+\phi \\\varepsilon (t_{0})=y_{0}-{\tilde {y}}_{0}=\varepsilon _{0}\end{cases}}}
f
(
t
,
y
)
−
f
(
t
,
y
~
)
=
f
y
′
(
t
,
y
~
~
)
(
y
−
y
~
)
⇒
{
ε
′
=
λ
ε
+
ϕ
,
λ
=
f
y
′
(
t
,
y
~
~
)
,
y
~
~
=
[
y
,
y
~
]
ε
(
t
0
)
=
ε
0
{\displaystyle f(t,y)-f(t,{\tilde {y}})=f'_{y}(t,{\tilde {\tilde {y}}})(y-{\tilde {y}})\Rightarrow {\begin{cases}\varepsilon '=\lambda \varepsilon +\phi ,\lambda =f'_{y}(t,{\tilde {\tilde {y}}}),{\tilde {\tilde {y}}}=\left[y,{\tilde {y}}\right]\\\varepsilon (t_{0})=\varepsilon _{0}\end{cases}}}
Получим линейное неоднородное уравнение
{
ε
1
′
=
ϕ
ε
(
t
0
)
=
0
{\displaystyle {\begin{cases}\varepsilon '_{1}=\phi \\\varepsilon (t_{0})=0\end{cases}}}
и
{
ε
2
′
=
λ
ε
2
ε
2
(
t
0
)
=
ε
0
⇒
ε
1
=
∫
t
0
t
ϕ
(
τ
)
d
τ
,
ε
2
=
ε
0
exp
{
∫
t
0
t
λ
(
τ
)
d
τ
}
{\displaystyle {\begin{cases}\varepsilon '_{2}=\lambda \varepsilon _{2}\\\varepsilon _{2}(t_{0})=\varepsilon _{0}\end{cases}}\Rightarrow \varepsilon _{1}=\int _{t_{0}}^{t}\phi (\tau )d\tau ,\varepsilon _{2}=\varepsilon _{0}\exp\{\int _{t_{0}}^{t}\lambda (\tau )d\tau \}}
Если
λ
(
t
)
>
0
⇒
{\displaystyle \lambda (t)>0\Rightarrow }
первое слагаемое
ε
2
→
[
t
→
∞
]
∞
{\displaystyle \varepsilon _{2}\rightarrow [t\to \infty ]{}\infty }
. Если
λ
(
t
)
<
0
⇒
ε
2
→
[
t
→
∞
]
0
{\displaystyle \lambda (t)<0\Rightarrow \varepsilon _{2}\rightarrow [t\to \infty ]{}0}
Оценка устойчивости ЗК
|
ε
(
t
)
|
≤
k
(
τ
)
(
|
ε
0
|
+
∫
t
0
τ
|
ϕ
(
t
)
|
d
t
)
,
k
(
τ
)
=
{
exp
{
∫
t
0
τ
λ
(
t
)
d
t
}
,
λ
(
t
)
≥
0
1
,
λ
(
t
)
<
0
{\displaystyle |\varepsilon (t)|\leq k(\tau )(|\varepsilon _{0}|+\int _{t_{0}}^{\tau }|\phi (t)|dt),k(\tau )={\begin{cases}\exp\{\int _{t_{0}}^{\tau }\lambda (t)dt\},&\lambda (t)\geq 0\\1,&\lambda (t)<0\end{cases}}}
Из этого неравенства вытекает устойчивость ЗК.
Рассмотрим возмущенную дискретную ЗК:
1
h
∑
j
=
0
k
y
~
i
+
1
−
j
λ
j
=
Φ
(
t
i
,
y
~
i
+
1
−
k
,
.
.
.
,
y
~
i
,
y
~
i
+
1
,
h
)
−
ϕ
{\displaystyle {\frac {1}{h}}\sum _{j=0}^{k}{\tilde {y}}_{i+1-j}\lambda _{j}=\Phi (t_{i},{\tilde {y}}_{i+1-k},...,{\tilde {y}}_{i},{\tilde {y}}_{i+1},h)-\phi }
ε
0
,
ε
1
,
.
.
.
,
ε
k
−
1
{\displaystyle \varepsilon _{0},\varepsilon _{1},...,\varepsilon _{k-1}}
- погрешности вычисления стартовых точек метода.
Определение. Будем называть метод (1) устойчивым, если:
max
0
≤
i
≤
n
|
y
i
−
y
~
i
|
≤
k
(
τ
)
max
0
≤
i
≤
k
−
1
(
|
ε
i
|
+
∑
i
+
k
N
−
1
h
|
ϕ
|
)
{\displaystyle \max {0\leq i\leq n}|y_{i}-{\tilde {y}}_{i}|\leq k(\tau )\max {0\leq i\leq k-1}(|\varepsilon _{i}|+\sum _{i+k}^{N-1}h|\phi |)}
Многошаговые методы. Методы Адамса
править
Порядок аппроксимации в построенных методах
править
Ψ
i
{\displaystyle \Psi _{i}}
- дискретная функция. Канонический вид:
y
i
+
1
−
y
i
h
=
1
2
f
i
+
1
2
f
i
+
1
{\displaystyle {\frac {y_{i+1}-y_{i}}{h}}={\frac {1}{2}}f_{i}+{\frac {1}{2}}f_{i+1}}
y
′
=
f
(
t
,
y
(
t
)
)
f
i
=
y
′
(
t
i
)
f
i
+
1
=
y
′
(
t
i
+
1
)
{\displaystyle y'=f(t,y(t))\quad f_{i}=y'(t_{i})\quad f_{i+1}=y'(t_{i+1})}
Ψ
i
=
y
(
t
i
+
1
)
−
y
(
t
i
)
h
−
1
2
y
′
(
t
i
)
−
1
2
y
′
(
t
i
+
1
)
=
y
(
t
i
)
+
h
y
′
(
t
i
)
h
2
2
y
″
(
t
i
)
+
h
3
6
y
‴
(
t
i
)
+
O
_
_
(
h
4
)
−
y
(
t
1
)
h
−
1
2
y
′
(
t
i
)
−
1
2
(
y
′
(
t
i
)
+
h
y
″
(
t
i
)
+
h
2
2
y
‴
(
t
i
)
+
O
_
_
(
h
3
)
)
=
h
2
y
″
(
t
i
)
+
h
6
y
‴
(
t
i
)
+
O
_
_
(
h
3
)
−
h
2
y
″
(
t
i
)
−
h
2
4
y
‴
(
t
i
)
+
O
_
_
(
h
3
)
=
−
h
2
12
y
‴
(
t
i
)
+
O
_
_
(
h
3
)
{\displaystyle \Psi _{i}={\frac {y(t_{i+1})-y(t_{i})}{h}}-{\frac {1}{2}}y'(t_{i})-{\frac {1}{2}}y'(t_{i+1})={\frac {y(t_{i})+hy'(t_{i}){\frac {h^{2}}{2}}y''(t_{i})+{\frac {h^{3}}{6}}y'''(t_{i})+{\underline {\underline {O}}}(h^{4})-y(t_{1})}{h}}-{\frac {1}{2}}y'(t_{i})-{\frac {1}{2}}(y'(t_{i})+hy''(t_{i})+{\frac {h^{2}}{2}}y'''(t_{i})+{\underline {\underline {O}}}(h^{3}))={\frac {h}{2}}y''(t_{i})+{\frac {h}{6}}y'''(t_{i})+{\underline {\underline {O}}}(h^{3})-{\frac {h}{2}}y''(t_{i})-{\frac {h^{2}}{4}}y'''(t_{i})+{\underline {\underline {O}}}(h^{3})=-{\frac {h^{2}}{12}}y'''(t_{i})+{\underline {\underline {O}}}(h^{3})}
Ψ
i
=
−
h
2
12
y
‴
(
t
i
)
,
Ψ
=
max
0
≤
i
≤
N
|
Ψ
i
|
≤
M
3
h
2
12
⇒
{\displaystyle \Psi _{i}=-{\frac {h^{2}}{12}}y'''(t_{i}),\Psi =\max {0\leq i\leq N}|\Psi _{i}|\leq {\frac {M_{3}h^{2}}{12}}\Rightarrow }
II аппроксимационный порядок
Если функция
y
(
t
)
{\displaystyle y(t)}
- аналитична и
k
{\displaystyle k}
раз дифференцируема, то
k
{\displaystyle k}
- шаговый метод А-Б и
k
−
1
{\displaystyle k-1}
шаговый метод А-М имеют
k
{\displaystyle k}
-й порядок аппроксимации по
h
{\displaystyle h}
.
Каноническая формула метода:
(
A
)
:
1
h
∑
j
=
0
k
α
j
y
i
+
1
−
j
=
Φ
(
t
i
,
y
i
+
1
−
k
,
.
.
.
,
y
i
,
y
i
+
1
)
{\displaystyle (A):{\frac {1}{h}}\sum _{j=0}^{k}\alpha _{j}y_{i+1-j}=\Phi (t_{i},y_{i+1-k},...,y_{i},y_{i+1})}
- дискретная ЗК. Стартовые точки
y
0
,
y
1
,
.
.
.
,
y
k
−
1
{\displaystyle y_{0},y_{1},...,y_{k-1}}
(
B
)
:
1
h
∑
j
=
0
k
α
j
y
i
+
1
−
j
∗
=
Φ
(
t
i
,
y
i
+
1
−
k
∗
,
.
.
.
,
y
i
∗
,
y
i
+
1
∗
)
{\displaystyle (B):{\frac {1}{h}}\sum _{j=0}^{k}\alpha _{j}y_{i+1-j}^{*}=\Phi (t_{i},y_{i+1-k}^{*},...,y_{i}^{*},y_{i+1}^{*})}
- дискретная ЗК с возмущенными данными. Стартовые точки
y
0
,
y
1
,
.
.
.
,
y
k
−
1
{\displaystyle y_{0},y_{1},...,y_{k-1}}
Определение. (Устойчивость
k
{\displaystyle k}
-шагового метода)
max
k
≤
i
≤
N
|
y
i
−
y
i
∗
|
≤
c
{
max
0
≤
i
≤
k
−
1
|
y
i
−
y
i
∗
|
+
max
|
∑
i
=
k
N
h
Ψ
i
|
}
{\displaystyle \max {k\leq i\leq N}|y_{i}-y_{i}^{*}|\leq c\{\max {0\leq i\leq k-1}|y_{i}-y_{i}^{*}|+\max |\sum _{i=k}^{N}h\Psi _{i}|\}}
Теорема. Если численный метод устойчив и имеет порядок аппроксимации
p
{\displaystyle p}
, то он сходится с
p
{\displaystyle p}
-м порядком точности по
h
{\displaystyle h}
.
Уравнения высокого порядка и системы ДУ
править
{
y
1
′
=
f
1
(
t
,
y
,
.
.
.
,
y
m
)
y
2
′
=
f
2
(
t
,
y
,
.
.
.
,
y
m
)
⋮
y
m
′
=
f
m
(
t
,
y
,
.
.
.
,
y
m
)
y
1
′
(
t
0
)
=
y
1
0
⋮
y
m
′
(
t
0
)
=
y
m
0
{\displaystyle {\begin{cases}y'_{1}=f_{1}(t,y,...,y_{m})\\y'_{2}=f_{2}(t,y,...,y_{m})\\\vdots \\y'_{m}=f_{m}(t,y,...,y_{m})\\y'_{1}(t_{0})=y_{1}^{0}\\\vdots \\y'_{m}(t_{0})=y_{m}^{0}\end{cases}}}
- нелинейная система ОДУ 1-го
y
=
(
y
1
⋮
y
m
)
y
0
=
(
y
1
0
⋮
y
m
0
)
F
=
(
f
1
(
t
,
y
1
,
.
.
.
,
t
m
)
⋮
f
1
(
t
,
y
1
,
.
.
.
,
t
m
)
)
{\displaystyle y={\begin{pmatrix}y_{1}\\\vdots \\y_{m}\end{pmatrix}}\quad y^{0}={\begin{pmatrix}y_{1}^{0}\\\vdots \\y_{m}^{0}\end{pmatrix}}\quad F={\begin{pmatrix}f_{1}(t,y_{1},...,t_{m})\\\vdots \\f_{1}(t,y_{1},...,t_{m})\end{pmatrix}}}
метод Эйлера:
{
y
′
=
F
(
t
,
y
)
y
(
t
0
)
=
y
0
{\displaystyle {\begin{cases}y'=F(t,y)\\y(t_{0})=y_{0}\end{cases}}}
y
i
+
1
=
y
i
+
h
F
(
t
i
,
y
i
)
{\displaystyle y_{i+1}=y_{i}+hF(t_{i},y_{i})}
{
a
m
y
(
m
)
(
t
)
+
a
m
−
1
y
(
m
−
1
)
+
.
.
.
+
a
1
y
′
(
t
)
+
a
0
y
(
t
)
=
f
(
t
)
y
(
t
0
)
=
y
0
⋮
y
(
m
−
1
)
(
t
0
)
=
y
0
(
m
−
1
)
{\displaystyle {\begin{cases}a_{m}y^{(m)}(t)+a_{m-1}y^{(m-1)}+...+a_{1}y'(t)+a_{0}y(t)=f(t)\\y(t_{0})=y_{0}\\\vdots \\y^{(m-1)}(t_{0})=y_{0}^{(m-1)}\end{cases}}}
y
1
(
t
)
=
y
(
t
)
;
y
2
(
t
)
=
y
′
(
t
)
=
;
.
.
.
;
y
m
(
t
)
=
y
(
m
−
1
)
(
t
)
{\displaystyle y_{1}(t)=y(t);y_{2}(t)=y'(t)=;...;y_{m}(t)=y^{(m-1)}(t)}
- замены
{
y
1
′
=
y
2
y
2
′
=
y
3
⋮
y
m
−
1
′
=
y
m
y
m
′
=
f
(
t
)
−
a
m
−
1
y
m
−
.
.
.
−
a
0
y
1
a
m
y
=
(
y
1
⋮
y
m
)
F
=
(
y
2
⋮
y
m
f
(
t
)
−
a
m
−
1
y
m
−
.
.
.
−
a
0
y
1
a
m
)
{\displaystyle {\begin{cases}y'_{1}=y_{2}\\y'_{2}=y_{3}\\\vdots \\y'_{m-1}=y_{m}\\y'_{m}={\frac {f(t)-a_{m-1}y_{m}-...-a_{0}y_{1}}{a_{m}}}\end{cases}}y={\begin{pmatrix}y_{1}\\\vdots \\y_{m}\end{pmatrix}}F={\begin{pmatrix}y_{2}\\\vdots \\y_{m}\\{\frac {f(t)-a_{m-1}y_{m}-...-a_{0}y_{1}}{a_{m}}}\end{pmatrix}}}
Устойчивость численных методов
править
1
h
∑
j
=
0
k
α
j
y
n
+
1
−
j
=
Φ
(
t
n
,
y
n
+
1
−
k
,
.
.
.
,
y
n
,
y
n
+
1
,
h
)
(
1
)
{\displaystyle {\frac {1}{h}}\sum _{j=0}^{k}\alpha _{j}y_{n+1-j}=\Phi (t_{n},y_{n+1-k},...,y_{n},y_{n+1},h)\quad (1)}
α
0
,
α
1
,
.
.
.
,
α
k
−
1
α
0
≠
0
{\displaystyle \quad \alpha _{0},\alpha _{1},...,\alpha _{k-1}\quad \alpha _{0}\neq 0}
Чтобы облегчить задачу будем исследовать устойчивость по начальным данным (без правой части)
y
0
,
y
1
,
.
.
.
,
y
k
−
1
{\displaystyle y_{0},y_{1},...,y_{k-1}}
Рассмотрим модельную задачу
y
′
=
0
{\displaystyle y'=0}
на нуль-устойчивость
1
h
∑
j
=
0
k
α
j
y
n
+
1
−
j
=
0
{\displaystyle {\frac {1}{h}}\sum _{j=0}^{k}\alpha _{j}y_{n+1-j}=0}
(
2
)
α
0
y
n
+
1
+
α
1
y
n
+
.
.
.
+
y
k
α
n
+
1
−
k
=
0
{\displaystyle (2)\quad \alpha _{0}y_{n+1}+\alpha _{1}y_{n}+...+y_{k}\alpha _{n+1-k}=0}
- конечно-разностное уравнение. Стартовые точки
y
n
+
1
−
k
.
.
.
y
n
+
1
{\displaystyle y_{n+1-k}...y_{n+1}}
Затем мы находим
y
0
,
y
1
,
.
.
.
,
y
k
−
1
,
y
k
,
y
k
+
1
,
.
.
.
,
y
N
{\displaystyle y_{0},y_{1},...,y_{k-1},y_{k},y_{k+1},...,y_{N}}
Возмущенная задача:
y
0
∗
.
.
.
y
k
−
1
∗
{\displaystyle y_{0}^{*}...y_{k-1}^{*}}
α
0
y
n
+
1
∗
+
α
1
y
n
∗
+
.
.
.
+
α
k
y
n
+
1
−
k
∗
=
0
{\displaystyle \alpha _{0}y_{n+1}^{*}+\alpha _{1}y_{n}^{*}+...+\alpha _{k}y_{n+1-k}^{*}=0}
y
i
∗
i
=
0
,
N
¯
{\displaystyle y_{i}^{*}\quad i={\overline {0,N}}}
- проверим на устойчивость
α
0
ε
n
+
1
+
α
1
ε
n
+
.
.
.
+
α
k
ε
n
+
1
−
k
=
0
{\displaystyle \alpha _{0}\varepsilon _{n+1}+\alpha _{1}\varepsilon _{n}+...+\alpha _{k}\varepsilon _{n+1-k}=0}
ε
n
=
ε
¯
q
n
{\displaystyle \varepsilon _{n}={\overline {\varepsilon }}q^{n}}
,
q
{\displaystyle q}
некоторое число.
ε
¯
{\displaystyle {\overline {\varepsilon }}}
- некоторая первоначальная погрешность.
Получим характеристическое уравнения
k
{\displaystyle k}
- шагового метода:
α
0
q
k
+
α
1
q
k
−
1
+
.
.
.
+
α
k
=
0
{\displaystyle \alpha _{0}q^{k}+\alpha _{1}q^{k-1}+...+\alpha _{k}=0}
P
(
q
)
=
∑
j
=
0
k
α
k
q
k
−
j
{\displaystyle P(q)=\sum _{j=0}^{k}\alpha _{k}q^{k-j}}
- характеристический полном
Уравнение
P
(
q
)
=
0
{\displaystyle P(q)=0}
имеет ровно
k
{\displaystyle k}
корней (кратных и простых). Пусть
q
1
{\displaystyle q_{1}}
- кратный корень кратности
r
{\displaystyle r}
, остальные простые:
q
=
c
1
(
1
)
q
1
n
+
c
1
(
2
)
n
q
n
+
.
.
.
+
c
1
(
r
−
1
)
n
r
−
1
q
1
n
+
c
2
q
2
n
+
.
.
.
+
c
m
q
m
n
(
3
)
{\displaystyle q=c_{1}^{(1)}q_{1}^{n}+c_{1}^{(2)}nq^{n}+...+c_{1}^{(r-1)}n^{r-1}q_{1}^{n}+c_{2}q_{2}^{n}+...+c_{m}q_{m}^{n}\quad (3)}
Общее решение уравнения:
|
q
i
|
<
1
,
n
→
∞
,
|
q
|
→
0
{\displaystyle |q_{i}|<1,n\to \infty ,|q|\to 0}
Определение. Будем говорить, что выполнено корневое условие, если все корни полнинома устойчивости по модулю меньше 1, или на гранирце круга
|
q
i
|
=
1
{\displaystyle |q_{i}|=1}
кратных корней нет
Теорема. Если выполнено корневое условие, то метод является устойчивым (называется
O
{\displaystyle O}
-устойчивым)
Утверждение. Все методы Рунге-Кутты и методы Адамаса являются нуль устойчивыми методами.
Р-К:
P
(
q
)
=
q
−
1
{\displaystyle P(q)=q-1}
А:
α
0
q
k
+
α
1
q
k
−
1
+
0
q
k
−
2
+
.
.
.
+
0
=
0
α
0
=
1
,
α
1
=
−
1
P
(
q
)
=
q
k
−
q
k
−
1
{\displaystyle \alpha _{0}q^{k}+\alpha _{1}q^{k-1}+0q^{k-2}+...+0=0\quad \alpha _{0}=1,\alpha _{1}=-1\quad P(q)=q^{k}-q^{k-1}}
Наличие
O
{\displaystyle O}
-устойчивости означает, что метод устойчив на конечном отрезке
[
t
0
,
T
]
{\displaystyle \left[t_{0},T\right]}
. Если
T
→
∞
{\displaystyle T\to \infty }
следует применять абсолютно устойчивые методы.
Модельная задача:
{
y
′
=
λ
y
λ
∈
C
y
(
t
0
)
=
y
0
R
e
λ
<
0
}
⇒
{\displaystyle \left.{\begin{cases}y'=\lambda y&\lambda \in C\\y(t_{0})=y_{0}&Re\lambda <0\end{cases}}\right\}\Rightarrow }
Ассимптотическая устойчивость
1
h
∑
j
=
0
k
α
j
y
n
+
1
−
j
=
λ
∑
j
=
0
k
β
j
(
λ
h
)
y
n
+
1
−
j
z
=
λ
h
{\displaystyle {\frac {1}{h}}\sum _{j=0}^{k}\alpha _{j}y_{n+1-j}=\lambda \sum _{j=0}^{k}\beta _{j}(\lambda h)y_{n+1-j}\qquad z=\lambda h}
∑
j
=
0
k
(
α
j
−
z
β
j
(
z
)
)
y
n
+
1
−
j
=
0
γ
=
α
j
−
z
β
j
(
z
)
{\displaystyle \sum _{j=0}^{k}(\alpha _{j}-z\beta _{j}(z))y_{n+1-j}=0\qquad \gamma =\alpha _{j}-z\beta _{j}(z)}
Для невозмущенной задачи:
∑
j
=
0
k
γ
j
y
n
+
1
−
j
=
0
y
0
,
.
.
.
,
y
k
−
1
{\displaystyle \sum _{j=0}^{k}\gamma _{j}y_{n+1-j}=0\qquad y_{0},...,y_{k-1}}
- стартовые точки
P
(
q
,
z
)
=
γ
0
q
k
+
γ
1
q
k
−
1
+
.
.
.
+
γ
k
=
0
{\displaystyle P(q,z)=\gamma _{0}q^{k}+\gamma _{1}q^{k-1}+...+\gamma _{k}=0}
. Ищем
q
1
,
q
2
,
.
.
.
,
q
n
{\displaystyle q_{1},q_{2},...,q_{n}}
(зависят от
z
{\displaystyle z}
)
Для явного метода Эйлера:
y
n
+
1
−
y
n
h
=
λ
y
n
{\displaystyle {\frac {y_{n+1}-y_{n}}{h}}=\lambda y_{n}}
ε
n
+
1
−
(
1
+
λ
h
)
ε
n
=
0
q
−
(
1
+
z
)
=
0
⇒
q
1
=
1
+
z
|
1
+
z
|
<
1
{\displaystyle \varepsilon _{n+1}-(1+\lambda h)\varepsilon _{n}=0\quad q-(1+z)=0\quad \Rightarrow \quad q_{1}=1+z\quad |1+z|<1}
−
1
<
1
+
z
<
1
⇒
z
∈
(
−
2
,
0
)
z
=
λ
h
;
λ
<
0
⇒
0
<
h
<
2
−
λ
{\displaystyle -1<1+z<1\quad \Rightarrow \quad z\in (-2,0)\quad z=\lambda h;\lambda <0\Rightarrow 0<h<{\frac {2}{-\lambda }}}
Метод называется устойчивым при ограничении на шаг
h
<
2
|
λ
|
{\displaystyle h<{\frac {2}{|\lambda |}}}
Для неявного:
y
n
+
1
−
y
n
h
=
λ
y
n
+
1
⇒
1
1
−
λ
h
=
q
⇒
|
1
−
λ
h
|
<
1
{\displaystyle {\frac {y_{n+1}-y_{n}}{h}}=\lambda y_{n+1}\Rightarrow {\frac {1}{1-\lambda h}}=q\Rightarrow |1-\lambda h|<1}
- верно т.к.
λ
<
0
{\displaystyle \lambda <0}
Ограничений на
h
{\displaystyle h}
нет. Область устойчивости:
|
1
−
z
|
<
1
{\displaystyle |1-z|<1}
Определение. Метод называется абсолютно устойчивым при конкретном значении
z
=
λ
h
{\displaystyle z=\lambda h}
если все его корни полинома устойчивости лежат внутри единичного круга, а на границе нет кратных корней. Множество точек
z
{\displaystyle z}
комплексной плоскости в которых выполнено корневое условие называется областью абсолютной устойчивости метода.
Определение. Численный метод называется
A
{\displaystyle A}
-устойчивым, если область абсолютой устойчивости включает в себя полу плоскость
R
e
z
<
0
{\displaystyle Rez<0}
Неявный метод Эйлера
A
{\displaystyle A}
-устойчив