Skip to content

Commit

Permalink
Fixes to Kepler
Browse files Browse the repository at this point in the history
  • Loading branch information
ChrisRackauckas committed May 19, 2024
1 parent 3ff5e35 commit b64455e
Show file tree
Hide file tree
Showing 2 changed files with 5 additions and 6 deletions.
2 changes: 1 addition & 1 deletion docs/Project.toml
Original file line number Diff line number Diff line change
Expand Up @@ -62,7 +62,7 @@ ModelingToolkit = "8"
ODEProblemLibrary = "0.1"
Optimization = "3"
OptimizationNLopt = "0.1, 0.2"
OrdinaryDiffEq = "6.31"
OrdinaryDiffEq = "6.76"
Plots = "1"
RecursiveArrayTools = "2, 3"
SDEProblemLibrary = "0.1"
Expand Down
9 changes: 4 additions & 5 deletions docs/src/examples/kepler_problem.md
Original file line number Diff line number Diff line change
Expand Up @@ -10,8 +10,7 @@ Also, we know that
$${\displaystyle {\frac {\mathrm {d} {\boldsymbol {p}}}{\mathrm {d} t}}=-{\frac {\partial {\mathcal {H}}}{\partial {\boldsymbol {q}}}}\quad ,\quad {\frac {\mathrm {d} {\boldsymbol {q}}}{\mathrm {d} t}}=+{\frac {\partial {\mathcal {H}}}{\partial {\boldsymbol {p}}}}}$$

```@example kepler
using OrdinaryDiffEq, LinearAlgebra, ForwardDiff, Plots;
gr();
using OrdinaryDiffEq, LinearAlgebra, ForwardDiff, Plots
H(q, p) = norm(p)^2 / 2 - inv(norm(q))
L(q, p) = q[1] * p[2] - p[1] * q[2]
Expand All @@ -30,12 +29,12 @@ sol = solve(prob, KahanLi6(), dt = 1 // 10);
Let's plot the orbit and check the energy and angular momentum variation. We know that energy and angular momentum should be constant, and they are also called first integrals.

```@example kepler
plot_orbit(sol) = plot(sol, vars = (3, 4), lab = "Orbit", title = "Kepler Problem Solution")
plot_orbit(sol) = plot(sol, idxs = (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),
plot(initial_first_integrals[1] .- map(u -> H(u.x[2], u.x[1]), sol.u),
lab = "Energy variation", title = "First Integrals")
plot!(initial_first_integrals[2] .- map(u -> L(u[2, :], u[1, :]), sol.u),
plot!(initial_first_integrals[2] .- map(u -> L(u.x[2], u.x[1]), sol.u),
lab = "Angular momentum variation")
end
analysis_plot(sol, H, L) = plot(plot_orbit(sol), plot_first_integrals(sol, H, L))
Expand Down

0 comments on commit b64455e

Please sign in to comment.