I’ve been working on a Python code to compute the orbit of a star around a black hole at the center of our Galaxy. The code includes an implementation of a gaseous disc with hydrodynamical drag force considered alongside the gravitational interaction. Initially, I implemented a 4th order Runge-Kutta method to solve the differential equations governing the motion. However, I wanted to enhance the code’s efficiency and accuracy by leveraging the solve_ivp function from the SciPy library.

According to my knowledge, the Runge-Kutta method works just fine. However, I was not able to recieve the right results with the solve_ivp “DOP853”.

Here is the code with DOP853 implementation:

################################################### Imports section #########################################################
import numpy as np
import matplotlib.pyplot as plt
from scipy.constants import G, parsec, year, c
from astropy.constants import M_sun
from scipy.integrate import solve_ivp


# Computing the velocity components of the gaseous disk
def gas_velocity_keplerian_disc(x, y):

    """
    Computes the gas velocity components in Keplerian motion around the black hole.
    Assumes the gas follows a circular orbit in the x-y plane.
    :param x, y: position of the star from the black hole in the x-y plane [m]
    :param m_bh: mass of the black hole [kg]
    :return: velocity components (vx_gas, vy_gas) in the x-y plane [m/s]
    """

    dist = np.sqrt(x**2 + y**2)

    # Magnitude of velocity for a gas particle in a circular orbit
    v_gas_mag = np.sqrt(G * m_bh / dist)

    # Components of velocity assuming the gas moves in a circular orbit in the x-y plane
    vx_gas = -y / dist * v_gas_mag
    vy_gas = x / dist * v_gas_mag

    return vx_gas, vy_gas

# Function to calculate acceleration
def acceleration(t, state, m, c_d, R, ro, h, R_disc):
    x, y, z, vx, vy, vz = state

    dist = np.sqrt(x**2 + y**2 + z**2)

    ax_grav = -G * m_bh * x / dist**3
    ay_grav = -G * m_bh * y / dist**3
    az_grav = -G * m_bh * z / dist**3

    # Check if z is close to +- h/2
    if -h/2 < z and z < h/2:
        # Hydrodynamical drag force components
        vx_gas, vy_gas = gas_velocity_keplerian_disc(x, y)

        vx_rel = vx - vx_gas
        vy_rel = vy - vy_gas
        vz_rel = vz

        v_rel_mag = np.sqrt(vx_rel**2 + vy_rel**2 + vz_rel**2)

        ax_drag = -0.5 * c_d * ro0 * np.pi * R**2 * v_rel_mag * vx_rel / m
        ay_drag = -0.5 * c_d * ro0 * np.pi * R**2 * v_rel_mag * vy_rel / m
        az_drag = -0.5 * c_d * ro0 * np.pi * R**2 * v_rel_mag * vz_rel / m
        print("Hydrodynamic drag force applied.")

    else:
        ax_drag, ay_drag, az_drag = 0.0, 0.0, 0.0
        print("Hydrodynamic drag force not applied.")

    # Total acceleration
    ax = ax_grav + ax_drag
    ay = ay_grav + ay_drag
    az = az_grav + az_drag

    return [vx, vy, vz, ax, ay, az]

# Event function to terminate integration when the distance from the black hole falls below 100 gravitational radii
def event(t, state, m, c_d, R, ro, h, R_disc):
    x, y, z, vx, vy, vz = state

    distance_from_bh = np.sqrt(x**2 + y**2 + z**2)
    return distance_from_bh - 100 * r_grav

event.terminal = True  # Terminate integration when the event is triggered
event.direction = -1   # Trigger the event only when crossing from positive to negative

############################################# Black hole's information section ##############################################
# Black hole info (in SI units)
m_bh = 4.31e6 * M_sun.value  # Mass of the black hole in kg
x_bh, y_bh, z_bh = 0.0, 0.0, 0.0
vx_bh, vy_bh, vz_bh = 0.0, 0.0, 0.0

############################################## Additional parameters definition #############################################
r_grav = (G * m_bh) / c**2
r_bondi = (2 * G * m_bh) / c**2

################################################# Initial conditions ########################################################
# Time range
tstart, tstop = 2000 * year, 2200 * year  # Integration time range in years
t_span = (tstart, tstop)
t_eval = np.linspace(tstart, tstop, 1000)

# Initial conditions of a star (in meters and km/s)
initial_state = np.array([-129598458359999.98, -117255748040000.0, 253025561560000.03,
                          -1325901.2, 723308.2, -449334.5])

#initial_state = np.array([1e14, 0.0, 0.0, 0.0, 1e6, 0.0])

# Disk information
c_d = 1
ro0 = 1e-14  #default 1e-11!!!
h = 2e20
R_disc = 1e20
# Shakura-Sunyaev disk parameters
alpha = 0.3  # typical value for the viscosity parameter
M_dot = 1e-6 * M_sun.value / year  # mass accretion rate per year

# Orbiting object info (dust particle)
# R = 696340e3  # Radius of star in meters
R = 2.0   # Radius of dust particle in meters
m = 0.028    # Mass of the dust particle in kg

# Solve the initial value problem with event function
sol = solve_ivp(acceleration, t_span, initial_state, method='DOP853',
                rtol=1e-10, atol=1e-10, args=(m, c_d, R, ro0, h, R_disc), events=event)

This is the expected result (my result from 4th order Runge Kutta code)

expected_rk4_result

And this is the result I’m getting with the DOP853 code:

dop853_result

Clearly, the spiraling is not physically correct.

Can anyone see an error in my implementation and can correct it? Thank you all!