Recursive integration physics problem

[DanielCar]

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. Let's say I'm dealing with a hydrogen atom, which consists of an electron and proton.

f=(k*q1*q2)/(r^2)
r approximately = r0 * v0*t * 0.5*a*t^2
f=ma
a=f/m

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,
Daniel

 

[robphy]

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.

 

[CarlB]

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:
http://arxiv.org/abs/quant-ph/9704021

Carl

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).

 

[DanielCar]

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.

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.

Daniel

 

[chronon]

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.

Daniel

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

 

[DanielCar]

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

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,
Daniel

 

[robphy]

This looks useful:
http://www.kopernekus.com/physics/orbitalmotion.asp

 

[einstone]

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

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.
Regards,
Einstone.

 

[chronon]

This looks useful:
http://www.kopernekus.com/physics/orbitalmotion.asp

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 d\theta is the important quantity. The error with the simple method is ~d\theta^2 while for 4th order Runge-Kutta it is ~d\theta^5. However, when the particles are close enough together we have d\theta~1, and so don't get any more accuracy

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,
Daniel

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.

 

[robphy]

My previous post was intended for the OP's "single particle in a central force" problem.

For n mutually interacting particles, these may be useful:
http://www.amara.com/papers/nbody.html
http://www.schlitt.net/xstar/n-body.pdf
http://www-2.cs.cmu.edu/~scandal/applets/bh.html

 

[ZapperZ]

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. Let's say I'm dealing with a hydrogen atom, which consists of an electron and proton.

f=(k*q1*q2)/(r^2)
r approximately = r0 * v0*t * 0.5*a*t^2
f=ma
a=f/m

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,
Daniel

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

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

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

a = \frac{F(x)}{m}
a(t) = \frac{dv}{dt} = \frac{F(x)}{m} (1)

Now, using the chain rule, we can write

\frac{dv}{dt} = \frac{dv}{dx} \frac{dx}{dt}

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

\frac{dv}{dx} v = \frac{F(x)}{m}

Rearranging, we obtain

v dv = \frac{F(x)}{m} dx

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.

Zz.