Задача Кеплера
Гамильтониан и момент импульса для задачи Кеплера имеют следующий вид:
Кроме того, известно следующее:
using OrdinaryDiffEq, LinearAlgebra, ForwardDiff, Plots;
gr();
H(q, p) = norm(p)^2 / 2 - inv(norm(q))
L(q, p) = q[1] * p[2] - p[1] * q[2]
pdot(dp, p, q, params, t) = ForwardDiff.gradient!(dp, q -> -H(q, p), q)
qdot(dq, p, q, params, t) = ForwardDiff.gradient!(dq, p -> H(q, p), p)
initial_position = [0.4, 0]
initial_velocity = [0.0, 2.0]
initial_cond = (initial_position, initial_velocity)
initial_first_integrals = (H(initial_cond...), L(initial_cond...))
tspan = (0, 20.0)
prob = DynamicalODEProblem(pdot, qdot, initial_velocity, initial_position, tspan)
sol = solve(prob, KahanLi6(), dt = 1 // 10);retcode: Success
Interpolation: 3rd order Hermite
t: 201-element Vector{Float64}:
0.0
0.1
0.2
0.30000000000000004
0.4
0.5
0.6
0.7
0.7999999999999999
0.8999999999999999
⋮
19.200000000000003
19.300000000000004
19.400000000000006
19.500000000000007
19.60000000000001
19.70000000000001
19.80000000000001
19.900000000000013
20.0
u: 201-element Vector{ArrayPartition{Float64, Tuple{Vector{Float64}, Vector{Float64}}}}:
([0.0, 2.0], [0.4, 0.0])
([-0.5830949354540153, 1.8556656829703986], [0.36982713146498514, 0.19503596514776078])
([-0.9788105843777312, 1.5274462532150213], [0.28987830863610903, 0.36495974735762693])
([-1.17547762665905, 1.1751394486895783], [0.18078065407309682, 0.4998457720618293])
([-1.2440239387295458, 0.8720450804540057], [0.05902925334751511, 0.6016956802132387])
([-1.2441259417439434, 0.6289994697149073], [-0.06577256855272472, 0.6762747102291482])
([-1.210142434136089, 0.4368770315976506], [-0.188677607179601, 0.7291942568591364])
([-1.159918613868923, 0.28408169071815415], [-0.30726896099260204, 0.7649583990935568])
([-1.1025329550493486, 0.16100716005909121], [-0.42042727561865095, 0.7869985179897889])
([-1.0426125487031446, 0.06047044972523817], [-0.5276934467175253, 0.7979089270264804])
⋮
([-1.2216434770974924, 1.0146166139269628], [0.12021680827049967, 0.5550113775144959])
([-1.2499499900381474, 0.7423750265723162], [-0.00391841642025021, 0.6423528468283755])
([-1.229831087369157, 0.5265058660314383], [-0.1281828192264333, 0.7053724817632377])
([-1.1861148292768198, 0.355578849211398], [-0.24911096207717606, 0.7491505605432118])
([-1.1314960903669982, 0.21881648422641656], [-0.3650489236780046, 0.7776241823022744])
([-1.0724336821492086, 0.10787192092145223], [-0.4752653800399116, 0.7937719633967358])
([-1.0122234000273318, 0.016617787590257133], [-0.5794991103107785, 0.7998530690972221])
([-0.9525349461453905, -0.059398567433266464], [-0.6777283550437248, 0.7976021348885333])
([-0.8941857396495508, -0.12343221924184411], [-0.7700512224747937, 0.7883718532084002])Давайте построим орбиту и проверим изменение энергии и момента импульса. Мы знаем, что энергия и момент импульса должны быть постоянными. Они также называются первыми интегралами.
plot_orbit(sol) = plot(sol, vars = (3, 4), lab = "Orbit", title = "Kepler Problem Solution")
function plot_first_integrals(sol, H, L)
plot(initial_first_integrals[1] .- map(u -> H(u[2, :], u[1, :]), sol.u),
lab = "Energy variation", title = "First Integrals")
plot!(initial_first_integrals[2] .- map(u -> L(u[2, :], u[1, :]), sol.u),
lab = "Angular momentum variation")
end
analysis_plot(sol, H, L) = plot(plot_orbit(sol), plot_first_integrals(sol, H, L))analysis_plot (generic function with 1 method)analysis_plot(sol, H, L)Давайте попробуем воспользоваться решателем Рунге-Кутты-Нюстрёма для решения этой задачи и проверим изменение первых интегралов.
sol2 = solve(prob, DPRKN6()) # В dt нет необходимости, потому что, в отличие от симплектических
# интеграторов, DPRKN6 является адаптивным
@show sol2.u |> length
analysis_plot(sol2, H, L)Затем попробуем решить ту же задачу с помощью решателя ERKN4, который предназначен для синусоидальных периодических функций.
sol3 = solve(prob, ERKN4()) # В dt нет необходимости, потому что, в отличие от симплектических
# интеграторов, ERKN4 является адаптивным
@show sol3.u |> length
analysis_plot(sol3, H, L)Как мы видим, ERKN4 плохо справляется с этой задачей, потому что она не синусоидальная.
Одним из преимуществ использования DynamicalODEProblem является возможность неявного преобразования задачи ОДУ второго порядка в нормальную систему ОДУ первого порядка, которая может решаться другими решателями ОДУ. В следующем примере давайте воспользуемся решателем Tsit5.
sol4 = solve(prob, Tsit5())
@show sol4.u |> length
analysis_plot(sol4, H, L)Примечание
Во всех решениях наблюдается дрейф, но у методов высокого порядка он меньше из-за их более высокой точности.
Заключение
Симплектический интегратор не всегда полностью сохраняет энергию, но она может восстановиться. Чтобы колебание энергии в конечном итоге восстановилось, симплектический интегратор должен иметь фиксированный временной шаг. Несмотря на изменение энергии, симплектический интегратор отлично сохраняет момент импульса.
Ни интегратор Рунге-Кутты-Нюстрёма, ни интегратор Рунге-Кутты не сохраняют ни энергию, ни момент импульса, и первые интегралы не имеют тенденции к восстановлению. Преимуществом интегратора Рунге-Кутты-Нюстрёма по сравнению с симплектическим интегратором является возможность адаптивности. Преимуществом интегратора Рунге-Кутты-Нюстрёма по сравнению с интегратором Рунге-Кутты является меньший объем вычислений функции на каждом шаге. Решатель ERKN4 лучше всего подходит для синусоидальных решений.
Множественная проекция
В этом примере нам известно, что энергия и момент импульса должны сохраняться. Этого можно достичь посредством проекции на многообразие. Как следует из названия, это процедура проецирования решения ОДУ на многообразие. Начнем с базового случая, когда проецирование на многообразие не применяется.
using DiffEqCallbacks
function plot_orbit2(sol)
plot(sol, vars = (1, 2), lab = "Orbit", title = "Kepler Problem Solution")
end
function plot_first_integrals2(sol, H, L)
plot(initial_first_integrals[1] .- map(u -> H(u[1:2], u[3:4]), sol.u),
lab = "Energy variation", title = "First Integrals")
plot!(initial_first_integrals[2] .- map(u -> L(u[1:2], u[3:4]), sol.u),
lab = "Angular momentum variation")
end
analysis_plot2(sol, H, L) = plot(plot_orbit2(sol), plot_first_integrals2(sol, H, L))
function hamiltonian(du, u, params, t)
q, p = u[1:2], u[3:4]
qdot(@view(du[1:2]), p, q, params, t)
pdot(@view(du[3:4]), p, q, params, t)
end
prob2 = ODEProblem(hamiltonian, [initial_position; initial_velocity], tspan)
sol_ = solve(prob2, RK4(), dt = 1 // 5, adaptive = false)
analysis_plot2(sol_, H, L)При отсутствии проекции на многообразие первые интегралы имеют существенные колебания.
function first_integrals_manifold(residual, u, p, t)
residual[1:2] .= initial_first_integrals[1] - H(u[1:2], u[3:4])
residual[3:4] .= initial_first_integrals[2] - L(u[1:2], u[3:4])
end
cb = ManifoldProjection(first_integrals_manifold)
sol5 = solve(prob2, RK4(), dt = 1 // 5, adaptive = false, callback = cb)
analysis_plot2(sol5, H, L)Как мы видим, благодаря проекции на многообразие изменение первых интегралов очень невелико, хотя мы используем решатель RK4, который не является симплектическим. Но подождите: что если выполнить проецирование только на многообразие сохранения энергии?
function energy_manifold(residual, u, p, t)
residual[1:2] .= initial_first_integrals[1] - H(u[1:2], u[3:4])
residual[3:4] .= 0
end
energy_cb = ManifoldProjection(energy_manifold)
sol6 = solve(prob2, RK4(), dt = 1 // 5, adaptive = false, callback = energy_cb)
analysis_plot2(sol6, H, L)Изменения энергии почти нет, но момент импульса изменяется существенно. А что если выполнить проецирование только на многообразие сохранения момента импульса?
function angular_manifold(residual, u, p, t)
residual[1:2] .= initial_first_integrals[2] - L(u[1:2], u[3:4])
residual[3:4] .= 0
end
angular_cb = ManifoldProjection(angular_manifold)
sol7 = solve(prob2, RK4(), dt = 1 // 5, adaptive = false, callback = angular_cb)
analysis_plot2(sol7, H, L)Результат ожидаемый.