Несколько вопросов по численному решению ОДУ

[Gickle]

День добрый. Имеются две учебные задачи по численному решению ОДУ.

1. Охлаждение плазмы.

Задача Коши:

$\frac{dT}{dt} = -\frac{2}{3 k_{B}}n_{H}\Lambda(T)\text{,} \quad T(0) = T_{init},$

где $\Lambda(T) = \begin{cases} \Lambda_{0} \left(\frac{T}{T_{0}}\right)^{\alpha}&\text{for $T \leq T_{0}$;}\\ \Lambda_{0} \left(\frac{T}{T_{0}}\right)^{\beta}&\text{for $T > T_{0}$.} \end{cases}$

На всякий случай перечислю константы всякие: $k_{B}$ - постоянная Больцмана, $n_H = 10^6 \,\text{$m^{-3}$}$, $T_0 = 20000 \,\text{$K$}$, $T_{init} = 10^7 \,\text{$K$}$, $\Lambda_0 = 10^{-35} \,\text{$J m^3 s^{-1}$}$ $\alpha = 10$, $\beta = -0.5$. Главное, что поведение функции $\Lambda(T)$ серьёзным образом меняется в точке $T_{0}$.

Сказано использовать метод Рунге-Кутты второго порядка. Итерационная схема при этом выглядит так:

$\\ k_1 = f(y_n, t_n)\\ k_2 = f(y_n + k_1 \Delta t_n, t_n)\\ y_{n+1} = y_n + \frac{k_1 + k_2}{2} \Delta t_n $

Считать нужно до момента, пока $T > T_{end}$, где $T_{end} = 6000 \,\text{$K$} $, а шаг по времени взять равным $\Delta t = 10^10 \text{$s$} $. В целом проблем особо никаких нет, всё хорошо.

Потом говорится, мол, поиграйте с шагом по времени. Вот тут начинаются проблемы: если поставить где-то в 10 раз больше шаг, то возникают проблемы из-за того, что $f(T,t)$ задана кусочно (в виде системы). То есть для вычисления где-то на шаге, скажем, в $30000 \,\text{$K$}$ мы используем ещё "здоровую" правую часть со степенью $\alpha = 10$, в результате чего из-за большого шага по времени "проскакиваем" область, в которой температура должна меняться по другому закону. Ну, лекарство от этого понятное - в области перед $T_0$ сделать шаг помельче просто, например.

Однако дальше просят подрубить к этому делу адаптивную регулировку размера шага, где критерием допустимости выступает абсолютная ошибка $\Delta T_{err}^{max} = 50 \,\text{$K$}$. И в конце, цитирую: "Confirm that your scheme is robust to large changes of the timestep size given as input for the first step". Правильно ли я понимаю, что от меня хотят, чтобы я показал, что схема устойчива относительно того, какой я размер шага на первой итерации выставлю? Но разве это имеет смысл вообще? Ведь размер шага на первой итерации не имеет никакого значения, поскольку регулировка его на старте всё равно изменит. Однако если понимать это в смысле "подтвердите, что при большом размере шага всё хорошо", то это не так, ибо см. выше. Что я неправильно понимаю?

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

1. Пусть у нас схема p-го порядка.

2. Для некоторого шага $\Delta t$:

$y_{A} - y(t_0 + \Delta t) = \alpha (\Delta t)^{p+1} + O((\Delta t)^{p+2})$

3. Делим шаг пополам. Делаем два шага:

$y_{B} - y(t_0 + \Delta t) = 2 \alpha \left(\frac{\Delta t}{2}\right)^{p+1} + O((\Delta t)^{p+2})$

4. Мы используем в качестве меры ошибки (поскольку истинная нам неизвестна) величину:

$\lvert y_{A} - y_{B} \rvert = \alpha \left(\Delta t \right)^{p+1} \left(1 - 2^{-p} \right)$

Вопрос: почему мы можем использовать такую величину? Мы предполагаем, что она $\sim \Delta y_{err}^{true}$? Ну, то есть порядка истинной ошибки, скажем, для случая A (не принципиально, как я понимаю) в этой точке. Или как?

2. Двойной маятник.

Лагранжиан этой системы:

$L = \frac{m_1}{2} \left(l_1 \dot\varphi_1 \right)^2 + \frac{m_2}{2} \left[ \left(l_1 \dot\varphi_1 \right)^2 + \left(l_2 \dot\varphi_2 \right)^2 + 2 l_1 l_2 \dot\varphi_1 \dot\varphi_2 \cos(\varphi_1 - \varphi_2) \right] - \\ - m_1 g l_1 (1 - \cos(\varphi_1)) - m_2 g \left[l_1 (1 - \cos(\varphi_1)) + l_2 (1 - \cos(\varphi_2)) \right]$

Вроде бы не соврал нигде (в любом случае ЛЛ-1, думаю, у всех имеется).

Нужно свести к системе уравнений первого порядка. Причём для этого почему-то настоятельно рекомендуется использоваться канонически сопряженные импульсы $q_i = \frac{\partial L}{\partial \dot\varphi_i}$. У меня в результате получился вот такой монстрик:

$$\begin{cases} \displaystyle \dot \varphi_1 = \frac{\displaystyle\frac{q_1}{l_1^2} - \frac{m_1 + m_2}{m_2} \frac{q_2}{l_1 l_2} \cos(\varphi_1 - \varphi_2)}{m_1 + m_2 \sin^2(\varphi_1 - \varphi_2)} \equiv f_1(\varphi_1, \varphi_2, q_1, q_2)\\ \displaystyle \dot \varphi_2 = \displaystyle \frac{1}{l_2 \left(m_1 + m_2 \sin^2(\varphi_1 - \varphi_2) \right)} \left\lbrace \frac{q_2}{l_2} \frac{m_1}{m_2} - \frac{q_1}{l_1} \right\rbrace \equiv f_2(\varphi_1, \varphi_2, q_1, q_2)\\ \dot q_1 = -m_2 l_1 l_2 f_1(...) f_2(...) \sin(\varphi_1 - \varphi_2) - (m_1 + m_2) g l_1 \sin(\varphi_1)\\ \dot q_2 = m_2 l_1 l_2 f_1(...) f_2(...) \sin(\varphi_1 - \varphi_2) - m_2 g l_2 \sin(\varphi_2)\end{cases} $$

Собственно говоря, вопроса два:

1. Почему так настоятельно рекомендуется использовать в качестве новых переменных канонически сопряженные импульсы? Везде видел, что используют просто $\dot \varphi_i$ и не парятся. Есть ли здесь какой-нибудь сакральный смысл за этим?
2. Может кто-нибудь подтвердить правильность полученного выражения? У меня, увы, ни Maple, ни чего бы то ни было ещё сейчас под рукой нет (а если бы и было что, то, если честно, я бы скорее ошибся там 100 раз).

Заранее спасибо!

(P.S.)

Перенос формулы на другую строку не работает с выключными формулами. Как в этом случае перенести на другую строку?

[Pavia]

Больше так не спрашивайте: 1 вопрос одна тема. Иначе будет каша.

Gickle в сообщении #1165095 писал(а):

Лагранжиан этой системы:

$L = \frac{m_1}{2} \left(l_1 \dot\varphi_1 \right)^2 + \frac{m_2}{2} \left[ \left(l_1 \dot\varphi_1 \right)^2 + \left(l_2 \dot\varphi_2 \right)^2 + 2 l_1 l_2 \dot\varphi_1 \dot\varphi_2 \cos(\varphi_1 - \varphi_2) \right] - \\ - m_1 g l_1 (1 - \cos(\varphi_1)) - m_2 g \left[l_1 (1 - \cos(\varphi_1)) + l_2 (1 - \cos(\varphi_2)) \right]$

Вроде бы не соврал нигде (в любом случае ЛЛ-1, думаю, у всех имеется).

В учебники формула по проще будет.

Обозначения неправильные. Импульс обозначается буковой $p$, а обобщённой координатой $q$.

Gickle в сообщении #1165095 писал(а):

Есть ли здесь какой-нибудь сакральный смысл за этим?

Да, избавится от косинусов в самом начале.

[Gickle]

Pavia в сообщении #1165413 писал(а):

Больше так не спрашивайте: 1 вопрос одна тема. Иначе будет каша.

Не хотел плодить темы. В следующий раз учту совет.

Pavia в сообщении #1165413 писал(а):

В учебники формула по проще будет.

Обозначения неправильные. Импульс обозначается буковой $p$, а обобщённой координатой $q$.

Тут я просто использовал номенклатуру, предложенную автором задачи. И лагранжиан там в таком виде записан, и канонически сопряжённый импульс $q$ обозначен. Разницы от этого всё равно никакой.

Pavia в сообщении #1165413 писал(а):

Gickle в сообщении #1165095 писал(а):

Есть ли здесь какой-нибудь сакральный смысл за этим?

Да, избавится от косинусов в самом начале.

Ну, это ведь скорее вопрос простоты вывода конечной системы уравнений и её результирующей громоздкости?

[Pavia]

Gickle в сообщении #1165565 писал(а):

Тут я просто использовал номенклатуру, предложенную автором задачи. И лагранжиан там в таком виде записан, и канонически сопряжённый импульс $q$ обозначен. Разницы от этого всё равно никакой.

Если все с крыши начнут прыгать и ты будешь? Нестандартные обозначения путают всех: и вас и тех, кто вас читает. Вас перестают понимать и никто такое проверять не будет.

Цитата:

Ну, это ведь скорее вопрос простоты вывода конечной системы уравнений и её результирующей громоздкости?

Так для этого математика и существует что-бы упрощать решения. Иначе-бы все так и считали на палочках как в 1 классе.

[AndrewSu]

Цитата:

Потом говорится, мол, поиграйте с шагом по времени

Цитата:

Правильно ли я понимаю, что от меня хотят, чтобы я показал, что схема устойчива относительно того, какой я размер шага на первой итерации выставлю? Но разве это имеет смысл вообще?

Для оценки устойчивости решения обычно проводят интегрирование два раза с шагом $t$ и $t/2$ и убеждаются в совпадении решений с заданной точностью. Если этого не происходит, то уменьшают шаг и решают еще раз.

Для ускорения вычислений придумали вложенные методы (https://en.wikipedia.org/wiki/List_of_Runge%E2%80%93Kutta_methods#Embedded_methods), которые позволяют контролировать ошибку на каждом шаге.

Ваша схема является частью вложенного метода Heun–Euler.