LEAPFROG метод и Ричард Фейнман
andreyilinskiyLEAPFROG ИНТЕГРИРОВАНИЕ и РИЧАРД ФЕЙНМАН¶
Введение¶
Одной из отличительных особенностей знаменитого курса общей физики "Феймановские лекции по физике" является широкий спектр освещаемых тем.
Так, например, в разделе, посвященном динамическим уравнениям Ньютона есть начала численных методов решения дифференциальных уравнений. Авторы предлагают численно решить задачу линейного маятника и задачу двух тел.
Начало стандартное - для вычислений координаты и скорости предлагается явный метод Эйлера (текст лекций отсюда):
Далее Фейнман предлагает для увеличения точности применить "little tricks" и немного изменить численную схему решения.
Так как скорость маятника меняется, то для шага по координате будет точнее взять скорость в середине интервала, а не в каком либо из его концов. Также и для шага по скорости лучше будет брать ускорение в середине временного шага.
Название метода в лекциях не приводится. Сейчас он называется "leapfrog method". Значения координаты и скорости при применении этого метода, как бы перепрыгивают друг через друга (координату мы получаем в моменты времени i, i+1, i+2, а скорость в i+1/2, i+3/2, i+5/2), что характерно для игры в чехарду. Русский перевод, который мне встретился - метод с перешагиванием.
Численные опыты¶
Интересно рассмотреть предлагаемый метод на нескольких численных экспериментах. Модельной задачей будут колебания нелинейного математического маятника на жестком стержне: \begin{equation*} \begin{cases} \varphi^\prime = \omega\\ \omega^\prime = -g/L\cdot \sin(\varphi) \end{cases} \end{equation*}
Точность¶
Решим задачу разными методами (для сравнения) на различных сетках.
# Функция для решения задачи разными методами на разных сетках
function solveODE(N, r, Tmax, u0, params, method)
Nr = N * r #количество интервалов на сгущенной сетке
h = Tmax / Nr #шаг интегрирования
g, L = params
phi = zeros(Nr + 1)
omega = zeros(Nr + 1)
phi[1], omega[1] = u0
for k in 2: Nr + 1
if method == "Euler"
phi[k] = phi[k-1] + h * omega[k-1]
omega[k] = omega[k-1] - h * g / L * sin( phi[k-1] )
elseif method == "RK2" #Метод Рунге-Кутты 2-ого порядка (метод средней точки)
phiMid = phi[k-1] + h/2 * omega[k-1]
omegaMid = omega[k-1] - h/2 * g / L * sin( phi[k-1] )
phi[k] = phi[k-1] + h * omegaMid
omega[k] = omega[k-1] - h * g / L * sin( phiMid )
elseif method == "leap"
if k == 2
omega[1] = u0[2] - h/2 * g / L * sin( u0[1] ) #для старта leapfrog метода необходимо найти скорость в центре первого шага интегрирования. Для этого делаем полушаг методом Эйлера
end
phi[k] = phi[k-1] + h * omega[k-1]
omega[k] = omega[k-1] - h * g / L * sin( phi[k] ) #эти значения соответствуют центру интервала от k-ого до k+1 узла
elseif method == "leap2" #переписанный leapfrog метод так, чтобы иметь скорость и угол в один и тот же момент времени
phi[k] = phi[k-1] + h * omega[k-1] - h^2/2 * g/L * sin( phi[k-1] )
omega[k] = omega[k-1] - h/2 * g / L * ( sin( phi[k-1] ) + sin( phi[k] ))
end
end
return phi[1: r: end], omega[1: r: end]
end;
Посмотрим что получилось. Построим зависимость угла от времени при решении задачи методом Рунге-Кутта второго порядка точности.
N = 200 #кол-во интервалов
r = 2 #коэф. сгущения сетки
Tmax = 100 #длительность симуляции
L = 9 #длина маятника
const g = 9.81 #ускорение свободного падения
# Начальные условия. Пусть угол будет очень большим, чтобы увидеть нелинейность колебаний
u0 = (0.9pi, 0)
phi, omega = solveODE(550, 1, Tmax, u0, (g, L), "RK2")
t = 0: Tmax/550: Tmax
plot(t, phi,
title = "Колебания нелинейного маятника",
xlabel = "Время, с",
ylabel = "Угол, рад",
legend = false,
marker = :circle,
markersize = 1)
Видно, что вершины графиков более плоские, чем у обычной синусоиды. Это говорит о нелинейности колебаний. Маятник почти "зависает" при углах близких к $\pi$.
По мере численного счета накапливается погрешность и в какой-то момент маятник начинает двигаться по кругу (маятник "набрал энергии" за счет погрешностей расчета).
Посмотрим как меняется погрешность в зависимости от величины шага при вычислениях разными методами. Оценивать погрешность будем по формуле Рунге-Ромберга.
ns = 18 #количество итераций сгущения сетки
solves = zeros(ns, N+1)
err = zeros(2, ns-1)
t = 0: Tmax/N: Tmax
methods = ("Euler", "RK2", "leap") #будем сравинвать эти три метода
p = plot()
for method in methods
for s = 0: ns-1
phi, omega = solveODE(N, r^s, Tmax, u0, (g, L), method)
solves[s+1, :] = phi #записываем приближенные решения на разных сетках в матрицу
if s > 0
err[1, s] = N*r^s
err[2, s] = norm(solves[s+1, :] - solves[s, :]) #норма разности между решениями на всё более сгущающейся сетке - оценка глобальной погрешности метода
end
end
plot!(err[1, :], err[2, :],
label = method,
title = "Погрешность в зависимости от количества интервалов",
xaxis = :log,
yaxis = :log,
lw=3,
marker = :circle,
markersize = 5,
xlims = (1e2, 1e8),
ylims = (1e-10, 1e4),
legend = :bottomright,
xlabel = "Количество интервалов",
ylabel = "Глобальная погрешность, рад")
end
display(p)
Комментарии к графику¶
- На линейных участках графики имеют разный наклон. Метод Эйлера имеет меньший наклон чем метод leapfrog и метод средней точки;
- Наклон графика показывает как быстро уменьшается погрешность при увеличении количества шагов интегрирования (уменьшении шага интегрирования);
- Погрешность по методу Эйлера при увеличении количества шагов в 2 раза уменьшается тоже в 2 раза (первый порядок точности);
- Погрешность по двум другим методам при увеличении количества шагов в 2 раза уменьшается в 4 раза. Это означает что предложенный leapfrog метод имеет второй порядок точности.
Метод leapfrog более точный. Для повышения точности расчетов он и вводился в лекциях.
Сохранение энергии¶
Системы, рассматриваемые в лекциях, имеют одну особенность - они консервативные, а значит в них действует закон сохранения энергии. Посчитаем энергию для данных, полученных методом RK2 и leapfrog. К сожалению, так как метод leapfrog в своей оригинальной записи не дает значений угла и скорости в один и тот же момент времени, его нельзя в лоб применить для расчета энергии системы. Но его можно переписать таким образом чтобы получить нужные данные в одни и те же моменты времени (в коде он назван как leap2).
phiRK2, omegaRK2 = solveODE(550, 1, Tmax, u0, (g, L), "RK2")
phiLeap, omegaLeap = solveODE(550, 1, Tmax, u0, (g, L), "leap2")
t = 0: Tmax/550: Tmax
energyRK2 = L^2 * omegaRK2.^2 / 2 - g*L*cos.(phiRK2)
energyLeap = L^2 * omegaLeap.^2 / 2 - g*L*cos.(phiLeap)
plot(t, energyRK2, lw = 3, label = "RK2")
plot!(t, energyLeap, label = "leapfrog",
title = "Энергия при разных числ. методах",
xlabel = "Время, с",
ylabel = "Энергия системы, Дж",
lw = 3)
Как видно энергия ведет себя различным образом, хотя решения получены методами одинакого порядка точности. И это свойство leapfrog метода в контексте выбранной задачи гораздо важнее и интереснее чем просто повышение точности.
Энергия колеблется оносительно истинного значения, но не смещается ни в плюс ни в минус. Такие интеграторы называются симплектическими и позволяют производить вычисления на больших временных промежутках. Для консервативных систем это чрезвычайно важно и, наверняка, появление именно этой численной схемы в лекциях - не случайно.
На анимации ниже видно, как фиолетовый шарик набирает энергию и начинает делать круги, в то время как синий (посчитанный с помощью leapfrog метода) выполняет корректные колебания раз за разом.