# How to add (energy) constraint when using NDSolve to Equation of Motion

To simplify my problem, I will try and solve the Equation of Motion for a particle in a 1D Harmonic Potential.

energy[x_, p_, m_, ω_] := p^2/(2 m) + (m ω^2)/2 x^2


I can define the following function to use NDSolve to solve for the Equation of Motion:

sol[x0_, p0_, m_, ω_, time_] := NDSolve[
{
x'[t] == p[t]/m,
p'[t] == -m ω^2 x[t],
x == x0,
p == p0
},
{x[t], p[t]},
{t, 0, time},
MaxSteps -> ∞
]


I can run the function "sol", however, if I plot how the energy changes in time:

    trajectory = sol[2, 1, 1, 1, 1000];

Plot[
Evaluate[
energy[x[t], p[t], 1, 1] /. trajectory
],
{t, 0, 1000}
] I understand in this case, the energy changes very slowly with time. However, I would still like to fix the energy. Naively, the first thing I tried was:

    solConstraint1[x0_, p0_, m_, ω_, time_] := NDSolve[
{
x'[t] == p[t]/m,
p'[t] == -m ω^2 x[t],
energy[x[t], p[t], m, ω] == energy[x0, p0, m, ω],
x == x0,
p == p0
},
{x[t], p[t]},
{t, 0, time},
MaxSteps -> ∞
]


"solConstraint1" doesn't run at all, then I came across the following Model Constrained Systems as DAEs

    solConstraint2[x0_, p0_, m_, ω_, time_] := NDSolve[
{
x'[t] == p[t]/m,
p'[t] == -m ω^2 x[t],
energy[x[t], p[t], m, ω] == energy[x0, p0, m, ω],
x == x0,
p == p0
},
{x[t], p[t]},
{t, 0, time},
Method -> {"IndexReduction" -> Automatic}
]


"solConstraint2" doesn't run as well. Is there something I am doing wrongly? How do I add an energy constraint to the equation of motion for NDSolve?

I do not know offhand what is a good way to recast as a DAE. One way to enforce the algebraic constraint, without getting an overdetermined (albeit consistent) system, is to use the projection method. I cribbed some of the submethod settings from advanced documentation in tutorial/NDSolveProjection.

sol2[x0_, p0_, m_, ω_, time_] :=
NDSolve[{x'[t] == p[t]/m, p'[t] == -m ω^2 x[t], x == x0,
p == p0}, {x[t], p[t]}, {t, 0, time}, MaxSteps -> ∞,
Method -> {"Projection",
Method -> {"SymplecticPartitionedRungeKutta",
"DifferenceOrder" -> 4, "PositionVariables" -> {x[t]}},
"Invariants" -> energy[x[t], p[t], m, ω]}, PrecisionGoal -> 12]


This seems to remove that drift off the constraint manifold and also to keep fluctiations to around 10^(-9). Could perhaps do better with more tuning of the various options and suboptions.

From the docs,

solConstraint2[x0_, p0_, m_, ω_, time_] :=
NDSolve[{x'[t] == p[t]/m, p'[t] == -m ω^2 x[t], x == x0, p == p0},
{x, p}, {t, 0, time},
Method -> {"TimeIntegration" -> {"Projection", "Invariants" -> energy[x[t], p[t], m, ω]}}]


This keeps the energy to within about 0.12 of its starting value. Not great, but it does not drift.

Note: Using the option WorkingPrecision -> 50 on the original sol also reigns in the drift and keeps it around 6 * 10^-10.