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

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!

• The 2nd ODE has a different power of $M$, so it has a completely different time scale (by $10^9$), which would cause trouble when trying to solve it numerically. Does the ODE solver give you any warnings at all? What ODE solver are you using? – Kirill Aug 1 '18 at 7:30
• Hi @Kirill. Online ODE's have not been able to solve this system. I do get a warning: ODEintWarning: Excess work done on this call (perhaps wrong Dfun type). Run with full_output = 1 to get quantitative information. warnings.warn(warning_msg, ODEintWarning). With this, I have tried to see what other differential type could be the problem, but have had no luck. – Nicole I. Aug 1 '18 at 14:16
• Since your ODE is poorly scaled, it can be tricky to find the right settings for the ODE solver. But could you not instead just solve the ODE $\frac{\mathrm{d}q}{\mathrm{d}e}=\cdots$, it even has a closed form? You could then recover $q(t)$ by solving another ODE for $t=t(e)$. – Kirill Aug 1 '18 at 14:46
• It might also help a numerical solver to rescale time, like $t=t'/M^a$ for some power $a$ of $M$ that will get the numerical solver to work. ODE solvers often make guesses about the typical timescales of the given ODEs and can choose default settings that might be inappropriate for ODEs with a wildly different scale. – Kirill Aug 1 '18 at 14:56
• Would M have anything to do with this though? It is a constant. I could try what you are saying, but with my error I feel like it is a stepping issue. See, my e value should only be between 0 and 1, which may be too small for this type of solver. With the error in mind, I found a similar problem on stack exchange (i'll add the link below) and I redid what they did but with my own equations and values. Unfortunately, the own solution code was littered with bugs I could not understand, but maybe it could help those experienced (like you) get a better idea of how to help me – Nicole I. Aug 1 '18 at 15:19

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