为什么算法在增加数字精度时会失败?我们如何降低算法对数字精度的敏感度?
Why the algorithm fails when increase the number precision? How can we decrease the sensitivity of the algorithm to the number precision?
我正在使用Newton Raphson + successive Substitute algorithm来进行flash calculation(化学过程模拟)。
当输入为低精度如0.1时,该算法可以很好地收敛,但当数字精度增加到0.11111或0.99999时。该算法不会收敛。
当我使用带有BFGS更新的准牛顿法时,同样的问题再次出现。我们如何降低代码对数值精度的敏感性?
下面是一个使用 matlab 求解 Rachford-Rice 方程的简单示例。当comp_overall=[0.9,1-0.9]时收敛的很好。但是,当数字精度增加到 [0.99999,1-0.99999] 时。它不会收敛。
K=[0.053154011443159 34.234731216532658],
comp_overall= [0.99999 1- 0.99999], phi=0.5; %initial values
epsilon = 1.0;
iter1 = 1;
while (epsilon >=1.e-05)
rc=0.0;
drc=0.0;
for i=1:2
% Rachford-Rice Equation
rc = comp_overall(i)*(K(i)-1.0)/(1.0+phi*(K(i)-1.0))+rc;
% Derivative
drc = comp_overall(i)*(K(i)-1.0)^2/(1.0+phiK(i)-1.0))^2+drc;
end
% Newton-Raphson
phi1 = phi +0.01 (rc / drc);
epsilon = abs( (phi1-phi)/phi );
% Convergence
phi = phi1;
iter1=iter1+1;
end
Newton-Raphson 方法依赖于函数在任意两个连续近似值之间是可微的。根据初始值的选择,z₁ = 0.99999
可能不是这种情况。让我们看一下 Rachford-Rice 函数的图形:
这个函数的根是φ₀ ≈ –0.0300781429
,最近的不连续点是–1/(K₂-1) ≈ –0.0300890052
。它们足够接近,足以让 Newton–Raphson 方法超调,跳过那个不连续点。
例如:
φ₁ = –0.025
f(φ₁) ≈ -0.9229770571
f'(φ₁) ≈ 1.2416569960
φ₂ = φ₁ + 0.01 * f(φ₁) / f'(φ₁) ≈ -0.0324334302
φ₂
位于不连续点的左侧,因此接下来的步骤将远离而不是朝向根。
φ₃ = -0.0358986759 < φ₂
可以做些什么:
- 当算法无法收敛时,以较小的步长重复。例如,从系数 0.01(现在)开始,每次失败后将其减小 10 次。
- 检测过冲。在每次迭代中检查当前近似值和前一个近似值之间是否存在不连续点 (
–1/(Kᵢ-1)
)。当它发生时,丢弃当前的近似值,减小系数并继续。
- 限制搜索范围。 [0, 1] 之外的解在物理上有意义吗?如果不是,您可以在近似值超出该范围后停止。
- 使用不同的方法。该函数在两个连续不连续点之间的任何区间上都是单调的,因此您可以对每个这样的区间执行二分查找。它将比 Newton–Raphson 方法更快、更稳健。
我正在使用Newton Raphson + successive Substitute algorithm来进行flash calculation(化学过程模拟)。
当输入为低精度如0.1时,该算法可以很好地收敛,但当数字精度增加到0.11111或0.99999时。该算法不会收敛。
当我使用带有BFGS更新的准牛顿法时,同样的问题再次出现。我们如何降低代码对数值精度的敏感性?
下面是一个使用 matlab 求解 Rachford-Rice 方程的简单示例。当comp_overall=[0.9,1-0.9]时收敛的很好。但是,当数字精度增加到 [0.99999,1-0.99999] 时。它不会收敛。
K=[0.053154011443159 34.234731216532658],
comp_overall= [0.99999 1- 0.99999], phi=0.5; %initial values
epsilon = 1.0;
iter1 = 1;
while (epsilon >=1.e-05)
rc=0.0;
drc=0.0;
for i=1:2
% Rachford-Rice Equation
rc = comp_overall(i)*(K(i)-1.0)/(1.0+phi*(K(i)-1.0))+rc;
% Derivative
drc = comp_overall(i)*(K(i)-1.0)^2/(1.0+phiK(i)-1.0))^2+drc;
end
% Newton-Raphson
phi1 = phi +0.01 (rc / drc);
epsilon = abs( (phi1-phi)/phi );
% Convergence
phi = phi1;
iter1=iter1+1;
end
Newton-Raphson 方法依赖于函数在任意两个连续近似值之间是可微的。根据初始值的选择,z₁ = 0.99999
可能不是这种情况。让我们看一下 Rachford-Rice 函数的图形:
这个函数的根是φ₀ ≈ –0.0300781429
,最近的不连续点是–1/(K₂-1) ≈ –0.0300890052
。它们足够接近,足以让 Newton–Raphson 方法超调,跳过那个不连续点。
例如:
φ₁ = –0.025
f(φ₁) ≈ -0.9229770571
f'(φ₁) ≈ 1.2416569960
φ₂ = φ₁ + 0.01 * f(φ₁) / f'(φ₁) ≈ -0.0324334302
φ₂
位于不连续点的左侧,因此接下来的步骤将远离而不是朝向根。φ₃ = -0.0358986759 < φ₂
可以做些什么:
- 当算法无法收敛时,以较小的步长重复。例如,从系数 0.01(现在)开始,每次失败后将其减小 10 次。
- 检测过冲。在每次迭代中检查当前近似值和前一个近似值之间是否存在不连续点 (
–1/(Kᵢ-1)
)。当它发生时,丢弃当前的近似值,减小系数并继续。 - 限制搜索范围。 [0, 1] 之外的解在物理上有意义吗?如果不是,您可以在近似值超出该范围后停止。
- 使用不同的方法。该函数在两个连续不连续点之间的任何区间上都是单调的,因此您可以对每个这样的区间执行二分查找。它将比 Newton–Raphson 方法更快、更稳健。