How do I Decrease Orbit Size in My Python Code?

I’m working on a project where my end goal is for the two orbiting stars to collide; I have some code from my class that I’ve been utilizing, and while I was successful in getting two objects to orbit one another/orbit the same point, I can’t get the radius of the orbit to decrease over time. No matter how I seem to manipulate my code, the orbit radius continues to increase. I’ve provided the base code I’ve been manipulating below.

G = 6.67e-11 #Actual gravitational constant
    dt = 1e-3
    mS = 1e6
    mp1 = 1e5
    mp2 = 1e5
    x1 = 1 #starting x position
    y1 = 0 #starting y position
    x2 = -1
    y2 = 0
    vx = 0 #starting x-velocity
    vy = 2 #starting y-velocity
    
    for param in ['figure.facecolor', 'axes.facecolor', 'savefig.facecolor']:
        plt.rcParams[param] = 'black' 

    timestep = 500 #timestep for video

    xlist1 = [x1] #x-position of BH 1
    ylist1 = [y1] #y-position of BH 1

    xlist2 = [x2] #x-position of BH 2
    ylist2 = [y2] #y-position of BH 2

    vxlist1 = [vx]
    vylist1 = [vy]

    vxlist2 = [vx]
    vylist2 = [vy]

    n = 10
    
    plt.style.use("cyberpunk")

    for param in ['figure.facecolor', 'axes.facecolor', 'savefig.facecolor']:
        plt.rcParams[param] = 'black' 
        
    tail_length = 100

    for t in range(1, timestep+1):

        for _ in range(n):
            r21 = x1**2 + y1**2
            x1 = xlist1[-1] + vxlist1[-1] * dt #updates x by taking the final value in the xlist array and adding the final value in the velocity (x direction) array multiplied by our time step
            y1 = ylist1[-1] + vylist1[-1] * dt
            vx = vxlist1[-1] - G*mS*mp1*xlist1[-1]/r21**1.5 * dt #gravitational force GM1M2/r^2
            vy = vylist1[-1] - G*mS*mp1*ylist1[-1]/r21**1.5 * dt #update velocity by taking the final value in the vy array and subtracting the gravitational force at the final location logged in the y position array
            xlist1.append(x1)
            ylist1.append(y1)
            vxlist1.append(vx)
            vylist1.append(vy)
        
            r22 = x2**2 + y2**2
            x2 = xlist2[-1] + vxlist2[-1] * dt #updates x by taking the final value in the xlist array and adding the final value in the velocity (x direction) array multiplied by our time step
            y2 = ylist2[-1] + vylist2[-1] * dt
            vx = vxlist2[-1] - G*mS*mp2*xlist2[-1]/r22**1.5 * dt #gravitational force GM1M2/r^2
            vy = vylist2[-1] - G*mS*mp2*ylist2[-1]/r22**1.5 * dt #update velocity by taking the final value in the vy array and subtracting the gravitational force at the final location logged in the y position array
            xlist2.append(x2)
            ylist2.append(y2)
            vxlist2.append(vx)
            vylist2.append(vy)

With this, I’ve changed velocity, increased and decreased the gravitational constant, and manually tried to decrease the x and y positions by subtracting 0.0001 from each value of the x and y array. At this point, I’m somewhat out of ideas and starting from scratch isn’t really an option, so any and all help is very much appreciated 🙂

>Solution :

If you want the two objects to spiral inward and eventually collide, you need to reduce their orbital radius over time.

In your current implementation, you have the line r21 = x1**2 + y1**2

This calculates distance between two objects (BH 1 and BH 2). However, you are not updating the orbital radius in each iteration. To make the objects spiral inward, you need modify the update equations for the positions and velocities to take into account the changing orbital radius.
Example:

for t in range(1, timestep+1):

    for _ in range(n):
        r21 = np.sqrt(x1**2 + y1**2)
        x1 = xlist1[-1] + vxlist1[-1] * dt
        y1 = ylist1[-1] + vylist1[-1] * dt
        vx = vxlist1[-1] - G * mS * mp1 * x1 / r21**3 * dt
        vy = vylist1[-1] - G * mS * mp1 * y1 / r21**3 * dt
        xlist1.append(x1)
        ylist1.append(y1)
        vxlist1.append(vx)
        vylist1.append(vy)
    
        r22 = np.sqrt(x2**2 + y2**2)
        x2 = xlist2[-1] + vxlist2[-1] * dt
        y2 = ylist2[-1] + vylist2[-1] * dt
        vx = vxlist2[-1] - G * mS * mp2 * x2 / r22**3 * dt
        vy = vylist2[-1] - G * mS * mp2 * y2 / r22**3 * dt
        xlist2.append(x2)
        ylist2.append(y2)
        vxlist2.append(vx)
        vylist2.append(vy)