Seeking Assistance with Numerical Solution of Spin-up Equation

Data Analysis
1

Hello Everyone, I’m trying to solve the following spin-up equation df/dt = C_GW*f^{11/3} numerically. I’m solving this equation using numerical methods, but my results don’t match with what I expect based on the analytical solutions. In particular, I expect to see a straight line when I plot the numerical solution in ‘log space,’ because the equation should follow a power law relationship, which looks like a straight line when plotted in log space. Can anyone help me understand why my numerical results aren’t matching with what I expect?

Here is the code I’m using to solve the equation;

import numpy as np
from scipy.integrate import solve_ivp
import matplotlib.pyplot as plt

# Constants
G = 4.5180 * 1e-30  # Units: Pc^3/M_sun sec^2 gravitational constant
c = 9.7221 * 1e-9    # Units: Pc/sec speed of light
xi = 0.58

# Define parameters function
def para(m1, m2):
    M = m1 + m2
    M_c = (m1 * m2) ** (3/5) / M ** (1/5)
    q = m2 / m1
    Lambda = np.sqrt(m1 / m2)
    r_6 = 1e-6
    rho_6 = (1.396 * (1e+13)) * m1 ** (3/4)
    C_GW = (96/5) * np.pi ** (8/3) * ((G * M_c / c ** 3) ** (5/3))
    return M_c, q, C_GW

# Define the ODE function
def df_dt1(t, f, C_GW):
    f_dot1 = C_GW * (f ** (11/3)) 
    return f_dot1
# Define parameters
n = 4
m1 = 1
m2 = np.logspace(-6, -3, n)

# Initialize lists to store results
M_c_val = []
q_val = []
C_GW_val = []


# Calculate C_GW 
for j in m2:
    M_c, q, C_GW, C_DF = para(m1, j)
    M_c_val.append(M_c)
    q_val.append(q)
    C_GW_val.append(C_GW)
  f_0 = [10]
t_span = [0,-1e11]


# Solve ODE
for i in range(len(q_val)):
    Sol_1 = solve_ivp(df_dt1, t_span, f_0, args=(C_GW_val[i],))  
    # Plot the solution
    plt.loglog(-Sol_1.t, Sol_1.y[0])
plt.xlabel('Time Left until Merger (Years)')

plt.ylabel('Frequency (f)')

plt.title('Frequency Evolution')

plt.legend(loc='center left', bbox_to_anchor=(1, 0.5))

plt.grid(True)

plt.show()

Here is my Analytical Solution

image
image1266×552 57 KB

2

Hi Charchit,

From a quick look, I think your solution seems OK. If you overplot the analytical solutions with your numerical results, you get something like this:

image
image730×454 42.2 KB

So the scaling seems fine. It’s just a matter of choosing correct initial conditions to completely compare.

The part of the code I changed is the following:

t_span = [-1e7,-1e13]

def f_analytic(fo, to, c_gw, t):
    return np.power(fo**(-8.3)+8/3*c_gw*(to-t), -3/8)

t_range = -np.logspace(7, 13, 50)

# Solve ODE
for i in range(len(q_val)):
    Sol_1 = solve_ivp(df_dt1, t_span, f_0, args=(C_GW_val[i],))  
    # Plot the solution
    plt.loglog(-Sol_1.t, Sol_1.y[0], label=f"q={q_val[i]}")
    plt.loglog(-t_range, f_analytic(f_0[0], 0, C_GW_val[i], t_range), ls='--', label=f"q_th={q_val[i]}")

plt.xlabel('Time Left until Merger (Years)')
plt.ylabel('Frequency (f)')
plt.title('Frequency Evolution')
plt.legend(loc='center left', bbox_to_anchor=(1, 0.5))
plt.grid(True)
plt.show()
3

Hi Marioskal, thanks for your suggestion.