重力

范围

此笔记本使用 （点质量动力学）PMD 类来模拟大质量物体之间的引力相互作用。

编码

%load_ext autoreload
%autoreload 2

The autoreload extension is already loaded. To reload it, use:
  %reload_ext autoreload

import numpy as np
import matplotlib.pyplot as plt
import matplotlib as mpl
from scipy import integrate, optimize, spatial
from matplotlib import animation, rc
import sympy as sp
from PMD import PMD, distances, MetaForce
sp.init_printing(use_latex = "mathjax")
rc('animation', html='html5')
%matplotlib nbagg

class Gravity(PMD):
    def __init__(self, G = 6.67e-11, **kwargs):
        """
        2D gravity
        """
        self.G  = G
        super().__init__(**kwargs)

    def derivative(self, X, t, cutoff_radius = 1.e-2):
        """
        Acceleration de chaque masse !
        """
        m, G = self.m, self.G
        n = len(m)
        P = X[:2 * n ].reshape(n ,2)
        V = X[ 2 * n:].reshape(n ,2)
        M = m * m[:, np.newaxis]
        D, R, U = distances(P)
        np.fill_diagonal(R, np.inf)
        if cutoff_radius > 0.: R = np.where(R > cutoff_radius, R, cutoff_radius)
        F =((G * M * R**-2)[:,:,np.newaxis] * U).sum(axis = 0)
        A = (F.T /m).T
        X2 = X.copy()
        X2[:2*n ] = V.flatten()
        X2[ 2*n:] = A.flatten()
        return X2

模拟一个高度对称的案例

# SETUP
G = 1.e03
nr = 4
nt = 4
nm = nr * nt + 1
m = np.ones(nm) * 4.e-3
m[0] = 1.

# CORRECTIONS
r = np.linspace(1., 2., nr)
theta = np.linspace(0., np.pi * 2, nt, endpoint = False)
R, Theta = np.meshgrid(r, theta)
r = np.concatenate([[0.],  R.flatten()])
theta = np.concatenate([[0.],  Theta.flatten()])

v = np.zeros_like(r)
v[1:] = (G * m[0] / r[1:])**.5 * .75 * np.random.normal(loc = 1., scale = .002, size = nm-1)
x  =   r * np.cos(theta)
y  =   r * np.sin(theta)
vx = - v * np.sin(theta)
vy =   v * np.cos(theta)


P = np.array([x,   y]).transpose()
V = np.array([vx, vy]).transpose()

vG = (V * m[:, np.newaxis]).sum(axis = 0) / m.sum()
V -= vG

s = Gravity(m = m, P = P, V = V, G = G, nk = 500)
dt = 1.e-3
nt = 50

pcolors = "r"
tcolors = "k"


fig = plt.figure()
ax = fig.add_subplot(1,1,1)
ax.set_aspect("equal")
margin = 1.
plt.axis([-2, 2, -2, 2])
plt.grid()
ax.axis("off")
points = []

msize = 10. * (s.m / s.m.max())**(1./6.)
for i in range(nm):
    plc = len(pcolors)
    pc = pcolors[i%plc]
    tlc = len(tcolors)
    tc = tcolors[i%tlc]
    trail, = ax.plot([], [], "-"+tc)
    point, = ax.plot([], [], "o"+pc, markersize = msize[i])
    points.append(point)
    points.append(trail)


def init():
    for i in range(2 * nm):
        points[i].set_data([], [])
    return points

def animate(i):
    s.solve(dt, nt)
    x, y = s.xy()
    for i in range(nm):
        points[2*i].set_data(x[i:i+1], y[i:i+1])
        xt, yt = s.trail(i)
        points[2*i+1].set_data(xt, yt)
    return points

anim = animation.FuncAnimation(fig, animate, init_func=init, frames=800, interval=20, blit=True)


plt.close()
anim
#plt.show()

<IPython.core.display.Javascript object>