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.

  • 2
    $\begingroup$ Works without a hitch for me. I get two nice-looking figures. Or is the problem that they are not what you expected? $\endgroup$ Nov 19, 2012 at 20:01
  • 2
    $\begingroup$ 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. $\endgroup$
    – chris
    Nov 19, 2012 at 20:02
  • $\begingroup$ @chris what do you mean by softening the force? $\endgroup$ Nov 19, 2012 at 20:42
  • $\begingroup$ @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. $\endgroup$ Nov 19, 2012 at 21:03
  • 3
    $\begingroup$ 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. $\endgroup$
    – Guillochon
    Nov 19, 2012 at 22:55


Your Answer

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

Browse other questions tagged or ask your own question.