如何在Fortran或MATLAB中处理BesselI(0,a)/BesselI(0,b)以避免NaN？

问题根源 #

直接计算大参数的修正第一类贝塞尔函数I₀(x)时，由于其指数增长特性，当x > ~709（双精度浮点数阈值）会溢出到inf，两个inf相除就得到NaN。而级数展开在大参数下收敛极慢甚至中途溢出，完全不适用。下面给出几种稳定计算比值的方法：

方法1：对数差分法 #

核心思路是计算ln(I₀(a)) - ln(I₀(b))，再取指数得到比值，避免直接计算超大的I₀(x)值。结合精确计算+渐近展开覆盖全参数范围：

MATLAB 实现

function logI0 = log_besseli0(x)
    if x < 20
        logI0 = log(besseli(0, x));
    else
        % 大参数渐近近似：ln(I₀(x)) ≈ x - 0.5*ln(2πx)
        logI0 = x - 0.5*log(2*pi*x);
    end
end

% 调用示例
a = 1200; b = 1150;
ratio = exp(log_besseli0(a) - log_besseli0(b));
disp(ratio);

Fortran 实现

function log_besseli0(x) result(logI0)
    implicit none
    double precision, intent(in) :: x
    double precision :: logI0
    double precision, parameter :: pi = 4.0d0 * atan(1.0d0)
    
    if (x < 20.0d0) then
        logI0 = log(besseli(0.0d0, x))
    else
        logI0 = x - 0.5d0 * log(2.0d0 * pi * x)
    end if
end function log_besseli0

program compute_ratio
    implicit none
    double precision :: a, b, ratio, log_besseli0
    a = 1200.0d0
    b = 1150.0d0
    ratio = exp(log_besseli0(a) - log_besseli0(b))
    print *, "Ratio: ", ratio
end program compute_ratio

方法2：渐近比值近似 #

当a和b都远大于1时，直接用渐近式计算比值，精度足够且效率极高：
I₀(a)/I₀(b) ≈ exp(a - b) * sqrt(b/a)

MATLAB 实现

function ratio = besseli0_ratio(a, b)
    if a < 20 && b < 20
        ratio = besseli(0,a)/besseli(0,b);
    else
        ratio = exp(a - b) * sqrt(b/a);
    end
end

Fortran 实现

function besseli0_ratio(a, b) result(ratio)
    implicit none
    double precision, intent(in) :: a, b
    double precision :: ratio
    
    if (a < 20.0d0 .and. b < 20.0d0) then
        ratio = besseli(0.0d0, a) / besseli(0.0d0, b)
    else
        ratio = exp(a - b) * sqrt(b / a)
    end if
end function besseli0_ratio

为什么级数展开行不通？ #

I₀(x)的级数形式为sum_{k=0}^∞ (x²/4)^k / (k!)²，当x很大时，前若干项会急剧增大到溢出阈值，后续项才开始衰减，不仅收敛速度极慢（需要数千项），还会在计算过程中直接溢出，完全不适合大参数场景。

内容的提问来源于stack exchange，提问作者user2665861