Notebook

SageMath para estudantes de física¶

Rogério T. C.

Aula V.1¶

Equação de Klein-Gordon no espaço-tempo de Schwarzschild¶

In [1]:

reset()
%display latex

Espaço-tempo¶

In [2]:

var('m', domain='positive')

Out[2]:

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

In [3]:

M = Manifold(4, 'M', structure='Lorentzian')
SD.<t, r, th, ph> = M.chart(r"t r:(0,+oo) th:(0,pi):\theta ph:(0,2*pi):\phi:periodic")

In [4]:

g = M.metric()
g[0, 0] = - (1 - 2*m/r)
g[1, 1] = 1/(1 - 2*m/r)
g[2, 2] = r^2
g[3, 3] = r^2*sin(th)^2
g.display()

Out[4]:

$\newcommand{\Bold}[1]{\mathbf{#1}}g = \left( \frac{2 \, m}{r} - 1 \right) \mathrm{d} t\otimes \mathrm{d} t + \left( -\frac{1}{\frac{2 \, m}{r} - 1} \right) \mathrm{d} r\otimes \mathrm{d} r + r^{2} \mathrm{d} {\theta}\otimes \mathrm{d} {\theta} + r^{2} \sin\left({\theta}\right)^{2} \mathrm{d} {\phi}\otimes \mathrm{d} {\phi}$

In [5]:

%time nabla = g.connection()

CPU times: user 5.25 s, sys: 160 ms, total: 5.41 s
Wall time: 4.51 s

Campo escalar¶

In [6]:

var ('omega ,k')
R = M.scalar_field(function('R')(r))
S = M.scalar_field(function('S')(th))

In [7]:

Phi = exp(-I*omega*t)*exp(I*k*ph)*R*S

In [8]:

Phi.set_name(name='Phi', latex_name=r'\Phi')

In [9]:

Phi.display()

Out[9]:

$\newcommand{\Bold}[1]{\mathbf{#1}}\begin{array}{llcl} \Phi:& M & \longrightarrow & \mathbb{R} \\ & \left(t, r, {\theta}, {\phi}\right) & \longmapsto & R\left(r\right) S\left({\theta}\right) e^{\left(i \, k {\phi} - i \, \omega t\right)} \end{array}$

Equação $\quad\nabla^\mu \nabla_\mu \Phi = 0$ ¶

In [10]:

%time KG = nabla(nabla(Phi).up(g)).trace()

CPU times: user 11 s, sys: 175 ms, total: 11.1 s
Wall time: 7.71 s

In [11]:

KG.display()

Out[11]:

$\newcommand{\Bold}[1]{\mathbf{#1}}\begin{array}{llcl} & M & \longrightarrow & \mathbb{R} \\ & \left(t, r, {\theta}, {\phi}\right) & \longmapsto & \frac{{\left({\left(2 \, m e^{\left(i \, k {\phi}\right)} - r e^{\left(i \, k {\phi}\right)}\right)} R\left(r\right) \cos\left({\theta}\right) \sin\left({\theta}\right) \frac{\partial\,S}{\partial {\theta}} + {\left(2 \, m e^{\left(i \, k {\phi}\right)} - r e^{\left(i \, k {\phi}\right)}\right)} R\left(r\right) \sin\left({\theta}\right)^{2} \frac{\partial^2\,S}{\partial {\theta} ^ 2} - {\left(\omega^{2} r^{3} R\left(r\right) e^{\left(i \, k {\phi}\right)} + 2 \, {\left(2 \, m^{2} e^{\left(i \, k {\phi}\right)} - 3 \, m r e^{\left(i \, k {\phi}\right)} + r^{2} e^{\left(i \, k {\phi}\right)}\right)} \frac{\partial\,R}{\partial r} + {\left(4 \, m^{2} r e^{\left(i \, k {\phi}\right)} - 4 \, m r^{2} e^{\left(i \, k {\phi}\right)} + r^{3} e^{\left(i \, k {\phi}\right)}\right)} \frac{\partial^2\,R}{\partial r ^ 2}\right)} S\left({\theta}\right) \sin\left({\theta}\right)^{2} - {\left(2 \, k^{2} m e^{\left(i \, k {\phi}\right)} - k^{2} r e^{\left(i \, k {\phi}\right)}\right)} R\left(r\right) S\left({\theta}\right)\right)} e^{\left(-i \, \omega t\right)}}{{\left(2 \, m r^{2} - r^{3}\right)} \sin\left({\theta}\right)^{2}} \end{array}$

In [12]:

KG = (KG/Phi).expr()

In [13]:

KG.expand()

Out[13]:

$\newcommand{\Bold}[1]{\mathbf{#1}}-\frac{\omega^{2} r^{3}}{2 \, m r^{2} - r^{3}} - \frac{4 \, m^{2} r \frac{\partial^{2}}{(\partial r)^{2}}R\left(r\right)}{{\left(2 \, m r^{2} - r^{3}\right)} R\left(r\right)} + \frac{4 \, m r^{2} \frac{\partial^{2}}{(\partial r)^{2}}R\left(r\right)}{{\left(2 \, m r^{2} - r^{3}\right)} R\left(r\right)} - \frac{r^{3} \frac{\partial^{2}}{(\partial r)^{2}}R\left(r\right)}{{\left(2 \, m r^{2} - r^{3}\right)} R\left(r\right)} - \frac{4 \, m^{2} \frac{\partial}{\partial r}R\left(r\right)}{{\left(2 \, m r^{2} - r^{3}\right)} R\left(r\right)} + \frac{6 \, m r \frac{\partial}{\partial r}R\left(r\right)}{{\left(2 \, m r^{2} - r^{3}\right)} R\left(r\right)} - \frac{2 \, r^{2} \frac{\partial}{\partial r}R\left(r\right)}{{\left(2 \, m r^{2} - r^{3}\right)} R\left(r\right)} - \frac{2 \, k^{2} m}{{\left(2 \, m r^{2} - r^{3}\right)} \sin\left({\theta}\right)^{2}} + \frac{k^{2} r}{{\left(2 \, m r^{2} - r^{3}\right)} \sin\left({\theta}\right)^{2}} + \frac{2 \, m \cos\left({\theta}\right) \frac{\partial}{\partial {\theta}}S\left({\theta}\right)}{{\left(2 \, m r^{2} - r^{3}\right)} S\left({\theta}\right) \sin\left({\theta}\right)} - \frac{r \cos\left({\theta}\right) \frac{\partial}{\partial {\theta}}S\left({\theta}\right)}{{\left(2 \, m r^{2} - r^{3}\right)} S\left({\theta}\right) \sin\left({\theta}\right)} + \frac{2 \, m \frac{\partial^{2}}{(\partial {\theta})^{2}}S\left({\theta}\right)}{{\left(2 \, m r^{2} - r^{3}\right)} S\left({\theta}\right)} - \frac{r \frac{\partial^{2}}{(\partial {\theta})^{2}}S\left({\theta}\right)}{{\left(2 \, m r^{2} - r^{3}\right)} S\left({\theta}\right)}$

In [14]:

KG = KG*(2*m*r^2-r^3)/(2*m-r)

In [15]:

KG.expand().combine()

Out[15]:

$\newcommand{\Bold}[1]{\mathbf{#1}}-\frac{\omega^{2} r^{3}}{2 \, m - r} - \frac{4 \, m^{2} r \frac{\partial^{2}}{(\partial r)^{2}}R\left(r\right) - 4 \, m r^{2} \frac{\partial^{2}}{(\partial r)^{2}}R\left(r\right) + r^{3} \frac{\partial^{2}}{(\partial r)^{2}}R\left(r\right) + 4 \, m^{2} \frac{\partial}{\partial r}R\left(r\right) - 6 \, m r \frac{\partial}{\partial r}R\left(r\right) + 2 \, r^{2} \frac{\partial}{\partial r}R\left(r\right)}{{\left(2 \, m - r\right)} R\left(r\right)} + \frac{2 \, m \frac{\partial^{2}}{(\partial {\theta})^{2}}S\left({\theta}\right) - r \frac{\partial^{2}}{(\partial {\theta})^{2}}S\left({\theta}\right)}{{\left(2 \, m - r\right)} S\left({\theta}\right)} - \frac{2 \, k^{2} m - k^{2} r}{{\left(2 \, m - r\right)} \sin\left({\theta}\right)^{2}} + \frac{2 \, m \cos\left({\theta}\right) \frac{\partial}{\partial {\theta}}S\left({\theta}\right) - r \cos\left({\theta}\right) \frac{\partial}{\partial {\theta}}S\left({\theta}\right)}{{\left(2 \, m - r\right)} S\left({\theta}\right) \sin\left({\theta}\right)}$

In [16]:

termos = KG.expand().combine().operands();termos

Out[16]:

$\newcommand{\Bold}[1]{\mathbf{#1}}\left[-\frac{\omega^{2} r^{3}}{2 \, m - r}, -\frac{4 \, m^{2} r \frac{\partial^{2}}{(\partial r)^{2}}R\left(r\right) - 4 \, m r^{2} \frac{\partial^{2}}{(\partial r)^{2}}R\left(r\right) + r^{3} \frac{\partial^{2}}{(\partial r)^{2}}R\left(r\right) + 4 \, m^{2} \frac{\partial}{\partial r}R\left(r\right) - 6 \, m r \frac{\partial}{\partial r}R\left(r\right) + 2 \, r^{2} \frac{\partial}{\partial r}R\left(r\right)}{{\left(2 \, m - r\right)} R\left(r\right)}, \frac{2 \, m \frac{\partial^{2}}{(\partial {\theta})^{2}}S\left({\theta}\right) - r \frac{\partial^{2}}{(\partial {\theta})^{2}}S\left({\theta}\right)}{{\left(2 \, m - r\right)} S\left({\theta}\right)}, -\frac{2 \, k^{2} m - k^{2} r}{{\left(2 \, m - r\right)} \sin\left({\theta}\right)^{2}}, \frac{2 \, m \cos\left({\theta}\right) \frac{\partial}{\partial {\theta}}S\left({\theta}\right) - r \cos\left({\theta}\right) \frac{\partial}{\partial {\theta}}S\left({\theta}\right)}{{\left(2 \, m - r\right)} S\left({\theta}\right) \sin\left({\theta}\right)}\right]$

In [17]:

termos[1]

Out[17]:

$\newcommand{\Bold}[1]{\mathbf{#1}}-\frac{4 \, m^{2} r \frac{\partial^{2}}{(\partial r)^{2}}R\left(r\right) - 4 \, m r^{2} \frac{\partial^{2}}{(\partial r)^{2}}R\left(r\right) + r^{3} \frac{\partial^{2}}{(\partial r)^{2}}R\left(r\right) + 4 \, m^{2} \frac{\partial}{\partial r}R\left(r\right) - 6 \, m r \frac{\partial}{\partial r}R\left(r\right) + 2 \, r^{2} \frac{\partial}{\partial r}R\left(r\right)}{{\left(2 \, m - r\right)} R\left(r\right)}$

In [18]:

gam = var('gamma')
KG_radial = sum(termos[:2]).canonicalize_radical() == gam
KG_angular = sum(termos[2:]).canonicalize_radical() == gam

In [19]:

KG_radial

Out[19]:

$\newcommand{\Bold}[1]{\mathbf{#1}}-\frac{\omega^{2} r^{3} R\left(r\right) + 2 \, {\left(2 \, m^{2} - 3 \, m r + r^{2}\right)} \frac{\partial}{\partial r}R\left(r\right) + {\left(4 \, m^{2} r - 4 \, m r^{2} + r^{3}\right)} \frac{\partial^{2}}{(\partial r)^{2}}R\left(r\right)}{{\left(2 \, m - r\right)} R\left(r\right)} = \gamma$

In [20]:

KG_radial = KG_radial*(2*m-r)*R.expr()-gam*(2*m-r)*R.expr();KG_radial

Out[20]:

$\newcommand{\Bold}[1]{\mathbf{#1}}-\omega^{2} r^{3} R\left(r\right) - \gamma {\left(2 \, m - r\right)} R\left(r\right) - 2 \, {\left(2 \, m^{2} - 3 \, m r + r^{2}\right)} \frac{\partial}{\partial r}R\left(r\right) - {\left(4 \, m^{2} r - 4 \, m r^{2} + r^{3}\right)} \frac{\partial^{2}}{(\partial r)^{2}}R\left(r\right) = 0$

In [21]:

KG_angular

Out[21]:

$\newcommand{\Bold}[1]{\mathbf{#1}}-\frac{k^{2} S\left({\theta}\right) - \cos\left({\theta}\right) \sin\left({\theta}\right) \frac{\partial}{\partial {\theta}}S\left({\theta}\right) - \sin\left({\theta}\right)^{2} \frac{\partial^{2}}{(\partial {\theta})^{2}}S\left({\theta}\right)}{S\left({\theta}\right) \sin\left({\theta}\right)^{2}} = \gamma$

In [22]:

desolve(KG_radial, R.expr(), ivar=r, algorithm='fricas')

Out[22]:

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

Solução numérica¶

In [23]:

DR = function('DR')(r)

In [24]:

KG_rad = [diff(R.expr(),r)-DR==0,KG_radial.subs({diff(R.expr(),r):DR(r), diff(R.expr(),r,2):diff(DR(r),r)})]

In [25]:

KG_rad

Out[25]:

$\newcommand{\Bold}[1]{\mathbf{#1}}\left[-{\rm DR}\left(r\right) + \frac{\partial}{\partial r}R\left(r\right) = 0, -\omega^{2} r^{3} R\left(r\right) - \gamma {\left(2 \, m - r\right)} R\left(r\right) - 2 \, {\left(2 \, m^{2} - 3 \, m r + r^{2}\right)} {\rm DR}\left(r\right) - {\left(4 \, m^{2} r - 4 \, m r^{2} + r^{3}\right)} \frac{\partial}{\partial r}{\rm DR}\left(r\right) = 0\right]$

In [26]:

KG_rad[1] = KG_rad[1].subs(m=0.1, omega=.2, gamma=2.0);KG_rad[1]

Out[26]:

$\newcommand{\Bold}[1]{\mathbf{#1}}-0.0400000000000000 \, r^{3} R\left(r\right) - 2 \, {\left(r^{2} - 0.300000000000000 \, r + 0.0200000000000000\right)} {\rm DR}\left(r\right) + 2.00000000000000 \, {\left(r - 0.200000000000000\right)} R\left(r\right) - {\left(r^{3} - 0.400000000000000 \, r^{2} + 0.0400000000000000 \, r\right)} \frac{\partial}{\partial r}{\rm DR}\left(r\right) = 0$

In [27]:

KG_rad[0] = solve(KG_rad[0],diff(R.expr(),r))[0].right()
KG_rad[1] = solve(KG_rad[1],diff(DR(r),r))[0].right()

In [28]:

KG_rad

Out[28]:

$\newcommand{\Bold}[1]{\mathbf{#1}}\left[{\rm DR}\left(r\right), -\frac{{\left(50 \, r^{2} - 15 \, r + 1\right)} {\rm DR}\left(r\right) + {\left(r^{3} - 50 \, r + 10\right)} R\left(r\right)}{25 \, r^{3} - 10 \, r^{2} + r}\right]$

In [29]:

var('x,y')

Out[29]:

$\newcommand{\Bold}[1]{\mathbf{#1}}\left(x, y\right)$

In [30]:

KG_rad[0] = KG_rad[0].subs({DR(r):y})
KG_rad[1] = KG_rad[1].subs({DR(r):y, R.expr():x})

In [31]:

KG_rad

Out[31]:

$\newcommand{\Bold}[1]{\mathbf{#1}}\left[y, -\frac{{\left(r^{3} - 50 \, r + 10\right)} x + {\left(50 \, r^{2} - 15 \, r + 1\right)} y}{25 \, r^{3} - 10 \, r^{2} + r}\right]$

In [32]:

radsol = desolve_system_rk4(KG_rad, [x,y], ics =[0.3 ,1 ,0.5], ivar=r, end_points=100, step=0.01)

In [33]:

len(radsol)

Out[33]:

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

In [34]:

points =[[ i , j ] for i ,j,k  in radsol ]

In [35]:

list_plot(points, axes_labels =['$r$' , '$R$'], color='darkred', legend_label = 'Solução numérica',
           gridlines=True, frame=True, axes=False)

Out[35]:

Variando $\omega$ entre 0,1 e 2,5¶

In [36]:

pl = []
for om in srange(0.1,2.6,.3):
    KG_rad = [diff(R.expr(),r)-DR==0,KG_radial.subs({diff(R.expr(),r):DR(r), diff(R.expr(),r,2):diff(DR(r),r)})]
    KG_rad[0] = solve(KG_rad[0],diff(R.expr(),r))[0].right().subs({DR(r):y})
    KG_rad[1] = solve(KG_rad[1].subs(m=0.1, omega=om, gamma=2.0),diff(DR(r),r))[0].right().subs({DR(r):y, R.expr():x})
    radsol = desolve_system_rk4 (KG_rad, [x,y], ics =[0.3 ,1 ,0.5], ivar=r, end_points=100, step=0.01)
    pl.append(list_plot([[i,j] for i,j,k in radsol], axes_labels =['$r$' , '$R$'], color='darkred', 
              legend_label = r'$\omega='+str(om.n(digits=3))+'$', gridlines=True, frame=True, axes=False))

In [37]:

ga = graphics_array(pl, nrows=3)
ga.show(figsize=[15,12])

In [ ]: