I tryin to plot the position of electron in the hidrogen atom by Coulomb's law, $F=K \frac{q_1 q_2}{r^3}$

This is mi code

    import numpy as np
import matplotlib.pyplot as plt

def coulomb(re, v0):

    #Set parameters:
    N = 500      
    dt = 2.2222222 / N  # Time Step:

    #Create an array, for all variables, of size N with all entries equal to zero:
    xe = np.zeros((N,))
    ye = np.zeros((N,))
    vxe = np.zeros((N,))
    vye = np.zeros((N,))

    # Initial Conditions:
    xe[0] = re                   # (x0 = r, y0 = 0) 
    vye[0] = v0                    #units in m/s^2

    #Implement Verlet Algorithm:
    for k in range(0, N-1):
        re = (xe[k]*2+ye[k]*2)*0.5
        vxe[k+1] = vxe[k] - ((mu * xe[k]) / (re**3)) * dt 
        xe [k+1] = xe[k] + vxe[k+1]*dt
        vye[k+1] = vye[k] - ((mu * ye[k]) / (re**3)) * dt 
        ye [k+1] = ye[k] + vye[k+1]*dt

    xi = plt.plot(xe, ye, 'go', markersize = 1)
    plt.plot(0,0,'yo')                  # yellow marker
    plt.plot(xe[0],0,'bo')  # dark blue marker
    plt.xlabel ('x')
    plt.ylabel ('y')

    return xi, xe, ye

# average distance electron-nucleus in meter
r = 5.1e-11
k = 8.9e9
mu = r*3 * 4 * k *np.pi*2  # coulomb parameter 

coulomb(r, np.sqrt(mu / r));

This works well, I've try to changing the value of Permittivity in l, but the graph becomes just 2 dots (the values of xe and vye remain constant) instead of a circular orbit and I don't know why this is? this happens

enter image description here

  • $\begingroup$ Where do you compute mu? It should be mu = k*el^2/me where el is the elementary charge and me the electron mass? What variable is or is influenced by "Permittivity in l"? $\endgroup$ Jan 30, 2021 at 11:09
  • $\begingroup$ oh sorry man, i just realize that i put the wrong code, i upload the correct version $\endgroup$
    – Rei D Gar
    Jan 30, 2021 at 18:21

1 Answer 1


The Kepler laws, or simple insertion of $z(t)=re^{iωt}$, tell that for a circular orbit of radius $r$ in the central field $\ddot z=-\mu\frac{z}{|z|^3}$ the angular speed $\omega$ is given by $$ r^3ω^2=μ. $$ That is how you get the initial velocity $(0,\sqrt{\mu/r})$ at the point $(r,0)$. The period of the orbit is $$ T=\frac{2\pi}{ω}=2π\sqrt{r^3/μ}. $$ Inserting $\mu=k·4\pi^2·r^3$ as coded gives $T=k^{-1/2}=1.06·10^{-5}$. Your sampling rate is just not high enough. Or you need to reduce the integration time to about $2.2·10^{-5}$ to integrate over two periods.

But most importantly, check your code and formula everywhere that you use the ** operator where you want to compute a power, at the moment there are many places where you have just one multiplication sign.

enter image description here

lines to repair:

re = (xe[k]**2+ye[k]**2)**0.5
mu = k * r**3 * 4 * np.pi**2  # coulomb parameter 
  • $\begingroup$ Tks bro, could you help me with that? I understood the theory, I tried to implement it but I got nothing. it's an assignment I have to pass the course :( $\endgroup$
    – Rei D Gar
    Jan 30, 2021 at 21:42
  • $\begingroup$ Hey man! I passed my assignment, I'm very thankfull for your help. $\endgroup$
    – Rei D Gar
    Feb 4, 2021 at 3:23

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.