Tell me more ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.
up vote 1 down vote favorite

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["DifferentialEquations`NDSolveUtilities`"];
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

Needs["DifferentialEquations`NDSolveUtilities`"];
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", 
   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.

share|improve this question
2  
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
2  
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
3  
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
1  
I really think your equations have an error in them. May be as Guillochon said, you have a sign problem. Becuase I tried many options, and still the solution blows up suddenly at around 20 seconds using your IC. If the physics is not correct, no matter how many options one tries, it will not help, and it is a waste of time to try to keep trying different software options to avoid this. If one is sure the physics is correct, then it is a different story. I'd say make your your ODE's are correct first. – Nasser Nov 20 '12 at 12:58
show 5 more comments

Know someone who can answer? Share a link to this question via email, Google+, Twitter, or Facebook.

Your Answer

 
discard

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

Browse other questions tagged or ask your own question.