迭代法是一种求解非线性方程根的方法, 它通过构造一个迭代过程, 将一个非线性方程转化为一个等价的不动点方程, 然后通过迭代逼近不动点, 从而得到非线性方程的根.

迭代法的基本思想是将隐式方程转化为显式的计算公式, 然后通过迭代, 求方程近似根.

问题描述

已知方程

$$\mathbf{f}(\mathbf{x})=0$$

在${\mathbb{R}}^n$的某个区域$D$上有根, 求方程的近似解$\mathbf{x}$.

算法描述

迭代法要求将上述方程化为如下形式:

$${\mathbf{x}}_{n+1}=\mathbf{\phi }({\mathbf{x}}_n)$$

则可取$\mathbf{\phi }$为迭代函数. 而迭代函数的构造方式并不是唯一的, 不同算法的构造方式不同, 下面将逐个进行介绍.

不动点迭代法

不动点迭代法是一种求解非线性方程的数值方法, 其基本思想是将原方程变形为关于$x$的迭代方程, 然后通过不断迭代求解$x$的值.

一维情形

假设有非线性方程

$$f(x)=0$$

则可以通过恒等变形改写为

$$x=\phi (x)$$

因此, 我们可以得到如下不动点迭代公式:

$$x_{n+1}=\phi (x_n)$$

故求解原方程的解就转化为求解函数$\phi (x)$的不动点. 设$x_0$为方程的一个近似解, 则利用上述不动点迭代公式即可逼近原方程的解.

多维情形

在一维情形的基础下, 进一步考虑非线性方程组

$$\mathbf{f}(\mathbf{x})=0$$

和一维情形类似, 可将方程变形为

$$\mathbf{x}=\mathbf{\phi }(\mathbf{x})$$

故有高维不动点迭代公式:

$${\mathbf{x}}_{n+1}=\mathbf{\phi }({\mathbf{x}}_n)$$

牛顿迭代法

牛顿迭代法(Newton-Raphson method)是求解方程近似根的一种方法, 其基本思想是利用泰勒公式将方程进行线性化, 然后通过不断迭代, 使得误差逐渐缩小, 最终得到近似解.

单根情形

首先, 将$f(x)$在$x=x_0$处进行泰勒展开, 得到

$$f(x)=f(x_0)+f^{\prime }(x_0)(x-x_0)+\frac{f^{\prime \prime }(x_0)}{2}(x-x_0{)}^2+\cdots$$

其中$f^{\prime }(x)$和$f^{\prime \prime }(x)$分别表示$f(x)$在$x$处的导数和二阶导数.

令

$$e(x)=\frac{f(x)}{f^{\prime }(x)}$$

则有

$$e(x)=-\frac{f(x_0)}{f^{\prime }(x_0)}-(x-x_0)-\frac{f^{\prime \prime }(x_0)}{2f^{\prime }(x_0)}(x-x_0{)}^2-\cdots$$

因此, 可以通过迭代公式

$$x_{n+1}=x_n-\frac{f(x_n)}{f^{\prime }(x_n)}$$

不断迭代求解, 直到$e(x)$的值足够小(即达到所需的精度), 即可得到方程的近似解.

值得注意的是, 牛顿迭代法的前提是方程的导数, 即$f^{\prime }(x)$在所求解的区间内不等于零, 否则会导致方法失效. 此外, 迭代过程中可能会出现震荡或收敛于非解的情况, 此时需要采取适当的策略进行优化和调整.

重根情形

若$x_0$为$f(x)$的根, 且其重数大于$1$, 此时$f(x_0)$和$f^{\prime }(x_0)$都为$0$, 这给牛顿法和后面的割线法都带来了问题, 因为在它们的公式分母中都包含导数或其估算值, 当解收敛到离根非常近时, 可能会导致除零错误.

对此, 我们可以令

$$u(x)=\frac{f(x)}{f^{\prime }(x)}$$

则可以证明: 若$x_0$为原方程的$m$重根, 则其为方程

$$u(x)=0$$

的单根, 从而对$u(x)=0$应用牛顿迭代法即可.

割线法

设$f:\mathbb{R}\to \mathbb{R}$, 但函数$f$的性质不好, 可能在某些地方不可导或导数难以计算. 这时我们就可以利用割线近似切线, 即

$$f^{\prime }(x_n)\approx \frac{f(x_n)-f(x_{n-1})}{x_n-x_{n-1}}$$

代入牛顿迭代公式, 就得到了割线法:

$$x_{k+1}=x_k-\frac{x_n-x_{n-1}}{f(x_n)-f(x_{n-1})}f(x_k)$$

抛物线法

抛物线法又称为米勒法, 是割线法的推广. 割线法是利用前两项构造线性插值, 而抛物线法则是利用前三项构造二次插值, 则我们可以求得如下迭代公式:

$$x_{k+1}=x_k-\frac{2f_k}{\omega +\mathrm{sgn}(\omega )\sqrt{{\omega }^2-4f_{k}f[x_k,x_{k-1},x_{k-2}]}}$$

其中,

$$\omega =\frac{f_k-f_{k-1}}{x_k-x_{k-1}}+\frac{x_k-x_{k-1}}{x_k-x_{k-2}}\left(\frac{f_k-f_{k-1}}{x_k-x_{k-1}}-\frac{f_{k-1}-f_{k-2}}{x_{k-1}-x_{k-2}}\right)$$

逆二次插值法

在抛物线法中, 当构造的抛物线不与$x$轴相交时, 可能会收敛到方程的复数根. 这一方面体现了抛物线法的一个优点: 实的初始值迭代求得方程的复数根; 但同时也有可能无法收敛到实数根. 使用逆二次插值法可以解决这个问题.

逆二次插值法又称为反抛物线法, 即使用以$y$为自变量的抛物线$x=g(y)$与横轴的交点逼近函数$f(x)$的实根. 可计算得迭代公式如下:

$$x_{k+1}=x_k-y_k\frac{x_k-x_{k-1}}{y_k-y_{k-1}}+\frac{y_k{y}_{k-1}}{y_k-y_{k-2}}\left(\frac{x_k-x_{k-1}}{y_k-y_{k-1}}-\frac{x_{k-1}-x_{k-2}}{y_{k-1}-y_{k-2}}\right)$$

算法实现

准备工作

由于牵涉到高维问题, 我们首先编写一个向量之间的度量函数:

double distance(const std::vector<double> &x, const std::vector<double> &y)
{size_t n(x.size());if (!n)return 0;if (n != y.size())return NAN;double r(0);auto i = x.cend(), j = y.cend();do{double t(*--i - *--j);r += t *= t;} while (--n);return sqrt(r);
}

由于抛物线法中涉及到符号函数, 另外编写符号函数sgn:

// 符号函数
template <class T>
inline int sgn(const T &x) noexcept
{if (x == 0)return 0;if (x > 0)return 1;return -1;
}

此外, 为了方便用户输入函数, 还要编写一个字符串替换的函数:

/** 替换子串* str   : 要替换的字符串* oldStr: 旧子串* newStr: 新子串* except: 排除的子串*/
std::string &substring_replace(std::string &str, const std::vector<std::string> &oldStr, const std::vector<std::string> &newStr, const std::vector<std::string> &except)
{if (oldStr.empty()){str.clear();return str;}// if(oldStr.size()!=newStr.size())//     if(oldStr.size()>newStr.size())//         oldStr.erase(oldStr.begin()+newStr.size(),oldStr.end());//     else//         newStr.erase(newStr.begin()+oldStr.size(),newStr.end());size_t n(oldStr.size() > newStr.size() ? newStr.size() : oldStr.size());QBitArray a(str.length()); // 这里用了Qt中的QBitArray, 使用bool数组也可以实现其功能for (const auto &i : except){size_t pos(0);while ((pos = str.find(i, pos)) != std::string::npos){a.fill(true, pos, pos + i.length());pos += i.length();}}size_t pos(0);for (size_t k(0); k < n; ++k){const std::string &pre = oldStr[k], &next = newStr[k];F:while ((pos = str.find(pre, pos)) != std::string::npos){size_t i(pre.length()), p(pos);doif (a.at(p++)){pos += pre.length();goto F;}while (--i);QBitArray t(a.size() + next.size() - pre.size());dot.setBit(i, a.at(i));while (++i < pos);p = pos + pre.length();str.replace(pos, pre.length(), next);if (i != (pos += next.length()))dot.setBit(i, false);while (++i != pos);while (p != a.size())t.setBit(i++, a.at(p++));a = t;}}return str;
}

以及浮点向量转字符串向量的函数:

// double向量转string向量
std::vector<std::string> vec_to_string(const std::vector<double> &x)
{std::vector<std::string> y(x.size());auto i = x.cbegin();auto k = y.begin();while (i < x.cend())*k++ = *i++;return y;
}

一般迭代法

为了方便用户输入, 我们首先引入字符串函数计算¹:

extern double calStr(string);
extern double calStr_x(string, const double &);
vector<string> func{"sqrt", "sin", "log", "exp", "pow", "abs", "cos", "tan", "asin", "acos", "atan"};

接着实现一般迭代法的单步迭代:

/** 迭代法* current_vec: 迭代向量* func       : 迭代函数* x          : 未知数** 返回(bool):*  true : 迭代失败*  false: 迭代成功*/
bool iterative_method(std::vector<double> &current_vec, const std::vector<std::string> &f, const std::vector<std::string> &x)
{std::vector<double> x0(current_vec);unsigned n(current_vec.size());do{std::string toCal(f[--n]);substring_replace(toCal, x, vec_to_string(current_vec), func);if (isnan(current_vec[n] = calStr(toCal))){current_vec = x0;return true;}} while (n);return false;
}

最终, 实现一般迭代法如下:

/** 迭代法* current_vec: 迭代向量* func       : 迭代函数* x          : 未知数* epsilon    : 迭代终止条件* mid        : 迭代中间结果导出** 返回(bool):*  true : 迭代失败*  false: 迭代成功*/
bool iterative_method(std::vector<double> &current_vec, const std::vector<std::string> &func, const std::vector<std::string> &x, double epsilon, QString &mid) // mid使用了QString, 可以对应转换为std::string或者std::vector<vector<double>>
{for (auto &i : x)mid.append(QString::fromStdString(i)).append(',');*(mid.end() - 1) = '\n';for (auto &i : current_vec)mid.append(QString::number(i, 'g', 14)).append(',');std::vector<double> x0(current_vec), x1(x0);if (iterative_method(current_vec, func, x)){current_vec = x0;return true;}*(mid.end() - 1) = '\n';for (auto &i : current_vec)mid.append(QString::number(i, 'g', 14)).append(',');while (distance(current_vec, x1) >= epsilon){x1 = current_vec;if (iterative_method(current_vec, func, x)){current_vec = x0;return true;}*(mid.end() - 1) = '\n';for (auto &i : current_vec)mid.append(QString::number(i, 'g', 14)).append(',');}mid.chop(1);return false;
}

割线法

/** 割线法* f   : 原方程函数* x1  : 第1个值* x2  : 第2个值* e   : 精度* path: 迭代保存路径*/
void secant_method(const function<double(double)> &f, double x1, double x2, double e = 1e-6, const QString &path = QStandardPaths::writableLocation(QStandardPaths::DesktopLocation))
{QFile file(path);if (!(file.open(QIODevice::WriteOnly)))throw "文件保存失败！";double f1(f(x1)), f2;if (isnan(f1))throw "区间内存在奇点！";file.write((to_string(x1) + '\n' + to_string(x2)).c_str());while (abs(x1 - x2) >= e){if (isnan(f2 = f(x2))){file.close();throw "区间内存在奇点！";}double t = x2 - (x2 - x1) * f2 / (f2 - f1);file.write(('\n' + to_string(t)).c_str());x1 = x2;x2 = t;f1 = f2;}file.close();
}

抛物线法

/** 抛物线法* f   : 原方程函数* x1  : 第1个值* x2  : 第2个值* x3  : 第3个值* e   : 精度* path: 迭代保存路径*/
void parabolic_method(const function<double(double)> &f, double x1, double x2, double x3, double e = 1e-6, const QString &path = QStandardPaths::writableLocation(QStandardPaths::DesktopLocation))
{QFile file(path);if (!(file.open(QIODevice::WriteOnly)))throw "文件保存失败！";double f1(f(x1)), f2(f(x2)), f3, omega0((f2 - f1) / (x2 - x1));if (isnan(f1) || isnan(f2)){file.close();throw "区间内存在奇点！";}file.write((to_string(x1) + '\n' + to_string(x2) + '\n' + to_string(x3)).c_str());while (abs(x1 - x2) >= e){if (isnan(f3 = f(x3))){file.close();throw "区间内存在奇点！";}double omega1 = (f3 - f2) / (x3 - x2), d123 = (omega1 - omega0) / (x3 - x1), omega = omega1 + (x3 - x2) * d123, delta = omega * omega - 4 * f3 * d123;if (delta < 0){file.close();throw "迭代出现复根！";}double t = x3 - 2 * f3 / (omega + sgn(omega) * sqrt(delta));file.write(('\n' + to_string(t)).c_str());x1 = x2;x2 = x3;x3 = t;f1 = f2;f2 = f3;omega0 = omega1;}file.close();
}

逆二次插值法

/** 逆二次插值法* f   : 原方程函数* x1  : 第1个值* x2  : 第2个值* x3  : 第3个值* e   : 精度* path: 迭代保存路径*/
void inverse_quadratic_interpolation(const function<double(double)> &f, double x1, double x2, double x3, double e = 1e-6, const QString &path = QStandardPaths::writableLocation(QStandardPaths::DesktopLocation))
{QFile file(path);if (!(file.open(QIODevice::WriteOnly)))throw "文件保存失败！";double f1(f(x1)), f2(f(x2)), f3, omega0((x2 - x1) / (f2 - f1));if (isnan(f1) || isnan(f2)){file.close();throw "区间内存在奇点！";}file.write((to_string(x1) + '\n' + to_string(x2) + '\n' + to_string(x3)).c_str());while (abs(x1 - x2) >= e){if (isnan(f3 = f(x3))){file.close();throw "区间内存在奇点！";}double omega1 = (x3 - x2) / (f3 - f2), t = x3 - omega1 * f3 + (omega1 - omega0) * f3 * f2 / (f3 - f1);file.write(('\n' + to_string(t)).c_str());x1 = x2;x2 = x3;x3 = t;f1 = f2;f2 = f3;omega0 = omega1;}file.close();
}

实例分析

例1

设有方程$x^3+2x^2+10x-20=0$.

首先计算出$f^{\prime }(x)=3x^2+4x+10$, 下面我们讨论选取不同初值时Newton迭代法的收敛速度.

分别取 x 0 ∈ { x ∈ Z : − 128 ⩽ x ⩽ 127 } x_0\in\{x\in\mathbb Z:-128\leqslant x\leqslant 127\} x0​∈{x∈Z:−128⩽x⩽127}, 并取终止标准为$|f(x_n)|<10^{-6}$, 计算得其误差与迭代次数如下所示:

观察可得, 在$x_0=1$附近时迭代次数较少, 收敛较快, 接着在$1$附近取更密集的初值, 计算结果如下所示:

利用卡当公式, 我们可以求得方程有唯一实根

$$x=-\frac{2}{3}-\frac{13\sqrt[3]{4}}{3\sqrt[3]{3\sqrt{3930}+176}}+\frac{1}{3}\sqrt[3]{6\sqrt{3930}+352}=1.36880\cdots$$

可见当初值越靠近根$x$, 迭代次数越少, 收敛速度越快.

例2

取初值$x_0=4,x_1=3.8$, 求方程$x^3-2x-5=0$在$[1,4]$上的根.

仍然取终止标准为$|f(x_n)|<10^{-6}$, 使用Newton迭代法可得误差为$-1.77636\times 10^{-5}$, 迭代次数为6次; 取终止标准为$|x_{n+1}-x_n|<10^{-6}$, 使用割线法可得误差为$2.37144\times 10^{-13}$, 迭代次数为8次.

---

1. 参见C++实现简单计算器(字符串替换). ↩︎