如何优化Float64精度下的表达式计算以避免负值问题？

我正在使用Julia语言，相信解决方案大多可跨语言迁移。我需要在Float64精度下对正数$T$计算如下表达式：$\sqrt{4e^{-\alpha T} - 3 - e^{-2\alpha T} + 2\alpha T}$。其中α为正常量，根号内的表达式理论上恒为正，但当$T$趋近于0时，因数值精度问题偶尔会得到负值。

使用Float64计算根号内表达式的结果如下：

julia> α=1/4.0; map( T -> ( 4 * exp(-α * T) - 3 - exp( -2α * T ) + 2α * T ), 10.0.^(-10:0) )
11-element Vector{Float64}:
 -4.137018548132909e-18
 -4.137018546840439e-17
  3.038735495922355e-17
 -2.512379594116203e-16
  4.113329634720811e-17
  1.8928899382176026e-16
  1.0552627836227929e-14
  1.041483522687403e-11
  1.0397158019953556e-8
  1.022361261644733e-5
  0.00867247257298609

而使用BigFloat计算相同值时则无此问题：

julia> α=BigFloat(1/4.0); map( T -> ( 4 * exp(-α * T) - 3 - exp( -2α * T ) + 2α * T ), 10.0.^(-10:0) )
11-element Vector{BigFloat}:
 1.041666666647135530517511139334179132983501942809441564565600382453160881163658e-32
 1.041666666471354361319426381902535606415006422809427598976663299435084708258933e-29
 1.041666664713541734328314928659467024540737963331449636847040682408936614887116e-26
 1.041666647135417166709396893449432138811473447225776037995767356684771923939405e-23
 1.041666471354189311710850303868877685043857627660615702953528423064833434922657e-20
 1.04166471354394556609940068924547895770393434711169898595204995967864381931136e-17
 1.041647135644529365261489400892467181239989440055364005317784883808499181692089e-14
 1.041471376951090709983860760560291162481306334994613272150860749219930444115473e-11
 1.039715818279496102769801850953898783936741961589734116598630369857383947507474e-08
 1.022361261666541821235659358423414459600515733306470783618814653113646037198644e-05
 0.008672472572986049376881532922102135745171026217378941282377306337925928800328869

请问是否可修改计算方式，在仍使用Float64的前提下提升精度并避免该负值问题？

解决方案 #

问题根源 #

当$T$趋近于0时，$x = \alpha T$趋近于0，直接计算根号内的表达式$4e^{-x} -3 -e^{-2x}+2x$会发生相近数相减抵消：多个高阶小量叠加后，低阶项相互抵消，剩余的四阶小量在Float64精度下容易因计算误差出现负值，导致根号无法计算。

通过泰勒展开可以看到，根号内的表达式本质是四阶小量：
$$f(x) = 4e^{-x} -3 -e^{-2x}+2x = \frac{x^4}{12} + \frac{x^6}{720} + O(x^8)$$

改进方案 #

采用分段计算策略：当$x$足够小时，用泰勒展开式近似计算；当$x$较大时，使用原表达式计算，同时加入非负保护确保根号输入合法。

Julia实现代码

function compute_expr(α::Float64, T::Float64)
    x = α * T
    # 阈值可根据实际需求调整，这里取1e-3平衡精度与效率
    if x < 1e-3
        # 泰勒展开到x^6项，足够覆盖Float64精度需求
        f = x^4 / 12 + x^6 / 720
    else
        f = 4 * exp(-x) - 3 - exp(-2x) + 2x
    end
    # 数值保险：确保根号输入非负
    f = max(f, 0.0)
    return sqrt(f)
end

测试对比

运行以下代码验证改进效果：

α = 1/4.0
Ts = 10.0.^(-10:0)

# 原方法（加非负保护）
original = map(T -> sqrt(max(4*exp(-α*T)-3-exp(-2α*T)+2α*T, 0.0)), Ts)
# 改进方法
improved = map(T -> compute_expr(α, T), Ts)
# BigFloat参考结果
α_big = BigFloat(1/4.0)
big_ref = map(T -> Float64(sqrt(BigFloat(4*exp(-α_big*T)-3-exp(-2α_big*T)+2α_big*T))), Ts)

# 打印对比
for (t, o, i, b) in zip(Ts, original, improved, big_ref)
    println("T = $t")
    println("原方法结果：$o")
    println("改进方法结果：$i")
    println("BigFloat参考：$b")
    println("---")
end

效果说明 #

• 当$T$很小时（如$10^{-10}$），改进方法通过泰勒展开直接计算四阶小量，避免了抵消误差，结果与BigFloat几乎一致，不会出现负值。
• 当$T$较大时，原表达式计算精度足够，无需额外处理。
• 该方案完全基于Float64计算，保持了原有的计算效率，且逻辑简单易迁移到其他语言。

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