[计算机数值分析]埃特金算法加速迭代法求根过程

Spring-_-Bear 的 CSDN 博客导航

---

由于 前有加速方案 需要提供迭代函数 $\phi (x)$ 的导数 ${\phi }^{\prime }(x)$ 而不便于实际应用。

若将迭代值 ${\overline{x}}_{k+1}$ = $\phi (x_k)$ 再迭代一次，得

$${\tilde{x}}_{k+1}=\phi ({\overline{x}}_{k+1})$$

由于

$$x^{\ast }-{\tilde{x}}_{k+1}\approx L(x^{\ast }-{\overline{x}}_{k+1})$$

将其与

$$x^{\ast }-{\overline{x}}_{k+1}\approx L(x^{\ast }-x_k)$$

联立，消去未知的 L，有

$$\frac{x^{\ast }-{\overline{x}}_{k+1}}{x^{\ast }-{\tilde{x}}_{k+1}}\approx \frac{x^{\ast }-x_k}{x^{\ast }-{\overline{x}}_{k+1}}$$

由此得

$$x^{\ast }\approx {\tilde{x}}_{k+1}-\frac{({\tilde{x}}_{k+1}-{\overline{x}}_{k+1}{)}^2}{{\tilde{x}}_{k+1}-2{\overline{x}}_{k+1}+x_k}$$

若以上式右端得出得结果作为新的改进值，则这样构造出得加速公式不再含有关于导数的信息，但它需要使用两次迭代值进行加工。其具体计算公式如下：

1. 迭代 ${\overline{x}}_{k+1}$ = $\phi (x_k)$
2. 迭代 ${\tilde{x}}_{k+1}=\phi ({\overline{x}}_{k+1})$
3. 改进 $x_{k+1}={\tilde{x}}_{k+1}-\frac{({\tilde{x}}_{k+1}-{\overline{x}}_{k+1}{)}^2}{{\tilde{x}}_{k+1}-2{\overline{x}}_{k+1}+x_k}$

上述方法称为埃特金（Aitken）加速方法。

使用埃特金加速算法解前例 求方程 $x^3-x-1=0$ 的唯一正根，如下图所示，埃特金算法的加速效果是显著的。

运行示例：

1. 未使用加速算法，需要迭代 9 次
   
2. 使用普通加速算法，需要迭代 4 次
   
3. 使用埃特金加速算法，仅需迭代 2 次
   

程序源码：

#include <iostream>
#include <cmath>

using namespace std;

/**
 * f(x) = (x+1)^{1/3)
 */
double f(double x)
{
    return pow(x + 1, 1.0 / 3);
}

int main(void)
{
    double x0;
    cout << "请输入迭代初值：";
    cin >> x0;

    double accrucy;
    cout << "请输入精度：";
    cin >> accrucy;

    int n;
    cout << "请输入您想要的最大迭代次数：";
    cin >> n;

    double x3;
    int count = 0;
    do
    {
        double x1 = f(x0);
        double x2 = f(x1);
        // x3 为利用埃特金加速公式处理后的近似解
        x3 = x2 - pow(x2 - x1, 2) / (x2 - 2 * x1 + x0);

        if (abs(x0 - x3) < accrucy)
        {
            cout << "近似解为：" << x3 << endl;
            break;
        }

        // 继续下一次迭代，直至找到符合精度要求的根或最大迭代次数用完
        cout << "第 " << ++count << " 次迭代！\t迭代函数的值为:" << x3 << "\t此次迭代精度为:" << abs(x3 - x0) << endl;
        double item = x0;
        x0 = x3;
        x3 = item;

        if (count > n)
        {
            cout << "迭代次数耗尽，迭代结束！未找到符合精度要求的根！！！" << endl;
            break;
        }

    } while (abs(x0 - x3) >= accrucy);

    return 0;
}