Многошаговые методы

Многошаговые методы Адамса

Так же как и неявные методы Рунге-Кутты, многошаговые методы хорошо решают жесткие задачи. Одно в силу барьера Далквиста, они никогда не будут иметь порядок аппроксимации больше 2.

В отличие от методов Рунге-Кутта, многошаговые методы основаны не на формуле Тейлора непосредственно, а на сеточной интерполяции.

Многошаговые методы построены на том, что для вычисления значения \(x_{k+1}\) применяются несколько предыдущих точек. При этом предыдущие несколько точек должны быть вычислены одношаговым методом (Эйлера, Рунге-Кутта).

Интерполяционный метод Адамса использует предыдущие значения и неизвестное значение.

\(y_{k+1}=y_{k}+h\left(f_{k+1}-\frac{\Delta f_{k}}{2}-\frac{\Delta^{2} f_{k-1}}{12}-\frac{\Delta^{3} f_{k-2}}{24}\right)\)

Экстраполяционный метод Адамса использует только предыдущие:

\(y_{k+1}=y_{k}+h\left(f_{k}+\frac{\Delta f_{k-1}}{2}+\frac{5 \Delta^{2} f_{k_{2}}}{12}+\frac{3 \Delta^{3} f_{k-3}}{8}\right)\)

Здесь

\[\begin{split} \begin{aligned} &\Delta f_{k}=f_{k+1}-f_{k}\\ &\Delta^{2} f_{k}=f_{k+2}-2 f_{k+1}+f_{k} \\ &\Delta^{3} f_{k}=f_{k+3}-3 f_{k+2}+3 f_{k+1}-f_{k} \end{aligned} \end{split}\]

и так далее.

Эти и другие варианты формул метода Адамса можно вывести самостоятельно, используя интерполяционный многочлен Ньютона.

Метод прогноза-коррекции Адамса

Многошаговые методы построены на том, что для вычисления значения \(x_{k+1}\) применяются несколько предыдущих точек. При этом предыдущие несколько точек должны быть вычислены одношаговым методом (Эйлера, Рунге-Кутта). Рассмотрим решение дифореренциального уравнения

\[ F\left(t, x, y, \ldots, y^{(n)}\right)=0 \]

на интервале \(\left[t_{i} ; t_{i+1}\right] .\)

Будем считать, что решение в точках \(t_{0}, t_{1}, t_{2}, \ldots, t_{i}\) уже найдено, и значения в этих точках будем использовать для нахождения значения \(x\left(t_{i+1}\right)\).

Проинтегрируем уравнение на интервале \(\left[t_{i} ; t_{i+1}\right]\) и получим соотношение:

\[ x\left(t_{i+1}\right)=x\left(t_{i}\right)+\int_{t_{i}}^{t_{i+1}} f(t, x(t)) d t \]

При вычислении интеграла вместо функции \(f(t, x(t))\) будем использовать интерполяционный полином Лагранжа, построенный по точкам \(\left(t_{i-3}, x_{i-3}\right),\left(t_{i-2}, x_{i-2}\right),\left(t_{i-1}, x_{i-1}\right),\left(t_{i}, x_{i}\right) .\)

Подставив полином Лагранжа в интеграл, получаем первое приближение (прогноз) \(\tilde{t}_{ i + 1 }\) для значения функции в точке \(t_{i+1}\)

\[\begin{split} \\ \tilde{x}_{i+1}=x_{i}+\frac{h}{24}\left(-9 f\left(t_{i-3}, x_{i-3}\right)+37 f\left(t_{i-2}, x_{i-2}\right)-59 f\left(t_{i-1}, x_{i-1}\right)+55 f\left(t_{i}, x_{i}\right)\right) . \end{split}\]

Как только \(\widetilde{x}_{i+1}\) вычислено, его можно использовать.

Следующий полином Лагранжа для функции \(f(t, x(t))\) построим по точкам \(\left(t_{i-2}, x_{i-2}\right),\left(t_{i-1}, x_{i-1}\right),\left(t_{i}, x_{i}\right)\) и новой точке \(\left(t_{i+1}, x_{i+1}\right)\), после чего подставляем его и получаем второе приближение (корректор):

\[ x_{i+1}=x_{i}+\frac{h}{24}\left(f\left(t_{i-2}, x_{i-2}\right)-5 f\left(t_{i-1}, x_{i-1}\right)+19 f\left(t_{i}, x_{i}\right)+9 f\left(t_{i+1}, \tilde{x}_{i+1}\right)\right) \]

Таким образом, для вычисления значения \(x\left(t_{i+1}\right)\) методом Адамса необходимо последовательно применять формулы прогноза и коррекции, а первые четыре точки можно получить методом Рунге-Кутта.

Метод Милна

Отличие метода Милна от метода Адамса состоит в том, что в качестве интерполяционного полинома используется полином Ньютона.

Подставив вместо функции \(f(t, x(t))\) интерполяционный полином Ньютона, построенный по точкам \(\left(t_{k-3}, x_{k-3}\right),\left(t_{k-2}, x_{k-2}\right),\left(t_{k-1}, x_{k-1}\right),\left(t_{k}, x_{k}\right)\) получаем первое приближение - прогноз Милна \(\widetilde{x}_{k+1}\) для значения функции в точке \(t_{k+1}\) :

\[ x_{k+1}=x_{k-3}+\frac{4 h}{3}\left(2 f\left(t_{k-2}, x_{k-2}\right)-f\left(t_{k-1}, x_{k-1}\right)+2 f\left(t_{k}, x_{k}\right)+2 f\left(t_{k}, \tilde{x}_{k}\right)\right) . \]

Следующий полином Ньютона для функции \(f(t, x(t))\) построим по точкам \(\left(t_{k-2}, x_{k-2}\right),\left(t_{k-1}, x_{k-1}\right),\left(t_{k}, x_{k}\right)\) и новой точке \(\left(t_{k+1}, x_{k+1}\right)\), после чего подставляем его и получаем второе приближение - корректор Милна:

\[ x_{k+1}=x_{k-1}+\frac{h}{3}\left(f\left(t_{k-1}, x_{k-1}\right)+4 f\left(t_{k}, x_{k}\right)+f\left(t_{k+1}, \tilde{x}_{k+1}\right)\right) \cdot \]

В методе Милна для вычисления значения \(x\left(t_{k+1}\right)\) также необходимо последовательно применять формулы прогноза и коррекции, а первые четыре точки можно получить методом Рунге-Кутта.

Существует модифицированный метод Милна. В нем сначала вычисляется первое приближение по формуле прогноза Милна, затем вычисляется управляющий параметр:

\[ m_{k+1}=\tilde{x}_{k+1}+\frac{28}{29}\left(x_{k}-\tilde{x}_{k}\right) \]

после чего вычисляется значение второго приближения - корректор Милна - по формуле:

\[ x_{k+1}=x_{k-1}+\frac{h}{3}\left(f\left(t_{k-1}, x_{k-1}\right)+4 f\left(t_{k}, x_{k}\right)+f\left(t_{k+1}, m_{k+1}\right)\right) \]

Анализ многошаговых методов на устойчивость

Пусть у нас есть произвольная многошаговая схема. Применим её для решения любимого модельного уравнения Далквиста \(\mathbf{x}^{\prime}=\lambda \mathbf{x}\).

В таком раскладе, многошаговая схема приводится к следующей форме:

\[ \alpha_s \mathbf{x}_{n+s}+\ldots+\alpha_0 \mathbf{x}_n=\mathbf{0} \]

При этом коэффициенты \(\alpha_i\) - явные функции от \(z=\lambda \cdot \tau\).

Представим решение как \(\mathbf{x}_i=\mathbf{z} \zeta^i\) и подставим его в общий вид :

\[ \mathbf{z} \zeta^n\left(\alpha_s \zeta^s+\ldots+\alpha_0\right)=\mathbf{0} \]

Отсюда \(\zeta\) должно удовлетворять так называемому характеристическому уравнению

\[ \alpha_s \zeta^s+\ldots+\alpha_0=0 \]

Исследование характеристического уравнения усложняется, когда у него имеются кратные корни кратности \(m\).

Таким образом, если характеристический многочлен имеет корни \(\zeta_1, \ldots, \zeta_l\) кратностей \(m_1, \ldots, m_l\) соответственно, то общее решение уравнения:

\[ \mathbf{x}_n=\mathbf{p}_1(n) \zeta_1^n+\ldots+\mathbf{p}_l(n) \zeta_l^n \]

где \(\mathbf{p}_j(n)\) - многочлены степеней \(\boldsymbol{m}_j-\mathbf{1}\), коэффициенты которых определяются из начальных условий. Поэтому для ограниченности решений \(\mathbf{x}_i\) необходимо потребовать:

1. Корни кратности более единицы \(\zeta\) лежат внутри единичной окружности \((|\zeta|<1)\),

2. Все корни на единичной окружности \((|\zeta|=1)\) должны быть простыми.

Определение. Многошаговые методы, удовлетворяющие этим требованиям, называются устойчивыми по Далквисту или D-устойчивыми.

Определение. Областью (абсолютной) устойчивости назовём множество \(S=\left \{z \in \mathbb{C}: \; \forall |\zeta_i(z)| \leq 1 \right \}\)

Пример исследования на устойчивость многошагового метода. (доделать)

Покажите, что метод ФДН 2-го порядка \(\frac{3 y_{l+1}-4 y_l+y_{l-1}}{2 h}=f\left(x_{l+1}, y_{l+1}\right)\) является \(A\)-устойчивым.

Решение.

Действуя, как и в случае методов Рунге-Кутты, получаем следующее уравнение для \(\Delta_l: \frac{3 \Delta_{l+1}-4 \Delta_l+\Delta_{l-1}}{2 h}=J \Delta_{l+1}\). Подставляя \(z=J h\), получим следующее линейное разностное уравнение: \(3 \Delta_{l+1}-4 \Delta_l+\Delta_{l-1}=2 z \Delta_{l+1}\). Его общее решение (при \(3-2 z \neq 0\), что выполняется, так как нам интересна область \(\operatorname{Re}(z) \leq 0\) ) равно \(\Delta_l=C_1 q_1^l+C_2 q_2^l\), если корни характеристического уравения \((3-2 z) q^2-4 q+1=0\) различны и \(\Delta_l=\left(C_1+C_2 l\right) q^l\), если корни кратные. Область устойчивости \(S\) (где решение не возрастает при больших \(l)\) - совокупность точек \(z\), для которых все решения \(\left|q_i(z)\right| \leq 1\) ( \(\left|q_i(z)\right|<1\) для кратных корней, что выполняется, т.к. \(q_1=q_2=1 / 2\) при \(z=-1 / 2)\). Вместо исследования неоднозначного отображения из комплексной плоскости \(z\) в комплексную плоскость \(q\) рассмотрим обратное (однозначное) отображение, задаваемое формулой \(z=\frac{3}{2}-\frac{2}{q}+\frac{1}{2 q^2}\). Внешность единичного круга в плоскости \(q\) целиком отображается в правую полуплоскость \(C^{+}\)плоскости \(z\). Действительно, при \(q=\rho e^{i \varphi}, \rho \geq 1\), \(\varphi \in[0,2 \pi]\) имеем

\[\begin{split} \begin{aligned} & \operatorname{Re}\left(z\left(\rho e^{i \varphi}\right)\right)=\operatorname{Re}\left(\frac{3}{2}-\frac{2}{\rho e^{i \varphi}}+\frac{1}{2 \rho^2 e^{2 i \varphi}}\right)=\frac{3}{2}-\frac{2}{\rho} \cos (\varphi)+ \\ & +\frac{1}{2 \rho}\left(2 \cos ^2(\varphi)-1\right)=\left(1-\frac{\cos (\varphi)}{\rho}\right)^2+\frac{1}{2}\left(1-\frac{1}{\rho^2}\right) \geq 0 . \end{aligned} \end{split}\]

Значит, точки из левой полуплоскости плоскости \(z\) при отображении в плоскость \(q\) могут перейти только во внутренность единичного круга либо его границу, и схема \(A\)-устойчива.