4.11 - 15 牛顿迭代法
Newton迭代法
泰勒级数： \(\displaystyle f(x)=\sum_{n=0}^{\infty}\frac{f^{(n)}(a)}{n!}(x-a)^n\)
泰勒展开公式： \(\displaystyle f(x)=\sum_{n=0}^n\frac{f^{(n)}(a)}{n!}(x-a)^n+R_n(x)\)
- 皮亚诺型余项： \(R_n(x)=o[(x-a)^n]\) ，即当 \(x\to a\) 时，余项为 \((x-a)^n\) 的高阶无穷小
- 拉格朗日型余项： \(R_n(x)=\displaystyle\frac{f^{(n+1)}(\theta)}{(n+1)!}(x-a)^{(n+1)},\ \theta\in(a,x)\)
- 积分型余项： \(R_{n}(x)=\displaystyle\int_{a}^{x}{\frac{f^{{(n+1)}}(t)}{n!}}(x-t)^{n}\,dt\)
原理：将非线性方程线性化——Taylor展开
取 \(x_0\) 作为初始近似值，将 \(f(x)\) 在 \(x_0\) 处做一阶Taylor展开： \[ \begin{array}{c}f(x)=f(x_0)+f'(x_0)(x-x_0)+\displaystyle{f''(\xi)}{2!}(x-x_0)^2,\quad \xi\in(x_0,x)\\ 0=f(x^*)\approx f(x_0)+f'(x_0)(x-x_0)\quad\Rightarrow\quad x^*\approx x_0-\displaystyle\frac{f(x_0)}{f'(x_0)}\\ \left\{\begin{array}{l}x_1=x_0-\displaystyle\frac{f(x_0)}{f'(x_0)}\\ x_{k+1}=x_k-\displaystyle\frac{f(x_k)}{f'(x_k)}\quad\leftarrow\ \textbf{Newton迭代公式}\end{array}\right.\end{array} \]
Newton迭代法的收敛性
定理：设 \(f\in C^2[a,b]\) （二阶连续可微），若 \(x^*\) 为 \(f(x)=0\) 在 \([a,b]\) 上的根，且 \(f'(x^*)\neq0\) ，则Newton迭代法是二阶收敛的，且有 \(\displaystyle\lim_{k\to\infty}\frac{x_{k+1}-x^*}{(x_k-x^*)^2}=\frac{f''(x^*)}{2f'(x^*)}\) 。
初值的选取：令 \(c=\displaystyle\frac{\max|f''(x)|}{2\min|f'(x)|}\) ，则有： \[ c|x_{k+1}-x^*|\le(c|x_{k}-x^*|)^2\le(c|x_{k-1}-x^*|)^4\le\dots\le\le(c|x_{k+1}-x^*|)^{2^{k+1}} \] 因此， \(c|x_0-x^*|=1\ \Rightarrow\ |x_0-x^*|\le\displaystyle\frac{2\min|f'(x)|}{\max|f''(x)|}\) 时，Newton迭代法收敛。
Newton下山法
调整 \(x_0\) 的选取来使得Newton迭代法满足收敛条件。
定义：对Newton迭代过程附加单调性要求： \(|f(x_{k+1})|<|x_k|\) ，满足该条件的Newton迭代法称为Newton下山法。
实现：若由 \(x_k\) 得到的 \(x_{k+1}\) 不能使得 \(|f|\) 减小，则在 \(x_k\) 和 \(x_{k+1}\) 之间找点 \(\overline{x_{k+1}}\) ，使得 \(|f(\overline{x_{k+1}})|<|f(x_k)|\) ， \[ \begin{aligned}\overline{x_{k+1}}&=\lambda x_{k+1}+(1-\lambda)x_k\\ &=\lambda[x_k-\frac{f(x_k)}{f'(x_k)}]+(1-\lambda)x_k\\ &=x_k-\lambda\frac{f(x_k)}{f'(x_k)}\end{aligned} \] 其中， \(\lambda=1\) 时即Newton迭代法，当 \(\lambda=1\) 效果不好时，将 \(\lambda\) 减半计算。
Newton迭代法的变形
Newton下山法中计算每次迭代都需要计算一阶导数，试图简化计算。
简化Newton迭代法
采用迭代格式： \(x_{k+1}=x_k-\displaystyle\frac{f(x_k)}{M},\ k=0,1,2,\dots\) ，即用常数来代替一阶导，通常取 \(M=f'(x_0)\) 。一般，简化Newton迭代法只具有线性收敛。
割线法
采用迭代格式： \(x_{k+1}=x_k-\displaystyle\frac{f(x_k)}{f(x_k)-f(x_{k-1})}\ k=0,1,2,3,\dots\) ，需要取两个初值 \(x_0,x_1\) 。收敛阶 \(p\approx1.618\) 。
求重根的Newton迭代法
重根：称 \(x^*\) 为方程 \(f(x)=0\) 的 \(m\) 重根时，是指 \(f(x)=(x-x^*)^mh(x)\) ，其中 \(h(x)\) 在 \(x=x^*\) 处连续且 \(h(x^*)\neq0\) ，若 \(h(x)\) 在 \(x^*\) 处充分可微，则 \(f(x^*)=f'(x^*)=\dots=f^{(m-1)}(x^*)=0,\ f^{(m)}(x^*)\neq0\) 。
由于 \([f(x)]^{\frac{1}{m}}=(x-x^*)[h(x)]^{\frac{1}{m}}\) ，知 \(x^*\) 恰是方程 \([f(x)]^{\frac{1}{m}}\) 的单根。应用Newton迭代法对该方程求解，得到： \[ \begin{aligned}x_{k+1}&=x_k-\frac{[f(x_k)]^{\frac{1}{m}}}{\frac{1}{m}[f(x_k)]^{\frac{1}{m}-1}f'(x_k)}\\ &=x_k-m\frac{f(x_k)}{f'(x_k)},\quad k=0,1,2,\dots\end{aligned} \] 也称之为带参数的Newton迭代法，求方程 \(f(x)=0\) 的 \(m\) 重根时具有平方收敛。
当 \(m\) 未知时：
根据函数 \(\displaystyle u(x)=\frac{f(x)}{f'(x)}=\frac{(x-x^*)h(x)}{mh(x)+(x-x^*)h'(x)}\) ，可见 \(x^*\) 恰是方程 \(u(x)\) 的单根。对之应用Newton迭代法有： \[ x_{k+1}=x_k-\frac{u(x_k)}{u'(x_k)}=x_k-\frac{f(x_k)f'(x_k)}{[f'(x_k)]^2-f(x_k)f''(x_k)},\quad k=0,1,2,\dots \] 在该迭代过程中，不需要知道根的重数，具有平方收敛。
代码：Newton迭代法（三种）
//
// Created by xa on 2021-03-04.
//
#include <iostream>
#include <functional>
#include <cmath>
//已知一阶导函数求单根
void iterate(std::function<double (double)> const& f, std::function<double (double)> const& f_1, double x_0);
//已知二阶导函数求重根（重数任意）
void iterate(std::function<double (double)> const& f, std::function<double (double)> const& f_1, std::function<double (double)> const& f_2, double x_0);
//未知导函数，有二初值，利用割线法求单根
void iterate(std::function<double (double)> const& f, double x_0, double x_1);
int main() {
auto f = [&](double x) -> double {
return pow(x,4) - 8.6*pow(x,3) - 35.51*pow(x,2) + 464.4*x - 998.46;
};
auto f_1 = [&](double x) -> double {
return 4*pow(x,3) - 3*8.6*pow(x,2) - 2*35.51*x + 464.4;
};
auto f_2 = [&](double x) -> double {
return 3*4*pow(x,2) - 2*3*8.6*x - 2*35.51;
};
iterate(f, f_1, 7);
iterate(f, f_1, f_2, 4);
iterate(f, 7, 8);
}
void iterate(std::function<double (double)> const& f, std::function<double (double)> const& f_1, double x_0) {
double x = x_0; double f_x = f(x); double f_1_x = f_1(x);
double next_x = x - f_x / f_1_x;
double lambda = 1;
int count = 0;
while (std::abs(next_x - x) >= 1e-6) {
++count;
x = next_x; f_x = f(x); f_1_x = f_1(x);
std::cout << "第" << count << "次：" << x << std::endl;
lambda = 1;
next_x = x - f_x / f_1_x;
while (std::abs(f(next_x)) > std::abs(f_x)) {
lambda /= 2;
next_x = x - lambda * f_x / f_1_x;
}
}
}
void iterate(std::function<double (double)> const& f, std::function<double (double)> const& f_1, std::function<double (double)> const& f_2, double x_0) {
double x = x_0; double f_x = f(x); double f_1_x = f_1(x); double f_2_x = f_2(x);
double next_x = x - f_x*f_1_x/(f_1_x*f_1_x-f_x*f_2_x);
double lambda = 1;
int count = 0;
while (std::abs(next_x - x) >= 1e-6) {
++count;
x = next_x; f_x = f(x); f_1_x = f_1(x); f_2_x=f_2(x);
std::cout << "第" << count << "次：" << x << std::endl;
lambda = 1;
next_x = x - f_x*f_1_x/(f_1_x*f_1_x-f_x*f_2_x);
while (std::abs(f(next_x)) > std::abs(f_x)) {
lambda /= 2;
next_x = x - lambda * f_x*f_1_x/(f_1_x*f_1_x-f_x*f_2_x);
std::cout << "OK";
}
}
}
void iterate(std::function<double (double)> const& f, double x_0, double x_1) {
double pre_x = x_0; double f_pre_x = f(pre_x);
double x = x_1; double f_x = f(x);
double g_x = f_x - f_pre_x;
double next_x = x - f_x / g_x;
double lambda = 1;
int count = 0;
while (std::abs(next_x - x) >= 1e-6) {
++count;
pre_x = x; x = next_x; f_pre_x = f(pre_x); f_x = f(x); g_x = f_x - f_pre_x;
std::cout << "第" << count << "次：" << x << std::endl;
lambda = 1;
next_x = x - f_x / g_x;
while (std::abs(f(next_x)) > std::abs(f_x)) {
lambda /= 2;
next_x = x - lambda * f_x / g_x;
}
}
}