Метод Рунге-Кутта для второго закона Ньютона

[Ulrih]

Второй закон Ньютона представляет собой дифференциальное уравнение второго порядка:

$$ \ddot x = - \frac{1}{m} \partial_x U(x) $$

Которую можно представить системой ДУ первого порядка:

$$ \dot P = - \frac{1}{m} \partial_x U(x) $$

$$ P = \dot Q $$

Численное решение в момент времени $t+dt$ можно найти, используя метод Рунге-Кутта:

$$ K_1 &= \begin{pmatrix} \dot{x}(t=0) \\ -\frac{1}{m} {\displaystyle \partial_x U[x(t=0)]} \end{pmatrix} = \begin{pmatrix} K_1^{(1)} \\ K_1^{(2)} \end{pmatrix} $$

$$ K_2 &= \begin{pmatrix} {\displaystyle \dot{x}(t=0)+K_1^{(2)}\frac{dt}{2} } \\ -\frac{1}{m} {\displaystyle \partial_x U\left[ x_0+K_1^{(1)}\frac{dt}{2}\right]}\end{pmatrix} = \begin{pmatrix} K_2^{(1)} \\ K_2^{(2)} \end{pmatrix} $$

$$ K_3 &= \begin{pmatrix} {\displaystyle \dot{x}(t=0)+K_2^{(2)}\frac{dt}{2} } \\ -\frac{1}{m} {\displaystyle \partial_x U\left[ x_0+K_2^{(1)}\frac{dt}{2}\right]} \end{pmatrix} = \begin{pmatrix} K_3^{(1)} \\ K_3^{(2)} \end{pmatrix} $$

$$ K_4 &= \begin{pmatrix} {\displaystyle \dot{x}(t=0)+K_4^{(2)} dt }\\ -\frac{1}{m} {\displaystyle \partial_x U\left[ x_0+K_3^{(1)} dt\right]} \end{pmatrix} = \begin{pmatrix} K_4^{(1)} \\ K_4^{(2)} \end{pmatrix} $$

Решением будет:

$$ Q(t+dt)=Q(t)+\left(K_1^{(1)}+2K_2^{(1)}+2K_3^{(1)}+K_4^{(1)} \right)dt/6 $$

$$ P(t+dt)=P(t)+\left(K_1^{(2)}+2K_2^{(2)}+2K_3^{(2)}+K_4^{(2)} \right)dt/6 $$

Правильно ли я понимаю, как работает этот метод?

[Yu_K]

http://alglib.sources.ru/diffequations/ ... ttasys.php

не нравится ваше $t=0$ в формулах - это только для начальной точки верно, ну и если перешли к $P$ и $Q$ - то перейдите полностью, зачем $x$ оставили в уравнениях.

[Ulrih]

Yu_K

Такую ссылку я тоже видел. Помогает она мало, т.к. я пишу программы в Fortran`e, а до С и Делфи я еще не дошел.
Формуды такие у меня тоже есть. Мне не понятно как именно здесь, в этом случае написать РК4 алгоритм... Интегрируем по времени, в правой части времени явно нет, правая часть от $y$ не зависит...

[Yu_K]

Лучше переписать так уравнение

$$ \ddot x = F(x) $$

и тогда его можно представить системой ДУ первого порядка:

$$ \dot x= V $$

$$ \dot V = F(x) $$

Тогда численное решение в момент времени $t+dt$ можно найти, используя метод Рунге-Кутта:

$$ K_1 &= \begin{pmatrix} \ V(t) \\ {\displaystyle F[x(t)]} \end{pmatrix} = \begin{pmatrix} K_1^{(1)} \\ K_1^{(2)} \end{pmatrix} $$

$$ K_2 &= \begin{pmatrix} \ V(t)+K_1^{(1)}/2 \\ {\displaystyle F[x(t)+K_1^{(2)}/2] } \end{pmatrix} = \begin{pmatrix} K_2^{(1)} \\ K_2^{(2)} \end{pmatrix} $$

и т.д.

Решением будет:

$$ x(t+dt)=x(t)+\left(K_1^{(1)}+2K_2^{(1)}+2K_3^{(1)}+K_4^{(1)} \right)dt/6 $$

$$ V(t+dt)=V(t)+\left(K_1^{(2)}+2K_2^{(2)}+2K_3^{(2)}+K_4^{(2)} \right)dt/6 $$

Нет никакой разницы на чем Вы программируете. Техническая часть формул не меняется - берутся значения правых частей в скорректированных точках по времени.

[Ulrih]

А как же размерности величин?

$$ [K_1^{(2)}] = L T^{-2} $$

в вектор столбце $K_2$ стоит

$$ x(t) + K_1^{(2)} $$

то есть метр + метр за секунду в квадрате...
вот что меня смущает...

[Yu_K]

Поторопился немного... показалось, что можно

$$dt$$

вынести за скобку, как в Вашем первом посте.

$$ K_1 &= \begin{pmatrix} \ V(t) dt\\ {\displaystyle F[x(t)]dt} \end{pmatrix} = \begin{pmatrix} K_1^{(1)} \\ K_1^{(2)} \end{pmatrix} $$

$$ K_2 &= \begin{pmatrix} \ [V(t)+K_1^{(2)}/2] dt\\ {\displaystyle [F[x(t)+K_1^{(1)}/2]]dt } \end{pmatrix} = \begin{pmatrix} K_2^{(1)} \\ K_2^{(2)} \end{pmatrix} $$

и т.д.

Решением будет:

$$ x(t+dt)=x(t)+\left(K_1^{(1)}+2K_2^{(1)}+2K_3^{(1)}+K_4^{(1)} \right)/6 $$

$$ V(t+dt)=V(t)+\left(K_1^{(2)}+2K_2^{(2)}+2K_3^{(2)}+K_4^{(2)} \right)/6 $$

Посмотрите правую часть как векторную функцию двух переменных

$$G(x,V)$$

- и делайте все как в формулах приведенных в http://alglib.sources.ru/diffequations/rungekuttasys.php - независимой переменной

$$t$$

в правой части не будет.

[Ulrih]

То есть получается, что коэффициенты $K_i$ зависят друг от друга.
Так, например, в $K_2^{1}$ входит $K_1^{2}$.
Это получается из-за того, что мы смотрим на это как на систему ДУ???

[ewert]

Yu_K писал(а):

метод Рунге-Кутта:

$$ K_1 &= \begin{pmatrix} \ V(t) \\ {\displaystyle F[x(t)]} \end{pmatrix} = \begin{pmatrix} K_1^{(1)} \\ K_1^{(2)} \end{pmatrix} $$

$$ K_2 &= \begin{pmatrix} \ V(t)+K_1^{(1)}/2 \\ {\displaystyle F[x(t)+K_1^{(2)}/2] } \end{pmatrix} = \begin{pmatrix} K_2^{(1)} \\ K_2^{(2)} \end{pmatrix} $$

и т.д.

$$ K_1 &= \begin{pmatrix} \ V(t) \\ {\displaystyle m^{-1}F[x(t)]} \end{pmatrix} = \begin{pmatrix} K_1^{(1)} \\ K_1^{(2)} \end{pmatrix} $$

$$ K_2 &= \begin{pmatrix} \ V(t)+K_1^{(2)}\cdot dt/2 \\ {\displaystyle m^{-1}F[x(t)+K_1^{(1)}\cdot dt/2] } \end{pmatrix} = \begin{pmatrix} K_2^{(1)} \\ K_2^{(2)} \end{pmatrix} $$

и т.д.

Ulrih писал(а):

То есть получается, что коэффициенты $K_i$ зависят друг от друга.
Так, например, в $K_2^{1}$ входит $K_1^{2}$.
Это получается из-за того, что мы смотрим на это как на систему ДУ???

Не "зависят", а получаются друг из друга рекуррентно. Это -- основная идея методов Рунге-Кутта, безразлично, идёт ли речь о системе или об одном скалярном уравнении.

[Ulrih]

При одном скалярном уравнение рекуррентность коэффициентов лежит на поверхности.
А вот с системой ДУ уже сложнее.
ewert

Цитата:

$$ K_1 &= \begin{pmatrix} \ V(t) \\ {\displaystyle m^{-1}F[x(t)]} \end{pmatrix} = \begin{pmatrix} K_1^{(1)} \\ K_1^{(2)} \end{pmatrix} $$

$$ K_2 &= \begin{pmatrix} \ V(t)+K_1^{(2)}\cdot dt/2 \\ {\displaystyle m^{-1}F[x(t)+K_1^{(1)}\cdot dt/2] } \end{pmatrix} = \begin{pmatrix} K_2^{(1)} \\ K_2^{(2)} \end{pmatrix} $$

и т.д.

вот это уже правильно?

[ewert]

Ничего сложного. Как сказал Yu_K,

Ulrih в сообщении #211119 писал(а):

При одном скалярном уравнение рекуррентность коэффициентов лежит на поверхности.
А вот с системой ДУ уже сложнее.

Ничего сложного. Как сказал Yu_K, те же соотношения, только векторные. Те же яйца,только в профиль.

Ulrih в сообщении #211119 писал(а):

вот это уже правильно?

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

[Yu_K]

Тольто, если в исходных уравнениях $m$ входило внурь функции $F$, почему оно оттуда выбралось вдруг?

http://i037.radikal.ru/0905/a5/efcfc5be6322.gif - вот как выглядят рекурсивные формулы для одного уравнения (см картинку в тексте ниже).

Вопрос - как можно просто получить разложение в точке $t+h$, и показать что там действительно 4-ый порядок аппроксимации? Техника получения формул описана у Березина Жидкова во втором томе. Однако вроде если взять уже готовые формулы для метода РК, то можно разложить в ряд с учетом вида уравнения и показать, что члены при соответствующих степенях $h$ одинаковы до $h^4$, но тогда рекурсивные производные придется вычислять...слишком громоздко...




http://s50.radikal.ru/i127/0905/37/e6633777c801.gif - а вот так выглядят конкретно рекурсивные формулы для системы подобной Вашей. Но понятно, что явно ими никто не пользуется - хотя для анализа они могут быть интересны.

[Ulrih]

Всем

Дорогие товарищи!
Я понял, что я делал не так.
Как обманывает нас подсознание, чем меньше $dt$, тем точнее результат. Это если мы смотрим в пределе и на бумажке.
Как только мы переходим к программированию численных методов, это нарушается до некоторой степени. Я прочел в хорошей книжке (Ильина В.А., Силаев П.К. - Численные методы для физиков-теоретиков (часть 2) 2004) надо быть предельно внимательным с условиями для шага по сетке.
Рунке-Кутта 4 порядка я наконец рассписал. Всем спасибо за это!

Теперь думаю о 5-м порядке и адаптивном шаге по сетке. Мне это кажется целесообразным, так как правая часть у меня с особенностями.

И подскажите на последок, правильно ли я понимаю Рунке-Кутта 5 порядка рекккурентные формулы для коэффициентов.

$$ K_1 = \begin{pmatrix} V(t) \\ F[x(t)] \end{pmatrix} = \begin{pmatrix} K_1^1 \\ K_1^2 \end{pmatrix} $$

$$ K_2 = \begin{pmatrix} V(t) + K_1^2 dt a_2 \\ F[x(t) + K_1^1 dt a_2] \end{pmatrix} =\begin{pmatrix} K_2^1 \\ K_2^2 \end{pmatrix} $$

$$ K_3 = \begin{pmatrix} V(t) + K_1^2 dt (a_3 - b_3) + K_2^2 dt b_3 \\ F[x(t) + K_1^1 dt (a_3 - b_3) + K_2^1 dt b_3] \end{pmatrix} = \begin{pmatrix} K_3^1 \\ K_3^2 \end{pmatrix} $$

и так далее... вплоть до $K_6$
где $a_i\,, b_i$ константы.

[Yu_K]

Хотелось бы увидеть пример, когда точности метода РК не хватает. И что значит не хватает - есть сравнение с точным решением? Может это связано с жесткостью системы? Тогда имеет смысл применять методы типа Розенброка. Если Вы работаете на фортране - то для жестких систем есть хороший фортрановский пакет Gear, на нем в свое время рассчитали полный набор реакций горения водорода и кислорода.

[Ulrih]

Yu_K
Нет. Сравнения с точным решением у меня нету. Маячком является закон сохранения энергии.

У меня есть 12 ДУ вида

$$ \ddot x_1 = -\frac{1}{m_1} \partial_x_1 U(x_1, y_1, z_1 ... z_4) $$

$$ ... $$

$$ \ddot z_4 = -\frac{1}{m_4} \partial_z_4 U(x_1, y_1, z_1... z_4) $$

То есть три частицы находятся в поле $ U(x_1, y_1, z_1... z_4)$. Массы частиц таковы: $m_1=m_2=m_3\,, m_4/m_1 \sim 14 $
Когда программа ищет решение по методу Рунге-Кутта 4 порядка, наблюдается очень сильная зависимость от величины $dt$ --- от шага интегрирования по времени; полная энергия очень быстро расходится (неустойчивое решение).
Сейчас считаю с достаточно малым шагом $dt$, но все равно имеет место отклонение от начального значения полной энергии... До этого я списывал не слишком большое отклонение от полной энергии на накопление локальной ошибки и как следствие большой величины глобальной на интервале $[t_1,t_2]$. Но даже с гораздо большим шагом --- все плохо.

Стоит сказать о правой части. Выражение $\partial_{\alpha} U(x_1, y_1, z_1... z_4) $ имеет разрыв при некоторой конфигурации частиц. То есть, функция $U$ не является гладкой функцией.

Вот код для РК4 из фортрана, может я какую ошибку проглядел.

do j = 1,4
do i = 1,3
   if (j/=4) then
   mass = mass1
      else
   mass = mass4
   end if

XXX(i,j,1) = veloc(i,j)*dt   ! K_1^1
x = ccoord
grad = derivate_of_direction(x,veloc) ;   du(i,j) = grad(i,j,numPEs)
VVV(i,j,1) = - du(i,j)/mass*dt   ! K_1^2

XXX(i,j,2) = (veloc(i,j) + VVV(i,j,1)/2.0d0)*dt   ! K_2^1
x = ccoord
x(i,j) = ccoord(i,j)+XXX(i,j,1)/2.0d0   ! x(t) + K_1^1/2
grad = derivate_of_direction(x,veloc) ;   du(i,j) = grad(i,j,numPEs)
VVV(i,j,2) = - du(i,j)/mass*dt   ! K_2^2

XXX(i,j,3) = (veloc(i,j) + VVV(i,j,2)/2.0d0)*dt   ! K_3^1
x = ccoord
x(i,j) = ccoord(i,j)+XXX(i,j,2)/2.0d0   ! x(t) + K_2^1/2
grad = derivate_of_direction(x,veloc) ;   du(i,j) = grad(i,j,numPEs)
VVV(i,j,3) = - du(i,j)/mass*dt   ! K_3^2

XXX(i,j,4) = (veloc(i,j) +VVV(i,j,3))*dt
x = ccoord
x(i,j) = ccoord(i,j)+XXX(i,j,3)*dt   ! x(t) + K_3^2
grad = derivate_of_direction(x,veloc) ;   du(i,j) = grad(i,j,numPEs)
VVV(i,j,4) = - du(i,j)/mass*dt

! \delta^j = K_1^j + 2*K_2^j + 2*K_3^j + K_4^j
! x(t+dt) = x(t) + \delta^1 / 6.0d0
! V(t+dt) = V(t) + \delta^2 / 6.0d0
RK4(i,j,1) = (XXX(i,j,1) + 2.0d0*XXX(i,j,2) + 2.0d0*XXX(i,j,3) + XXX(i,j,4))/6.0d0 ! coord
RK4(i,j,2) = (VVV(i,j,1) + 2.0d0*VVV(i,j,2) + 2.0d0*VVV(i,j,3) + VVV(i,j,4))/6.0d0 ! veloc

end do   !for i
end do   !for j

ccoord = ccoord + RK4(:,:,1)
veloc = veloc + RK4(:,:,2)

[Yu_K]

Ulrih писал(а):

Когда программа ищет решение по методу Рунге-Кутта 4 порядка, наблюдается очень сильная зависимость от величины $dt$ --- от шага интегрирования по времени; полная энергия очень быстро расходится (неустойчивое решение).

Однако бред какой-то. Как это энергия расходится? И что значит зависимость от шага, если есть два разных расчета с разными шагами - то они не накладываются друг на друга? А какое поле, какое взаимодействие между частицами?

Возьмите решение задачи трех тел с гравитационным взаимодействием - точное и проверьте свою программу. Например такое решение http://www.youtube.com/watch?gl=RU&v=NliTsLmkDKo. Или четыре гравитационных одинаковых тела поставьте в углы квадрата и закрутите симметрично вокруг центра масс.