Computando Integrais de Sommerfeld

Integração numérica para calcular integrais de Sommerfeld.

Introdução

Em várias aplicações de eletromagnetismo, especialmente aqueles envolvendo dipolos harmônicos no tempo num plano horizontal estratificado, surgem Integrais de Sommerfeld (em coordenadas cilíndricas) da forma

\[ \int_0^\infty \tilde G(k_\rho ; z|z') J_\nu(k_\rho \rho) k_\rho \, dk_\rho, \quad \nu=0,1,2 \]

onde \(\tilde G\) é a função de Green no domínio espectral, \(z\) e \(z'\) são as coordenadas verticais (perpendicular à estratificação) dos pontos de observação e da fonte, respectivamente, \(\rho\) é a distância radial entre tais pontos, \(J_\nu\) é função de Bessel de primeira espécie, \(k_\rho\) é o número de onda da componente radial dela. Esta integral é, frequentemente, difícil de resolver por não ter solução em fórmula fechada, ser altamente oscilatória, singular e divergente. Para resolver estas integrais, revisaremos o artigo de Krzysztof A. Michalski e Mosig (2016).

A variável \(k_\rho\) é complexa e há pontos sobre o eixo real (ou próximos dele) nos quais o integrando é singular. Portanto, faz-se a integração utilizando um caminho que evite tais pontos, como ilustrado abaixo. A altura (parte imaginária) desse caminho deve ser limitada de forma a evitar que a função de Bessel exploda (tenha valores muito elevados), pois isto pode levar a erros de cancelamento catastrófico.

using Plots
x = range(0, 6, length=100)
xd = 0
y = [xd < xi < pi + xd ? sin(xi - xd) : 0 for xi in x]
plot(x,y, ylims=(-0.2,1.1), xlabel="Re(kρ)", ylabel="Im(kρ)", label="caminho de integração")

Alternativamente, pode-se identificar os pontos no eixo \(\Re(k_\rho)\) que tornam o integrando singular e fazer a integração em intervalos com alguma técnica capaz de lidar com singularidades nos pontos finais e iniciais como, por exemplo, a quadratura de Gauss-Jacobi ou a fórmula tanh-sinh de exponencial dupla (DE). Contudo, computar a cauda da integral (de \(a\) até \(\infty\)) é mais complicado devido à natureza oscilatória de \(\tilde G(k_\rho ; z|z')\), que tem a seguinte forma:

\[ \tilde G(k_\rho ; z|z') \sim \frac{e^{-k_\rho \zeta}}{k_\rho^\delta} [C + \mathcal{O}(k_\rho^{-1})], \quad k_\rho \to \infty \]

onde \(C\) é uma constante, \(\zeta > 0\) é um parâmetro relacionado a \(z\) e \(z'\) (por exemplo, \(\zeta = |z - z'|\) em um meio homogêneo) e \(\mathcal{O}(k_\rho^{-1})\) é um valor assintótico.

Aceleração de Sequências por Extrapolação

Para computar a cauda das integrais de Sommerfeld, é interessante acelerar a convergência de uma série da forma

\[ S_n = \sum_{i=0}^n u_i, \quad n\to \infty \]

A sequência \(S_n\) pode ser particionada como \(S_n = S + r_n\), onde \(r_n\) é o erro de truncamento dela, também chamado de resto. Busca-se uma transformada \(\{S_n\} \to \{S_n'\}\) tal que a convergência seja mais rápida, isto é,

\[ \frac{|r_n'|}{|r_n|} = \mathcal{O}(x_n^{-\mu}), \quad \mu > 0, \quad n \to \infty \]

A convergência destas séries pode ser caracterizada pelo limite

\[ \lambda = \lim_{n\to\infty} \frac{r_{n+1}}{r_n} = \lim_{n\to\infty} \frac{u_{n+1}}{u_n} \]

• Se \(|\lambda| > 1\), a sequência diverge;
• Se \(|\lambda| = 1\), a convergência é linear;
• Se \(|\lambda| < 1\), a convergência é logarítmica;
• Se \(|\lambda| = 0\), a convergência é hiperlinear;

O procedimento de extrapolação mais geral assume que o resto tem uma expansão assintótica

\[ r_n \sim \sum_{j=0}^\infty c_j \, \psi_j(n), \quad n \to \infty \]

onde \(c_j\) são coeficientes desconhecidos e \(\psi_j(n)\) são funções arbitrárias, mas conhecidas de \(n\), com a propriedade de \(\psi_{j+1}=o(\psi_j(n))\). Para muitas sequências encontradas na prática, o resto tem uma expansão do tipo Poincaré da forma

\[ r_n \sim \omega_n \sum_{j=0}^\infty c_j \, x_n^{-j}, \quad n \to \infty \]

onde \(\omega_n\) são as estimativas do resto, que proveem a informação do comportamento de \(r_n\) para \(n\) grande, e \(\{x_n\}\) é uma sequência auxiliar de pontos de interpolação tal que \(x_{n+1}>x_n\). Os melhores aceleradores de convergência são obtidos tanto numericamente, quanto analiticamente, baseados no comportamento assintótico conhecido da sequência.

O método das médias ponderadas e o algorítmo de Mosig-Michalski

Usando médias ponderadas, podemos expressar a soma parcial por

\[ S_n' = \frac{W_n S_n + W_{n+1} S_{n+1}}{W_n + W_{n+1}}, \]

onde \(W_n\) é o peso associado a \(S_n\). Usando o coeficiente de proporção dos pesos \(\eta_n\), esta expressão pode ser reformulada em

\[ \eta_n = -\frac{W_n}{W_{n+1}} \\ S_n' = \frac{S_{n+1} - \eta_n S_n}{1 - \eta_n} \]

Assumindo que os restos \(r_n\) tem expanssão assintótica, sabedo que

\[ \frac{r_n'}{r_n} = \frac{r_{n+1} / r_n - \eta_n}{1 - \eta_n}, \]

então

\[ \frac{r_{n+1}}{r_n} - \frac{\omega_{n+1}}{\omega_n}[1 + \mathcal{O}(x_n^{-2})], \quad n \to \infty, \]

o que sugere escolher a proporção dos pesos como

\[ \eta_n = \frac{\omega_{n+1}}{\omega_n}. \]

Com esta escolha de \(\eta_n\), e assumindo que o denominador das soma parcial \(S_n'\) tem o comportamento assintótio \[ 1 - \eta_n = \mathcal{O}(x_n^{-\sigma}), \quad n \to \infty, \]

onde \(\sigma < 2\). A condição

\[ \frac{|r_n'|}{|r_n|} = \mathcal{O}(x_n^{-\mu}, \quad \mu > 0, n \to \infty \]

é satisfeita para \(\mu = 2 - \sigma\). Tem-se \(\sigma=0\) para todas as sequências que não tem convergência logarítmica e \(\sigma=1\) neste último caso. Assim sendo, chega-se ao esquema recursivo

\[ S_n^{(k+1)} = \frac{S_{n+1}^{(k)} - \eta_n{(k)} S_n{(k)}}{1 - \eta_n{(k)}}, \quad n,k \ge 0, \]

com valor inicial \(S_n^{(0)} = S_n\) e coeficientes

\[ \eta_n^{(k)} = \eta_n \left ( \frac{x_n}{x_{n+1}} \right )^{\mu k} = \frac{\eta_n}{(1 + \Delta x_n / x_n)^{\mu k}} \approx \frac{\eta_n}{(1 + \mu k \Delta x_n / x_n)}, \]

onde a segunda expressão é assintótica para valores grandes de \(n\).

Integrais de Cauda Oscilatórias

Para integração das caudas das integrais de Sommerfeld, combinam-se métodos de quadraturas progressivas, como a regra DE, e os métodos de extrapolação discutidos acima.

Considere uma integral \(S\) feita sobre o eixo real com um integrando \(f(\xi)\) oscilatório, podendo ter convergência lenta ou ser divergente. Calcula-se \(S\) como o limite de uma sequência de somas parciais, conforme a seguir, sendo \(\xi_{-1} = a\). Os pontos \(\xi_i\) selecionados sãos os zeros de \(f(\xi)\).

\[ S = \int_ a^\infty f(\xi)\, d\xi = \lim_{n \to \infty} \sum_{i=0}^n \int_{\xi_{i-1}}^{\xi_i} f(\xi)\, d\xi \]

Esta metodologia de particionar a integral e aplicar extrapolação à série resultante é chamada de Método de Partição-Extrapolação (PE).

Primeiro, investiga-se o comportamento dos restos \(r_n\) para fazer uma extrapolação eficiente.

\[ r_n = S_n - S = - \int_{\xi_n}^\infty f(\xi) \, d\xi \]

Assumindo \(\xi > b > a\), expressa-se \(f(\xi) \approx g(\xi)\, p(\xi)\), onde \(g(\xi)\) tem uma expanssão assintótica

\[ g(\xi) \sim \frac{e^{-\xi \, \zeta}}{\xi^\alpha} \sum_{j = 0}^\infty a_j \xi^{-j}, \quad \xi \to \infty, \]

com o parâmetro \(\zeta \ge 0\), enquanto \(p(\xi)\) é periódica com um meio período \(q\) tal que

\[ p(\xi + q) = - p(\xi) \]

Desta forma, podemos assumir pontos equidistantes \[ \xi_k = b + k q, \quad k = 0, 1, \dots, \]

onde \(b\) é o primeiro zero de \(p(\xi)\). A partir disto, deduz-se

\[ r_n \sim -(-1)^n e^{-nq\zeta} \sum^\infty_{j=0} a_j \int_b^\infty \frac{p(t) e^{-t \zeta}}{(t + nq)^{\alpha + j}} dt. \]

Fazendo uma expansão em série de Taylor em \(t=b\) do denominador do integrando acima, chega-se a uma expansão do tipo Poincaré \[ r_n \sim \omega_n \sum^\infty_{j=0} c_j \xi_n^{-j}, \]

onde \(c_j\) são \(n\) coeficientes independentes e \(\omega_n\) são as estimativas do resto, dadas por \[ \omega_n = (-1)^{n+1} \frac{e^{-nq\zeta}}{\xi_n^\alpha}. \]

Os pontos \(\xi_n\) podem ser substituídos pelos pontos de interpolação \(x_n\) com fator de deslocamento \(\beta = b / q\). Usando as duas expressões acima para computar o limite \(\lambda\), obtendo-se \(\lambda = -e^{-q\zeta}\), o que indica que a sequência \(\{S_n\}\) é alternada e linearmente convergente, portanto todos os métodos do tipo Levin e o algorítimo de Shanks-Wynn são aplicáveis para acelerar a convergência. Usando a estimativa do resto \(\omega_n\) em fórmula fechada, deduzida acima, no algorítimo de Levin-Sidi e Mosig-Michalski, deduz-se a variante-a destes métodos. O parâmetro \(\Omega_k\) no algorítimo de Mosig-Michalski se torna \[ \Omega_k = \frac{\omega_k}{\omega_{k-1}} = -e^{-q\zeta} \left ( \frac{x_{k-1}}{x_k} \right )^\alpha, \quad k \ge 1, \]

sendo \(\Omega_0 = 0\) e os parâmetros do algorítimo \(\sigma = 0\) e \(\mu = 2\).

Como exemplo, considere a integral \[ I_\nu(a,z,\rho) = \int_a^\infty e^{-\xi z} J_\nu (\xi \rho) \xi^\nu d\xi \] \[ I_\nu(a,z,\rho) = \frac{(2 \rho)^\nu \Gamma(\nu + 1/2)}{(z^2 + \rho^2)^{\nu + 1/2} \sqrt{\pi}} - \int_0^a e^{-\xi z} J_\nu (\xi \rho) \xi^\nu d\xi. \]

Para \(\nu = 0,1\) e \(2\), as integrais \(I_\nu(a,z,\rho)\) se assemelham à forna estática das integrais de cauda que aparecem no cômputdo de campos eletromagnéticos in meios estratificados. Com \(z=0\) e \(\nu =1,2\), as integrais do lado esquerdo divergem, mas são definidas no sentido da Soma de Abel (é possível redefinir o somatório como uma série de potências convergente).

Nota: Krzysztof A. Michalski e Mosig (2016) apresentam a fórmula acima usando um duplo fatorial, válida para valores inteiros de \(\nu\). Eu escolhi apresentar a fórmula mais geral usando a Função Gama, assim como apresentado por Gradshteyn e Ryzhik (2007) (página 1089, equação 6.623.1), válida para qualquer valor de \(\nu\).

Tendo em vista a forma assintótica da função de Bessel \[ J_\nu (\xi \rho) \sim \sqrt{\frac{2}{\pi \xi \rho}} \cos(\xi \rho - \nu \pi/2 - \pi/4),\quad z\to\infty, \]

a integral \(I_\nu(a,z,\rho)\) se encaixa no escopo do método PE com \(\zeta=z\), \(\alpha=1/2-\nu\), e \(q=\pi/\rho\). A escolha natural de pontos de quebra \(\xi_i\) para as integrais parciais são, neste caso, os zeros consecutivos de \(J_\nu (\xi \rho)\) (Partição Longman). Contudo, os zeros extados da função de Bessel não tem fórmula fechada, necessitando de processos iterativos. Por isto, outras alternativas para escolher \(\xi_i\) foram propostas, tal como a partição de Lyness, que usa zeros assintóticos baseados na expansão assintótica da função de Bessel.

Há, também, a partição de Sidi, que usa nós equidistantes separados pelo meio período \(q\). A partição de Sidi é a mais simples, porém não confiável para certos valores de \(\rho\) para os quais as somas parciais deixam de ser estritamente alternadas, o que leva a uma degradação do desempenho do método. Contudo, Krzysztof A. Michalski e Mosig (2016) afirmam que este problema não aparece se o limite inferior de integração \(a\) for escolhido como o primeiro zero da função de Bessel que seja maior do que o \(a\) original.

Código

O código a seguir, assim como versões em Fortran e Python, estão disponíveis no Github.

using SpecialFunctions
using DoubleExponentialFormulas

Aceleração de Sequência

"""
Sequence acceleation of an infinite series by extrapolation via the Mosig–Michalski algorithm (t-transformation weighted averages).

# Arguments
- partial_sums: Vector of partial sums of the sequence to be accelerated, where the n-th entry should be the partial sum up to index n.
- u: Transformation parameter controlling the extrapolation's sensitivity (1 for logarithmic convergence, 2 otherwise).

# Returns
- extrapolated_sum: The extrapolated sum, providing an improved estimate for the series true value beyond available terms.
"""
function mosig_michalski(partial_sums, u = 1)
    kmax = length(partial_sums) - 1

    beta = 1
    X = (0:kmax) .+ beta  # equidistant nodes

    # Transformed sequence values at each extrapolation order k
    R = deepcopy(partial_sums)

    # remainder estimates
    wk0 = 1  # w_{k-1}
    wk1 = 1  # w_{k}

    local extrapolated_sum
    for k = 0:(kmax)
        if k >= 1
            wk1 = partial_sums[k + 1] - partial_sums[k]
            Gk = wk1 / wk0
        else
            Gk = 0
        end
        for j = 1:k
            d = X[k - j + 2] - X[k - j + 1]
            eta = Gk / (1 + u * (j - 1) * d / X[k - j + 1])
            R[k - j + 1] = (R[k - j + 2] - eta * R[k - j + 1]) / (1 - eta)
        end
        extrapolated_sum = R[1]
        wk0 = wk1
    end
    return extrapolated_sum
end

Teste da extrapolação com as seguintes sequências, usando \(N=20\) somas parciais.

\[ S_1 = \sum_{i=0}^{\infty} \frac{(-1)^{i}}{\sqrt{i + 1}} \approx 0.604898643421630 \]

\[ S_2 = \sum_{i=0}^{\infty} \frac{(4/5)^{i+1}}{i + 1} = \ln(5) \]

\[ S_3 = \sum_{i=0}^{\infty} \frac{1}{(i + 1)^2} = \frac{\pi^2}{6} \]

let N = 20, u = 1
    partial_sum(n::Int, series_term) = cumsum([series_term(i) for i = 0:(n - 1)])
    
    series_term1(i::Int) = (-1)^i / sqrt(i + 1)
    extrapolated_sum = mosig_michalski(partial_sum(N, series_term1), u)
    @show isapprox(extrapolated_sum, 0.604898643421630)

    series_term2(i::Int) = (4/5)^(i+1) / (i + 1)
    extrapolated_sum = mosig_michalski(partial_sum(N, series_term2), u)
    @show isapprox(extrapolated_sum, log(5))

    series_term3(i::Int) = 1 / (i + 1)^2
    extrapolated_sum = mosig_michalski(partial_sum(N, series_term3), u)
    @show isapprox(extrapolated_sum, pi^2 / 6; atol = 1e-2)
end

isapprox(extrapolated_sum, 0.60489864342163) = true
isapprox(extrapolated_sum, log(5)) = true
isapprox(extrapolated_sum, pi ^ 2 / 6; atol = 0.01) = true

true

Partição-Extrapolação

"""
Numerical integragion of Sommerfeld integral tails with Partition-Extrapolation
via the Mosig-Michalski method.

The expected form of the integral is the following, where ``\\tau`` is the
radius of the circular region that contains nonzero current, ``z`` and ``z'``
are the perpendicular distances to the medium stratification of the field
and the source, respectively, and ``\\rho`` is the lateral distance between
those points.

# Parameters
- fun: The integrand function ``f(x)`` to be integrated over each partition subinterval.
- a: The lower limit of integration (typically the start of the tail, after any head subtraction).
- q: Partition step size; the width of each subinterval (typically matched to Bessel oscillation period or chosen for balance between cost and convergence).
- z: Tail asymptotic decay parameter (often related to physical geometry or stratification; used in remainder estimate).
- alpha: Decay exponent in analytic remainder estimate for the extrapolation transformation.
- tol: Desired relative tolerance for extrapolation stopping criterion.
- kmax: Maximum number of partitions/subintervals for extrapolation.
- u: Transformation parameter controlling the extrapolation's sensitivity (1 for logarithmic convergence, 2 otherwise).

# Returns
- val: Accelerated estimate of the semi-infinite tail integral.
- error_estimate: Estimate of extrapolation absolute error.
"""
function part_extrap(fun::Function, a, q, z, alpha, tol = 1e-6, kmax = 10, u = 1)
    X = zeros(kmax + 2)
    X[1] = a
    partial_sums = zeros(ComplexF64, kmax + 2)
    partial_errors = similar(partial_sums)
    for k = 1 : kmax+1
        X[k + 1] = X[k] + q
        i, e = quadde(fun, X[k], X[k + 1])
        partial_sums[k + 1] = i + partial_sums[k]
        partial_errors[k + 1] = e + partial_errors[k]
    end

    R = partial_sums
    local old0 = 0.0
    local old1 = 0.0
    local val = 0.0
    local error_estimate = 0.0
    for k = 0 : kmax + 1
        if k > 0
            # analytical remainder estimates
            Gk = -exp(-q * z) * (X[k] / X[k + 1])^alpha
        else
            Gk = 0
        end

        # Mosig-Michalski extrapolation
        for j = 1:k
            d = X[k - j + 2] - X[k - j + 1]
            eta = Gk / (1 + u * (j - 1) * d / X[k - j + 1])
            R[k - j + 1] = (R[k - j + 2] - eta * R[k - j + 1]) / (1 - eta)
        end
        val = R[1]

        error_estimate = max(abs(val - old0), abs(val - old1))
        if k > 1 && abs(error_estimate) < tol * abs(val)
            break
        end
        old0 = old1
        old1 = val
    end
    return val, error_estimate
end

Teste

Testamos com a seguinte integral:

\[ I = \int_{a}^{\infty} e^{-x z} \, J_\nu(x p) \, x^\nu \]

com parâmetros \(\nu = 2, p = 1, z = 0, a = 5.13562\)

begin
    v = 2
    p = 1
    a = 5.13562
    q = pi / p
    z = 0
    tol = 1e-9
    alpha = 1/2 - v
    kmax = 10
    u = 2

    fun = x -> exp(-x * z) * besselj(v, x * p) * x^v

    y1 = (2 * p)^v * gamma(v + 1/2) / (sqrt(pi) * (z^2 + p^2)^(v + 1/2))
    y2, E = quadde(fun, 0, a)
    Sn = y1 - y2  # "exact" value of the integral tail

    val, E = part_extrap(fun, a, q, z, alpha, tol, kmax, u)
    @show abs(val - Sn) < tol
end

abs(val - Sn) < tol = true

true

Conclusão

Quando o interesse é computar campos nas zonas distantes (acima de algumas centenas de vezes o comprimento de onda), outros métodos, como o de Descida Mais Íngreme, são mais apropriados do que os apresentados acima; vide K. A. Michalski (1985).

Referências

Gradshteyn, I. S., e I. M. Ryzhik. 2007. Table of Integrals, Series, and Products. Editado por Alan Jeffrey e Daniel Zwillinger. 7.ª ed. New York: Academic Press.

Michalski, K. A. 1985. «On the Efficient Evaluation of Integrals Arising in the Sommerfeld Halfspace Problem». IEE Proceedings H (Microwaves, Antennas and Propagation) 132 (5): 312–18.

Michalski, Krzysztof A., e Juan R. Mosig. 2016. «Efficient Computation of Sommerfeld Integral Tails – Methods and Algorithms». Journal of Electromagnetic Waves and Applications 30 (3): 281–317. https://doi.org/10.1080/09205071.2015.1129915.