Тест 3-го порядка против симплектического интегратора 4-го порядка со странным результатом

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

Тогда я подумал, что было бы неплохо продемонстрировать влияние разных временных шагов на глобальную точность методов с разными порядками, и я написал и запустил скрипт Python / Pylab для этого. Для сравнения я выбрал:

( leap2 ) Пример 2-го порядка Википедии, с которым я знаком, хотя я знаю его под именем leapfrog ,
( ruth3 ) Симплектический интегратор Рут 3-го порядка ,
( ruth4 ) Симплектический интегратор 4-го порядка Рут .

Странно то, что какой бы временной шаг я ни выбрал, метод Рут 3-го порядка кажется более точным в моем тесте, чем метод Рут 4-го порядка, даже на порядок.

Поэтому мой вопрос: что я здесь делаю не так? Подробности ниже.

Методы раскрывают свою силу в системах с разделимыми гамильтонианами, то есть теми, которые можно записать как

H (q, p) = T (p) + V (q)

$H(q,p) = T(p) + V(q)$ где

q

$q$ содержит все координаты положения,

p

$p$ содержит сопряженные импульсы,

T

$T$ представляет кинетические энергия и

V

$V$ потенциальная энергия.

В нашей настройке мы можем нормализовать силы и импульсы по массам, к которым они применяются. Таким образом, силы превращаются в ускорения, а импульсы превращаются в скорости.

Симплектические интеграторы поставляются со специальными (заданными, постоянными) коэффициентами, которые я обозначу $a_1,\ldots,a_n$ и $b_1,\ldots,b_n$ . С этими коэффициентами один шаг для развития системы от времени $t$ до времени $t+\delta t$ принимает вид

Для $i=1,\ldots,n$ :
1. Вычислить вектор $g$ всех ускорений, учитывая вектор $q$ всех положений
2. Изменить вектор $v$ всех скоростей на $b_i\,g\,\delta t$
3. Изменить вектор $q$ всех позиций на $a_i\,v\,\delta t$

Мудрость теперь заключается в коэффициентах. Это

\begin{aligned} [\begin{matrix} a_{1} & a_{2} \\ b_{1} & b_{2} \end{matrix}] & = [\begin{matrix} \frac{1}{2} & \frac{1}{2} \\ 0 & 1 \end{matrix}] & (leap2) \\ [\begin{matrix} a_{1} & a_{2} & a_{3} \\ b_{1} & b_{2} & b_{3} \end{matrix}] & = [\begin{matrix} \frac{2}{3} & - \frac{2}{3} & 1 \\ \frac{7}{24} & \frac{3}{4} & - \frac{1}{24} \end{matrix}] & (ruth3) \\ [\begin{matrix} a_{1} & a_{2} & a_{3} & a_{4} \\ b_{1} & b_{2} & b_{3} & b_{4} \end{matrix}] & = \frac{1}{2 - \sqrt[3]{2}} [\begin{matrix} \frac{1}{2} & \frac{1 - \sqrt[3]{2}}{2} & \frac{1 - \sqrt[3]{2}}{2} & \frac{1}{2} \\ 0 & 1 & - \sqrt[3]{2} & 1 \end{matrix}] & (ruth4) \end{aligned}

$\begin{align} \begin{bmatrix} a_1 & a_2 \\ b_1 & b_2 \end{bmatrix} &= \begin{bmatrix} \frac{1}{2} & \frac{1}{2} \\ 0 & 1 \end{bmatrix} &&\textsf{(leap2)} \\ \begin{bmatrix} a_1 & a_2 & a_3 \\ b_1 & b_2 & b_3 \end{bmatrix} &= \begin{bmatrix} \frac{2}{3} & -\frac{2}{3} & 1 \\ \frac{7}{24} & \frac{3}{4} & -\frac{1}{24} \end{bmatrix} &&\textsf{(ruth3)} \\ \begin{bmatrix} a_1 & a_2 & a_3 & a_4 \\ b_1 & b_2 & b_3 & b_4 \end{bmatrix} &= \frac{1}{2-\sqrt[3]{2}} \begin{bmatrix} \frac{1}{2} & \frac{1-\sqrt[3]{2}}{2} & \frac{1-\sqrt[3]{2}}{2} & \frac{1}{2} \\ 0 & 1 & -\sqrt[3]{2} & 1 \end{bmatrix} &&\textsf{(ruth4)} \end{align}$

\begin{aligned} y^{″} + y & = 0 & y (0) & = 1 & y^{'} (0) & = 0 \end{aligned}

$\begin{align} y''+y &= 0 & y(0) &= 1 & y'(0) &= 0 \end{align}$

∴ (y (t), y^{'} (t)) = (\cos t, - \sin t)

$\therefore\quad \left(y(t),y'(t)\right) = (\cos t, -\sin t)$

(q, v)

$(q,v)$

(y, y^{'})

$(y,y')$

Я интегрировал IVP с вышеупомянутыми методами над с размером шага с целым числом выбранным где-то между и . Принимая во внимание скорость прыжка 2 , я утроил для этого метода. Затем я нанес на график полученные кривые в фазовом пространстве и увеличил масштаб в где в идеале кривые должны снова прийти к . $t\in[0,2\pi]$ $\delta t=\frac{2\pi}{N}$ $N$ $10$ $100$ $N$ $(1,0)$ $t=2\pi$

Вот графики и увеличения для и : $N=12$ $N=36$

Для , leap2 с размером шага случается , чтобы прибыть ближе к дому , чем ruth4 с размером шага . Для , ruth4 победы над leap2 . Однако ruth3 с тем же размером шага, что и ruth4 , прибывает гораздо ближе к дому, чем остальные, во всех настройках, которые я тестировал до сих пор. $N=12$ $\frac{2\pi}{3N}$ $\frac{2\pi}{N}$ $N=36$

Вот скрипт Python / Pylab:

import numpy as np
import matplotlib.pyplot as plt

def symplectic_integrate_step(qvt0, accel, dt, coeffs):
    q,v,t = qvt0
    for ai,bi in coeffs.T:
        v += bi * accel(q,v,t) * dt
        q += ai * v * dt
        t += ai * dt
    return q,v,t

def symplectic_integrate(qvt0, accel, t, coeffs):
    q = np.empty_like(t)
    v = np.empty_like(t)
    qvt = qvt0
    q[0] = qvt[0]
    v[0] = qvt[1]
    for i in xrange(1, len(t)):
        qvt = symplectic_integrate_step(qvt, accel, t[i]-t[i-1], coeffs)
        q[i] = qvt[0]
        v[i] = qvt[1]
    return q,v

c = np.math.pow(2.0, 1.0/3.0)
ruth4 = np.array([[0.5, 0.5*(1.0-c), 0.5*(1.0-c), 0.5],
                  [0.0,         1.0,          -c, 1.0]]) / (2.0 - c)
ruth3 = np.array([[2.0/3.0, -2.0/3.0, 1.0], [7.0/24.0, 0.75, -1.0/24.0]])
leap2 = np.array([[0.5, 0.5], [0.0, 1.0]])

accel = lambda q,v,t: -q
qvt0 = (1.0, 0.0, 0.0)
tmax = 2.0 * np.math.pi
N = 36

fig, ax = plt.subplots(1, figsize=(6, 6))
ax.axis([-1.3, 1.3, -1.3, 1.3])
ax.set_aspect('equal')
ax.set_title(r"Phase plot $(y(t),y'(t))$")
ax.grid(True)
t = np.linspace(0.0, tmax, 3*N+1)
q,v = symplectic_integrate(qvt0, accel, t, leap2)
ax.plot(q, v, label='leap2 (%d steps)' % (3*N), color='black')
t = np.linspace(0.0, tmax, N+1)
q,v = symplectic_integrate(qvt0, accel, t, ruth3)
ax.plot(q, v, label='ruth3 (%d steps)' % N, color='red')
q,v = symplectic_integrate(qvt0, accel, t, ruth4)
ax.plot(q, v, label='ruth4 (%d steps)' % N, color='blue')
ax.legend(loc='center')
fig.show()

Я уже проверил на простые ошибки:

Нет опечатка в Википедии. Я проверил ссылки, в частности ( 1 , 2 , 3 ).
Я правильно понял последовательность коэффициентов. Если вы сравните с порядком в Википедии, обратите внимание, что последовательность операций оператора работает справа налево. Моя нумерация совпадает с Candy / Rozmus . И если я все же попробую другой порядок, результаты ухудшатся.

Мои подозрения:

Неправильный порядок шага: может быть, схема 3-го порядка Рут имеет несколько меньшие подразумеваемые константы, и если бы размер шага был сделан очень маленьким, то метод 4-го порядка победил бы? Но я даже попробовал , и метод 3-го порядка все еще лучше. $N=360$
Неправильный тест: Что-то особенное в моем тесте позволяет методу Рут третьего порядка вести себя как метод высшего порядка?

python time-integration ccorn
источник

Можете ли вы дать числовые значения ошибок? Сложно сказать по сюжету. Как ошибки масштабируются при изменении ? Масштабируются ли они, как ожидается, от заказов методов? Обычно для проверки этого на графике журнала регистрируются ошибки против

N

$N$

N

$N$

Кирилл

@ Кирилл: Работаем над этим. Скоро отредактирую.

кукуруза

Одна вещь, в которой я подозреваю, это выбор линейных прав: ошибки усечения методов обычно зависят от некоторой старшей производной от прав, поэтому, если все старшие производные от правых чисел исчезают, вы можете наблюдать странное поведение сходимости. Вероятно, стоит попробовать более необычные rhs.

Кирилл

Ответы:

Следуя совету Кирилла , я выполнил тест с из списка приблизительно геометрически растущих значений, и для каждого вычислил ошибку как где представляет приближение получены путем численного интегрирования. Вот результат на графике log-log: $N$ $N$

ϵ (N) = ‖ \tilde{z} (2 π) - \tilde{z} (0) ‖_{2} where \tilde{z} (t) = (\tilde{y} (t), {\tilde{y}}^{'} (t))

$\epsilon(N) = \|\tilde z(2\pi) - \tilde z(0)\|_2 \quad\text{where}\quad \tilde z(t) = \left(\tilde{y}(t),\tilde y'(t)\right)$

\tilde{z}

$\tilde z$

Таким образом, ruth3 действительно имеет тот же порядок что и ruth4 для этого теста, и подразумевает, что константы имеют только величину. $4$ $\frac{1}{100}$

Интересно. Мне придется продолжить расследование, возможно, попробовать другие тесты.

Кстати, вот дополнения к скрипту Python для графика ошибок:

def int_error(qvt0, accel, qvt1, Ns, coeffs):
    e = np.empty((len(Ns),))
    for i,N in enumerate(Ns):
        t = np.linspace(qvt0[2], qvt1[2], N+1)
        q,v = symplectic_integrate(qvt0, accel, t, coeffs)
        e[i] = np.math.sqrt((q[-1]-qvt1[0])**2+(v[-1]-qvt1[1])**2)
    return e

qvt1 = (1.0, 0.0, tmax)
Ns = [12,16,20,24,32,40,48,64,80,96,128,160,192,
      256,320,384,512,640,768,1024,1280,1536,2048,2560,3072]

fig, ax = plt.subplots(1)
ax.set_xscale('log')
ax.set_xlabel(r"$N$")
ax.set_yscale('log')
ax.set_ylabel(r"$\|z(2\pi)-z(0)\|$")
ax.set_title(r"Error after 1 period vs #steps")
ax.grid(True)
e = int_error(qvt0, accel, qvt1, Ns, leap2)
ax.plot(Ns, e, label='leap2', color='black')
e = int_error(qvt0, accel, qvt1, Ns, ruth3)
ax.plot(Ns, e, label='ruth3', color='red')
e = int_error(qvt0, accel, qvt1, Ns, ruth4)
ax.plot(Ns, e, label='ruth4', color='blue')
ax.legend(loc='upper right')
fig.show()

ccorn
источник

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

Кирилл

@Kirill: Это является ответом. У ruth3 действительно есть более высокий порядок и меньшие постоянные здесь. Обнаружен из-за вашего предложения составить график ошибок журнала. Поэтому на вопрос дан ответ, и я решительно не буду менять суть вопроса после того, как на него дан ответ, даже если ответ был составлен мной.

кукуруза

Тем не менее, я был бы рад принять дальнейший анализ. (Ответы на вопросы принимаются автоматически, но я полагаю, что это можно изменить.)

ccorn

Я немного посмотрел на это и не смог найти убедительного объяснения. Эта сходимость 4-го порядка по ruth3 имеет много общего с начальными условиями: попробуйте установить и не интегрировать в течение полного периода (ни полупериода). Одна вещь, которая может легко произойти, состоит в том, что ошибка усечения имеет некоторый компонент с «нулевым средним», который компенсируется при интеграции в полный период. Я также попробовал чтобы проверить, имеет ли это какое-либо отношение к имеющему мало высоких производных, но в моих тестах кажется, что начальные условия и периодичность больше связаны с этим.

p_{0} \neq 0

$p_0\neq0$

V (q) = 1 / q + \log q

$V(q)=1/q+\log q$

V

$V$

Кирилл

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

Крис

Построение ошибки , по всему интервалу, масштабированное по степени размера шага, заданного ожидаемым порядком, дает графики $\ddot q=-q$ $q(0)=1,\dot q(0)=0$

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

В методе Рут 3-го порядка старший коэффициент компонента представляется равным нулю, видимая предельная кривая близка или равна горизонтальной оси. Видимые графики четко показывают преобладание члена ошибки 4-го порядка. Масштабирование для ошибки 4-го порядка дает график, аналогичный другим. $p$

Как можно видеть, во всех 3 случаях коэффициент ошибки старшего порядка в компоненте равен нулю при после полного периода. При объединении ошибок обоих компонентов поведение компонента таким образом, доминирует, ложно создавая впечатление метода 4-го порядка на графике журнала. $q$ $t=2\pi$ $p$

Максимальный коэффициент в компоненте может быть найден около . График журнала должен отображать правильные глобальные порядки ошибок. $q$ $3\pi/2$

То, что исчезновение члена ошибки 3-й степени в Ruth3p является артефактом простоты линейного уравнения, показывает нелинейный пример , с соответствующими участками $\ddot q=-\sin(q)$ $q(0)=1.3,~\dot q(0)=0$

Лутц Леманн
источник