O métoddo de disparo é uma técnica numérica para resolver equações diferenciais com condições de contorno. Embora nem sempre seja o mais estável, sua implementação relativamente simples o torna uma ferramenta didática valiosa para problemas de autovalores em mecânica quântica, como as equações de Schrödinger.

1. Oscilador Harmônico Unidimensional

Equação Adimensional

A equação de Schrödinger para o oscilador harmônico, após a transformação t = x/x₀ e E' = E/(ℏω/2), torna-se:

Código de Implementação

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint
from scipy.optimize import fsolve

# Equação do oscilador harmônico: Phi''(t) = (t^2 - E' ) Phi(t)
def osc_deriv(y, t, energy_param):
    return [y[1], (t * t - energy_param) * y[0]]

boundary_max = 10.0
num_points = 1000

def osc_shoot(energy_val):
    time_array = np.linspace(boundary_max, 0, num_points)
    # Condições de contorno em t=boundary_max: Phi~1e-10, Phi'=0
    initial_conditions = [1E-10, 0.0]
    solution = odeint(osc_deriv, initial_conditions, time_array, args=(energy_val,))
    # Para paridade ímpar, o alvo é Phi(0) ~ 0
    return solution[-1, 0]

# Encontrar o autovalor para o primeiro estado excitado ímpar
eigen_value = fsolve(osc_shoot, [-0.5])
print(f"Autovalor E' encontrado: {eigen_value[0]:.8f}")

# Calcular e plotar a autofunção correspondente
time_plot = np.linspace(boundary_max, 0, num_points)
wave_function = odeint(osc_deriv, [1E-10, 0], time_plot, args=(eigen_value[0],))

plt.figure()
plt.plot(time_plot, wave_function[:, 0])
plt.xlabel('t')
plt.ylabel('Phi(t)')
plt.title(f'Oscilador Harmônico - Estado de Paridade Ímpar (E\'={eigen_value[0]:.2f})')
plt.grid(True)
plt.show()

A estratégia consiste em integrar numericamente a partir de um ponto no inifnito (t=10) onde a função de onda deve decair, na direção da origem. Ajusta-se o parâmetro E' até que a condição de contorno na origem seja satisfeita (Φ(0)≈0 para paridade ímpar). Integrra da direita para a esquerda suprime a solução divergente não física (∼e+t²/2).

2. Átomo de Hidrogênio

Equação Radial

Após separação de variáveis e transformação para a variável adimensional ρ = κr, a equação radial para o potencial de Coulomb resulta em:

Código de Implementação

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint
from scipy.optimize import fsolve

# Equação radial do Hidrogênio: u''(ρ) = (1 - 2*ρ0/ρ + l(l+1)/ρ²) * u(ρ)
def h_deriv(y, rho, rho0_param, angular_momentum):
    return [y[1], (1.0 - 2.0 * rho0_param / rho + angular_momentum * (angular_momentum + 1) / rho**2) * y[0]]

rho_limit = 50.0
num_grid = 1000
l_val = 2  # Momento angular

def h_shoot(rho0_trial):
    # Integra de rho=rho_limit para rho~0
    rho_array = np.linspace(rho_limit, 1e-3, num_grid)
    # Condições iniciais: u(rho_limit) pequena, u'(rho_limit)=0
    sol = odeint(h_deriv, [1e-10, 0.0], rho_array, args=(rho0_trial, l_val))
    # Condição de contorno: u(rho->0) deve ser zero (R(r) finita)
    return sol[-1, 0]

# Encontrar rho0, que define o número quântico principal n
rho0_found = fsolve(h_shoot, [2.5])
n_principal = rho0_found[0]
energy_ev = -13.6 / n_principal**2
print(f"rho0 = {n_principal:.8f}, n ~ {n_principal:.2f}, E = {energy_ev:.4f} eV")

# Gerar o gráfico da autofunção radial
rho_plot = np.linspace(rho_limit, 1e-3, num_grid)
u_solution = odeint(h_deriv, [1e-10, 0.0], rho_plot, args=(rho0_found[0], l_val))

plt.figure()
plt.plot(rho_plot, u_solution[:, 0])
plt.xlabel('ρ (unidades adimensionais)')
plt.ylabel('u(ρ)')
plt.title(f'Átomo de Hidrogênio: Função Radial para l={l_val}')
plt.grid(True)
plt.show()

3. Potencial de Woods-Saxon

Modelo e Parâmetros

Para núcleos, um potencial mais realista é o de Woods-Saxon. A forma central é:

Código de Implementação

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint
from scipy.optimize import fsolve

# Parâmetros do potencial de Woods-Saxon
V0_central = -53.0
r0_central = 1.25
a0_central = 0.65
V0_so = 22.0
r0_so = 1.25
a0_so = 0.65

def woods_saxon_potential(r, l_quantum, j_total, mass_number):
    R_central = r0_central * mass_number**(1.0/3.0)
    R_so = r0_so * mass_number**(1.0/3.0)
    
    # Termo central
    exp_c = np.exp((r - R_central) / a0_central)
    f_central = 1.0 / (1.0 + exp_c)
    V_central = V0_central * f_central
    
    # Termo spin-órbita (aproximação simplificada)
    exp_so = np.exp((r - R_so) / a0_so)
    denominator = (1.0 + exp_so)**2
    derivative_term = - (1.0 / a0_so) * exp_so / denominator
    V_so = V0_so * r0_so**2 * (1.0 / r) * derivative_term * (j_total*(j_total+1) - l_quantum*(l_quantum+1) - 0.75) / 2.0
    
    return V_central + V_so

def ws_deriv(y, rho, k_wave, l_quantum, j_total, mass_number):
    # rho é a variável adimensional, r = rho / k_wave
    r_physical = rho / k_wave
    potential = woods_saxon_potential(r_physical, l_quantum, j_total, mass_number)
    # fator 19.4: ℏ²/(2m) em MeV·fm² ≈ 19.4
    term = 1.0 + potential / (19.4 * k_wave**2) + l_quantum*(l_quantum+1) / rho**2
    return [y[1], term * y[0]]

def ws_shoot(k_trial, l_quantum, j_total, mass_number):
    rho_max = 20.0
    rho_array = np.linspace(rho_max, 0.02, 1000)
    sol = odeint(ws_deriv, [1e-10, 0.0], rho_array, args=(k_trial, l_quantum, j_total, mass_number))
    return sol[-1, 0]

def find_eigenvalue(k_guess, l_quantum, j_total, mass_number):
    return fsolve(ws_shoot, k_guess, args=(l_quantum, j_total, mass_number))

# Exemplo: Splitting dos orbitais g para diferentes números de massa
mass_numbers = np.arange(60, 201, 5)
l_orbital = 4

energies_j92 = []
energies_j72 = []

for A in mass_numbers:
    # Encontrar k para j=9/2
    k92 = find_eigenvalue(0.8, l_orbital, 4.5, A)
    E92 = -19.4 * k92**2
    energies_j92.append(E92[0])
    
    # Encontrar k para j=7/2
    k72 = find_eigenvalue(0.7, l_orbital, 3.5, A)
    E72 = -19.4 * k72**2
    energies_j72.append(E72[0])

plt.figure(figsize=(8, 5))
plt.plot(mass_numbers, energies_j92, 'o-', label='0g₉/₂')
plt.plot(mass_numbers, energies_j72, 's-', label='0g₇/₂')
plt.xlabel('Número de Massa (A)')
plt.ylabel('Energia de Ligação (MeV)')
plt.title('Splitting Spin-Órbita no Potencial de Woods-Saxon')
plt.legend()
plt.grid(True)
plt.show()

O procedimento numérico é análogo ao do hidrogênio, mas com um potencial mais complexo e a busca por parâmetros de onda k. A dependência com o número de massa A revela o splitting spin-órbita característico da estrutura nuclear.