# Coulomb/kepler potential dynamics. NDSolve breaks

I've been trying to do a simple dynamics in coulomb potential (electron(s) around a nucleus). My equations break down. I think it's because of 1/0.

is there a way to make it work?

this is what I have so far:

Needs["DifferentialEquationsNDSolveUtilities"];
eqs = {{Derivative[q][T] == p[T],
Derivative[p][T] == -q[T]/(4 Pi Abs[ q[T]]^3)}, {q == 2,
p == 0.1}};
vars = {q[T], p[T]};
time = {T, 0, 20};
step = 1/25;
solee = NDSolve[eqs, vars, time, Method -> "ExplicitEuler",
StartingStepSize -> step, MaxSteps -> Infinity];
ParametricPlot[Evaluate[vars /. First[solee]], Evaluate[time],
PlotPoints -> 100]
Plot[Evaluate[vars /. First[solee]], Evaluate[time]]


I've used the documentation in tutorial/NDSolveSPRK

Update (1)

In the second plot, you can see that the trajectories start to jump. It's an unphysical behavior. A correct solution would be periodic or pseudo-periodic orbits.

Plot[Evaluate[{p[T]^2/2 - 1/(4 Pi q[T]), p[T]^2/2, -(1/(4 Pi q[T]))} /. First[solee]], Evaluate[time]]


Plots total, kinetic and potential energy. the total energy (blue line) should be a straight horizontal line.

Update (2)

if I comment out 4 Pi, I get the correct result

Needs["DifferentialEquationsNDSolveUtilities"];
eqs = {{Derivative[q][T] == p[T],
Derivative[p][T] == -q[T]/(*4 Pi *) (
Norm[q[T]]^3)}, {q == {1, 0.1, 0.1}, p == {0.1, 1, 0.1}}};
vars = {q[T], p[T]};
time = {T, 0, 20};
step = 0.01;
solee = NDSolve[eqs, vars, time(*, Method->"ExplicitEuler",
StartingStepSize->step,MaxSteps->Infinity*)];
ParametricPlot[Evaluate[{{p[T], q[T]}} /. First[solee]],
Evaluate[time], PlotPoints -> 100]
Plot[Evaluate[{vars} /. First[solee]], Evaluate[time]]
Plot[Evaluate[{Norm[q[T]], q[T]} /. First[solee]], Evaluate[time]]
Plot[Evaluate[{Norm[p[T]]^2/2 - 1/(*4 Pi*) Norm[q[T]], Norm[p[T]]^2/
2, -(1/(*4 Pi*) Norm[q[T]])} /. First[solee]], Evaluate[time]]


So it seems that the problem might be with initial conditions.

• Works without a hitch for me. I get two nice-looking figures. Or is the problem that they are not what you expected? Nov 19, 2012 at 20:01
• Two comments: i) the real problem is 3 dimensional; here you treat it as 1 dimensional; ii) you need to soften your force so that when particles go through each other the force does not become infinite. Nov 19, 2012 at 20:02
• @chris what do you mean by softening the force? Nov 19, 2012 at 20:42
• @SjoerdC.deVries at the time just before 20, the second curve starts jumping. if you increase time (for example to 200), the motion doesn't look like an orbit. Nov 19, 2012 at 21:03
• I would recommend using "StiffnessSwitching" as the method, and to soften the potential (as suggested by others) by just adding a small constant to the denominator: 1/(r + dr). Also, make sure you switch the sign of the force when it crosses the origin in 1D. Nov 19, 2012 at 22:55