非线性方程(组)：一维非线性方程（一）二分法、不动点迭代、牛顿法 [MATLAB]

1. 二分法(Bisection)

1) 原理

　　【介值定理】 对于连续的一元非线性函数，若其在两个点的取值异号，则在两点间必定存在零点。

　　【迭代流程】 若左右两端取值不同，则取其中点，求其函数值，取中点和与中点取值异号的端点构成新的区间（其中必有零点）。进行下一次迭代。

2) 实现二分求根算法

　　使用MATLAB实现二分法代码如下。捕捉异常主要是为了在无法进行二分法的区间内发生输出zeropt为空的错误。

function [ zeropt ] = bisection( func, left, right, prec )
%　　二分法找非线性函数零点
%   输入4参数：函数句柄func，上/下界left/right，要求绝对精度prec
%   返回1参数：零点
leftVal = func(left);
rightVal = func(right);
if leftVal*rightVal >= 0        % 捕捉异常：若上下界处符号非相反
    error('Function vals at the boundaries are of the same sign, bisection unable to continue!');
end
flag = 0;
while (right - left) > prec     % 当区间长度大于精度时
    mid = (left + right)/2;
    midVal = func(mid);
    if midVal == 0              % 若中点值恰好为零则直接得到该值
        zeropt = mid;
        flag = 1;
        break;
    end
    if midVal*leftVal < 0       % 否则找到取值与中点取值异号的端点，将该区间作为新的查找区间
        right = mid;
    else
        left = mid;
        leftVal = midVal;
    end
end
if flag == 0                    % 若由于区间长度满足精度而退出循环，则取区间中点为零点
    zeropt = (left + right)/2;
end

　

　　对于任意满足要求的函数（连续）和区间（两端函数值异号），二分法保证能够找到一个满足要求的解。比如输入函数句柄，并执行调用二分法函数，就得到如下结果：

func = @(x)x^3 - 20*x^2 - 25*x + 500;
bisection(func, 0, 15, 0.0001)
ans = 5.0000

　　即算出该函数在(0, 15)区间上的一个零点是5.0000.

 

3) 二分法的迭代步复杂度、收敛速度及其他性质

　　【迭代步复杂度】一次求中点，一次函数值求值（由x求y为一次求值）。注意由于每次迭代时只有中点值是没算过的，其他两个值之前都算过可以继续使用，故每次迭代仅一次函数值求值。

　　【收敛速度】最大可预期的误差上限为区间长度，而每次迭代恒定地将区间长度减半。$\Rightarrow$ 线性收敛且收敛常数为0.5( $r=1, C=0.5$ ).

　　【二分法的优点】1）安全可控：区间总是减半缩小，在区间缩小过程中总是至少有一个零点存在于区间中；2）迭代复杂度很低，完全不考虑函数的性质，仅仅使用一次函数值求值。之后在其他更高效的方法之中，二分法经常被作为起安全保护作用的算法，避免其他算法迭代产生的点位置过于离谱。

　　【二分法的缺点】1）完全没有利用函数的性质！太可惜了，放弃了很重要的信息！2）收敛率仅为1，收敛速度缓慢。事实上是几种迭代方法中收敛最慢的两种之一；3）需要一个两端异号的区间作为开头，那如果是一个仅仅在中间较小的一段区域小于零，两侧均为正的函数呢？等到人工搜寻到异号的区间几乎都已经找到零点了！

 

2. 不动点迭代(Fixed Point Iteration)

1) 原理

　　【不动点】 $x=g(x)$ 的解称为函数 $g(x)$ 的不动点，几何图像为函数图像与正比例函数 $y=x$ 的交点。

　　【不动点迭代】 迭代形式为：$x_{k+1}=g(x_k)$ 。设函数的不动点为 $x^*$ ，当 $|g'(x^*)|<1$ 时，在不动点附近的存在一个领域，使得从领域内的某个值开始的不动点迭代收敛，收敛到不动点。通过 $\delta x_{k+1}\approx g'(x^*) \delta x_{k}$ 式，当导数绝对值小于1时，迭代后误差的线性主部较迭代前误差小，以导数绝对值为收敛常数线性收敛。唯一不同的是若导数值为零，此时该迭代至少二次收敛。用不动点迭代解非线性方程的核心在于将非线性方程转化为一个收敛的不动点迭代。

　　例如，求解一元二次方程：$f(x)=x^2-x-2$ 时，可以分别将其转化为解等价的不动点迭代：$g(x)=x^2-2,x=g(x)$或$g(x)=\sqrt{x+2},x=g(x)$ 。但是，前者在2附近发散，只有后者收敛。

 

2) 实现不动点迭代算法

　　不动点迭代的算法可谓没有任何技术含量，但构造不动点迭代倒是颇有些技术含量。以下为MATLAB代码，其中iteration只是用于标识迭代次数

function [ zeropt, iteration ] = fixed_point( func, seed, prec )
prev = seed;
next = func(prev);
iteration = 1;
while abs(next - prev) >= prec && iteration < 500
    prev = next;
    next = func(prev);
    iteration = iteration + 1;
end
zeropt = next;
if iteration >= 500
    error('Method fails to converge.');
end
end

　　

　　以上文的两个不动点函数为例：$g(x)=x^2-2, g(x)=\sqrt{x+2}$ ，均满足 $g(2)=2$ ，即2均为不动点。分别从偏离2一定距离的点开始进行试验：

func1 = @(x)x^2 - 2;
func2 = @(x)sqrt(x+2);

% 命令行输入
[zero, iteration] = fixed_point(func1, 1, 0.0001)
% 命令行输出
zero = -1, iteration = 2

% 命令行输入
[zero, iteration] = fixed_point(func1, 1.9, 0.0001)
% 命令行输出
错误使用 fixed_point (line 12)
Method fails to converge.

% 命令行输入
[zero, iteration] = fixed_point(func1, 2.1, 0.0001)
% 命令行输出
zero = Inf, iteration = 13

% 命令行输入
[zero, iteration] = fixed_point(func2, 1.9, 0.0001)
% 命令行输出
zero = 2.0000, iteration = 6

% 命令行输入
[zero, iteration] = fixed_point(func2, 2.1, 0.0001)
% 命令行输出
zero = 2.0000, iteration = 6

　　　　显示出前一个函数func1由于2附近不动点迭代不收敛，会出现：跳跃到另一个零点；超过迭代步数；上溢出 的情形。而func2在2附近开始进行不动点迭代总能快速达到2.

3) 不动点迭代的迭代步复杂度、收敛速度及其他性质

　　【迭代步复杂度】 一次函数值求值

　　【收敛速度】 如上所述，若导数为零则至少二次收敛，若导数非零则以导数绝对值为收敛常数线性收敛。

　　【优点】 1）迭代复杂度很低；2）若已知迭代函数，则实现起来没有任何技术含量，只需要不断调用函数即可。

　　【缺点】 1）收敛速度很低，除非导数恰好为零，否则收敛率仅为1；2）对于一个非线性方程，必须先找到它对应的一个收敛的不动点迭代，否则一切都没用；3）必须从距离真解足够近开始迭代才有效！不过这也不只是不动点迭代的问题。

 

3. 牛顿法(Newton)

1) 原理

　　设函数 $f(x)$ 零点的真解 $x^*$ ，即 $f(x^*)=0$，则在其一个领域内应有：$$f(x^*)\approx f(x^*-h)+f'(x^*-h)h=0$$　　那么假设函数表达式已知，而又已知当前的值x在真解附近，为了估计真解的值就可以利用上式反解出 $x^*$：$$x^*=x+h=x-\frac{f(x)}{f'(x)}$$　　若经过该计算，新的值距离真值更接近了，就可以把该过程写成不断迭代的形式来获得更好的精度，即：$$ x_{k+1}=x_k-\frac{f(x_k)}{f'(x_k)}$$　　牛顿法的几何图像为用当前点的切线与x轴的交点估计零点的位置。

2) 牛顿法的实现

　　在一般的迭代方法中，若要确定是否已经达到收敛，可以采用相邻两点间隔是否小于要求精度的方法，虽然这样的方法并不严格。在牛顿法的实现中，目前使用了待定精度（手动控制精度）的方法。另外，牛顿法中如果出现零斜率的情形（虽然在绝大多数情况下都不会出现）如何处理也是一个问题，在这里采用了“向右移动一个单位”选取下一个点的比较随意的方法。

function [ zeropt ] = nonlinNewton( func, prec, seed )
%	牛顿法求零点，输入3个参数：函数句柄func，精度prec，起始点seed；返回零点
syms f x;
f(x) = func(x);             % 生成符号函数类型f
fdiff(x) = diff(f);         % 生成f的符号函数类型导数fdiff
while fdiff(seed) == 0
    seed = seed + 1;
end
prev = seed;
next = prev - double(f(prev)/fdiff(prev));
while abs(next - prev) >= prec  % 执行迭代，直到精度符合要求
    prev = next;
    if fdiff(prev) == 0 % 若斜率为零则用任意算法取下一个点
        next = prev + 1;
    else
        next = prev - double(f(prev)/fdiff(prev));
    end
end
zeropt = next;
end

 

　　同样用测试不动点迭代方法的函数（原函数）$f(x)=x^2-x-2$ 进行测试：

func = @(x)x^2-x-2;

% 命令行输入
nonlinNewton(func, 0.0001, 1)
% 命令行输出
ans = 2.0000

%命令行输入
nonlinNewton(func, 0.0001, 3)
% 命令行输出
ans = 2.0000

　　

3) 牛顿法的迭代步复杂度、收敛速度及其他性质

　　【迭代步复杂度】 一次原函数的求值，一次导函数的求值，以及不超过三次的四则运算。需要注意的是，牛顿法需要求导函数，并且要求可以得到导函数在任意一点的取值，这本身将耗费很大的时间代价。以上实现的算法只适用于已知解析式（符号表达式）、可以求函数的导函数的情形。

　　【收敛速度】 考虑牛顿法和不动点迭代的关系，可知牛顿法的迭代公式等价于不动点迭代：$g(x)=x-f(x)/f'(x)$ 。按照不动点迭代收敛速度的分析方法，可以将 $g(x)$ 求导，结果得到 $g'(x)=f(x)f''(x)/(f'(x))^2$ ，在 $f(x)=0,f'(x)\neq 0$ 的条件下，必然有 $g'(x)=0$，因此在多数情况下牛顿法都是二次收敛的。而在零点处 $f'(x)=0$ 的情形不是别的，正是多重根的情形，该情形下问题本身的条件就属病态。

　　【优点】 收敛速度比较快，对于单根均满足二次收敛。

　　【缺点】 1）需要求导数，时间代价大；2）和不动点迭代类似，也需要选取合适的起点，这本身就是一件非常微妙玄学的事情。到了高维情况下，起点的选择还会显得玄学得多。

　　牛顿法用一条切线进行估值。可以想见：对于凹凸性不发生变化且存在零点的函数，牛顿法从任意起点都一定收敛。另外，对于线性函数，牛顿法一步迭代立即收敛，因为这步迭代形成的切线就是函数本身了。

　　最后，如果func函数只能得到它的数值输入输出而没有显式的符号表达式，如何得到可以确定它在任意一点导数值的函数呢？思路：先进行数值微分 $\rightarrow$ 转化为数值求导 $\rightarrow$ 对结果进行插值形成完整的函数。