根求解数量增加时精度丢失的原因排查

问题背景 #

使用二分法（while循环实现）求解方程根时，求解100个根时精度可达250多位小数，但求解1000个根时，最后一个根的精度仅剩6位小数。相关代码及结果如下：

代码实现 #

setprecision(BigFloat, 250; base=10)

Mp = BigFloat("250.0")
L = BigFloat("1.0")
rp = BigFloat("1.0")


exp = Mp - 10
pk = 50
h = 1/10
ts = -2*pk
te = -h
v0 = -pk
vmax = 0
u0 = 0
umax = pk
k1 = (vmax -v0)/h
k2 = (umax -u0)/h
 
max = 100.0
function arctanh(x)
    if abs2(x) < 1
        res = BigFloat(log( ( 1 + x ) / (1 - x)  ))
    else
        res = BigFloat(sign(x) * log1p( ( 2 * abs(x) ) / (1 - abs2(x)) ))
    end
    return 0.5*res
end

function tor_2(r)
    a = BigFloat( rp / r )
    b = BigFloat(arctanh(a))
    res = BigFloat(L^2 * (-b / rp))
    return res
end

function find_root(f, a, b, precision)
    oldm = big(a)
    n = big(a)
    p = big(b)
    m = (n + p) / big(2)
    epsilon = big(10)^(-precision)

    while abs(f(m)) > epsilon && oldm != m
        oldm = m
        val = f(m)
        if val > 0
            p = m
        elseif val < 0
            n = m
        else
            break
        end
        m = (n + p) / BigFloat(2.0)
    end
    return m
 
end
t = @elapsed begin
    tort = [find_root(r -> tor_2(r) - test, rp + 10.0^-exp, max, Mp) for test in ts:h:te]
end

结果对比 #

• 正确结果（tort[1]）：

1.00000000000000000000000000000000000000000000000000000000000000000000000000000000000000276779305347347506129736291395816937080609516467895441878785070622487206090198597561760179491324750990568982322199122021923608490714617064493293276254845275543133336

• 错误结果（tort[1000]）：

10.033311132253989056880435010549127161920612410233012969744805223401820078850954716601912612057050747701277915837004398085186404284469580227204776685726634389654992756053090443362119682617587346589647165675974434032301224533470334376859854777986199219

• 预期结果：

10.03331113225398961014527049285149912673425410710831002117675758378789785273643535498548789545079535319530347609313690968614200882388184338974490564270302639426889010992841078245555898544027132450710912524700458841030045198039559947676313086937548943

精度丢失原因 #

• tor_2(r)的数值稳定性随r增大下降：当r很大时，a = rp/r趋近于0，此时用log((1+x)/(1-x))计算arctanh(x)会引入较大相对误差——因为(1+x)/(1-x)接近1，log函数在1附近的输入微小变化会被放大，导致b的精度损失，进而影响tor_2(r)的计算精度，让二分法的收敛判断出现偏差。
• 二分法收敛条件不严谨：当前仅用abs(f(m)) > epsilon判断是否收敛，但当f(r)在根附近斜率极小时（tor_2(r)的导数L²/(rp² - r²)绝对值趋近于0），即使f(m)的绝对值小于阈值，m与真实根的偏差也可能很大，无法保证精度。
• 测试序列的累积误差：生成test序列时，ts:h:te使用默认浮点数而非BigFloat，多次累加h会产生累积误差，最后几个test值本身精度不足，直接影响根的求解结果。
• 迭代提前终止：循环终止条件包含oldm != m，当BigFloat的精度用尽时，m不再变化，即使abs(f(m))仍未达标，循环也会停止。对于大r场景，需要更多迭代次数才能收敛到指定精度，却可能提前终止。

解决建议 #

• 优化arctanh计算：当x趋近于0时，改用泰勒展开式arctanh(x) = x + x^3/3 + x^5/5 + ...计算，直到新增项小于指定精度，避免log函数的误差。
• 改进收敛条件：同时监测m的变化量，比如添加abs(m - oldm) < epsilon/10作为终止条件，确保m的精度满足要求。
• 用BigFloat生成测试序列：将ts、h、te全部转换为BigFloat，比如ts = BigFloat(-2*pk)，h = BigFloat(1)/10，避免累加时的精度损失。
• 调整精度与迭代次数：针对大r场景，适当提高BigFloat的精度上限，或者增加迭代次数，确保收敛到足够精确的根。

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