pr4

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

# Визначення параметрів
A = 10
m = 10
r1 = 1
r2 = 1.1
r3 = 4

# Обчислення констант B, C та D
B = (4/3) * A * (1/r1**4 + 1/r2**4 + 1/r3**4)
C = 2 * A * (1/(r1**4*r2**4) + 1/(r2**4*r3**4) + 1/(r1**4*r3**4))
D = 4 * A / (r1**4 * r2**4 * r3**4)

# Визначення функції для системи диференціальних рівнянь
def model(Y, t):
    return [Y[1]/m, 16*A/(Y[0]**17 - 12*B/(Y[0]**13) + 8*C/(Y[0]**9) - 4*D/(Y[0]**5))]

# Визначення часового інтервалу
tspan = np.linspace(0, 100, 1000)  # Змініть кількість точок за потреби

# Визначення початкових умов
q0_values = np.linspace(0.1, 2, 20)
p0_values = np.zeros_like(q0_values)

# Ініціалізація масивів для збереження результатів
q_values = []
p_values = []

# Розв'язання диференціальних рівнянь для різних початкових умов
for q0, p0 in zip(q0_values, p0_values):
    Y = odeint(model, [q0, p0], tspan)

    # Додавання результатів до списків
    q_values.append(Y[:, 0])
    p_values.append(Y[:, 1])

# Побудова фазового портрету
plt.figure()
for q, p in zip(q_values, p_values):
    plt.plot(q, p)

plt.xlabel('Розтяг (q)')
plt.ylabel('Імпульс (p)')
plt.title('Фазовий портрет осцилятора')
plt.grid(True)
plt.show()