Podúloha 1 – Riešenie cez klasickú mechaniku

Obrázok 1: Zakreslenie síl pôsobiacich na tyčku.

Najprv si zakreslíme všetky sily, ktoré pôsobia na tyčku (Obrázok 1). Na tyčku pôsobí 5 síl: tiažová sila $\vec G$ (veľkosť $m g$), normálová sila od podlahy $\vec F_1$, normálová sila od steny $\vec F_2$, trecia sila od podlahy $\vec F_{f,1}$ (veľkosť $f_1 F_1$) a trecia sila od steny $\vec F_{f,2}$ (veľkosť $f_2 F_2$). Môžeme si všimnúť, že 2 premenné ($F_1$ a $F_2$) sú neznáme. Ich veľkosť musí byť taká, aby sa tyčka neprepadla cez stenu alebo podlahu. Môžeme si ešte všimnúť z geometrie obrázka, že ťažisko (bod $T$) sa pri šmýkaní tyčky pohybuje po kružnici s polomerom $L/2$ okolo rohu steny $(0,0)$.

Z klasickej mechaniky vieme, že pre popis pohybu tuhého telesa nám postačí poznať aká celková sila $\vec F_T$ a celkový moment síl $\vec M_T$ pôsobia na teleso.

Pre náš prípad dostávame

$$ \vec F_T = \sum_i \vec F_i = \vec G + \vec F_1 + \vec F_2 + \vec F_{f,1} + \vec F_{f,2} $$ $$ = ( -f_1 F_1 + F_2, - mg + F_1 + f_2 F_2, 0) \, , $$ $$ \vec M_T = \sum_i \vec r_i \times \vec F_i = \vec 0 \times \vec G + \vec r_{T1} \times\vec F_1 + \vec r_{T2} \times \vec F_2 + \vec r_{T1} \times \vec F_{f,1} + \vec r_{T2} \times \vec F_{f,2} $$ $$ = \left(0,0, \frac{L}{2} \left( F_1 \cos \alpha - f_1 F_1 \sin \alpha - F_2 \sin \alpha - f_2 F_2 \cos \alpha \right) \right) \, . $$

Celková sila $\vec F_T$ a celkový moment síl $\vec M_T$ spôsobujú zrýchlenie telese $\vec a_T$ a uhlové zrýchlenie $\varepsilon_T$ podľa Newtonového zákona

$$ \vec F_T = m \vec a_T \, , \quad \quad \vec M_T = I \varepsilon_T \vec e_M \, , $$

kde $I$ je moment zotrvačnosti telesa okolo rotačnej osi, ktorá je určená jednotkovým vektorom v smere celkového momentu síl $\vec e_M$.

Pre tyčku hmotnosti $m$ a dĺžky $L$ platí, že jej moment zotrvačnosti okolo ťažiska je $\frac{1}{12} m L^2$. Pre zrýchlenie ťažiska $\vec a_T$ a uhlové zrýchlenie $\varepsilon_T$ dostaneme teda

$$ \vec a_T = \frac{1}{m} \vec {F_T} = \left( -f_1 \frac{F_1}{m} + \frac{F_2}{m}, - mg + \frac{F_1}{m} + f_2 \frac{F_2}{m}, 0 \right) \, , $$

$$ \varepsilon_T = \frac{1}{ \frac{1}{12} m L^2 } \left| \vec M_T \right| = \frac{6}{L} \left( \frac{F_1}{m} \cos \alpha - f_1 \frac{F_1}{m} \sin \alpha - \frac{F_2}{m} \sin \alpha - f_2 \frac{F_2}{m} \cos \alpha \right) \, . $$

Zatiaľ sme dostali iba všeobecné pohybové rovnice pre dané sily $F_1$ a $F_2$. My však máme 2 väzobné podmienky: bod 1 nesmie prepadnúť cez podlahu a bod 2 nesmie prepadnúť cez stenu. Musíme teda najprv vypočítať pohybové rovnice pre tieto body. Keď už poznáme $\vec a_T$ a $\vec \varepsilon_T$, tak vieme vypočítať zrýchlenie každého bodu tuhého telesa pomocou vzťahu

$$ \vec a_A = \vec a_T + \vec \varepsilon_T \times \vec r_{TA} + \vec \omega \times \left( \vec \omega \times \vec r_{TA} \right) \, , $$

kde prvý člen zodpovedá celkovému zrýchleniu celého telesa, druhý člen zodpovedá uhlovému zrýchleniu a tretí člen zodpovedá odstredivému zrýchleniu.

Pre bod 1 ($\vec r_{T1} = \left( \frac{L}{2} \cos \alpha, - \frac{L}{2} \sin \alpha, 0 \right)$) a bod 2 ($\vec r_{T1} = \left( - \frac{L}{2} \cos \alpha, \frac{L}{2} \sin \alpha, 0 \right)$) dostávame:

$$ \vec a_1 = \left( -f_1 \frac{F_1}{m} + \frac{F_2}{m}, - mg + \frac{F_1}{m} + f_2 \frac{F_2}{m}, 0 \right) + \left( +\frac{L}{2} \varepsilon_T \sin \alpha , +\frac{L}{2} \varepsilon_T \cos \alpha, 0 \right) + \left( - \frac{\omega^2 L}{2} \cos \alpha , +\frac{\omega^2 L}{2} \sin \alpha, 0 \right) \, , $$

$$ \vec a_2 = \left( -f_1 \frac{F_1}{m} + \frac{F_2}{m}, - mg + \frac{F_1}{m} + f_2 \frac{F_2}{m}, 0 \right) + \left( - \frac{L}{2} \varepsilon_T \sin \alpha , - \frac{L}{2} \varepsilon_T \cos \alpha, 0 \right) + \left( + \frac{\omega^2 L}{2} \cos \alpha , - \frac{\omega^2 L}{2} \sin \alpha, 0 \right) \, . $$

Aby sme zabránili prepadnutiu cez podlahu, tak v každom momente predpokladáme $a_{1,y} = 0$. Podobne na stene dostávame podmienku $a_{2,x} = 0$.¹ Dostaneme teda dve rovnice (2 podmienky) o dvoch neznámych ($F_1$ a $F_2$). Aby sme si zjednodušili zápis, tak si zadefinujeme nasledujúce skratky: $a_1 := \frac{F_1}{m}$, $a_2 := \frac{F_2}{m}$, $s := \sin \alpha$, $c:= \cos \alpha$.

Podmienky $a_{1,y} = 0$ a $a_{2,x} = 0$ sa dajú potom zapísať ako

$$ a_1 (1 + 3 c^2 - 3 f_1 s c) + a_2 (f_2 - 3 s c - 3 f_2 c^2) = g - \frac{\omega^2 L}{2} s $$

$$ a_1 (-f_1 - 3 s c + 3 f_1 s^2) + a_2 (1 + 3 s^2 + 3 f_2 s c) = - \frac{\omega^2 L}{2} c \, . $$

Ide o jednoduchú lineárnu sústavu rovníc, takže môžeme použiť maticový zápis

$$ \begin{pmatrix} A & B \\ C & D \end{pmatrix} \begin{pmatrix} a_1 \\ a_2 \end{pmatrix} = \begin{pmatrix} g - \frac{\omega^2 L}{2} s \\ - \frac{\omega^2 L}{2} c \end{pmatrix} \, , $$

kde koeficienty $A$, $B$, $C$ a $D$ zodpovedajú vyššie uvedením zátvorkám. Nájsť $a_1$ a $a_2$ znamená vyriešiť maticovú rovnicu. Rovnicu stačí zľava vynásobiť inverznou maticou a dostaneme

$$ \begin{pmatrix} a_1 \\ a_2 \end{pmatrix} = \begin{pmatrix} A & B \\ C & D \end{pmatrix} ^{-1} \begin{pmatrix} g - \frac{\omega^2 L}{2} s \\ - \frac{\omega^2 L}{2} c \end{pmatrix} \, . $$

Pre inverznú maticu veľkosti 2x2 máme vzťah

$$ \begin{pmatrix} a_1 \\ a_2 \end{pmatrix} = \frac{1}{AD - BC} \begin{pmatrix} D & -B \\ -C & A \end{pmatrix} \begin{pmatrix} g - \frac{\omega^2 L}{2} s \\ - \frac{\omega^2 L}{2} c \end{pmatrix} \, . $$

Keď už poznáme premenné $a_1$ a $a_2$, tak poznáme aj aké normálové sily $F_1$ a $F_2$ pôsobia a poznáme aj pohybové rovnice pre pohyb ťažiska.

$$ a_{T,x} = -f_1 a_1 + a_2 $$ $$ = \frac{ g s \left(-6 f_1 s + 3 c - 3 f_1 f_2 c \right) + \omega^2 L \left( - 2 c + f_1 f_2 c + 3 f_1 s \right) }{4 - 2 f_1 f_2 + 6 (f_2 -f_1) s c } $$ $$ a_{T,y} = -g + a_1 + f_2 a_2 $$ $$ = \frac{ - g c \left(-6 f_1 s + 3 c - 3 f_1 f_2 c \right) + \omega^2 L \left( - 2 s + f_1 f_2 s - 3 f_2 c \right) }{4 - 2 f_1 f_2 + 6 (f_2 -f_1) s c } $$

Síce poznáme celkové zrýchlenie ťažiska v smeroch $x$ a $y$, nás bude zaujímať rozloženie zrýchlenia podľa iných smerov. Z obrázka vieme, že ťažisko sa bude pohybovať po kružnici. Môžeme teda povedať, že zrýchlenie v normálovom smere (v smere $(0,0) - T$) bude odstredivé zrýchlenie a zrýchlenie v tangenciálnom smere (kolmo na normálový smer v bode $T$) bude určovať uhlové zrýchlenie. Nesmieme zabudnúť, že tieto uhlové a odstredivé zrýchlenia sú iné od uhlového a odstredivého zrýchlenia tuhého telesa, keďže ide o zrýchlenia iba bodového objektu – ťažiska. Rozložiť vektor na normálový a tangenciálny smer sa dá urobiť aj pomocou matíc.

$$ \begin{pmatrix} a_{T,N} \\ a_{T,T} \end{pmatrix} = \begin{pmatrix} +c & +s \\ -s & +c \end{pmatrix} \begin{pmatrix} a_{T,x} \\ a_{T,y} \end{pmatrix} $$

Nás zaujíma tangenciálne zrýchlenie $a_{T,T}$, pre ktoré dostaneme vzťah $$ a_{T,T} = \frac{ - 3 g \left[-2 f_1 \sin \alpha + (1 - f_1 f_2) \cos \alpha \right] - 3 \omega^2 L \left( f_1 \sin^2 \alpha + f_2 \cos^2 \alpha \right) }{4 - 2 f_1 f_2 + 6 (f_2 -f_1) \sin \alpha \cos \alpha } \, . $$ Toto tangenciálne zrýchlenie je priamo spojené uhlovým zrýchlením uhla $\alpha$ (cez faktor $\frac{L}{2}$). Ak všetko prevedieme do jednej súradnice $\alpha$, tak dostaneme diferenciálnu rovnicu druhého rádu $$ \ddot \alpha = \frac{ - \frac{3 g}{L} \left[-2 f_1 \sin \alpha + (1 - f_1 f_2) \cos \alpha \right] - 3 \dot \alpha^2 \left( f_1 \sin^2 \alpha + f_2 \cos^2 \alpha \right) }{2 - f_1 f_2 + 3 (f_2 -f_1) \sin \alpha \cos \alpha } \, . $$

Podúloha 1 – Riešenie cez parameter

Za toto riešenie ďakujeme Richardovi Dudekovi. Najprv si uvedomíme, že sklz tyčky po stene sa správa rovnako, ako keby oba konce tyčky boli pripevnené ku podlahe a stene. V takom prípade celý pohyb tyčky v každom momente je jednoznačne určeným jedným parametrom – uhlom $\alpha$. Keďže sa ťažisko pohybuje po kružnici, tak jeho poloha je jasne definovaná ako

$$ x_T (t) = \frac{L}{2} \cos \alpha (t) \, , $$ $$ y_T (t) = \frac{L}{2} \sin \alpha (t) \, . $$

Tak ako pri predchádzajúcom riešení nájdeme celkovú silu a moment síl

$$ \vec F_T = ( -f_1 F_1 + F_2, - mg + F_1 + f_2 F_2, 0) \, , $$ $$ \vec M_T = \left(0,0, \frac{L}{2} \left( F_1 \cos \alpha - f_1 F_1 \sin \alpha - F_2 \sin \alpha - f_2 F_2 \cos \alpha \right) \right) \, . $$

Dostaneme 3 rovnice pre pohyb: 2 rovnice pre pohyb ťažiska (v smere $x$ a $y$) a 1 rovnicu pre rotačný pohyb.

$$ \ddot x_{T} = -f_1 \frac{F_1}{m} + \frac{F_2}{m} \, , $$ $$ \ddot y_{T} = -g + \frac{F_1}{m} + f_2 \frac{F_2}{m} \, , $$ $$ \varepsilon = \frac{6}{L} \left( \frac{F_1}{m} \cos \alpha - f_1 \frac{F_1}{m} \sin \alpha - \frac{F_2}{m} \sin \alpha - f_2 \frac{F_2}{m} \cos \alpha \right) \, . $$

Tieto 3 zrýchlenia sú viazané parametrom $\alpha$. Zrýchlenie ťažiska dostaneme dvoma deriváciami vzťahu pre polohu

$$\dot x_T (t) = - \frac{L}{2} \sin \alpha \, \dot \alpha \, , $$ $$ \ddot x_T (t) = \frac{L}{2} \left( - \sin \alpha \, \ddot \alpha - \cos \alpha \, \dot \alpha^2 \right) \, , $$ $$ \dot y_T (t) = \frac{L}{2} \cos \alpha \, \dot \alpha \, , $$ $$ \ddot y_T (t) = \frac{L}{2} \left( \cos \alpha \, \ddot \alpha - \sin \alpha \, \dot \alpha^2 \right) \, . $$

Vzťah medzi uhlovým zrýchlením $\varepsilon$ a uhlom $\alpha$ je cez druhú časovú deriváciu s opačným znamienkom (lebo pri pozitívnom smere momentu, uhol $\alpha$ bude klesať a naopak).

$$ \varepsilon = - \ddot \alpha \, . $$

Po dosadení dostaneme 3 rovnice o 3 neznámych ($F_1$, $F_2$ a $\ddot \alpha$).

$$ \frac{L}{2} \left( - \sin \alpha \, \ddot \alpha - \cos \alpha \, \dot \alpha^2 \right) = -f_1 \frac{F_1}{m} + \frac{F_2}{m} \, , $$ $$ \frac{L}{2} \left( \cos \alpha \, \ddot \alpha - \sin \alpha \, \dot \alpha^2 \right) = -g + \frac{F_1}{m} + f_2 \frac{F_2}{m} \, , $$ $$ - \ddot \alpha = \frac{6}{L} \left( \frac{F_1}{m} \cos \alpha - f_1 \frac{F_1}{m} \sin \alpha - \frac{F_2}{m} \sin \alpha - f_2 \frac{F_2}{m} \cos \alpha \right) \, . $$

Pre priehľadnosť ich vieme napísať v maticovom tvare

$$ \begin{pmatrix} f_1 & - 1 & - \frac{L}{2} \sin \alpha \\ -1 & - f_2 & + \frac{L}{2} \cos \alpha \\ - \cos \alpha + f_1 \sin \alpha & \sin \alpha + f_2 \cos \alpha & - \frac{L}{6} \end{pmatrix} \begin{pmatrix} \frac{F_1}{m} \\ \frac{F_2}{m} \\ \ddot \alpha \end{pmatrix} = \begin{pmatrix} \frac{L}{2} \cos \alpha \, \dot \alpha^2 \\ -g +\frac{L}{2} \sin \alpha \, \dot \alpha^2 \\ 0 \end{pmatrix} \, . $$

Riešením sústavy rovníc nájdeme

$$ \begin{pmatrix} \frac{F_1}{m} \\ \frac{F_2}{m} \\ \ddot \alpha \end{pmatrix} = \begin{pmatrix} \frac{ g \left[ 1 + 3 f_2 \sin \alpha \cos \alpha + 3 \sin^2 \alpha \right] + L \dot \alpha^2 \left( - 2 \sin \alpha - f_2 \cos \alpha \right) }{ 2 \left[2 - f_1 f_2 + 3 (f_2 -f_1) \sin \alpha \cos \alpha \right] } \\ \frac{ g \left[ f_1 - 3 f_1 \sin^2 \alpha + 3 \sin \alpha \cos \alpha \right] + L \dot \alpha^2 \left( - 2 \cos \alpha + f_1 \sin \alpha \right) }{ 2 \left[2 - f_1 f_2 + 3 (f_2 -f_1) \sin \alpha \cos \alpha \right] } \\ \frac{ - \frac{3 g}{L} \left[-2 f_1 \sin \alpha + (1 - f_1 f_2) \cos \alpha \right] - 3 \dot \alpha^2 \left( f_1 \sin^2 \alpha + f_2 \cos^2 \alpha \right) }{2 - f_1 f_2 + 3 (f_2 -f_1) \sin \alpha \cos \alpha } \end{pmatrix} \, . $$

Podúloha 2

Aby sme dostali rovnicu pohybu kyvadla musíme dosiahnuť rovnicu v tvare

$$ \ddot \phi = - K \sin \phi \, , $$

kde $K$ je nejaká vlastná konštanta kyvadla a uhol $\phi$ môže byť v nejakom lineárnom vzťahu voči uhlu $\alpha$.

Keď sa pozrieme na rovnicu pohybu pre $\alpha$, tak pokiaľ menovateľ nebude konštantný, tak budeme mať stále komplikovanú závislosť. Menovateľ je konštantný iba za podmienky $f_1 = f_2 = f$. Vtedy sa rovnica zjednoduší na tvar

$$ \ddot \alpha = - \frac{3 g (1 + f^2)}{L (2 - f^2 )} \left[- \frac{2 f}{1 + f^2} \sin \alpha + \frac{1 - f^2}{1 + f^2} \cos \alpha \right] - \frac{3 f}{2 - f^2 } \dot \alpha^2 \, . $$

Faktor $(1+f^2)$ sme vybrali dopredu zámerne. Keď zadefinujeme novú konštantu $\tan \varphi := \frac{2 f}{1 - f^2}$, tak platí $\sin \varphi = \frac{2 f}{1 + f^2}$ a $\cos \varphi = \frac{1- f^2}{1 + f^2}$ a rovnica sa zjednoduší na tvar

$$ \ddot \alpha = - \frac{3 g (1 + f^2)}{L (2 - f^2 )} \cos \left( \alpha + \varphi \right) - \frac{3 f}{2 - f^2 } \dot \alpha^2 \, , $$

čo sa dá interpretovať ako kyvadlo s konštantou $K = \frac{3 g (1 + f^2)}{L (2 - f^2 )}$, stabilnou polohou $\alpha_0 = - \varphi- \frac{\pi}{2}$ a s odporom vzduchu, keďže odpor vzduchu závisí od druhej mocniny rýchlosti.

Ak chceme mať kyvadlo bez odporu, tak člen s druhou mocninou rýchlosti musí zmiznúť. To nastane iba pre $f_1 = f_2 = f = 0$.

$$ \ddot \alpha = - \frac{3 g }{2 L } \cos \alpha \, . $$

Podúloha 3

Predtým, ako numericky spočítame čas pádu tyčky, tak musíme najprv spomenúť 2 dôležité aspekty. Prvý aspekt, musíme najprv zistiť, pri akých uhloch sa tyčka vôbec zošmykne (empirickým pozorovaním sme zistili, že niekedy sa tyčka nezošmykne). Druhý aspekt, musíme si overiť, že v našej simulácii sú normálové sily kladné počas celého času.

Na nájdenie kritického uhla $\alpha_k$ si postačíme s už nájdenými vzťahmi. Predpokladajme, že dynamické a statické koeficienty trenia sú rovnaké. Potom hraničný uhol nastane, keď rovnica pre $\ddot \alpha$ nám dá nulové zrýchlenie a stabilnú polohu. V stabilnej polohe nemáme rýchlosť, preto člen $\dot \alpha = 0$ vypadne. Po zopár algebraických úpravách dostaneme

$$ \tan \alpha_k = \frac{1 - f_1 f_2}{2 f_1} \, . $$

Vieme preto, že počiatočný uhol musí byť menší

$$ \alpha (t = 0) < \alpha_k = \arctan \frac{1 - f_1 f_2}{2 f_1} \, . $$

Môžeme si všimnúť, že v prípade $f_1 f_2 > 1$ sa tyčka nikdy nezošmykne a je držaná stenou a podlahou v stabilnej polohe (predstavte si stenu, podlahu a tyčku z gumeného materiálu).

Pokiaľ ide o druhý aspekt, v každom čase simulácie vieme vypočítať normálové sily a overiť, či sú kladné. Môžeme si však všimnúť, že ak zanedbáme členy s $\dot \alpha^2$ (napríklad pri začiatku pádu), tak obe normálové sily budú vždy kladné.² Preto môžeme očakávať, že sa tyčka odlepí od steny/podlahy až v prípade, ak bude dostatočne rýchlo rotovať.

Ďalej si môžeme všimnúť, že pre parametre, keď tyčka môže spadnúť ($f_1 f_2< 1$) je menovateľ kladný a teda nám nehrozí žiadne nefyzikálne stúpanie uhla $\alpha$.

Príklad numerickej simulácie, ktorá využíva Eulerovu integračnú metódu a je implementovaná v jazyku Python, je uvedený nižšie. V skratke, simulácia zadefinuje potrebné parametre systému ($g$, $L$, $f_1$, $f_2$), počiatočné parametre ($\alpha(t = 0)$ a $\dot \alpha(t = 0) = 0$), skontrolujeme, či sa tyčka šmýka, nastavíme si polia, kam si budeme ukladať premenné, vykonáme slučku while dokiaľ uhol neklesne pod nulu, vypočítame v každom čase zmenu uhla, rýchlosti uhla a normálové sily ($F_1$ a $F_2$) a nakoniec výsledky zobrazíme v grafe.

Napríklad pre parametre $L = \SI{1}{\metre}$, $f_1 = 0,2$, $f_2 = 0,1$, $\alpha_0 = \frac{\pi}{3}$ dostaneme čas $\SI{0.685}{\second}$. Môžeme si však všimnúť, že normálová sila $N_2$ tesne pred koncom klesne ku záporným hodnotám. Takže tyčka v spodnej časti trajektórie sa pohybuje dostatočnou rýchlosťou, aby sa ‘’odrazila’’ od podlahy. Pri skúšaní iných parametrov dostávame podobné správanie. Kvalitatívne sa pri páde potenciálna energia tyčky nedokonale (kvôli treniu) premení na rotačnú kinetickú energiu až do bodu, keď na tyčku prestane pôsobiť normálová sila od podlahy. Takže naša lenivá tyčka si tesne pred koncom vyskočí.

Závislosť uhla $\alpha$ (modrá), normálovej sily $F_1$ (oranžová) a $F_2$ (zelená) od času. Nulový uhol je vyznačený prerušovanou červenou čiarou.

import numpy as np
import matplotlib.pyplot as plt

# Constants
g = 9.81  # gravity, m/s^2
l = 1.0  # length, m
f1 = 0.2  # friction coefficient
f2 = 0.1  # friction coefficient

# Initial conditions
alpha_0 = np.pi / 3  # initial angle in radians
dot_alpha_0 = 0  # initial angular velocity

# Check if the initial angle is below the critical angle
critical_angle = np.arctan((1 - f1 * f2) / (2 * f1))
if alpha_0 > critical_angle:
    raise ValueError(
        f"Initial angle {alpha_0} is not below the critical angle {critical_angle:.4f} radians."
    )
# Time step and maximum time
dt = 0.001  # time step, seconds
max_time = 2  # maximum simulation time, seconds

# Initialize variables
alpha = alpha_0
dot_alpha = dot_alpha_0
time = 0
final_time = None

# Lists to store results for plotting
times = []
alphas = []
F1_values = []
F2_values = []

# Euler integration loop
while alpha > 0 and time <= max_time:
    # Calculate second derivative using the given formula
    d2alpha = (
        -(3 * g / l)
        * (
            -2 * f1 * np.sin(alpha)
            + (1 - f1 * f2) * np.cos(alpha)
            + (l / g)
            * dot_alpha**2
            * (f2 * np.cos(alpha) ** 2 + f1 * np.sin(alpha) ** 2)
        )
        / (2 - f1 * f2 + 3 * (f2 - f1) * np.sin(alpha) * np.cos(alpha))
    )

    # Calculate N1 and N2
    F1 = (
        g * (1 + 3 * f2 * np.sin(alpha) * np.cos(alpha) + 3 * np.sin(alpha) ** 2)
        + l * dot_alpha**2 * (-2 * np.sin(alpha) - f2 * np.cos(alpha))
    ) / (2 * (2 - f1 * f2 + 3 * (f2 - f1) * np.sin(alpha) * np.cos(alpha)))

    F2 = (
        g * (f1 - 3 * f1 * np.sin(alpha) ** 2 + 3 * np.sin(alpha) * np.cos(alpha))
        + l * dot_alpha**2 * (-2 * np.cos(alpha) + f1 * np.sin(alpha))
    ) / (2 * (2 - f1 * f2 + 3 * (f2 - f1) * np.sin(alpha) * np.cos(alpha)))

    # Update angular velocity and angle using Euler method
    dot_alpha += d2alpha * dt
    alpha += dot_alpha * dt

    # Append current values to lists
    times.append(time)
    alphas.append(alpha)
    F1_values.append(F1)
    F2_values.append(F2)

    # Increment time
    time += dt

    # Record the final time when the angle becomes 0
    if alpha <= 0 and final_time is None:
        final_time = time

# Print the final time
if final_time is not None:
    print(f"The angle became 0 at time {final_time:.4f} seconds.")
else:
    print("The angle did not reach 0 within the simulation time.")

# Plot the results
plt.figure(figsize=(12, 6))
plt.plot(times, alphas, label="Angle \u03b1 (radians)")
plt.plot(times, F1_values, label="F1")
plt.plot(times, F2_values, label="F2")
plt.axhline(0, color="red", linestyle="--", label="\u03b1 = 0")
plt.title("Time Dependence of Angle \u03b1 and Functions F1, F2")
plt.xlabel("Time (s)")
plt.ylabel("Values")
plt.legend()
plt.grid()
plt.show()