Занятие №12:

Центральные силовые поля

Постановка задачи

НО в общем виде сила это вектор!

F_G

Одномерная формулировка

=
\dfrac{GMm}{r^2}
F_K
=
\dfrac{kq_1q_2}{r^2}

Одномерные задачи

M, Q
m, q
F_G
F_K
r
O

Двумерная формулировка

F_G
=
\dfrac{GMm}{r^2}
|
\dfrac{\vec{r}}{r}
\cdot
r
=
|\vec{r}|
\vec{F}_G
\vec{F}_G
=
\dfrac{GMm}{r^2} \cdot \dfrac{\vec{r}}{r}
=
\dfrac{GMm}{r^3}\cdot \vec{r}
r
=
\sqrt{x^2 + y^2}
=
\left(x^2 + y^2\right)^{\frac{1}{2}}
\vec{F}_K
=
\dfrac{kq_1q_2}{r^3}\cdot \vec{r}

Аналогично тут и далее 

Модель центрального силового поля

\begin{cases} F_{Gx} = - \dfrac{GMm}{\left(x^2 + y^2\right)^{\frac{3}{2}}}\cdot x \\ \\ F_{Gy} = - \dfrac{GMm}{\left(x^2 + y^2\right)^{\frac{3}{2}}} \cdot y \end{cases}

Центральное поле

M, Q
\vec{r}
x
x
y
y
M
\ll
m
Q
\ll
q
m, q
\vec{F}_G
=
\{F_{Gx}, F_{Gy}\}

Дифференциальная задача

\vec{F}_G
=
\{ma_{x}, ma_{y}\}
\vec{F}_P
=
m\vec{a}
=
\begin{cases} m \cdot \dfrac{d^2x}{dt^2} = - \dfrac{GMm}{\left(x^2 + y^2\right)^{\frac{3}{2}}}\cdot x \\ \\ m \cdot \dfrac{d^2y}{dt^2} = - \dfrac{GMm}{\left(x^2 + y^2\right)^{\frac{3}{2}}} \cdot y \end{cases}
\begin{cases} a_{x} = \dfrac{d^2x}{dt^2} \\ \\ a_{y} = \dfrac{d^2y}{dt^2} \end{cases}
\begin{cases} \dfrac{dx}{dt} = v_x \\ \\ \dfrac{dv_x}{dt} = - \dfrac{GM}{\left(x^2 + y^2\right)^{\frac{3}{2}}} \cdot x \\ \dfrac{dy}{dt} = v_y \\ \\ \dfrac{dv_y}{dt} = - \dfrac{GM}{\left(x^2 + y^2\right)^{\frac{3}{2}}}\cdot y \\ \end{cases}

Дифференциальная задача

Начальные условия

y
M, Q
x
\vec{v}
\vec{v}
\vec{v}
\vec{v}
x_0
=
149\cdot10^9 \, \text{м}
y_0
=
0 \, \text{м}
v_{x_0}
=
0 \, \text{м/с}
v_{y_0}
=
30 000 \, \text{м/с}
x_0
=
0 \, \text{м}
y_0
=
-149\cdot10^9 \, \text{м}
v_{x_0}
=
30 000 \, \text{м/с}
v_{y_0}
=
0 \, \text{м/с}
x_0
=
-149\cdot10^9 \, \text{м}
y_0
=
0 \, \text{м}
v_{x_0}
=
0 \, \text{м/с}
v_{y_0}
=
-30 000 \, \text{м/с}
x_0
=
0 \, \text{м}
y_0
=
149\cdot10^9 \, \text{м}
v_{x_0}
=
-30000 \, \text{м/с}
v_{y_0}
=
0 \, \text{м/с}
earth
x

Алгоритм решения

import numpy as np
from scipy.integrate import odeint
import matplotlib.pyplot as plt
from matplotlib.animation import FuncAnimation

# Определяем переменную величину
frames = 365
seconds_in_year = 365 * 24 * 60 * 60
years = 1
t = np.linspace(0, years*seconds_in_year, frames)


# Определяем функцию для системы диф. уравнений
def move_func(s, t):
    x, v_x, y, v_y = s
    
    dxdt = v_x
    dv_xdt = - G * m * x / (x**2 + y**2)**1.5
    dydt = v_y
    dv_ydt = - G * m * y / (x**2 + y**2)**1.5
    
    return dxdt, dv_xdt, dydt, dv_ydt
  
# Определяем начальные значения и параметры
G = 6.67 * 10**(-11)
m = 1.98 * 10**(30)

x0 = 149 * 10**9
v_x0 = 0
y0 = 0
v_y0 = 30000

s0 = (x0, v_x0, y0, v_y0)
# Решаем систему диф. уравнений
def solve_func(i, key):
    sol = odeint(move_func, s0, t)
    if key == 'point':
        x = sol[i, 0]
        y = sol[i, 2]
    else:
        x = sol[:i, 0]
        y = sol[:i, 2]
    return x, y
# Строим решение в виде графика и анимируем
fig, ax = plt.subplots()

ball, = plt.plot([], [], 'o', color='b')
ball_line, = plt.plot([], [], '-', color='b')
plt.plot([0], [0], 'o', color='y', ms=20)

def animate(i):
    ball.set_data(solve_func(i, 'point'))
    ball_line.set_data(solve_func(i, 'line'))
    
ani = FuncAnimation(fig,
                    animate,
                    frames=frames,
                    interval=30)
plt.axis('equal')
edge = 2*x0
ax.set_xlim(-edge, edge)
ax.set_ylim(-edge, edge)

ani.save('earth_sun.gif')

Несколько тел в центральном поле

Начальные условия

y
M, Q
x
\vec{v}
\vec{v}
earth
mercury
import numpy as np
from scipy.integrate import odeint
import matplotlib.pyplot as plt
from matplotlib.animation import FuncAnimation

# Определяем переменную величину
frames = 365
seconds_in_year = 365 * 24 * 60 * 60
years = 1
t = np.linspace(0, years*seconds_in_year, frames)
# Определяем функцию для системы диф. уравнений
def move_func(s, t):
    (x1, v_x1, y1, v_y1,
     x2, v_x2, y2, v_y2) = s

    dxdt1 = v_x1
    dv_xdt1 = - G * m * x1 / (x1**2 + y1**2)**1.5
    dydt1 = v_y1
    dv_ydt1 = - G * m * y1 / (x1**2 + y1**2)**1.5

    dxdt2 = v_x2
    dv_xdt2 = - G * m * x2 / (x2**2 + y2**2)**1.5
    dydt2 = v_y2
    dv_ydt2 = - G * m * y2 / (x2**2 + y2**2)**1.5

    return (dxdt1, dv_xdt1, dydt1, dv_ydt1,
            dxdt2, dv_xdt2, dydt2, dv_ydt2)
# Определяем начальные значения и параметры
G = 6.67 * 10**(-11)
m = 1.98 * 10**(30)

x10 = 149 * 10**9
v_x10 = 0
y10 = 0
v_y10 = 30000

x20 = 0
v_x20 = -47360
y20 = 0.387 * 149 * 10**9
v_y20 = 0

s0 = (x10, v_x10, y10, v_y10,
      x20, v_x20, y20, v_y20)
# Решаем систему диф. уравнений
def solve_func(i, key):
    sol = odeint(move_func, s0, t)
    if key == 'point':
        x1 = sol[i, 0]
        y1 = sol[i, 2]
        x2 = sol[i, 4]
        y2 = sol[i, 6]
    else:
        x1 = sol[:i, 0]
        y1 = sol[:i, 2]
        x2 = sol[:i, 4]
        y2 = sol[:i, 6]
    return ((x1, y1), (x2, y2))
# Строим решение в виде графика и анимируем
fig, ax = plt.subplots()

ball1, = plt.plot([], [], 'o', color='b')
ball_line1, = plt.plot([], [], '-', color='b')

ball2, = plt.plot([], [], 'o', color='r')
ball_line2, = plt.plot([], [], '-', color='r')

plt.plot([0], [0], 'o', color='y', ms=20)

def animate(i):
    ball1.set_data(solve_func(i, 'point')[0])
    ball_line1.set_data(solve_func(i, 'line')[0])

    ball2.set_data(solve_func(i, 'point')[1])
    ball_line2.set_data(solve_func(i, 'line')[1])
ani = FuncAnimation(fig,
                    animate,
                    frames=frames,
                    interval=30)

plt.axis('equal')
edge = 2*x10
ax.set_xlim(-edge, edge)
ax.set_ylim(-edge, edge)

plt.show()

Спасибо за понимание!

Made with Slides.com