Mathematica Stack Exchange is a question and answer site for users of Mathematica. Join them; it only takes a minute:

Sign up
Here's how it works:
  1. Anybody can ask a question
  2. Anybody can answer
  3. The best answers are voted up and rise to the top

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:

eqs = {{Derivative[1][q][T] == p[T], 
    Derivative[1][p][T] == -q[T]/(4 Pi Abs[ q[T]]^3)}, {q[0] == 2, 
    p[0] == 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

eqs = {{Derivative[1][q][T] == p[T], 
    Derivative[1][p][T] == -q[T]/(*4 Pi *) (
     Norm[q[T]]^3)}, {q[0] == {1, 0.1, 0.1}, p[0] == {0.1, 1, 0.1}}};
vars = {q[T], p[T]};
time = {T, 0, 20};
step = 0.01;
solee = NDSolve[eqs, vars, time(*, Method->"ExplicitEuler", 
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.

share|improve this question
Works without a hitch for me. I get two nice-looking figures. Or is the problem that they are not what you expected? – Sjoerd C. de Vries Nov 19 '12 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. – chris Nov 19 '12 at 20:02
@chris what do you mean by softening the force? – kirill_igum Nov 19 '12 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. – kirill_igum Nov 19 '12 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. – Guillochon Nov 19 '12 at 22:55

Your Answer


By posting your answer, you agree to the privacy policy and terms of service.

Browse other questions tagged or ask your own question.