# Разработать программу, для решения численным методом гравитационной задачи N тел.

Эволюция системы N материальных точек описывается следующей системой уравнений:
\begin{equation}
    \begin{cases}
        \frac{dr_i}{dt} = v_i\\
        \frac{dv_i}{dt} = \displaystyle\sum_{i \ne j}^{N} Gm_j\frac{r_j - r_i}{|r_j - r_i|^3}
    \end{cases}
\end{equation}


```python
import numpy as np
import time
from scipy.constants import g
import matplotlib.pyplot as plt
from matplotlib.animation import ArtistAnimation
from IPython.display import HTML
```

#### Функция `getAcc` находит проекции ускорения на оси x и y. 


```python
def getAcc(pos, mass, G, softening):
    x = pos[:, 0:1]
    y = pos[:, 1:2]

    dx = x.T - x
    dy = y.T - y

    inv_r3 = (dx ** 2 + dy ** 2 + softening ** 2) ** (-1.5)

    ax = G * (dx * inv_r3) @ mass
    ay = G * (dy * inv_r3) @ mass

    a = np.hstack((ax, ay))

    return a
```

### Задание начальных условий:


```python
t = 0
t_end = 10.0
dt = 0.01
softening = 0.1
field_size = 500
```


```python
n = int(input("Введите количество частиц: "))
des_1 = int(input("Введите 1, если масса различна, иначе 2: "))
des_2 = int(input("Введите 1, если первоначальная скорость = 0, иначе 2: "))

```

    Введите количество частиц: 50
    Введите 1, если масса различна, иначе 2: 1
    Введите 1, если первоначальная скорость = 0, иначе 2: 2



```python
if des_1 == 1:
    mass = np.random.random((n, 1)).astype(float) * 10 ** 4
else:
    mass = np.ones((n, 1)) * 2 * 10**4

if des_2 == 1:
    vel = np.zeros((n, 2))
else:
    vel = np.random.randn(n, 2)

pos = np.random.randn(n, 2) * 100
```


```python
acc = getAcc(pos, mass, g, softening)
```


```python
color = []
pos_x = {}
pos_y = {}
for i in range(n):
    pos_x[i] = []
    pos_y[i] = []
    color.append((
        np.random.random(),
        np.random.random(),
        np.random.random()
    ))
```

### Построение модели

#### Переменная `step_count` содержит количетсво шагов моделирования. 


```python
step_count = int(np.ceil(t_end / dt))
```


```python
fig, ax = plt.subplots()

ax.set_xlim((-field_size, field_size))
ax.set_ylim((-field_size, field_size))

plt.close()
frames = []
```

### Осуществление вычислений

Код просчитывает ускорение и скорость для всех заданных точек.

А также сохраняет кадры для анимации.


```python
for i in range(step_count):
    vel += acc * dt / 2.0
    pos += vel * dt
    line = []
    for j in range(len(pos)):
        pos_x[j].append(pos[j][0])
        pos_y[j].append(pos[j][1])
        temp, = ax.plot(pos_x[j], pos_y[j], color=color[j], linewidth=1)
        line.append(temp,)
    frames.append(line)
    acc = getAcc(pos, mass, g, softening)
    vel += acc * dt / 2.0
    t += dt
```

Создать анимацию:


```python
anim = ArtistAnimation(
    fig,
    frames,
    interval=60,
    blit=True,
    repeat=True
)
```

Отобразить анимацию.


```python
HTML(anim.to_html5_video())
```








0fa42ebeb8e23ace335a26934207c0072894a2fd
ipynb
Jupyter Notebook
individual tasks/N Body simulation.ipynb
HAXF13D/TRO-lab2
e57f0375fb935e7a9c1b42c966ebba194023ee68
MIT
individual tasks/N Body simulation.ipynb
HAXF13D/TRO-lab2
e57f0375fb935e7a9c1b42c966ebba194023ee68
MIT
individual tasks/N Body simulation.ipynb
HAXF13D/TRO-lab2
e57f0375fb935e7a9c1b42c966ebba194023ee68
MIT
Qwen/Qwen-72B
1. YES
2. YES