Solving for a set of coupled ODEs to get correct variable values

Question

My question is about how I can solve a coupled system of ODE's, and print out the variables in a plot.

I am solving for an q value and an e value, seen in this set of coupled ODE's below: $$ \begin{aligned} \frac{dq}{dt} &= \frac{48}{5\pi M^2}(2\pi Mq)^{11/3}e\frac{1+\frac{73}{24}e^2+\frac{37}{96}e^4}{(1-e^2)^{7/2}},\\ \frac{de}{dt} &= -\frac{304}{15M}(2\pi Mq)^{8/3}e\frac{1+\frac{121}{304}e^2}{(1-e^2)^{5/2}}. \end{aligned}$$

With my initial value for $e$ is 0.95 an my initial $q$ is 1, where constant $M$ = 1e9.

For my results, I want to plot q vs. e. I expect the e value to go to zero while q increases.

def q_e(x,t):
    # M, a constant
    M = 1e9

    # q and e input assignment
    q = x[0]
    e = x[1]

    # Eq.
    dqdt = (48/(5*math.pi*M**2))*(2*math.pi*M*q)**(11/3)*((1+(73/24)*e**2+(37/96)*e**4)/(1-e**2)**(7/2))
    dedt = -1*(304/(15*M))*(2*math.pi*M*q)**(8/3)*e*((1+(121/304)*e**2)/(1-e**2)**(5/2))
    return [dfdt,dedt]

x0 = [1,0.95]   # initial values q = 1, e = 0.99
t = np.linspace(0,100)
y = odeint(q_e,x0,t)

q = x[:,0]      # output in sep. columns 
e = x[:,1]

plt.plot(q,e)
#plt.semilogy(q,e) # test

I assumed that if I printed my q and e values, I could see where I got wrong. So I tried it,

print(y)

and got an output of

[[  1.00000000e+00   9.50000000e-01]
[  6.29120765e+05   3.44703571e-05]
[  1.39158341e+02   4.86785907e+02]
[  1.39158341e+02   4.86793360e+02]

Where the [ 1.39158341e+02 4.86793360e+02] repeats itself for many columns after for the length of $t$. It is somewhat positive to notice that the q and e values of [ 6.29120765e+05 3.44703571e-05] are normal.

I have tried messing around with my code, putting the e value instead of time t, but nothing is working. Looking for pointers and help, thanks!

Answer 1

The function $q(e)$ satisfies a first order linear ODE $$ \frac{\mathrm{d}q}{\mathrm{d}e} = \frac{111 e^4+876 e^2+288}{(e^2-1) (121 e^2+304)} q(e), $$ which can be solved very easily by using an integrating factor: $$ q(e) = C_1 \exp\left( \tfrac{111}{121} e-3\operatorname{arctanh}e+\tfrac{10440}{1331 \sqrt{19}} \operatorname{arctan}\left(\tfrac{11}{4 \sqrt{19}} e\right) \right). $$ For the initial conditions you gave this gives $$ q|_{t=\infty} = q(0)/q(e_0) = 38.561177784533462867. $$

To follow up on my comments, this is what rescaling does to the numerical solution (I added a factor of $e$ to dqdt to match the formula).

from numpy import *
from scipy.integrate import odeint, solve_ivp

M = 1e9
scale = 1e22

def rhs(t, u):
  q, e = u[0], u[1]
  dqdt = (48/(5*math.pi*M**2))*(2*math.pi*M*q)**(11/3)*((1+(73/24)*e**2+(37/96)*e**4)/(1-e**2)**(7/2)) * e
  dedt = -1*(304/(15*M))*(2*math.pi*M*q)**(8/3)*e*((1+(121/304)*e**2)/(1-e**2)**(5/2))
  return [dqdt/scale, dedt/scale]

ic = [1.0, 0.95]

sol = solve_ivp(rhs, (0.0, 10.0), ic, method="LSODA", atol=0, rtol=1e-8)
print("success: ", sol.success)
print("final t: ", sol.t[-1])
print("final y: ", sol.y[:,-1])
print("y error: ", abs(sol.y[0,-1] - 38.561177784533462867))

With the rescaling, it succeeds with no warning. Since you say you are interested in plotting $q(e)$, it's only necessary to increase $t$ until $e(t)$ becomes really small, and that $t$ is substantially smaller than $100$. Output:

success:  True
final t:  10.0
final y:  [3.85611807e+01 7.80753973e-19]
y error:  2.9253639937110165e-06