lr3(1)

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

omega0 = 1.02
a = 0.0001
b = 0.2

def system(x, t, a, b):
    dxdt = x[1]
    d2xdt2 = -omega0**2 * x[0] - b * x[1] + a
    return [dxdt, d2xdt2]

def initial_conditions(a, b, top=True):
    x0 = [0, 0]
    if top:
        x0[1] = np.sqrt(abs(a)) / omega0
    else:
        x0[1] = -np.sqrt(abs(a)) / omega0
    return x0

t = np.linspace(0, 30, 1000)

x0_positive = initial_conditions(a, b, top=True)
solution_positive = odeint(system, x0_positive, t, args=(a, b))

V_x = 0.5 * omega0**2 * solution_positive[:, 0]**2 + a * solution_positive[:, 0]

plt.figure(figsize=(8, 6))
plt.plot(solution_positive[:, 0], V_x, label='V(x)')
plt.xlabel('x')
plt.ylabel('V(x)')
plt.title('Залежність потенціальної енергії V(x) від x')
plt.legend()
plt.grid()
plt.show()

plt.figure(figsize=(8, 6))
plt.plot(solution_positive[:, 0], solution_positive[:, 1], label='a > 0 (dx/dt < 0)')
plt.scatter([x0_positive[0]], [x0_positive[1]], color='blue', marker='>', label='Start', s=100)
plt.xlabel('x')
plt.ylabel('dx/dt')
plt.title('Фазовий портрет')
plt.legend()
plt.grid()
plt.show()