Notebook

SageMath para físicos¶

Rogério T. C.

Aula II.1¶

O problema de Kepler¶

Na mecânica clássica, o problema de Kepler é um caso especial do problema de dois corpos, no qual os dois corpos interagem por uma força central que varia em intensidade com o inverso do quadrado da distância $\mathbf{r}$ entre eles. O problema consiste em encontrar a posição dos dois corpos ao longo do tempo, dadas suas massas, posições e velocidades.

Vamos considerar o caso em que a posição de uma das massas coincide com a do centro de massa (como no sistema sol+planeta). Por se tratar um problema de força central, o movimento é confinado em um plano. Isso nos permite descrevê-lo em coordenadas polares.

O potencial e a energia cinética são dados respectivamente por

$V(r) = -\frac{GMm}{r} \quad \text{e} \quad T = \frac{1}{2}m\left(\dot{r}^2+r^2\dot{\phi}^2\right).$

In [1]:

%display latex
reset()

Variáveis dinâmicas

In [2]:

var('t m G M', domain='positive')
phi(t) = function('phi')(t)
r(t) = function('r')(t)

In [3]:

assume(r(t)>0)

In [4]:

assumptions()

Out[4]:

$\newcommand{\Bold}[1]{\mathbf{#1}}\left[t > 0, m > 0, G > 0, M > 0, r\left(t\right) > 0\right]$

Vetor posição

In [5]:

R = vector([r(t)*cos(phi(t)),r(t)*sin(phi(t))]);R

Out[5]:

$\newcommand{\Bold}[1]{\mathbf{#1}}\left(\cos\left(\phi\left(t\right)\right) r\left(t\right),\,r\left(t\right) \sin\left(\phi\left(t\right)\right)\right)$

In [6]:

R.norm().simplify_trig().canonicalize_radical()

Out[6]:

$\newcommand{\Bold}[1]{\mathbf{#1}}r\left(t\right)$

Potencial

In [7]:

V = -G*M*m/(r(t));V

Out[7]:

$\newcommand{\Bold}[1]{\mathbf{#1}}-\frac{G M m}{r\left(t\right)}$

Energia cinética

In [8]:

T = m/2*diff(R,t)*diff(R,t);T

Out[8]:

$\newcommand{\Bold}[1]{\mathbf{#1}}\frac{1}{2} \, {\left(\cos\left(\phi\left(t\right)\right) r\left(t\right) \frac{\partial}{\partial t}\phi\left(t\right) + \sin\left(\phi\left(t\right)\right) \frac{\partial}{\partial t}r\left(t\right)\right)}^{2} m + \frac{1}{2} \, {\left(r\left(t\right) \sin\left(\phi\left(t\right)\right) \frac{\partial}{\partial t}\phi\left(t\right) - \cos\left(\phi\left(t\right)\right) \frac{\partial}{\partial t}r\left(t\right)\right)}^{2} m$

In [9]:

T = T.simplify_trig();T

Out[9]:

$\newcommand{\Bold}[1]{\mathbf{#1}}\frac{1}{2} \, m r\left(t\right)^{2} \frac{\partial}{\partial t}\phi\left(t\right)^{2} + \frac{1}{2} \, m \frac{\partial}{\partial t}r\left(t\right)^{2}$

Lagrangiana

In [10]:

L = T-V;L

Out[10]:

$\newcommand{\Bold}[1]{\mathbf{#1}}\frac{1}{2} \, m r\left(t\right)^{2} \frac{\partial}{\partial t}\phi\left(t\right)^{2} + \frac{1}{2} \, m \frac{\partial}{\partial t}r\left(t\right)^{2} + \frac{G M m}{r\left(t\right)}$

Equações de Euler-Lagrange

$\displaystyle {\frac {\partial L}{\partial q_{i}}}(t,{\boldsymbol {q}}(t),{\dot {\boldsymbol {q}}}(t))-{\frac {\mathrm {d} }{\mathrm {d} t}}{\frac {\partial L}{\partial {\dot {q}}_{i}}}(t,{\boldsymbol {q}}(t),{\dot {\boldsymbol {q}}}(t))=0$

$\frac {\partial L}{\partial q_{i}}(t,{\boldsymbol {q}}(t),{\dot {\boldsymbol {q}}}(t))=0 \Rightarrow {\frac {\partial L}{\partial {\dot {q}}_{i}}}(t,{\boldsymbol {q}}(t),{\dot {\boldsymbol {q}}}(t))=C$

Usaremos as quantidades conservadas no lugar das equações de Euler-Lagrange

Variáveis simbólicas para energia e momento angular

In [11]:

var('E', domain='real')
var('l', latex_name=r'\ell', domain='positive')

Out[11]:

$\newcommand{\Bold}[1]{\mathbf{#1}}{\ell}$

Energia conservada

In [12]:

eq_E = E == T+V;eq_E

Out[12]:

$\newcommand{\Bold}[1]{\mathbf{#1}}E = \frac{1}{2} \, m r\left(t\right)^{2} \frac{\partial}{\partial t}\phi\left(t\right)^{2} + \frac{1}{2} \, m \frac{\partial}{\partial t}r\left(t\right)^{2} - \frac{G M m}{r\left(t\right)}$

In [13]:

def formal_derivative(f, x):
    r"""
    The formal derivative of `f` with respect to the symbolic function `x`.
    """
    tempX = SR.symbol()
    return f.subs({x: tempX}).diff(tempX).subs({tempX: x})

In [14]:

formal_derivative(L,phi)

Out[14]:

$\newcommand{\Bold}[1]{\mathbf{#1}}0$

Momento angular conservado

In [15]:

eq_phi = formal_derivative(L,diff(phi,t)) ==l;eq_phi

Out[15]:

$\newcommand{\Bold}[1]{\mathbf{#1}}m r\left(t\right)^{2} \frac{\partial}{\partial t}\phi\left(t\right) = {\ell}$

In [16]:

eq_phi1 = solve(eq_phi,diff(phi(t),t))[0];eq_phi1

Out[16]:

$\newcommand{\Bold}[1]{\mathbf{#1}}\frac{\partial}{\partial t}\phi\left(t\right) = \frac{{\ell}}{m r\left(t\right)^{2}}$

Substituindo o momento angular na energia conservada

In [17]:

eq_E1 = eq_E.subs(eq_phi1);eq_E1

Out[17]:

$\newcommand{\Bold}[1]{\mathbf{#1}}E = \frac{1}{2} \, m \frac{\partial}{\partial t}r\left(t\right)^{2} - \frac{G M m}{r\left(t\right)} + \frac{{\ell}^{2}}{2 \, m r\left(t\right)^{2}}$

In [18]:

eq_r = solve(eq_E1,diff(r(t),t))[1];eq_r

Out[18]:

$\newcommand{\Bold}[1]{\mathbf{#1}}\frac{\partial}{\partial t}r\left(t\right) = \frac{\sqrt{2 \, G M m^{2} r\left(t\right) + 2 \, E m r\left(t\right)^{2} - {\ell}^{2}}}{m r\left(t\right)}$

Pela regra da cadeia,

$\frac{d\phi}{dr}=\frac{\frac{d\phi}{dt}}{\frac{dr}{dt}}$

$\displaystyle \phi = \int \frac{d\phi}{dt}\left(\frac{dr}{dt}\right)^{-1}dr$

Usaremos a variável simbólica $\varphi$ no lugar da função $\phi$ e uma variável simbólica para a função $r(t)$ .

In [19]:

vr = var('vr', latex_name='r', domain='positive')
vp = var('vp', latex_name=r'\varphi')

In [20]:

vr

Out[20]:

$\newcommand{\Bold}[1]{\mathbf{#1}}{r}$

$\int \frac{d\phi}{dt}\left(\frac{dr}{dt}\right)^{-1}dr$

In [21]:

phrt = integral((eq_phi1.rhs()/eq_r.rhs()).subs({r(t):vr}),vr, hold=True);phrt

Out[21]:

$\newcommand{\Bold}[1]{\mathbf{#1}}\int \frac{{\ell}}{\sqrt{2 \, G M m^{2} {r} + 2 \, E m {r}^{2} - {\ell}^{2}} {r}}\,{d {r}}$

Calculando a integral

In [22]:

assume(G^2*M^2*m^3+2*E*l^2>0)
phrt.unhold()

Out[22]:

$\newcommand{\Bold}[1]{\mathbf{#1}}-\arcsin\left(-\frac{G M m^{2}}{\sqrt{G^{2} M^{2} m^{4} + 2 \, E {\ell}^{2} m}} + \frac{{\ell}^{2}}{\sqrt{G^{2} M^{2} m^{4} + 2 \, E {\ell}^{2} m} {r}}\right)$

Solução em termos de $\varphi$

In [23]:

eq_vp = vp ==phrt.unhold();eq_vp

Out[23]:

$\newcommand{\Bold}[1]{\mathbf{#1}}{\varphi} = -\arcsin\left(-\frac{G M m^{2}}{\sqrt{G^{2} M^{2} m^{4} + 2 \, E {\ell}^{2} m}} + \frac{{\ell}^{2}}{\sqrt{G^{2} M^{2} m^{4} + 2 \, E {\ell}^{2} m} {r}}\right)$

In [24]:

sol_vp = sin(eq_vp.lhs()) == sin(eq_vp.rhs());sol_vp

Out[24]:

$\newcommand{\Bold}[1]{\mathbf{#1}}\sin\left({\varphi}\right) = \sin\left(-\arcsin\left(-\frac{G M m^{2}}{\sqrt{G^{2} M^{2} m^{4} + 2 \, E {\ell}^{2} m}} + \frac{{\ell}^{2}}{\sqrt{G^{2} M^{2} m^{4} + 2 \, E {\ell}^{2} m} {r}}\right)\right)$

In [25]:

sol_vp = sol_vp.trig_simplify();sol_vp

Out[25]:

$\newcommand{\Bold}[1]{\mathbf{#1}}\sin\left({\varphi}\right) = \frac{G M m^{2} {r} - {\ell}^{2}}{\sqrt{G^{2} M^{2} m^{4} + 2 \, E {\ell}^{2} m} {r}}$

Isolando o $r$

In [26]:

sol_r = solve(sol_vp,vr)[0];sol_r

Out[26]:

$\newcommand{\Bold}[1]{\mathbf{#1}}{r} = \frac{{\ell}^{2}}{G M m^{2} - \sqrt{G^{2} M^{2} m^{4} + 2 \, E {\ell}^{2} m} \sin\left({\varphi}\right)}$

Solução de livro texto em termos da excentricidade $\epsilon$ da órbita

In [27]:

ep = var('epsilon', domain='real');ep

Out[27]:

$\newcommand{\Bold}[1]{\mathbf{#1}}\epsilon$

In [28]:

eq_ep = ep == sqrt(1+2*l^2*E/(m*(G*m*M)^2));eq_ep

Out[28]:

$\newcommand{\Bold}[1]{\mathbf{#1}}\epsilon = \sqrt{\frac{2 \, E {\ell}^{2}}{G^{2} M^{2} m^{3}} + 1}$

Para incluir a excentricidade na solução, isolamos $E$ na expressão da excentricidade

In [29]:

assume(ep>0)
sub_E = solve(eq_ep,E)[0]
forget(ep>0)
sub_E

Out[29]:

$\newcommand{\Bold}[1]{\mathbf{#1}}E = \frac{{\left(G^{2} M^{2} \epsilon^{2} - G^{2} M^{2}\right)} m^{3}}{2 \, {\ell}^{2}}$

e substituimos na solução

In [30]:

sol_r = sol_r.subs(sub_E).canonicalize_radical().factor();sol_r

Out[30]:

$\newcommand{\Bold}[1]{\mathbf{#1}}{r} = -\frac{{\ell}^{2}}{{\left(\epsilon \sin\left({\varphi}\right) - 1\right)} G M m^{2}}$

Faremos $\varphi \to \varphi - \pi/2$ pra deixar com a cara mais usual nos livros texto.

In [31]:

sol_r = sol_r.subs(vp=vp-pi/2).expand_trig();sol_r

Out[31]:

$\newcommand{\Bold}[1]{\mathbf{#1}}{r} = \frac{{\ell}^{2}}{{\left(\epsilon \cos\left({\varphi}\right) + 1\right)} G M m^{2}}$

Deixando os dois lados sem dimensão

In [32]:

sol_r = sol_r*G*M*m^2/l^2;sol_r

Out[32]:

$\newcommand{\Bold}[1]{\mathbf{#1}}\frac{G M m^{2} {r}}{{\ell}^{2}} = \frac{1}{\epsilon \cos\left({\varphi}\right) + 1}$

Vejamos os gráficos das órbitas em coordenadas polares (polar=True)

In [33]:

orbita_0 = plot(sol_r.rhs().subs(ep==0),(vp,0,2*pi),polar=True,aspect_ratio=1, color='red', axes=False, frame=True, title=r'Órbita circular ($\epsilon=0$)')

In [34]:

orbita_0

Out[34]:

In [35]:

centro = circle((0,0),.01, fill=True, color='black', axes=False)

In [36]:

centro+orbita_0

Out[36]:

In [37]:

orbita_1 = plot(sol_r.rhs().subs(ep==1),(vp,-.07*pi,.07*pi), polar=True, color='red', axes=False, frame=True, title=r'Órbita parabólica ($\epsilon=1$)')

In [38]:

orbita_1

Out[38]:

In [39]:

orbita_m = plot(sol_r.rhs().subs(ep==.8),(vp,-pi,pi),aspect_ratio=1, polar=True, color='red', axes=False, frame=True, title=r'Órbita elíptica ($\epsilon<1$)')

In [40]:

centro + orbita_m

Out[40]:

In [41]:

orbita_M = plot(sol_r.rhs().subs(ep==2),(vp,-sqrt(2)/3,sqrt(2)/3), polar=True, color='red', axes=False, frame=True, title=r'Órbita hiperbólica ($\epsilon<1$)')

In [42]:

orbita_M

Out[42]:

In [43]:

orbita_0+centro+orbita_m+orbita_M+orbita_1

Out[43]:

Calculando a Constante na Lei de Deslocamento de Wien¶

A lei do deslocamento de Wien é a lei da física que relaciona o comprimento de onda do pico de emissão de radiação eletromagnética de corpo negro e sua temperatura:

$\lambda _{\text{max}}={\frac {b}{T}}$ onde

$\displaystyle \lambda _{\text{max}}$ é o comprimento de onda do pico de intensidade da radiação eletromagnética;
$T\,$ é a temperatura absoluta do corpo negro;
$b\,$ é a constante de proporcionalidade, chamada constante de dispersão de Wien.

Vamos calcular a constante $b$ usando a lei de Planck para a radiação de um corpo negro

$u(\lambda ,T)={2hc^{2} \over \lambda ^{5}}{1 \over e^{hc/\lambda kT}-1}.$

Constantes		Unidade
$h$	constante de Planck	J $\cdot$ s
$c$	velocidade da luz no vácuo	m $\cdot$ s $^{-1}$
$e$	número de Euler	sem dimensão
$k$	constante de Boltzmann	J $\cdot$ K $^{-1}$

In [44]:

var('h c k', domain='positive')
la = var('la', latex_name = r'\lambda', domain='positive')
lam = var('lam', latex_name = r'\lambda_{\text{max}}', domain='positive')
u(la ,T) = 2*h*c^2/(la^5)*1/(e^(h*c/(k*la*T))-1)

In [45]:

u(la,T)

Out[45]:

$\newcommand{\Bold}[1]{\mathbf{#1}}\frac{2 \, c^{2} h}{{\lambda}^{5} {\left(e^{\left(\frac{c h}{T k {\lambda}}\right)} - 1\right)}}$

In [46]:

u0 = diff(u(lam,T),lam)==0;u0

Out[46]:

$\newcommand{\Bold}[1]{\mathbf{#1}}-\frac{10 \, c^{2} h}{{\lambda_{\text{max}}}^{6} {\left(e^{\left(\frac{c h}{T k {\lambda_{\text{max}}}}\right)} - 1\right)}} + \frac{2 \, c^{3} h^{2} e^{\left(\frac{c h}{T k {\lambda_{\text{max}}}}\right)}}{T k {\lambda_{\text{max}}}^{7} {\left(e^{\left(\frac{c h}{T k {\lambda_{\text{max}}}}\right)} - 1\right)}^{2}} = 0$

In [47]:

solve(u0,lam)

Out[47]:

$\newcommand{\Bold}[1]{\mathbf{#1}}\left[{\lambda_{\text{max}}} = \frac{c h e^{\left(\frac{c h}{T k {\lambda_{\text{max}}}}\right)}}{5 \, {\left(T k e^{\left(\frac{c h}{T k {\lambda_{\text{max}}}}\right)} - T k\right)}}\right]$

In [48]:

u1 = (u0*lam^6).expand();u1

Out[48]:

$\newcommand{\Bold}[1]{\mathbf{#1}}-\frac{10 \, c^{2} h}{e^{\left(\frac{c h}{T k {\lambda_{\text{max}}}}\right)} - 1} + \frac{2 \, c^{3} h^{2} e^{\left(\frac{c h}{T k {\lambda_{\text{max}}}}\right)}}{T k {\lambda_{\text{max}}} {\left(e^{\left(\frac{2 \, c h}{T k {\lambda_{\text{max}}}}\right)} - 2 \, e^{\left(\frac{c h}{T k {\lambda_{\text{max}}}}\right)} + 1\right)}} = 0$

In [49]:

eqx = x == c*h/(lam*T*k);eqx

Out[49]:

$\newcommand{\Bold}[1]{\mathbf{#1}}x = \frac{c h}{T k {\lambda_{\text{max}}}}$

In [50]:

eql = solve(eqx,lam)[0];eql

Out[50]:

$\newcommand{\Bold}[1]{\mathbf{#1}}{\lambda_{\text{max}}} = \frac{c h}{T k x}$

In [51]:

u2 = u1.subs(eql);u2

Out[51]:

$\newcommand{\Bold}[1]{\mathbf{#1}}\frac{2 \, c^{2} h x e^{x}}{e^{\left(2 \, x\right)} - 2 \, e^{x} + 1} - \frac{10 \, c^{2} h}{e^{x} - 1} = 0$

In [52]:

eqx2 = solve(u2,x)[0];eqx2

Out[52]:

$\newcommand{\Bold}[1]{\mathbf{#1}}x = 5 \, {\left(e^{x} - 1\right)} e^{\left(-x\right)}$

Exercício¶

Encontre uma solução numérica para a equação

$x = 5 \, {\left(e^{x} - 1\right)} e^{\left(-x\right)}$

** find_root(equação, min, max)

In [53]:

plot(eqx2,x,(x,0,10))

Out[53]:

In [54]:

eqx3 = x == find_root(eqx2,4,6);eqx3

Out[54]:

$\newcommand{\Bold}[1]{\mathbf{#1}}x = 4.96511423174416$

In [55]:

eql2 = (eql*T).subs(eqx3);eql2

Out[55]:

$\newcommand{\Bold}[1]{\mathbf{#1}}T {\lambda_{\text{max}}} = \frac{0.20140523527264692 \, c h}{k}$

In [56]:

from scipy.constants import c as nc
from scipy.constants import h as nh
from scipy.constants import k as nk

In [57]:

constantes = {c:nc, k:nk, h:nh}

In [58]:

constantes

Out[58]:

$\newcommand{\Bold}[1]{\mathbf{#1}}\left\{c : 299792458.0, k : 1.380649 \times 10^{-23}, h : 6.62607015 \times 10^{-34}\right\}$

In [59]:

eql2.subs(constantes)

Out[59]:

$\newcommand{\Bold}[1]{\mathbf{#1}}T {\lambda_{\text{max}}} = 0.0028977719551852403$

Logo,

$b \approx 2,89777\, \text{mm}\cdot\text{K}$

In [60]:

u(500/(10^9),4500).subs(constantes)

Out[60]:

$\newcommand{\Bold}[1]{\mathbf{#1}}6377929856000.104$

In [61]:

plot(u(la/(10^9),4500).subs(constantes),0,3000, frame=True, axes=False)

Out[61]:

In [62]:

plot(u(la/(10^9),4500).subs(constantes),0,3000, color=colors.darkgreen,axes=False, frame=True, gridlines=True, axes_labels=[r'$\lambda\, (nm)$',r'$u(\lambda)$'], legend_label='$T=4500$')

Out[62]:

In [63]:

lista_u = [u(la/(10^9),k).subs(constantes) for k in [3500, 4000, 4500, 5000, 5500]]

In [64]:

lista_u

Out[64]:

$\newcommand{\Bold}[1]{\mathbf{#1}}\left[\frac{1.1910429723971882 \times 10^{29}}{{\lambda}^{5} {\left(e^{\frac{4110.791078582668}{{\lambda}}} - 1\right)}}, \frac{1.1910429723971882 \times 10^{29}}{{\lambda}^{5} {\left(e^{\frac{3596.9421937598345}{{\lambda}}} - 1\right)}}, \frac{1.1910429723971882 \times 10^{29}}{{\lambda}^{5} {\left(e^{\frac{3197.2819500087417}{{\lambda}}} - 1\right)}}, \frac{1.1910429723971882 \times 10^{29}}{{\lambda}^{5} {\left(e^{\frac{2877.553755007867}{{\lambda}}} - 1\right)}}, \frac{1.1910429723971882 \times 10^{29}}{{\lambda}^{5} {\left(e^{\frac{2615.9579590980616}{{\lambda}}} - 1\right)}}\right]$

In [65]:

plot(lista_u,0,3000, axes=False, frame=True, gridlines=True, axes_labels=[r'$\lambda\, (nm)$',r'$u(\lambda)$'])

Out[65]: