5
$\begingroup$

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[0] == x0,
   p[0] == 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}
    ]

enter image description here

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[0] == x0,
       p[0] == 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[0] == x0,
       p[0] == 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?

$\endgroup$
3
$\begingroup$

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[0] == x0, 
   p[0] == 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.

$\endgroup$
3
$\begingroup$

From the docs,

solConstraint2[x0_, p0_, m_, ω_, time_] := 
 NDSolve[{x'[t] == p[t]/m, p'[t] == -m ω^2 x[t], x[0] == x0, p[0] == 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.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.