一般来说,牛顿法的收敛性依赖于初值 x 0 x_{0} x0 的选取,如果 x 0 x_{0} x0 偏离方程的正解根 x ∗ x^{*} x∗ 较远,则牛顿法可能发散。
例:用牛顿法求方程 x 3 − x − 1 = 0 x^{3}-x-1=0 x3−x−1=0 在 x = 1.5 x = 1.5 x=1.5 附近的一个根。
解:取迭代初值 x 0 = 1.5 x₀=1.5 x0=1.5,用牛顿公式
x k + 1 = x k − x k 3 − x k − 1 3 x k 2 − 1 x_{k+1} = x_{k} - \frac{x^{3}_{k}-x_{k}-1}{3x^{2}_{k}-1} xk+1=xk−3xk2−1xk3−xk−1
计算结果如下:
x
1
=
1.34783
,
x
2
=
1.32520
,
x
3
=
1.32472
x₁ = 1.34783, x₂ = 1.32520, x₃ = 1.32472
x1=1.34783,x2=1.32520,x3=1.32472
其中
x
3
x_{3}
x3 的每一位数字都是有效数字。
但是,如果改用 x 0 = 0.6 x₀ = 0.6 x0=0.6 作为初值,则按照牛顿公式迭代一次得 x 1 = 17.9 x₁ = 17.9 x1=17.9,这个结果反而比 x 0 x₀ x0 更偏离了所求的根 x ∗ x^{*} x∗。
为了防止迭代发散,通常对迭代过程再附加一项要求,即保证函数值单调下降,满足这项要求的算法称为下山法。
将牛顿法与下山法结合起来使用,即在下山法保证函数值稳定下降的前提下,用牛顿法加快收敛速度。
为此,将牛顿法的上一次迭代初值与本次迭代值进行适当加权平均后作为新的改进值,我们有理由相信,加权后的改进值是比未改进前的近似根更好的的根,它更好地逼近了真实解。即采用新的迭代公式
x k + 1 = x k − λ f ( x k ) f ′ ( x k ) x_{k+1}=x_{k}-\lambda\frac{f(x_{k})}{f'(x_{k})} xk+1=xk−λf′(xk)f(xk)
其中 λ 称为下山因子,适当选择下山因子可以使函数满足单调下降条件。下山因子的选取是一个逐步探索的过程,从 λ = 1 λ = 1 λ=1 开始反复将因子 λ 的值减半进行试算,一旦单调性条件成立,则 “下山成功”,否则 “下山失败”。若 “下山失败”,这时则需另选初值 x 0 x₀ x0 重算。
如前有例子描述,当选取初值 x 0 = 0.6 x₀ = 0.6 x0=0.6 时,如果取下山因子 λ = 1 32 λ=\frac{1}{32} λ=321,则可求出 x 1 = 1.140625 x₁ = 1.140625 x1=1.140625,这个结果纠正了原有的严重偏差。
运行示例:
程序源码:
#include <iostream>
#include <cmath>
using namespace std;
/**
* 原函数 f(x)
*/
double f(double x)
{
return pow(x, 3) - x - 1;
}
/**
* 导函数 f1(x) = f'(x)
*/
double f1(double x)
{
return 3 * pow(x, 2) - 1;
}
int main(void)
{
double x0;
cout << "请输入迭代初值:";
cin >> x0;
double accuracy;
cout << "请输入精度:";
cin >> accuracy;
int N;
cout << "请输入最大迭代次数:";
cin >> N;
double x1;
int count = 0;
do
{
if (f1(x0) == 0)
{
cout << "在 x0 附近 f(x) 的一阶导数值为 0,不适用牛顿法求方程的根!" << endl;
break;
}
double L = 1;
x1 = x0 - L * f(x0) / f1(x0);
// 函数值未单调下降,控制下山因子的大小使得单调性条件得到满足,若不满足,则换初值重算
while (abs(f(x1)) - abs(f(x0)) >= 0)
{
L /= 2;
// 应用下山因子对迭代值进行改进
x1 = x0 - L * f(x0) / f1(x0);
}
cout << "\n下山因子 L = " << L << "第 " << count << " 次迭代,方程的近似根为:" << x1 << endl;
double temp = x1;
x1 = x0;
x0 = temp;
if (++count > N)
{
cout << "达到允许的最大迭代次数!迭代结束!" << endl;
break;
}
} while ((abs(x1 - x0) > accuracy));
return 0;
}