О методе Адамса-Бэшфорта-Мултона

[maisvendoo]

Итак, разобрался вроде с построением интерполяционных полиномов, написал небольшой код на Maple позволяющий получить явный и неявный метод любого порядка.

Вопросы такого плана, может они покажутся глупыми (но сделайте скидку - две недели пялюсь в колонки цифр без особого успеха), тем не менее.

1. На каком шаге алгоритма необходимо сдвигать значения предвычисленных производных $f_{n-N},\cdots, f_{n-3}, f_{n-1}, f_{n}$ - после прогноза или уже после коррекции?
2. Необходимы ли итерации на этапе коррекции для достижения запданной точности?

P.S.: Кстати, если интересно, вот упомянутый код

restart;
with(LinearAlgebra):

# Явный метод 
p := 4;

X := Vector(p):
Y := Vector(p):

for i from 1 to p do
  X[p - i + 1] := x - (i-1)*h;
  Y[p - i + 1] := f[k - i + 1];
end do:

poly := interp(X, Y, u):
y[k+1] = y[k] + factor(int(poly, u=x..x + h));

# Неявный метод
for i from 1 to p do
 X[p - i + 1] := x + (2 - i)*h;
 Y[p - i + 1] := f[k + (2 - i)];
end do: 

poly := interp(X, Y, u):
y[k+1] = y[k] + factor(int(poly, u=x..x + h));

[maisvendoo]

Думаю необходимы некоторые пояснения. В своей работе я использую предиктор-корректор метод - предсказание по формуле Адамса-Бэшфорта, коррекция по Адамсу-Мултону. Алгоритм таков

0. Делаю 3 шага методом Рунге-Кутты 4-го порядка для определения $\vec{f}_{k-3}, \vec{f}_{k-2}, \vec{f}_{k-1}$
1. На этапе предсказания вычисляю производную $\vec{f}_k\left(t_k,\vec{y}_k\right)$ и вектор
$\vec{p}_k = \vec{y}_k + \frac{h}{24}\left(55\vec{f}_k -59\vec{f}_{k-1}+37\vec{f}_{k-2}-9\vec{f}_{k-3} \right) $
2. Вычисляю производную $\vec{f}_{k+1}\left(t_{k+1}, \vec{p}_k\right)$
3. Выполняю коррекцию
$\vec{y}_k^{*} = \vec{y}_k + \frac{h}{24}\left(9\vec{f}_{k+1} +19\vec{f}_{k}-5\vec{f}_{k-1}+\vec{f}_{k-2} \right) $
4. Проверяю точность
$\max{ \left|\vec{y}_k^* - \vec{p}_k\right| < \varepsilon}$
5. Если точность не достигнута, возвращаюсь к шагу 2.
6. Если достигнута точность, то
$\vec{y}_{k+1} = \vec{y}_k^*$
и сдвигаю производные
$\vec{f}_{k-1} = \vec{f}_{k}$
$\vec{f}_{k-2} = \vec{f}_{k-1}$
$\vec{f}_{k-3} = \vec{f}_{k-2}$

Есть сомнения что делаю всё верно...

Данное применяю к жесткой задаче (дифф. уравнения движения поезда с учетом продольных динамических реакций). Не лучше ли для жесткой задачи построить неявную схему на основе формулы из шага 3?

[ShMaxG]

А что у Вас означают $k$ и $i$? Что такое $f_k(t_i,y_i)$?

-- Вт дек 10, 2013 14:25:49 --

P.S. Не знаю, насколько для жестких систем хорош Адамс , но может формулы дифференцирования назад больше подойдут (ФДН, BDF)?

[maisvendoo]

Исправился. Вместо i надо k. Индексы массивов в голове сидят.
На удивление метод Гира не дает такого результата как адамс применительно к моей проблеме

$\vec{f}_{k}\left(t_k,\vec{y}_{k}\right)$ - вектор правых частей системы ОДУ, в моем случае существенно нелинейная вектор-функция

[ShMaxG]

Вроде верно. А в чем конкретно сомнения?

maisvendoo в сообщении #798632 писал(а):

Не лучше ли для жесткой задачи построить неявную схему на основе формулы из шага 3?

Ну стройте, только придумайте сначала начальное приближение для соответствующей итеративной процедуры. Использование же решения явного метода кажется вполне естественным...

Для жесткого метода еще можно метод Рунге-Кутты выбрать неявным и поточнее.

maisvendoo в сообщении #798632 писал(а):

5. Если точность не достигнута, возвращяюсь к шагу 2.

К шагу 2, но только $p_k$ там надо заменить на $y^{*}_k$.

[maisvendoo]

Да, конечно выполняю указанное Вами присваивание.
Сомнения были в верности алгоритма.

К неявным схемам идти не хочется из-за необходимости вычислять матрицу Якоби. В правую часть входит сила, вида

где T - усилие в межвагонной связи; s - деформация связи. Стрелки на графике указывают на направление скорости деформации.
В принципе это допускает получение матрицы Якоби в аналитическом виде. Получается следующая матрица
$J = \left(\begin{matrix}O & E \\ \cfrac{\partial \vec{a}}{\partial \vec{s}} & \cfrac{\partial \vec{a}}{\partial \vec{v}}\end{matrix}\right)$
где O, E - соответственно нулевая и единичная матрицы размерности n/2 (n - размерность системы ОДУ); $\vec{a}$ - вектор ускорений подвижных единиц. Однако даже если вместо указанных производных от ускорения поставить нули, метод Гира и неявный Эйлер находят решение, причем проходящее по энергетическому балансу на основе теоремы об изменении кинетической энергии.

Решение неявным методом расходится при попадании на участок LN указанной характеристики, а Адамс дает вполне приемлемые результаты (правда с некоторыми оговорками)

P. S.: Приятно удивлен возможностью использовать LaTeX на форуме )