Órbitas e Trajetórias: Equações Numéricas na Física

Em 1705, o astrônomo inglês Edmond Halley cruzou registros de cometas observados em 1531, 1607 e 1682 e chegou a uma conclusão ousada: os três podiam ser o mesmo objeto, retornando ao Sol a cada cerca de 76 anos. Usando as leis de Kepler — publicadas quase um século antes — e a mecânica recém-formalizada por Newton, Halley calculou a órbita e previu o retorno do cometa para 1758. Ele não viveu para conferir: morreu em 1742. Mas o cometa voltou, pontualmente, dezesseis anos depois — e ganhou o nome de quem acertou a previsão (veja o verbete Halley’s Comet).

A órbita do Cometa Halley é uma elipse extremamente alongada: excentricidade \(e \approx 0{,}967\) (a da Terra é \(e \approx 0{,}017\), quase um círculo perfeito). Para saber onde o cometa está numa data qualquer, é preciso resolver a equação de Kepler:

onde \(M\) é a anomalia média (proporcional ao tempo decorrido) e \(E\) é a anomalia excêntrica — o ângulo que efetivamente localiza o cometa na elipse. Essa equação mistura \(E\) livre com \(\sin(E)\): não existe manipulação algébrica capaz de isolar \(E\), e nenhuma fórmula fechada resolve o problema em geral. Chegou a hora de resolvê-la de verdade, com números.

Por Que Órbitas e Trajetórias Importam?

#

• Navegação aeroespacial: software de bordo de satélites e sondas resolve a equação de Kepler (ou variantes dela) continuamente, para saber a posição exata a cada instante.
• Balística: calcular o alcance de um projétil sob resistência do ar real — não o modelo ideal sem atrito do ensino médio — é um problema clássico de engenharia militar e, mais tarde, de esportes como golfe, beisebol e atletismo.
• Determinação de órbitas: da mecânica celeste clássica ao rastreamento moderno de satélites e detritos espaciais, o mesmo tipo de equação aparece.
• Motores de física em jogos e simulações: qualquer simulação de corpo em queda livre com arrasto, ou em órbita, resolve por baixo dos panos o mesmo tipo de problema deste artigo.

Os métodos usados aqui — newton e brentq — são exatamente os apresentados em Resolvendo Equações Numericamente com SciPy, o primeiro artigo deste bloco. Lá, foram aplicados a exemplos didáticos, escolhidos para ilustrar o comportamento de cada método. Aqui, os mesmos métodos resolvem dois problemas físicos completos, com números reais.

Resolvendo a Equação de Kepler

#

scipy.optimize.newton resolve a equação de Kepler rapidamente quando fornecemos a derivada analítica, \(1 - e\cos(E)\) — convergência quadrática, a mesma vantagem já vista no artigo anterior.

import numpy as np
from scipy import optimize

def kepler_residual(E, M, ecc):
    return E - ecc * np.sin(E) - M

def kepler_residual_prime(E, M, ecc):
    return 1 - ecc * np.cos(E)

ecc = 0.967  # excentricidade do Cometa Halley
M_val = 1.0
E_sol = optimize.newton(
    kepler_residual, x0=M_val, fprime=kepler_residual_prime, args=(M_val, ecc)
)
print("E (rad):", E_sol)
print("checagem do resíduo:", kepler_residual(E_sol, M_val, ecc))

E (rad): 1.9114369764896801
checagem do resíduo: 0.0

O Que É, Afinal, a Anomalia Excêntrica?

#

Até aqui \(E\) foi só “o ângulo que localiza o cometa” — o suficiente para resolver a equação, mas não para entender o que se resolveu. E as fórmulas que vêm a seguir, \(x = a(\cos E - e)\) e \(y = a\sqrt{1-e^2}\sin E\), parecem tiradas do nada. Vale gastar um parágrafo com a construção geométrica, porque ela torna as duas fórmulas óbvias.

Desenhe, em volta da elipse, o círculo auxiliar: o círculo de raio \(a\) (o semi-eixo maior) centrado no centro da elipse. Agora, dado um ponto \(P\) do planeta sobre a elipse, suba na vertical até encontrar o círculo, num ponto \(Q\).

A figura a seguir monta essa construção numa órbita de excentricidade moderada (\(e = 0{,}6\), escolhida para o desenho ficar legível): o círculo auxiliar tracejado, a projeção vertical de \(Q\) para \(P\), e os dois ângulos marcados em seus respectivos vértices — \(E\) no centro, \(\nu\) no foco.

Com isso no lugar, a posição cartesiana do cometa a partir de \(E\) (com o Sol na origem, semi-eixo maior \(a\) normalizado e excentricidade \(e\)) é \(x = a(\cos E - e)\), \(y = a\sqrt{1-e^2}\sin E\). O gráfico a seguir mostra dois painéis lado a lado: à esquerda, a própria equação de Kepler como um problema de raiz — a curva \(g(E) = E - e\sin(E)\) cruzando a reta horizontal \(y=M\); à direita, o que isso significa fisicamente — a órbita elíptica, o Sol num dos focos, e o cometa na posição correspondente ao \(E\) encontrado.

import matplotlib.pyplot as plt
import numpy as np
from scipy import optimize


def kepler_residual(E, M, ecc):
    return E - ecc * np.sin(E) - M


def kepler_residual_prime(E, M, ecc):
    return 1 - ecc * np.cos(E)


ecc = 0.967  # excentricidade do Cometa Halley
M_val = 1.0
E_sol = optimize.newton(
    kepler_residual, x0=M_val, fprime=kepler_residual_prime, args=(M_val, ecc)
)

a_orbit = 1.0  # semi-eixo maior normalizado

E_range = np.linspace(0, 2 * np.pi, 400)
g_E = E_range - ecc * np.sin(E_range)

fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(13, 6))

# painel esquerdo: a equação de Kepler como raiz
ax1.plot(E_range, g_E, color="C0", label=r"$g(E) = E - e\sin(E)$")
ax1.axhline(M_val, color="C1", linestyle="--", label=f"M = {M_val}")
ax1.scatter([E_sol], [M_val], color="C2", zorder=3, s=80, label=f"E ≈ {E_sol:.3f} rad")
ax1.set_xlabel("E (rad)")
ax1.set_ylabel("g(E)")
ax1.set_title("Equação de Kepler como problema de raiz")
ax1.legend()

# painel direito: a órbita
E_orbit = np.linspace(0, 2 * np.pi, 400)
x_orbit = a_orbit * (np.cos(E_orbit) - ecc)
y_orbit = a_orbit * np.sqrt(1 - ecc**2) * np.sin(E_orbit)
x_planeta = a_orbit * (np.cos(E_sol) - ecc)
y_planeta = a_orbit * np.sqrt(1 - ecc**2) * np.sin(E_sol)

ax2.plot(x_orbit, y_orbit, color="C0")
ax2.scatter([0], [0], color="C1", s=200, marker="*", zorder=3, label="Sol (foco)")
ax2.scatter([x_planeta], [y_planeta], color="C2", s=100, zorder=3, label="cometa")
ax2.plot([0, x_planeta], [0, y_planeta], color="gray", linestyle=":", linewidth=1)
ax2.set_aspect("equal")
ax2.set_xlabel("x")
ax2.set_ylabel("y")
ax2.set_title(f"Órbita elíptica (e = {ecc})")
ax2.legend()

A 2ª Lei de Kepler, Finalmente Visível

#

Até aqui, usamos a órbita elíptica só como forma geométrica, para localizar um único ponto nela. Falta mostrar a outra metade do que Kepler descobriu: não é só a forma da órbita que é diferente de um círculo — a velocidade do cometa ao longo dela também varia, de um jeito específico:

Até aqui, resolvemos a equação de Kepler para um único instante. Mas a anomalia média \(M\) foi definida, por construção, para crescer linearmente com o tempo — é exatamente por isso que ela é fácil de calcular a partir de qualquer data. Isso significa que amostrar \(M\) em passos iguais equivale a amostrar o tempo em passos iguais. Resolvendo a equação de Kepler para vários \(M\) igualmente espaçados ao longo de uma órbita completa, obtemos as posições do cometa em intervalos de tempo iguais — e a 2ª lei fica visível: os pontos deveriam se aglomerar perto do periélio e ficar espaçados perto do afélio.

import matplotlib.pyplot as plt
import numpy as np
from scipy import optimize


def kepler_residual(E, M, ecc):
    return E - ecc * np.sin(E) - M


def kepler_residual_prime(E, M, ecc):
    return 1 - ecc * np.cos(E)


ecc = 0.967  # excentricidade do Cometa Halley
a_orbit = 1.0  # semi-eixo maior normalizado

E_orbit = np.linspace(0, 2 * np.pi, 400)
x_orbit = a_orbit * (np.cos(E_orbit) - ecc)
y_orbit = a_orbit * np.sqrt(1 - ecc**2) * np.sin(E_orbit)


def posicoes_orbita(M_amostras, ecc):
    E_amostras = np.array([
        optimize.newton(kepler_residual, x0=M_i, fprime=kepler_residual_prime, args=(M_i, ecc))
        for M_i in M_amostras
    ])
    x = a_orbit * (np.cos(E_amostras) - ecc)
    y = a_orbit * np.sqrt(1 - ecc**2) * np.sin(E_amostras)
    return x, y


n_pontos = 16
M_amostras = np.linspace(0, 2 * np.pi, n_pontos, endpoint=False)
x_amostras, y_amostras = posicoes_orbita(M_amostras, ecc)

fig, ax = plt.subplots(figsize=(7, 7.5))
ax.plot(x_orbit, y_orbit, color="C0", zorder=1)
ax.scatter([0], [0], color="C1", s=200, marker="*", zorder=3, label="Sol (foco)")
ax.scatter(x_amostras, y_amostras, color="C2", s=40, zorder=3, label="posições (Δt iguais)")
for i in range(n_pontos):
    ax.plot([0, x_amostras[i]], [0, y_amostras[i]], color="gray", linewidth=0.7, zorder=0)
ax.set_aspect("equal")
ax.set_xlabel("x")
ax.set_ylabel("y")
ax.set_title("2ª lei de Kepler: pontos mais densos perto do periélio")
ax.legend(loc="upper center", bbox_to_anchor=(0.5, -0.12), ncol=2)

Visualmente já dá para notar a aglomeração perto do periélio (a região da elipse próxima ao Sol, à direita no gráfico).

Mas essa figura mostra os raios, e a lei fala de áreas. Vale, então, desenhar as áreas de fato: três setores correspondentes ao mesmo intervalo de tempo (um décimo do período, cada), tomados em regiões bem diferentes da órbita.

import numpy as np
from scipy import optimize


def kepler_residual(E, M, ecc):
    return E - ecc * np.sin(E) - M


def kepler_residual_prime(E, M, ecc):
    return 1 - ecc * np.cos(E)


ecc = 0.967  # excentricidade do Cometa Halley
a_orbit = 1.0  # semi-eixo maior normalizado


def posicoes_orbita(M_amostras, ecc):
    E_amostras = np.array([
        optimize.newton(kepler_residual, x0=M_i, fprime=kepler_residual_prime, args=(M_i, ecc))
        for M_i in M_amostras
    ])
    x = a_orbit * (np.cos(E_amostras) - ecc)
    y = a_orbit * np.sqrt(1 - ecc**2) * np.sin(E_amostras)
    return x, y


delta_M = 2 * np.pi / 10  # um décimo do período
setores = {"perto do periélio": 0.0,
           "no meio do caminho": 2 * np.pi * 0.20,
           "perto do afélio": 2 * np.pi * 0.45}

for rotulo, M0 in setores.items():
    M_setor = np.linspace(M0, M0 + delta_M, 300)
    xs_s, ys_s = posicoes_orbita(M_setor, ecc)
    xs = np.concatenate(([0.0], xs_s))   # fecha o polígono no Sol
    ys = np.concatenate(([0.0], ys_s))
    area = 0.5 * abs(np.dot(xs, np.roll(ys, -1)) - np.dot(ys, np.roll(xs, -1)))
    print(f"{rotulo:>20}: área = {area:.5f}")

   perto do periélio: área = 0.08002
  no meio do caminho: área = 0.08004
     perto do afélio: área = 0.08004

Iguais até a quarta casa decimal. E o mais impressionante é o contraste entre os formatos:

Perto do periélio, o cometa percorre um arco curto mas a distância ao Sol é pequena — o setor é gordo e baixo. Perto do afélio, ele mal se move, mas o raio é enorme — o setor vira uma lasca comprida e finíssima. Os dois efeitos se compensam exatamente, e é isso que a 2ª lei afirma.

Medir essas áreas com poucas amostras, porém, reserva uma pegadinha.

import numpy as np
from scipy import optimize


def kepler_residual(E, M, ecc):
    return E - ecc * np.sin(E) - M


def kepler_residual_prime(E, M, ecc):
    return 1 - ecc * np.cos(E)


ecc = 0.967  # excentricidade do Cometa Halley
a_orbit = 1.0  # semi-eixo maior normalizado


def posicoes_orbita(M_amostras, ecc):
    E_amostras = np.array([
        optimize.newton(kepler_residual, x0=M_i, fprime=kepler_residual_prime, args=(M_i, ecc))
        for M_i in M_amostras
    ])
    x = a_orbit * (np.cos(E_amostras) - ecc)
    y = a_orbit * np.sqrt(1 - ecc**2) * np.sin(E_amostras)
    return x, y


def areas_das_fatias(x, y):
    x_prox, y_prox = np.roll(x, -1), np.roll(y, -1)
    return 0.5 * np.abs(x * y_prox - x_prox * y)


for n_teste in [16, 100, 1000, 5000]:
    M_teste = np.linspace(0, 2 * np.pi, n_teste, endpoint=False)
    x_teste, y_teste = posicoes_orbita(M_teste, ecc)
    areas_teste = areas_das_fatias(x_teste, y_teste)
    var_rel = (areas_teste.max() - areas_teste.min()) / areas_teste.mean()
    print(f"n = {n_teste:5d} pontos: variação relativa das áreas = {var_rel:.2%}")

n =    16 pontos: variação relativa das áreas = 104.29%
n =   100 pontos: variação relativa das áreas = 69.60%
n =  1000 pontos: variação relativa das áreas = 12.45%
n =  5000 pontos: variação relativa das áreas = 0.72%

Perto do Periélio Extremo: Newton Fica Mais Lento

#

Newton com derivada converge rápido para excentricidades moderadas, mas a convergência piora conforme \(e \to 1\) — órbitas cada vez mais alongadas, o regime dos cometas de longo período. Resolvendo a mesma equação de Kepler para \(M=1\) rad, variando apenas a excentricidade, e contando quantas iterações Newton precisa até o resíduo cair abaixo de \(10^{-12}\):

import matplotlib.pyplot as plt
import numpy as np


def kepler_residual(E, M, ecc):
    return E - ecc * np.sin(E) - M


def kepler_residual_prime(E, M, ecc):
    return 1 - ecc * np.cos(E)


def newton_conta_iteracoes(f, fp, x0, args, tol=1e-12, max_iter=100):
    x = x0
    for i in range(1, max_iter + 1):
        fx = f(x, *args)
        if abs(fx) < tol:
            return i
        x = x - fx / fp(x, *args)
    return max_iter  # não convergiu nesse limite


ecc = 0.967  # excentricidade do Cometa Halley
M_val = 1.0

eccs = np.linspace(0.1, 0.99, 40)
iteracoes = [
    newton_conta_iteracoes(kepler_residual, kepler_residual_prime, M_val, args=(M_val, e_i))
    for e_i in eccs
]

iter_halley = newton_conta_iteracoes(
    kepler_residual, kepler_residual_prime, M_val, args=(M_val, ecc)
)
print(f"Halley (e={ecc}): {iter_halley} iterações")

fig, ax = plt.subplots(figsize=(8, 5))
ax.plot(eccs, iteracoes, "o-", color="C0")
ax.scatter([ecc], [iter_halley], color="C1", zorder=3, s=80, label=f"Halley (e={ecc})")
ax.legend()
ax.set_xlabel("excentricidade (e)")
ax.set_ylabel("iterações até convergir")
ax.set_title("Newton fica mais lento perto de órbitas quase parabólicas")

Halley (e=0.967): 7 iterações

A relação não é suave: o número de iterações sobe em degraus (de 4 para 5 em \(e \approx 0{,}23\), de 5 para 6 em \(e \approx 0{,}55\), de 6 para 7 em \(e \approx 0{,}92\)), porque é uma contagem inteira de passos até cruzar a tolerância — não uma quantidade contínua. Mesmo assim, a tendência é clara: o Cometa Halley, com \(e \approx 0{,}967\), já está no patamar mais lento observado nessa faixa. Órbitas ainda mais alongadas — cometas de período muito longo, com \(e\) acima de \(0{,}99\) — exigiriam ainda mais iterações, ou um chute inicial melhor que \(M\) para acelerar a convergência.

Por Que Kepler Precisou de Marte, e Não da Terra

#

O mesmo código responde uma pergunta histórica que costuma passar batido: por que Kepler levou anos analisando dados de Marte, e não do planeta que tinha logo abaixo dos pés?

import numpy as np
from scipy import optimize


def kepler_residual(E, M, ecc):
    return E - ecc * np.sin(E) - M


def kepler_residual_prime(E, M, ecc):
    return 1 - ecc * np.cos(E)


def newton_conta_iteracoes(f, fp, x0, args, tol=1e-12, max_iter=100):
    x = x0
    for i in range(1, max_iter + 1):
        fx = f(x, *args)
        if abs(fx) < tol:
            return i
        x = x - fx / fp(x, *args)
    return max_iter  # não convergiu nesse limite


for nome, e_i in [("Terra", 0.0167), ("Marte", 0.0934), ("Halley", 0.967)]:
    E_i = optimize.newton(kepler_residual, x0=1.0, fprime=kepler_residual_prime,
                          args=(1.0, e_i))
    iteracoes_i = newton_conta_iteracoes(kepler_residual, kepler_residual_prime,
                                         1.0, args=(1.0, e_i))
    print(f"{nome:>7}  e={e_i:.4f}  b/a={np.sqrt(1 - e_i**2):.4f}  "
          f"E={E_i:.4f} rad  ({iteracoes_i} iterações)")

  Terra  e=0.0167  b/a=0.9999  E=1.0142 rad  (3 iterações)
  Marte  e=0.0934  b/a=0.9956  E=1.0825 rad  (4 iterações)
 Halley  e=0.9670  b/a=0.2548  E=1.9114 rad  (7 iterações)

A órbita da Terra tem \(b/a = 0{,}9999\): ela difere de um círculo perfeito em 0,014%. Nenhum instrumento do século XVI detectaria esse achatamento — desenhada em papel, a órbita da Terra é um círculo. Marte, com \(e = 0{,}0934\), ainda é quase circular (\(b/a = 0{,}9956\)), mas tem uma vantagem decisiva: como o Sol está no foco, e não no centro, a distância Sol-Marte varia \(\pm 9{,}3\%\) ao longo do ano marciano — contra \(\pm 1{,}7\%\) na Terra.

Era essa variação, e não o achatamento, que aparecia nos dados de Tycho Brahe. Kepler passou a maior parte de uma década no que chamou de “guerra contra Marte” justamente porque Marte era o planeta com dados bons cuja órbita destoava o suficiente de um círculo para denunciar a elipse. Com a Terra, ou com Vênus (\(e = 0{,}007\)), a descoberta simplesmente não estaria nos números.

Repare também na coluna de iterações: 3 para a Terra, 7 para o Halley. É a mesma tendência que a próxima seção mede em detalhe.

Projétil com Resistência do Ar

#

O segundo problema deixa a mecânica celeste e volta à Terra: um projétil lançado com velocidade inicial \(v_0\) e ângulo \(\theta\), sob arrasto proporcional à velocidade (\(F_{arrasto} = -kmv\)). A altura em função do tempo tem termo exponencial:

Achar o tempo de voo — o instante em que a altura volta a zero — é resolver uma equação transcendental em \(t\).

E há uma armadilha antes disso, exatamente do tipo que o artigo anterior ensinou a procurar: a altura vale zero em dois instantes — no lançamento e na queda. \(t = 0\) é uma raiz legítima da equação, e sem cuidado é ela que o solver encontra. Por isso o bracket usado a seguir começa em \(0{,}1\) s, e não em zero: ele exclui deliberadamente a raiz trivial.

import numpy as np
from scipy import optimize


g = 9.81
k_drag = 0.2  # coeficiente de arrasto / massa
v0 = 30.0  # m/s
angle_deg = 45.0
v0x = v0 * np.cos(np.radians(angle_deg))
v0y = v0 * np.sin(np.radians(angle_deg))


def altura(t):
    return (v0y + g / k_drag) / k_drag * (1 - np.exp(-k_drag * t)) - g * t / k_drag


t_voo = optimize.brentq(altura, 0.1, 20)
print("tempo de voo (com resistência):", t_voo, "s")

tempo de voo (com resistência): 3.838446938469727 s

O gráfico a seguir mostra a altura ao longo do tempo, com a raiz — o tempo de voo — marcada onde a curva volta a cruzar zero.

import matplotlib.pyplot as plt
import numpy as np
from scipy import optimize


g = 9.81
k_drag = 0.2  # coeficiente de arrasto / massa
v0 = 30.0  # m/s
angle_deg = 45.0
v0x = v0 * np.cos(np.radians(angle_deg))
v0y = v0 * np.sin(np.radians(angle_deg))


def altura(t):
    return (v0y + g / k_drag) / k_drag * (1 - np.exp(-k_drag * t)) - g * t / k_drag


t_voo = optimize.brentq(altura, 0.1, 20)

t_range = np.linspace(0, t_voo * 1.2, 300)
fig, ax = plt.subplots(figsize=(8, 5))
ax.plot(t_range, altura(t_range), color="C0")
ax.axhline(0, color="gray", linewidth=1, zorder=0)
ax.scatter([t_voo], [0], color="C2", zorder=3, label="tempo de voo")
ax.legend()
ax.set_xlabel("t (s)")
ax.set_ylabel("altura (m)")
ax.set_title("Projétil com resistência do ar linear")

Quanto a Resistência do Ar Realmente Muda a Trajetória?

#

O caso ideal (sem arrasto) tem trajetória e tempo de voo em forma fechada: \(t_{voo}^{ideal} = 2v_{0y}/g\). Sobrepondo as duas trajetórias completas — posição \(x(t)\) contra \(y(t)\), não só altura contra tempo — fica visível o quanto o arrasto reduz o alcance.

import matplotlib.pyplot as plt
import numpy as np
from scipy import optimize


g = 9.81
k_drag = 0.2  # coeficiente de arrasto / massa
v0 = 30.0  # m/s
angle_deg = 45.0
v0x = v0 * np.cos(np.radians(angle_deg))
v0y = v0 * np.sin(np.radians(angle_deg))


def altura(t):
    return (v0y + g / k_drag) / k_drag * (1 - np.exp(-k_drag * t)) - g * t / k_drag


t_voo = optimize.brentq(altura, 0.1, 20)

t_voo_ideal = 2 * v0y / g
t_range_ideal = np.linspace(0, t_voo_ideal, 300)
x_ideal = v0x * t_range_ideal
y_ideal = v0y * t_range_ideal - 0.5 * g * t_range_ideal**2

t_range_resist = np.linspace(0, t_voo, 300)
x_resist = v0x / k_drag * (1 - np.exp(-k_drag * t_range_resist))
y_resist = altura(t_range_resist)

fig, ax = plt.subplots(figsize=(9, 5))
ax.plot(x_ideal, y_ideal, color="C0", label="sem resistência do ar (fórmula fechada)")
ax.plot(x_resist, y_resist, color="C1", label="com resistência do ar (numérico)")
ax.scatter([x_ideal[-1]], [0], color="C0", zorder=3)
ax.scatter([x_resist[-1]], [0], color="C1", zorder=3)
ax.axhline(0, color="gray", linewidth=1, zorder=0)
ax.set_ylim(0, y_ideal.max() * 1.4)  # espaço no topo para a legenda não sobrepor as curvas
ax.legend(loc="upper right")
ax.set_xlabel("x (m)")
ax.set_ylabel("y (m)")
ax.set_title("Trajetória: com × sem resistência do ar")

import numpy as np
from scipy import optimize


g = 9.81
k_drag = 0.2  # coeficiente de arrasto / massa
v0 = 30.0  # m/s
angle_deg = 45.0
v0x = v0 * np.cos(np.radians(angle_deg))
v0y = v0 * np.sin(np.radians(angle_deg))


def altura(t):
    return (v0y + g / k_drag) / k_drag * (1 - np.exp(-k_drag * t)) - g * t / k_drag


t_voo = optimize.brentq(altura, 0.1, 20)

t_voo_ideal = 2 * v0y / g
x_alcance_ideal = v0x * t_voo_ideal
x_alcance_resist = v0x / k_drag * (1 - np.exp(-k_drag * t_voo))

print(f"alcance ideal: {x_alcance_ideal:.1f} m (t_voo = {t_voo_ideal:.2f} s)")
print(f"alcance com resistência: {x_alcance_resist:.1f} m (t_voo = {t_voo:.2f} s)")

alcance ideal: 91.7 m (t_voo = 4.32 s)
alcance com resistência: 56.8 m (t_voo = 3.84 s)

Com \(k = 0{,}2\), o arrasto derruba o alcance em cerca de 38% — de \(91{,}7\) m para \(56{,}8\) m — e reduz o tempo de voo em meio segundo. A lição não é o número exato (que depende do modelo linear, como avisado acima), e sim a ordem de grandeza do erro: ignorar a resistência do ar não produz um desvio de alguns por cento, produz um desvio de dezenas por cento. É por isso que a parábola do ensino médio, por mais elegante que seja, não serve para nenhuma aplicação prática de balística — e por que vale a pena aprender a resolver o caso resistido, mesmo quando ele custa uma raiz numérica em vez de uma fórmula.

Quanto o Coeficiente de Arrasto Importa?

#

k_drag = 0{,}2 foi um valor único, escolhido para o exemplo acima. Esse coeficiente depende de forma do projétil, área frontal, densidade do ar — e, crucialmente, da massa: o \(k\) que aparece nas fórmulas é a razão entre o coeficiente de arrasto bruto e a massa do projétil, \(k = b/m\), com unidade de s⁻¹. É por isso que uma bola de tênis e uma bala de canhão do mesmo tamanho se comportam de formas tão diferentes: mesma área frontal, mesmo \(b\), massas muito distintas — e portanto \(k\) muito distintos. Repetindo o mesmo brentq para uma faixa de valores de \(k\), dá para ver como o alcance e o tempo de voo se comportam sem precisar de nenhuma ferramenta nova.

import matplotlib.pyplot as plt
import numpy as np
from scipy import optimize


g = 9.81
k_drag = 0.2  # coeficiente de arrasto / massa
v0 = 30.0  # m/s
angle_deg = 45.0
v0x = v0 * np.cos(np.radians(angle_deg))
v0y = v0 * np.sin(np.radians(angle_deg))


def altura(t):
    return (v0y + g / k_drag) / k_drag * (1 - np.exp(-k_drag * t)) - g * t / k_drag


t_voo = optimize.brentq(altura, 0.1, 20)

x_alcance_resist = v0x / k_drag * (1 - np.exp(-k_drag * t_voo))


def altura_com_k(t, k):
    return (v0y + g / k) / k * (1 - np.exp(-k * t)) - g * t / k


def alcance_com_k(k):
    t_voo_k = optimize.brentq(altura_com_k, 0.1, 20, args=(k,))
    x_voo_k = v0x / k * (1 - np.exp(-k * t_voo_k))
    return x_voo_k, t_voo_k


k_range = np.linspace(0.02, 0.6, 40)
resultados = [alcance_com_k(k_i) for k_i in k_range]
alcances = [r[0] for r in resultados]
tempos_voo = [r[1] for r in resultados]

fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(12, 5))
ax1.plot(k_range, alcances, color="C0")
ax1.scatter([k_drag], [x_alcance_resist], color="C2", zorder=3, label=f"k = {k_drag}")
ax1.set_xlabel("coeficiente de arrasto (k)")
ax1.set_ylabel("alcance (m)")
ax1.set_title("Alcance x arrasto")
ax1.legend()

ax2.plot(k_range, tempos_voo, color="C1")
ax2.scatter([k_drag], [t_voo], color="C2", zorder=3, label=f"k = {k_drag}")
ax2.set_xlabel("coeficiente de arrasto (k)")
ax2.set_ylabel("tempo de voo (s)")
ax2.set_title("Tempo de voo x arrasto")
ax2.legend()

Nos extremos testados, o alcance cai de \(86{,}7\) m (\(k = 0{,}02\), arrasto quase desprezível) para \(30{,}5\) m (\(k = 0{,}6\), arrasto forte) — menos de um terço. A queda é suave e monótona: não há surpresa qualitativa aqui, ao contrário do que aconteceu com a excentricidade de Kepler (onde apareceram degraus). Faz sentido: mais arrasto sempre significa menos alcance, sem comportamento especial de limiar.

Tabela-Resumo

#

Próximos Passos

#

Os dois problemas deste artigo vieram da física — um da mecânica celeste, outro da balística — e cada um usou um único cenário numérico real (o Cometa Halley, um lançamento a \(45°\)) em vez de exemplos didáticos isolados. Mas a física não é o único domínio onde equações sem forma fechada aparecem o tempo todo: a química e a engenharia têm as suas próprias, com armadilhas diferentes — cúbicas com múltiplas raízes físicas, aproximações de sala de aula que deixam de valer, e equações implícitas que engenheiros resolviam à mão décadas antes de existir scipy.optimize.

No próximo artigo, Da Química à Engenharia, resolvemos três dessas equações reais — o volume molar de um gás real, o pH de um ácido fraco e o fator de atrito de Colebrook — todas com as mesmas ferramentas apresentadas neste bloco.