Úloha 7

Na hmotný bod o hmotnosti $m=1$ kg působí síla, která ho přitahuje do počátku soustavy souřadnic, a současně síla, která ho odpuzuje. Přitažlivá síla $F_{+}$ závisí na vzdálenosti od počátku jako $F_{+}(r)=h/r^2$, odpudivá síla $F_{-}$ klesá se vzdáleností jako $F_{-}(r)=k/r^6$. Při pohybu je hmotný bod zpomalován tlumící silou $F_o$ úměrnou rychlosti $lv$ jako $F_o(v)=lv$. Najděte rovnovážnou polohu $r_0$ hmotného bodu (tj. vzdálenost od středu ve které je výslednice sil působících na hmotný bod nulová) pokud je $h = 1$ Nm${^2}$, $k = 10$ ${\rm Nm^6}$ a $l = 0.1$ Ns/m. Nakreslete trajektorii hmotného bodu, který se v čase $t=0$ nachází v klidu ve vzdálenosti $r = 1$ m od počátku.

import numpy as np
import matplotlib.pyplot as plt

síla $F_+$

def Fp(r):
    return -1/r**2

síla $F_-$

def Fm(r):
    return 10/r**6

síla $F_o$

def Fo(v):
    return -0.1*v

r_plot=np.arange(-5,5,0.001)

Nakresleme si nejdříve jak vypadá průběh síly $F_+$ a $F_-$

plt.plot(r_plot,Fp(r_plot),c='blue',label="$F_{+}$")
plt.plot(r_plot,Fm(r_plot),c='red',label="$F_{-}$")
plt.ylim(-2,2)
plt.xlim(-5,5)
plt.plot([-5,5],[0,0],ls='dashed',c='grey')
plt.plot([0,0],[-5,5],ls='dashed',c='grey')
plt.legend()
plt.show()

počáteční podmínky

m=1 #hmotnost
dt=0.001 #krok v čase
tmax=400 #max. čas
t=np.arange(0,tmax,dt) #pole časů
n=np.size(t)
x=np.empty(n) #pole souřadnic
vx=np.empty(n) #pole rychlostí
x[0]=1 #v čase t = 0 je x = 0
r=np.abs(x[0]) #vzdálenost od počátku
vx[0]=(Fp(r)+Fm(r))/m*dt/2 #v čase t = 0 je vx = 0, jako první prvek pole je rychlost v čase dt/2

numerické řešení pohybové rovnice
$\ddot{x}=-\frac{h}{x^2}+\frac{k}{x^6}-l\dot{x}\frac{\dot{x}}{|\dot{x}|}$

for i in range(n-1):
    x[i+1]=x[i]+vx[i]*dt
    r=np.abs(x[i])
    v=np.abs(vx[i])         
    vx[i+1]=vx[i]+(Fp(r)+Fm(r)+Fo(v)*vx[i]/np.abs(vx[i]))/m*dt

nakreslíme is časový průběh souřadnice

plt.plot(t,x,c='red')
r_equilibrium=(10)**(1/4)
plt.xlim(0,tmax)
plt.plot([0,400],[r_equilibrium,r_equilibrium],ls='dashed',c='grey')
plt.xlabel('čas (s)')
plt.ylabel('poloha (m)')
plt.show()

a časový průběh rychlosti

plt.plot(t,vx,c='blue')
plt.xlim(0,tmax)
plt.plot([0,400],[0,0],ls='dashed',c='grey')
plt.xlabel('čas (s)')
plt.ylabel('rychlost (m/s)')
plt.show()