BESSELK

What it does

BESSELK(x, n)

BESSELK(x, n) returns K_n(x), the modified Bessel function of the second kind, of integer order n at the real point x > 0. It is the second solution of the modified Bessel equation

x² w″ + x w′ − (x² + n²) w  =  0

the one that decays as x grows, where BESSELI is the one that stays finite at the origin. Microsoft's page defines it through the J and Y Bessel functions at an imaginary argument, in an equation given as an image. The standard form of that relation is

K_n(x)  =  (π/2) · i^(n+1) · [ J_n(i·x) + i · Y_n(i·x) ]

For real x > 0 a more useful form is the integral

K_n(x)  =  ∫₀^∞  e^(−x cosh t) cosh(n t) dt

which makes three things visible at once: K_n(x) is positive, it falls as x grows, and it rises as the order rises.

• Decay. For large x, K_n(x) ≈ √(π/(2x)) · e^(−x). As with I_n, every order looks alike at that scale. K_0(100) is already about 4.7E-45.
• The singularity at 0. K_0(x) ≈ −ln x as x → 0, and K_n(x) ≈ ½ (n−1)! (2/x)^n for n ≥ 1. The function goes to infinity at the origin, faster at higher order, and has no value there.
• Negative x. K_n has a branch point at 0 and is not real on the negative axis. There is no parity rule of the kind I_n has, which is why this function has a smaller domain than BESSELI.

Arguments

x must be positive, and the documentation does not say so. Microsoft's remarks cover only a nonnumeric argument and n < 0. The recorded answers show #NUM! on all 24 calls with x ≤ 0 and a non-negative order, including x = 0 at orders 0, 1, 2 and 10 (Excel 16.0 build 20430, Windows 64-bit, CPU not recorded).

Truncation and the negative test. As with BESSELI, the sign test sees n before truncation. All 51 recorded calls with x > 0 and n strictly between −1 and 0 returned #NUM!, although truncation would give order 0. Positive fractions truncate as documented, and BESSELK(x, 0.9) reads the same as BESSELK(x, 0) at the grid points where both were asked.

BESSELK is not a function of real order, even though K_ν is perfectly well defined for real ν. BESSELK(x, 1.5) is K_1(x), not the closed-form K_(3/2)(x).

The slots are numeric. No reading of text, a logical, an empty cell or an array in either slot is on record for this function; see What has not been checked.

Try it

Quirks

Result and edge cases

Returns Number. Unless stated otherwise, the behaviour below is read from the 1,236 recorded answers on build 20430.

• The documented example BESSELK(1.5, 1) is printed by Microsoft as 0.277387804. The true value is 0.2773878005, so the printed value differs from it in the ninth significant digit. This is the function's normal accuracy, not an error in the documentation.
• Small x, ordinary order behaves correctly and the numbers grow quickly. The smallest positive x read is 0.001, where order 10 is about 1.9E38, the largest value in the recording.
• Large x decays toward the bottom of the range. The largest x read is 100, where order 0 is about 4.7E-45. In the recording no call with x > 0 and a non-negative order raised an error, and neither end of the range was reached.
• Overflow and underflow. K_n is not representable for high order at small x, or at x beyond about 740, where K_0 rounds to 0 below the smallest subnormal. How Excel reports either has not been read. BESSELI raises #NUM! well before its own limit, so a silent 0 cannot be assumed here.

Errors

As documented on Microsoft's page:

As recorded, in addition to the documented #NUM! for n ≤ −1 (build 20430):

The documentation states neither row. The first is forced by the mathematics and is what any reader would expect; it is still a gap in the page. The second is a choice, and an implementation that truncates first gets it wrong.

How to implement

Numerical notes

Low-order polynomials, by design. The notes below record that BESSELK is correctly rounded on 0 of 1,155 recorded rows, with a largest relative error of 1.7E-7. Treat it as a seven-digit function: the same order of accuracy as BESSELI, and for the same reason.

The form the recorded answers follow. The notes do not give a staging for this function. The published candidates are the polynomials of Abramowitz & Stegun 9.8.5 to 9.8.8: for x ≤ 2 a polynomial in (x/2)² plus a ln(x/2) term multiplied by the polynomial I_0 or I_1, and for x > 2 a polynomial in 2/x multiplied by e^(−x)/√x. Orders 2 and above then step upward from K_0 and K_1. While writing this page that recipe was transcribed in plain binary64 and compared with the recorded rows: it gave the same double on 1,141 of the 1,155 numeric rows and came within 10 units on the other 14 (Excel 16.0 build 20430, Windows 64-bit, CPU not recorded). That is strong evidence for the form, and it is a spot check on the same rows, not a held-out measurement. It says nothing about the units in the remaining 14, which is where a staging difference such as the x87 storage the notes describe for BESSELI would show.

Upward is the stable direction. K_n grows with n, so stepping up amplifies nothing dangerous. There is no starting point to choose and no analogue of BESSELI's large-x defect. The relative error at orders 2 and above in the recording comes from K_0 and K_1 and thins out as the order rises: the worst row at each order falls from about 1.7E-7 at order 1 to between 1E-8 and 2.4E-8 at orders 12 and above.

Where the largest errors sit. The worst rows in the recording are at order 1, and then order 0, with x near 38 to 40. That is a swing of the large-x polynomial's own error curve, with nothing to suggest a defect.

There is a seam at x = 2. Evaluating the printed coefficients on both sides of 2 gives a jump of about 6.7E-8 relative in K_0 and 1.1E-8 in K_1. That is arithmetic on the published coefficients: no recorded row sits within 0.0003 of 2. The matching step in BESSELI is visible in Excel's answers at adjacent doubles, so there is every reason to expect this one too.

Rust code

62 besselk "BESSELK" (2, 2) [V "x", V "n"] |s, a, sv| { bessel(s, a, k::bessel::besselk_kernel, false) }

fn ln<const S: bool>(x: D<S>) -> D<S> {
    D::new(excel_log(x.0))
}

fn i0<const S: bool>(x: D<S>) -> D<S> {
    let ax = x.abs();
    if ax.0 < 3.75 {
        let t = x / D(3.75);
        horner(&I0_SMALL, t * t)
    } else {
        exp(ax) / ax.sqrt() * horner(&I0_LARGE, D(3.75) / ax)
    }
}

fn k0<const S: bool>(x: D<S>) -> D<S> {
    if x.0 <= 2.0 {
        let h = x / D(2.0);
        -(ln(h) * i0(x)) + horner(&K0_SMALL, x * x / D(4.0))
    } else {
        exp(-x) / x.sqrt() * horner(&K0_LARGE, D(2.0) / x)
    }
}

fn k1<const S: bool>(x: D<S>) -> D<S> {
    if x.0 <= 2.0 {
        let h = x / D(2.0);
        ln(h) * i1(x) + D(1.0) / x * horner(&K1_SMALL, x * x / D(4.0))
    } else {
        exp(-x) / x.sqrt() * horner(&K1_LARGE, D(2.0) / x)
    }
}

/// `BESSELK(x, n)`.
pub fn besselk_kernel(x: f64, n: f64) -> Result<f64, WorksheetErrorCode> {
    if !x.is_finite() { return Err(WorksheetErrorCode::Num); }
    let n = order(n)?;
    let x: D<true> = D::new(x);
    if !(x.0 > 0.0) {
        return Err(WorksheetErrorCode::Num);
    }
    let value = match n {
        0 => k0(x),
        1 => k1(x),
        _ => forward(x, n, k0(x), k1(x), true),
    };
    finish(Ok(D(value.0)))
}

No implementation has been published with a Handbook test suite for this function.

Version differences

The recorded environments below distinguish a reading from a same-input comparison. Counts describe captured inputs in each environment; overlapping captures and comparisons are not totals of unique worksheet calls.

The record contains no same-input comparison between these environments. Reading different inputs on different builds does not establish either a version change or cross-build agreement.

Recorded environments and comparison details.

Alternatives

Relationships

• BESSELI — the growing solution of the same equation. The two are tied by the exact cross-product I_n(x)·K_(n+1)(x) + I_(n+1)(x)·K_n(x) = 1/x, a useful check that needs no oracle. In the form described under Numerical notes, K_0 and K_1 use I_0 and I_1 for x ≤ 2, so the two functions are not independent there. Their accuracy note is shared.
• The recurrence. K_(n+1)(x) = K_(n−1)(x) + (2n/x)·K_n(x). Every term is positive, which is why stepping upward in n is safe here and unsafe for I_n.
• Where readers meet it. K_0 is the kernel of the two-dimensional screened Poisson equation and the steady-state well function of a leaky aquifer, and K_1 appears in the density of the normal-inverse-Gaussian distribution.

Named variants

No named variants are recorded yet.