Skip to content

numeric_integral.md

Latest commit

 

History

History
395 lines (298 loc) · 18.8 KB

numeric_integral.md

File metadata and controls

395 lines (298 loc) · 18.8 KB

Appendix B of 10.1109/TPWRD.2005.855448

A novel numerical integration technique is applied for the calculation of the infinite integrals of (1) and (2). It is based on the presence of the terms $e^{-(h_i+h_j)u}$ and $e^{-2 h_i u}$ in the integrand and on the observation that the integral of (2) has as a factor the term $cos(uy)$. Therefore, some of the roots of the integrand are known. Thus, when the horizontal distance $y$ between the two conductors is not zero, the integral in the right part of (2) can be calculated as follows.

The Gauss–Legendre method [14] is applied to calculate the integral between zero and the first root of $cos(uy)$, as this method is capable of handling the initial steep descent of the integrand. Then the Lobatto rule [10] is applied in the intervals between subsequent roots of $cos(uy)$. The procedure is repeated until convergence is achieved, controlled by a user defined tolerance, set to $10^{-9}$ in this paper.

For the calculation of the integral in (1) or the integral in (2) in the case where the horizontal distance $y$ is zero, the procedure is described below.

The Gauss–Legendre method is applied in order to find the integral at an initial interval between zero and $u_0 = 4 \times 10^{-4}/(h_i+h_j)$ in (2) or $u_0 = 4 \times 10^{-4}/(2h_i)$ in (1). Then the shifted Gauss-Laguerre method [14] is used for the evaluation of the rest of the integral. The procedure is repeated iteratively. In each iteration, the initial interval is bisected and the use of the Gauss–Legendre method is extended by $10 u_0$ intervals to the right of $u_0$. Convergence is achieved when the absolute difference between two succedent values of calculated integral is less than the predefined tolerance.

In all cases examined, convergence is achieved after 3–4 iterations.

A detailed analysis of the numerical algorithm is provided below.

1. Initial computation using Gauss-Legendre quadrature

Paper’s description:

"The Gauss–Legendre method is applied to calculate the integral between zero and the first root... as this method is capable of handling the initial steep descent of the integrand."

Code mapping:

This is handled by the first Gauss-Legendre integral computation in the function:

Gu_1 = compute_legendre_integral(u(1), height, s_legendre, w_legendre, permittivity_layers, omega);
Gu_1: Gauss-Legendre quadrature
function G = compute_legendre_integral(interval, height, points, weights, permittivity, omega)
    % Compute integral using Gauss-Legendre quadrature
    G = 0;
    for w = 1:length(weights)
        u_new = ((interval * (points(w) - 1)) / 2) + interval;
        F_u = nakagawa(u_new, omega, permittivity);
        G = G + weights(w) * exp(-2 * height * u_new) * F_u;
    end
    G = G * interval / 2;
end

This term evaluates the integral over the interval $[0, u_0]$ using Gauss-Legendre quadrature, which is suited for finite intervals. The general description of the computation is as follows:

  • Compute the initial segment of the integral from 0 to the first integration step size u(1), where the function typically has a steep descent (as the paper mentions).
  • compute_legendre_integral implements Gauss-Legendre quadrature, which uses roots of Legendre polynomials and weights for high accuracy in this segment.
  • Gauss-Legendre quadrature is well-suited for this, as it accurately handles functions with rapid changes in this region.
  • It divides $[0, u_0]$ into sub-intervals using the Legendre points and weights, which optimally sample the steep descent.

The first interval $u_0$ is estimated as:

u(1) = 4e-4 / (2 * height);

This scales the step size by the conductor height, ensuring it is small enough to capture the steep descent behavior, captured by the exponential term within the integration loop:

G = G + weights(w) * exp(-2 * height * u_new) * F_u;
  • The Gauss-Legendre method integrates over finite intervals by summing weighted function evaluations at specific quadrature nodes $u_{\text{new}}$ within the interval.
  • Each node $u_{\text{new}}$ contributes to the integral based on the value of the integrand at that point.
  • Since $u_{\text{new}}$ varies for each quadrature point, the exponential decay must be applied individually at each node.

Integral formulation

For a finite interval $[a, b]$, the integral is approximated as:

$$\int_a^b e^{-2 \cdot \text{height} \cdot u} f(u) \, du \approx \sum_{w=1}^N w_w \cdot e^{-2 \cdot \text{height} \cdot u_w} \cdot f(u_w)$$

Here, $w_w$ are the weights and $u_w$ are the quadrature nodes.

2. Evaluation using Laguerre quadrature

Paper’s description:

"Then the shifted Gauss-Laguerre method is used for the evaluation of the rest of the integral."

This term evaluates the semi-infinite part of the integral starting from $u_0$. Specifically, it uses Gauss-Laguerre quadrature, which is designed for intervals of the form $[u_0, \infty)$.

Code mapping:

This is handled by the Laguerre quadrature computation:

Gu_2 = compute_laguerre_integral(u(1), height, s_laguerre, w_laguerre, permittivity_layers);
Gu_2: Gauss-Laguerre Quadrature
function G = compute_laguerre_integral(u_start, height, points, weights, permittivity)
    % Compute integral using Laguerre quadrature
    G = 0;
    for v = 1:length(weights)
        u_new = points(v) / (2 * height) + u_start;
        F_u = nakagawa(u_new, 2 * pi * permittivity, permittivity);
        G = G + weights(v) * F_u;
    end
    G = G * exp(-2 * height * u_start) / (2 * height);
end

The general description of the computation is as follows:

  • Compute the contribution of the integral for the semi-infinite part of the domain using Gauss-Laguerre quadrature, i.e. the part of the integral extending to infinity, starting at $u_0$.
  • Laguerre quadrature naturally handles the exponential decay of terms like $e^{-2 h \cdot u}$, which dominate in the semi-infinite domain.
  • Gauss-Laguerre is well-suited for semi-infinite intervals. The weights $w_\text{laguerre}$ and nodes $s_\text{laguerre}$ are designed to approximate integrals weighted by exponential decay terms $e^{-x}$. It assumes the exponential decay $e^{-x}$ is already factored into the quadrature weights and nodes.
  • In this case, $e^{-2 \cdot \text{height} \cdot u_{\text{new}}}$ is already accounted for by the Laguerre quadrature weights (inherently designed for $e^{-x}$), so only $e^{-2 \cdot \text{height} \cdot u_{\text{start}}}$ needs to be applied after the summation.

Integral formulation

For the semi-infinite integral, the shifted domain becomes:

$$\int_{u_{\text{start}}}^\infty e^{-2 \cdot \text{height} \cdot u} f(u) \, du = e^{-2 \cdot \text{height} \cdot u_{\text{start}}} \cdot \int_0^\infty e^{-x} f(x + u_{\text{start}}) \, dx$$

Gauss-Laguerre handles $\int_0^\infty e^{-x} f(x + u_{\text{start}}) dx$ directly, so the exponential $e^{-2 \cdot \text{height} \cdot u_{\text{start}}}$ is applied once, outside the loop.

The compute_laguerre_integral function implements:

u_new = points(v) / (2 * height) + u_start;
  • $u_{\text{start}} = u(1)$ ensures a shift to match the end of the Gauss-Legendre interval.
  • $F_u$ is computed via the Nakagawa function, which is the main part of the integrand.
function F = nakagawa(u, omega, permittivity)
    % Compute the Nakagawa function
    mu_0 = 4 * pi * 1e-7;
    epsilon_0 = 8.854187817e-12;

    % Free-space propagation constant
    k_0 = omega * sqrt(mu_0 * epsilon_0);

    % Return Nakagawa function value
    F = 1 / (u + sqrt(u^2 + 1j * omega * mu_0 * (1 / permittivity)));
end

The nodes and weights of Laguerre quadrature are tailored for integrals involving exponentials like $e^{-x}$. The interval $[u_0, \infty)$ is effectively compressed into a manageable summation by Laguerre weights and nodes.

Hybrid quadrature method

The paper explains why two methods are used:

  1. Gauss-Legendre for $[0, u_0]$:
  • This interval contains a steep descent of the integrand, which Gauss-Legendre handles efficiently for short, finite intervals.
  • The sharp initial decay means that uniform quadrature (e.g., trapezoidal) would need many points to achieve similar accuracy.
  1. Gauss-Laguerre for $[u_0, \infty)$:
  • This region is dominated by the exponential decay $e^{-2 h \cdot u}$, making Laguerre quadrature the natural choice for efficient evaluation.
  • Using Gauss-Legendre here would require transforming the infinite domain or truncating it, which is less efficient.

It is noted that Gu_1 and Gu_2 do not integrate over the same interval. They are complementary contributions to the overall integral:

  • Gu_1: Covers $[0, u_0]$.
  • Gu_2: Covers $[u_0, \infty)$.

How Gauss-Laguerre handles semi-infinite domains

Gauss-Laguerre quadrature is specifically designed for semi-infinite integrals of the form:

$$\int_{0}^{\infty} e^{-x} f(x) \, dx$$

To compute an integral over $[u_0, \infty)$, the Laguerre transformation shifts the semi-infinite domain $[0, \infty)$ to start from $u_0$. This shift is performed in the code using the term:

u_new = points(v) / (2 * height) + u_start;

Step 1: Original semi-infinite interval The base Laguerre quadrature nodes (points) and weights (weights) are computed for the standard interval $[0, \infty)$, using the formula:

$$\int_{0}^{\infty} e^{-u} f(u) \, du \approx \sum_{v=1}^n w_v f(x_v)$$

Where:

  • $x_v$: Gauss-Laguerre nodes.
  • $w_v$: Gauss-Laguerre weights.

Step 2: Shifting the interval to start at $u_0$ To extend the integration domain to $[u_0, \infty)$, the variable $u$ is shifted by adding $u_0$:

$u_{\text{new}} = u_0 + \frac{x_v}{2 \cdot h}$

From the code:

u_new = points(v) / (2 * height) + u_start;
  • u_start: Acts as the starting point of the shifted interval $u_0$.
  • points(v) / (2 * height): Scales the Laguerre nodes appropriately for the integration range.

This transformation moves the integration domain from $[0, \infty)$ to $[u_0, \infty)$.

Step 3: How the exponential decay is handled

Laguerre quadrature implicitly accounts for the exponential decay $e^{-u}$. When $u_{\text{new}} = u_0 + x_v / (2 \cdot h)$, the decay factor $e^{-2 h \cdot u_{\text{new}}}$ is naturally incorporated into the computation. The scaling factor $1 / (2 \cdot h)$ ensures proper handling of the semi-infinite range.

From the code:

G = G + weights(v) * F_u; % Sum the weighted function values

In the computation of Gu_2, the Laguerre transformation shifts the effective interval to $[u_0, \infty)$ without explicitly stating so. Instead, this happens "automatically" due to:

  1. The use of shifted $u_{\text{new}} = u_0 + x_v / (2 \cdot h)$.
  2. The exponential decay factor $e^{-u}$, inherent in Laguerre quadrature.

Thus, even though Gu_2 is computed using Gauss-Laguerre points and weights, it effectively integrates from $u_0$ to infinity. This clever use of Gauss-Laguerre quadrature avoids having to directly handle the infinite upper bound $\infty$. The semi-infinite domain $[u_0, \infty)$ is "transformed" into a summation over Laguerre nodes and weights, with the scaling and shifting done via:

u_new = points(v) / (2 * height) + u_start;

In conclusion, Gu_2 does indeed compute the contribution from $[u_0, \infty)$, even though it’s not immediately obvious from the code.

3. Iterative extension of integration

Paper’s description:

"The procedure is repeated iteratively. In each iteration, the initial interval is bisected and the use of the Gauss–Legendre method is extended by intervals to the right of $u_0$."

Code mapping:

This is implemented in the iterative loop of the code:

for i = 2:iteration_limit
    % Update step size
    step_factor = 1 / (2^(i - 1));
    u(i) = 10 * u(i - 1);

    % Gauss-Legendre segments
    Gu_3 = compute_legendre_integral(u(1) * step_factor, height, s_legendre, w_legendre, permittivity_layers, omega);
    Gu_4 = compute_legendre_integral(u(i) - u(i - 1), height, s_legendre, w_legendre, permittivity_layers, omega);

    % Laguerre contribution
    Gu_6 = compute_laguerre_integral(u(i), height, s_laguerre, w_laguerre, permittivity_layers);

    % Total impedance
    dZ(i) = Gu_3 + Gu_4 + Gu_6;

    % Check convergence
    if abs(dZ(i) - dZ_previous) <= tolerance
        break;
    end

    % Update for next iteration
    dZ_previous = dZ(i);
end

The general description of the computation is as follows:

  • Extends the integration by progressively adding contributions from new intervals, refining the solution iteratively.
  • The step factor halves the segment size at each iteration, improving resolution.
  • u(i) defines the new integration bound, increasing logarithmically to cover the semi-infinite domain.

The main terms in the algorithm are:

  • Gu_3: Computes a small segment at the leftmost interval $[0, u_0/2]$ using Gauss-Legendre.

  • Gu_4: Computes the interval between the two successive bounds $[u(i-1), u(i)]$ using Gauss-Legendre.

  • Gu_6: Adds the Laguerre contribution for the new interval $[u(i), \infty)$.

  • Convergence check: The code compares successive impedance estimates $dZ(i)$ to $dZ(i-1)$ using:

if abs(dZ(i) - dZ_previous) <= tolerance

The loop stops when the change falls below the predefined tolerance.

4. Final computation of impedance

Paper’s description:

"Convergence is achieved when the absolute difference between two successive values of the calculated integral is less than the predefined tolerance."

Code mapping:

This is reflected in the final line of the main function:

d_Z = 1j * omega * mu_0 * dZ(i) / pi;

Once convergence is achieved, the final self-impedance $\Delta Z$ is computed as: $\Delta Z = \frac{j \cdot \omega \cdot \mu_0 \cdot dZ(i)}{\pi}$ Here, dZ(i) represents the converged value of the integral.

Summary: how each part matches the paper

Step in the paper Code implementation
Initial steep descent handled by Gauss-Legendre Gu_1 = compute_legendre_integral(u(1), ...)
Shifted Gauss-Laguerre for semi-infinite domain Gu_2 = compute_laguerre_integral(u(1), ...)
Iterative refinement via step bisection Iterative for loop, halving the step size: step_factor = 1 / (2^(i-1)), adding contributions with Gu_3, Gu_4, and Gu_6.
Convergence check using tolerance if abs(dZ(i) - dZ_previous) <= tolerance, breaking the loop when the result stabilizes.
Final impedance computed after convergence d_Z = 1j * omega * mu_0 * dZ(i) / pi, assembling the result using the final converged value of the total integral.

Complete code implementation

function d_Z = method_self_aeras_1_strwma(frequency, permittivity_layers, height)
    % Numerically evaluate self-impedance integrals using hybrid quadrature
    % frequency: frequency of operation [Hz]
    % permittivity_layers: permittivity values of the layers
    % height: distance from conductor to ground [m]

    % --- Constants ---
    tolerance = 1e-10;    % Convergence tolerance
    mu_0 = 4 * pi * 1e-7; % Magnetic permeability
    omega = 2 * pi * frequency;

    % Load Gauss-Legendre and Laguerre points/weights
    [s_legendre, w_legendre] = load_quadrature('ah16.txt');
    [s_laguerre, w_laguerre] = load_quadrature('ah35.txt');

    % --- Initialization ---
    u(1) = 4e-4 / (2 * height);
    dZ_previous = 0;
    iteration_limit = 100; % Fail-safe iteration cap

    % Compute initial segments
    Gu_1 = compute_legendre_integral(u(1), height, s_legendre, w_legendre, permittivity_layers, omega);
    Gu_2 = compute_laguerre_integral(u(1), height, s_laguerre, w_laguerre, permittivity_layers);
    dZ(1) = Gu_1 + Gu_2;

    % Iterative computation
    for i = 2:iteration_limit
        % Update step size
        step_factor = 1 / (2^(i - 1));
        u(i) = 10 * u(i - 1);

        % Gauss-Legendre segments
        Gu_3 = compute_legendre_integral(u(1) * step_factor, height, s_legendre, w_legendre, permittivity_layers, omega);
        Gu_4 = compute_legendre_integral(u(i) - u(i - 1), height, s_legendre, w_legendre, permittivity_layers, omega);

        % Laguerre contribution
        Gu_6 = compute_laguerre_integral(u(i), height, s_laguerre, w_laguerre, permittivity_layers);

        % Total impedance
        dZ(i) = Gu_3 + Gu_4 + Gu_6;

        % Check convergence
        if abs(dZ(i) - dZ_previous) <= tolerance
            break;
        end

        % Update for next iteration
        dZ_previous = dZ(i);
    end

    % Final result
    d_Z = 1j * omega * mu_0 * dZ(i) / pi;
end

% --- Helper Functions ---
function [points, weights] = load_quadrature(filename)
    % Load quadrature points and weights from file
    data = load(filename);
    points = data(:, 1);
    weights = data(:, 2);
end

function G = compute_legendre_integral(interval, height, points, weights, permittivity, omega)
    % Compute integral using Gauss-Legendre quadrature
    G = 0;
    for w = 1:length(weights)
        u_new = ((interval * (points(w) - 1)) / 2) + interval;
        F_u = nakagawa(u_new, omega, permittivity);
        G = G + weights(w) * exp(-2 * height * u_new) * F_u;
    end
    G = G * interval / 2;
end

function G = compute_laguerre_integral(u_start, height, points, weights, permittivity)
    % Compute integral using Laguerre quadrature
    G = 0;
    for v = 1:length(weights)
        u_new = points(v) / (2 * height) + u_start;
        F_u = nakagawa(u_new, 2 * pi * permittivity, permittivity);
        G = G + weights(v) * F_u;
    end
    G = G * exp(-2 * height * u_start) / (2 * height);
end

function F = nakagawa(u, omega, permittivity)
    % Compute the Nakagawa function
    mu_0 = 4 * pi * 1e-7;
    epsilon_0 = 8.854187817e-12;

    % Free-space propagation constant
    k_0 = omega * sqrt(mu_0 * epsilon_0);

    % Return Nakagawa function value
    F = 1 / (u + sqrt(u^2 + 1j * omega * mu_0 * (1 / permittivity)));
end