“勾三股四弦五”是毕达哥拉斯定理的一个特例。人们在研究三体问题的时候,尝试把三个星体摆放到边长分别是 3、4、5 的直角三角形的顶点之上,三个星体的质量分别是对边的长度,初始速度都为零,然后看系统轨道的演化,出乎人们的意料,轨道呈现出非常复杂的图案。
import numpy as np
import wq.core.physics.unit.au as au
from math import sqrt
from wq.core.math.ode import verlet as solver
from wq.core.physics.nbody.body3p import accelerationOf
accel = accelerationOf(au, 5.0, 3.0, 4.0)
step = solver(accel)
time = 0
x = np.array([0.0, 0.0, 0.0, 4.0, 3.0, 0.0])
v = np.array([0.0, 0.0, 0.0, 0.0, 0.0, 0.0])
cx = 4.0 * 3.0 / (5.0 + 3.0 + 4.0)
cy = 3.0 * 4.0 / (5.0 + 3.0 + 4.0)
K = 0.0
U = - au.G * (5.0 * 4.0 / 3.0 + 3.0 * 5.0 / 4.0 + 3.0 * 4.0 / 5.0)
E = K + U
def evolve(tao):
global time, x, v, K, U, E
time, x, v = step(time, x, v, tao)
x1 = x[0]
y1 = x[1]
x2 = x[2]
y2 = x[3]
x3 = x[4]
y3 = x[5]
vx1 = v[0]
vy1 = v[1]
vx2 = v[2]
vy2 = v[3]
vx3 = v[4]
vy3 = v[5]
r12 = sqrt((x1 - x2) * (x1 - x2) + (y1 - y2) * (y1 - y2))
r13 = sqrt((x1 - x3) * (x1 - x3) + (y1 - y3) * (y1 - y3))
r23 = sqrt((x2 - x3) * (x2 - x3) + (y2 - y3) * (y2 - y3))
v1 = sqrt(vx1 * vx1 + vy1 * vy1)
v2 = sqrt(vx2 * vx2 + vy2 * vy2)
v3 = sqrt(vx3 * vx3 + vy3 * vy3)
K = (5.0 * (vx1 * vx1 + vy1 * vy1) + 3.0 * (vx2 * vx2 + vy2 * vy2) + 4.0 * (vx3 * vx3 + vy3 * vy3)) /2
U = - au.G * (5.0 * 3.0 / r12 + 5.0 * 4.0 / r13 + 3.0 * 4.0 / r23)
E = K + U
dx = (5.0 * x1 + 3.0 * x2 + 4.0 * x3) / (5.0 + 3.0 + 4.0) - cx
dy = (5.0 * y1 + 3.0 * y2 + 4.0 * y3) / (5.0 + 3.0 + 4.0) - cy
x[0] = x[0] - dx
x[1] = x[1] - dy
x[2] = x[2] - dx
x[3] = x[3] - dy
x[4] = x[4] - dx
x[5] = x[5] - dy
return min(r12, r13, r23), max(v1, v2, v3)
%matplotlib inline
import matplotlib.pyplot as plt
serial_time = [time]
serial_x1 = [0]
serial_y1 = [0]
serial_x2 = [0]
serial_y2 = [4]
serial_x3 = [3]
serial_y3 = [0]
serial_energy = [E]
dt = 0.1
for i in range(6):
for j in range(5000):
minDist, maxVelo = evolve(dt)
dt = min(minDist / maxVelo / 100, 0.1)
x1 = x[0]
y1 = x[1]
x2 = x[2]
y2 = x[3]
x3 = x[4]
y3 = x[5]
serial_time.append(time)
serial_x1.append(x1)
serial_y1.append(y1)
serial_x2.append(x2)
serial_y2.append(y2)
serial_x3.append(x3)
serial_y3.append(y3)
serial_energy.append(E)
plt.subplot(321+i)
plt.scatter(serial_x1, serial_y1, facecolor='r', edgecolor='none', marker='.')
plt.scatter(serial_x2, serial_y2, facecolor='g', edgecolor='none', marker='.')
plt.scatter(serial_x3, serial_y3, facecolor='b', edgecolor='none', marker='.')
serial_time = [serial_time[-1]]
serial_x1 = [serial_x1[-1]]
serial_y1 = [serial_y1[-1]]
serial_x2 = [serial_x2[-1]]
serial_y2 = [serial_y2[-1]]
serial_x3 = [serial_x3[-1]]
serial_y3 = [serial_y3[-1]]
serial_energy = [serial_energy[-1]]
plt.show()
import base64
from tempfile import NamedTemporaryFile
VIDEO_TAG = """<video controls>
<source src="data:video/mp4;base64,{0}" type="video/mp4">
Your browser does not support the video tag.
</video>"""
def anim_to_html(anim):
if not hasattr(anim, '_encoded_video'):
with NamedTemporaryFile(suffix='.mp4') as f:
anim.save(f.name, fps=200, extra_args=['-vcodec', 'libx264'])
video = open(f.name, "rb").read()
anim._encoded_video = base64.b64encode(video)
return VIDEO_TAG.format(anim._encoded_video)
from IPython.display import HTML
def display_animation(anim):
plt.close(anim._fig)
return HTML(anim_to_html(anim))
import matplotlib.pyplot as plt
# First set up the figure, the axis, and the plot element we want to animate
fig = plt.figure()
ax = plt.axes(xlim=(-3, 5), ylim=(-2, 8))
liner, = ax.plot([], [], marker='.', color='r', lw=2)
lineg, = ax.plot([], [], marker='.', color='g', lw=2)
lineb, = ax.plot([], [], marker='.', color='b', lw=2)
# initialization function: plot the background of each frame
def init():
liner.set_data([0], [0])
lineg.set_data([0], [4])
lineb.set_data([3], [0])
return [liner, lineg, lineb]
# animation function. This is called sequentially
dt = 0.1
def animate(i):
global dt
minDist, maxVelo = evolve(dt)
dt = min(minDist / maxVelo / 100, 0.1)
x1 = x[0]
y1 = x[1]
x2 = x[2]
y2 = x[3]
x3 = x[4]
y3 = x[5]
liner.set_data([x1], [y1])
lineg.set_data([x2], [y2])
lineb.set_data([x3], [y3])
return [liner, lineg, lineb]
%matplotlib inline
from matplotlib import animation
# call the animator. blit=True means only re-draw the parts that have changed.
anim = animation.FuncAnimation(fig, animate, init_func=init,
frames=30000, interval=5, blit=True)
# call our new function to display the animation
display_animation(anim)