
    Dismiss Notice
Dismiss Notice
Join Physics Forums Today!

Recursive integration physics problem

  1. Jul 12, 2005 #1
    I'm trying for fun to do a subatomic particle simulation. I'm having trouble doing something that should be simple which is calculating electron position given the time. Lets say I'm dealing with a hydrogen atom, which consists of an electron and proton.

    r approximately = r0 * v0*t * 0.5*a*t^2

    So how do I calculate the force if it depends on the distance (radius) and the distance depends on the force?

    More stuff you already know:
    integrate(acceleration with respect to time)=velocity
    integrate(velocity with respect to time)=distance (radius)

    Looks like this is the same problem one would have trying to figure out when a rock drop from far above the Sun would hit it.

    Let me know if you have any suggestions on how to find the function I'm looking for.

    Thanks in advance,
  2. jcsd
  3. Jul 12, 2005 #2


    User Avatar
    Science Advisor
    Homework Helper
    Gold Member

    In your force equation, r is the distance between two positions, one being the position (x_p, y_p, z_p) of the proton and the other being the position (x_e,y_e, z_e) of the electron. r is simply the distance using the pythagorean theorem r=\sqrt( (x_p-x_e)^2 + (y_p-y_e)^2 + (z_p-z_e)^2 ). Then you have F.

    Now, given F, you can get a, then, by integration, the "new" v and "new" r. (You probably want to do this vectorially, if you want to model some kind of orbit.) With that new r, you can go back and find F... etc.

    By the way, your expression for r is not the same as the r in your force equation. Your expression for r assumes a constant-wrt-time acceleration, which is valid only for a small interval around a large value of r.
  4. Jul 13, 2005 #3


    User Avatar
    Science Advisor
    Homework Helper

    Dear DanielCar,

    It's a basic belief of modern physics that it is impossible to specify the position and momentum of an electron at the same time. This is the Heisenberg uncertainty principle and the simulation you're writing is in violation of it.

    There is a version of quantum mechanics due to the famous physicist D. Bohm which allows simultaneous exact positions and momenta. It does this by defining a quantum potential, and then calculating, as you are, the classical trajectory for the electron with respect to this quantum potential. The quantum potential is generated by the electron itself, more or less, and it corresponds to the usual quantum wave function. (Actually something like the phase of the wave function.)

    If you work the problem in Bohmian mechanics, which appears to me to be the achieving a realistic simulation that contains both velocity and position, you will find that in the 1S wave states the electron is stationary, and in states with nonzero orbital angular momentum, the electron makes circles centered on the spin axis (at various offsets from the nucleus).

    If you're interested, it's probably worthwhile to pick up one of the several texts by Bohm or others, such as "The Undivided Universe" by Bohm and Hiley, or the many many papers that one can find by searching for "bohmian mechanics".

    For example, the formulas needed for the simulation are given in section 2 of this paper by E. Deotto and G.C. Ghirardi:


    Also, I should mention that various authors have given different interpretations of spin in Bohmian mechanics and that some of these will modify the velocity calculation (though all of these various interpretations will give the same results as standard quantum mechanics).
  5. Jul 13, 2005 #4
    I'm thinking the answer lies with in the two-body problem: http://scienceworld.wolfram.com/physics/Two-BodyProblem.html
    So I need to integrate the velocity equation with respect to dt to come up with a distance equation. Let me know if someone has already done something like this.

  6. Jul 13, 2005 #5
    Yes, ever since Kepler people have been working on the two body problem, but the trouble is that there's no closed form for the angle as a function of time. I suggest you look into the theory of orbital mechanics. Try a google search for "Kepler's equation" and "true anomaly"

    I've been writing a simulation of the N-body problem which works by numerical integration, i.e. take a small time step dt and repeat v:=v+f(x).dt x:=x+v.dt . This has a problem when the particles get close together, pick up a lot of velocity but don't lose it again on the next time step because they are too far apart
  7. Jul 13, 2005 #6
    Interesting. My simulation has had the same issue on occasion. It is caused by the following: Body A approaches Body B. Assume it is on a direct path for simplicity in explanation. The calculation for 1/r^2 fails if body A passes body B, because the force will invert its direction once the bodies pass each other and the simple dt calculation doesn't account for this.
    So in summary:
    1. Two bodies close to each other, forces approach infinity
    2. Once bodies pass it other force inverts and tends to cancel out
    3. Simple dt numerical integration doesn't account for this and just takes the huge value in one direction.

    This can be solved by using paths such as eliptical orbits or by some type of collision detection and adjusting the time quantum and or force accordingly.

    Let me know if this helps,
  8. Jul 13, 2005 #7


    User Avatar
    Science Advisor
    Homework Helper
    Gold Member

  9. Jul 14, 2005 #8
    The above formula applies only to uniformly accelerating bodies.
    When acceleration & distance are interdependent , the corresponding differential equation needs to be solved. You'll find these tricks in any resource
    on the motions of planets.
    But accelerating charges radiate & one must use Maxwell's equations
    to deal with charged particles.
  10. Jul 15, 2005 #9
    Thanks for the link. It suggests the 4th order Runge-Kutta method, but I had actually already tried that, and it doesn't help very much. It is useful for a more accurate overall calculation (e.g. making sure that the ellipse "joins up"). It looks like [tex]d\theta[/tex] is the important quantity. The error with the simple method is ~[tex]d\theta^2[/tex] while for 4th order Runge-Kutta it is ~[tex]d\theta^5[/tex]. However, when the particles are close enough together we have [tex]d\theta[/tex]~1, and so don't get any more accuracy

    As you have found using elliptical orbits isn't that simple, and I'm reluctant to have different time quanta for different bodies in the simulation. I think the easiest way might be to ensure conservation of energy for the close 2-particle systems. I've also been looking at approximating the motion by a parabolic orbit (which is simpler than the elliptical case) when the particles are close enough together.
  11. Jul 15, 2005 #10


    User Avatar
    Science Advisor
    Homework Helper
    Gold Member

  12. Jul 15, 2005 #11


    User Avatar
    Staff Emeritus
    Science Advisor
    Education Advisor
    2016 Award

    I don't know, guys. Maybe all of you understood more about what's going on than what I've read. However, NEGLECTING the issue of the validity of the model to simulate "subatomic particle", the central question that is being asked remains, I think, unanswered, which is how does one find the solution to the equation of motion for a force that has an explicit position dependence, but not time dependence, i.e. F(x), let's say. I believe this is what is being asked when the OP asked


    Note that if F is given as F(t), how to get the solution to the velocity and position is trivial and we can all go home. However, if F is given as F(x) instead, then what?

    The trick here is to use what we all learn in calculus classes and use the chain rule. We all know that

    [tex]a = \frac{F(x)}{m}[/tex]
    [tex]a(t) = \frac{dv}{dt} = \frac{F(x)}{m}[/tex] (1)

    Now, using the chain rule, we can write

    [tex]\frac{dv}{dt} = \frac{dv}{dx} \frac{dx}{dt}[/tex]

    But that last fraction dx/dt is just v. Thus, when we substitute this into (1), we get

    [tex]\frac{dv}{dx} v = \frac{F(x)}{m}[/tex]

    Rearranging, we obtain

    [tex]v dv = \frac{F(x)}{m} dx[/tex]

    Each side now has the "right" integration variable, and v can be solved (assuming you know the boundary conditions).

    Do the same thing to solve for x.


Have something to add?

Similar Discussions: Recursive integration physics problem
  1. Integrals, Physics (Replies: 3)

  2. Integration problem (Replies: 3)