# Phase portrait of non-linear system of ode (Triple Galaxy System)

Context: Hello, I am an undergrad student of physics self-studying Python Programming. I am trying to find the value of H_lam for which limit cycles corresponding to the triple galaxy system governed by the equations given here occur

Code

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

# Define function
def TGS(t, z):
x = np.empty(3); u = np.empty(3)
y = np.empty(3); v = np.empty(3)
x[0]= z[8]; x[1]= z[0]; x[2]= z[4]
y[0]= z[9]; y[1]= z[1]; y[2]= z[5]
u[0]= z[10]; u[1]= z[2]; u[2]= z[6]
v[0]= z[11]; v[1]= z[3]; v[2]= z[7]
G = 6.67e-17; m= [2e+43,2e+43,2e+43]
R = []
for i in [0,1,2]:
R1=[]
for j in [0,1,2]:
R1.append(np.sqrt((x[j]-x[i])**2+(y[j]-y[i])**2))
R.append(R1)
return [u[1], sum((G*m[j]*(x[j]-x[1])/((R[1][j])**3)) for j in [0,2])+(H_lam**2)*x[1],\
u[2], sum((G*m[j]*(x[j]-x[2])/((R[2][j])**3)) for j in [0,1])+(H_lam**2)*x[2],\
v[1], sum((G*m[j]*(y[j]-y[1])/((R[1][j])**3)) for j in [0,2])+(H_lam**2)*y[1],\
v[2], sum((G*m[j]*(y[j]-y[2])/((R[2][j])**3)) for j in [0,1])+(H_lam**2)*y[2]]

# Set limits for evaluation
a, b = 0,1000

# Define H_lams and set the stage for solving the Diff eq
H_lams = [61]
t = np.linspace(a, b, 50000)

# Solve the Diff eq
for H_lam in H_lams:
sol = solve_ivp(tbs, [a, b], [0.15, -0.4, 0, -90, 0, 0.8, 0, 0], t_eval=t)
plt.plot(sol.y[0],sol.y[2], "-")

# Make a little extra room for legend
plt.xlim([-1e+10,1e+10])
plt.axis("equal")
plt.legend([f"$$\H_lam={m}$$" for m in H_lams])
plt.axes().set_aspect(1)


My thoughts and the problem: x, y are the x and y coordinates while X and Y are the velocities along the x-axis and y-axis respectively. Note that I wanted to keep the first galaxy at origin with zero velocity. The masses for each body, in this case, are equal but I would like to alter them later.

Now, this code that I have written does give me a plot, but I am not sure if it is truly the phase portrait I am looking for. Can someone please suggest a way to verify it?

Moreover, is there a way to make this code more efficient?

Edit: The plot I am getting is

• I would not use lowercase for position and uppercase for velocity. You’re going to confuse yourself eventually. I’d suggest (x,y,z) and (u,v,w), conventional names that are clearer. Or use X and U vectors or tuples that are clear and can be better optimized by python. – Bill Barth Nov 10 '20 at 18:38
• Thanks a lot for your response. I have changed all the X's as u's and Y's as v's. But is the code correct? Because the only plot I am getting is a straight line for every value of H_lam. I have posted the link to the image in the main post. Thank you – Nirmal Padwal Nov 11 '20 at 8:46
• Welcome to scicomp! A good check for plausibility in your case could be the preservation of momentum and kinetic energy. You could write two simple methods calculating these and plotting them alongside your phase portrait. For your initial conditions you should take care that your total momentum is zero, otherwise the bodies will drift away in your coordinate system. – MPIchael Nov 12 '20 at 15:34
• You made some efforts towards the barycentric frame, but the derivatives function does not seem to be complete, receiving 12 dimensions and returning 8. Also, you need to zero out the mass center and total impulse. The numbers you use in constants, masses and initial conditions seem not very related to each other. With G and the masses you seem to use meters, seconds and kilograms as base units, but then the initial conditions for positions and velocities are ridiculously low. If you use parsec or similar as length unit, I think the gravitational constant would be different. – Lutz Lehmann Nov 23 '20 at 13:30