alpha(t) for ring oscillations https://phys.pro/problems/881/

881.py

import numpy as np
import matplotlib.pyplot as plt
from matplotlib import rc, rcParams
from scipy.integrate import solve_ivp

from matplotlib.ticker import FuncFormatter
from matplotlib.ticker import MultipleLocator

# Plot setup {{{
lw = 2
rcParams['font.size'] = 45
rcParams['text.usetex'] = True
rcParams['font.sans-serif'] = ['Computer Modern']
rc('text.latex', preamble=r'\usepackage[T2A]{fontenc} \usepackage[utf8]{inputenc} \usepackage{amsmath} \usepackage{bm}')

rc('axes', linewidth=1.8)
rc('xtick.major', size=4, width=2.0)
rc('ytick.major', size=4, width=2.0)
fig, ax = plt.subplots(figsize=(12,7))
fig.subplots_adjust(left=.20, bottom=.20, right=.90, top=.70)

# Set comma as a decimal separator (bad)
#import locale
#locale.setlocale(locale.LC_NUMERIC, 'ru_RU')
#rcParams['axes.formatter.use_locale'] = True
def format_decimal(value, tick_number):
    value = round(value, 2) # Fixes 0,60000000000001 instead of 0,6
    return '$' + str(value).replace('.', '{,}') + '$'
#ax.yaxis.set_major_formatter(FuncFormatter(format_decimal))
#ax.xaxis.set_major_formatter(FuncFormatter(format_decimal))
def format_bold(value, tick_number):
    value = round(value)
    return r'$\bm{' + str(value).replace('.', '{,}') + '}$'
ax.yaxis.set_major_formatter(FuncFormatter(format_bold))
ax.xaxis.set_major_formatter(FuncFormatter(format_bold))
#}}}

# Set parameters
m = 0.100 # kg
M = 0.200 # kg
R = 0.5   # m
g = 10    # m/s²

w0 = np.sqrt(m*g / (2*M*R))
T = 2*np.pi/w0

# ä = - mȧ²sina / (M+2m·sin²(a/2)) - mg sina / 2R(M+2m·sin²(a/2))

def system_of_equations(t, y):
    a, a_dot = y
    a_double_dot = (
        -m*a_dot**2*np.sin(a)/(M+2*m*np.sin(a/2)**2)
        -m*g*np.sin(a)/(2*R*(M+2*m*np.sin(a/2)**2))
    )
    return [a_dot, a_double_dot]

# Initial parameters
a0     = np.radians(30)  # Initial angle, rad
a_dot0 = 0    # Initial angular velocity, rad/s

# Set time interval
t_min = 0   # s
#t_max = 10  # s
t_max = 3 * T
times = np.linspace(t_min, t_max, 10000)

# Solve differential equation
solution = solve_ivp(system_of_equations,
                     (t_min, t_max), 
                     [a0, a_dot0],
                     t_eval=times
                     )

# Get solution
t = solution.t
a_values = solution.y[0]


# Plot numerical solution
ax.plot( t, np.degrees(a_values),
        label=r'$\textrm{\bf Численное решение}$',
        color='k',
        lw=lw*1.8,
        #ls = (0, (13, 3.5))
        ls = (0, (1, 1))
        )
ax.set_xlabel(r'$\bm{t}, \textrm{\bf с}$')
ax.set_ylabel(r'$\bm{\alpha, {}^\circ}$')


# Plot simplified solution a = a0·cos(ωt) for comparison
alphas = a0 * np.cos (w0*times)
ax.plot(times, np.degrees(alphas),
        label=r'$\bm{\alpha = \alpha_0\cdot\cos\omega_0 t}$',
        lw=lw*2, color='k')

#plt.figtext(0.77, 0.92, r'$m={}\;\textrm{{кг}}$'.format(m).replace('.', '{,}') )
#plt.figtext(0.77, 0.85, r'$M={}\;\textrm{{кг}}$'.format(M).replace('.', '{,}') )
#plt.figtext(0.77, 0.78, r'$R={}\;\textrm{{м}}$'.format(R).replace('.', '{,}') )
#
#plt.figtext(0.54, 0.89, r'$\omega_0 = \sqrt{\dfrac{m}{2M}\cdot\dfrac{g}{R}}$')
#plt.figtext(0.54, 0.78, r'$\alpha_0={}^\circ$'.format(round(np.degrees(a0),1)).replace('.', '{,}'))

ax.grid()
ax.legend(loc='upper center', bbox_to_anchor=(0.525, 1.67),
          ncol=1, fancybox=True, shadow=True)
ax.yaxis.set_major_locator(plt.MaxNLocator(5))
ax.xaxis.set_major_locator(plt.MaxNLocator(6))


plt.show()