# Seeking Assistance with Numerical Solution of Spin-up Equation

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

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:

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()
``````

Hi Marioskal, thanks for your suggestion.