非线性方程(组)：一维非线性方程（二）插值迭代方法 [MATLAB]

　　通常而言，方程没有可以广泛求解的silver bullet，可是有几类方程的求解方法已经很是清晰确凿了，好比线性方程、二次方程或一次分式。一次方程能够直接经过四则运算反解出答案，二次方程的求根公式也给出了只须要四则运算和开根号的符号表达式。而一次分式的分子即为一次函数。更多的方程并无普适的符号表达式，但经过用便于求零点的函数模仿、代替之也能够估计零点的位置。插值方法能够实现这一思路。

　　插值迭代方法包括割线法、二次插值法等多项式插值方法，反插法以及线性分式插值法等等，其核心是用几个点及其函数值信息，经过插值产生一条容易计算零点的函数图线，而后用插值函数的零点来估计原函数的零点，不断迭代以达到足够的精度。每一个算法的大体思路均相同，再也不分别阐述具体原理。

1. 割线法(Secant Method)

　　用一次函数模拟原函数，用该一次函数的零点做为原函数零点的下一个估计。起始须要两个点。迭代式：$$x_{k+1}=x_k-f(x_k)\frac{x_k-x_{k-1}}{f(x_k)-f(x_{k-1})}$$　　迭代式和牛顿法的迭代式相似。实际上，割线法正是用两点的割线斜率代替了牛顿法中使用的切线斜率（即导数）：$f'(x_k)\approx [f(x_k)-f(x_{k-1})]/(x_k-x_{k-1})$ .

function [ zeropt, iteration ] = SecantMethod( func, start0, start1, prec )
%　　割线法求零点，函数句柄func，两个起点start0,start1，绝对精度prec
%   返回：零点zeropt，迭代步数iteration
prev = start0;
current = start1;
next = current - func(current)*(current - prev)/(func(current) - func(prev));
iteration = 0;
while abs(next - current) > prec && iteration < 500
    prev = current;
    current = next;
    next = current - func(current)*(current - prev)/(func(current) - func(prev));
    iteration = iteration + 1;
end
if iteration >= 500
    error('Method fails to converge');
end
zeropt = next;
end

　　进行试验：

% 用二次函数x^2-x-2
func = @(x)x^2-x-2;
[zero, iter] = SecantMethod(func, 6, 10, 0.0001)
% 输出
zero = 2.0000, iter = 7

% 输入
[zero, iter] = SecantMethod(func, -3, -9, 0.0001)
% 输出
zero = -1.0000, iter = 6

　　【迭代步复杂度】两次函数值求值+固定次数的四则运算（5）。

　　【收敛速度】超线性收敛，$r\approx 1.618$ 

　　【优点】1）每迭代步复杂度低，收敛速度中等；2）不用求导，只须要进行函数求值操做

　　【劣势】收敛速度不够快。

 

2. 二次插值法(Muller's Method)

　　用二次函数模拟原函数，将二次函数的根做为下一个零点估计值。每次迭代删去最老的一个点，利用包括新点的连续三个点进行插值。起始须要三个点。

　　二次插值面临的问题是，二次函数并不必定有实根。事实上，二次插值法的过程当中彻底不排斥出现复根，并且理论上依然能够收敛。至于对于二次方程的两个根，选取哪个，主要是在迭表明达式中为了减免减法所带来的有效数字损失而定的。

　　二次插值法相对复杂，其图像因为复根的出现也不直观，可是其收敛率却达到了：$r\approx 1.839$. Muller证实了任意m次多项式插值方法的收敛率r知足方程：$$r^{m+1}-r^m-...-r^2-r-1=0, \quad 或\quad r^{m+1}=\sum\limits_{k=0}^mr^k$$

 

3. 二次反插法

　　一样用二次函数模拟原函数，但反插法使用的是旋转90°的抛物线，即x关于y的二次函数。这样能够同时解决没有实根和选择恐惧症的问题，由于x关于y的二次函数与横轴一定有且只有一个交点。

　　根据Lagrange插值的方法，假设已有点 $(a, f_a),(b,f_b),(c, f_c)$ ，他们的二次反插应当为：$$x=a\frac{(y-f_b)(y-f_c)}{(f_a-f_b)(f_a-f_c)}+b\frac{(y-f_a)(y-f_c)}{(f_b-f_a)(f_b-f_c)}+c\frac{(y-f_a)(y-f_b)}{(f_c-f_b)(f_c-f_a}$$　　如今下一个点的纵坐标是零（做为零点的估计），则应当有：$$x=-\frac{af_bf_c(f_b-f_c)+bf_af_c(f_c-f_a)+cf_af_b(f_a-f_b)}{(f_a-f_b)(f_b-f_c)(f_c-f_a)}$$

 

function [ zeropt, iteration ] = InveQuadra( func, start0, start1, start2, prec )
%   二次反插求零点。输入函数句柄func，三个起始点start0-start2，要求精度prec
%   返回两个值：零点zeropt，迭代步数iteration
a = start0; b = start1; c = start2;
fa = func(a); fb = func(b); fc = func(c);
diffab = fa - fb; diffbc = fb - fc; diffca = fc - fa;
next = a*fb*fc/(diffab*diffca) + b*fa*fc/(diffab*diffbc) + c*fa*fb/(diffbc*diffca);
next = - next;
iteration = 0;
while abs(next - c) > prec && iteration < 500
    a = b; b = c; c = next;
    fa = func(a); fb = func(b); fc = func(c);
    diffab = fa - fb; diffbc = fb - fc; diffca = fc - fa;
    next = a*fb*fc/(diffab*diffca) + b*fa*fc/(diffab*diffbc) + c*fa*fb/(diffbc*diffca);
    next = - next;
    iteration = iteration + 1;
end
if iteration >= 500
    error('Method fails to converge.');
end
zeropt = next;
end

 

　　用一样的函数进行试验：

func = @(x)x^2-x-2;
% 输入
[zero, iter] = InveQuadra(func, -3, -9, -7,  0.0001)
% 输出
zero = -1.0000, iter = 5

% 输入
[zero, iter] = InveQuadra(func, 31, 16, 67,  0.0001)
% 输出
zero = 2.0000, iter = 8

func = @(x)x^3 - 20*x^2 - 25*x + 500;
% 输入
[zero, iter] = InveQuadra(func, -10, 10, -80,  0.0001)
% 输出
zero = 20.0000, iter = 4

　　能够看出，当应用于一样函数相似地靠近某一真值的种子时，二次反插比割线法收敛更快一点点。

　　【迭代步复杂度】须要三次函数求值+若干次四则运算。考虑到须要插值的是二次曲线，四则运算的数量显著高于割线法。

　　【收敛速度】和二次插值相同，二次反插的收敛率也是 $r\approx 1.839$.

　　【优点】1）不用求导。若是用牛顿法同样的思路，造成二次曲线不只要求一阶导，还要求二阶导！2）避免了二次插值的复根问题和选根问题；3）$r\approx 1.839$ 已经比较使人满意。

　　【劣势】计算复杂度相对于割线法来讲较大，且须要三个起始点。

 

 

4. 线性分式插值法

　　线性分式插值法使用形如 $\phi(x)=\frac{x-u}{vx-w}$ 的形式来插值以前的迭代点，并将插值函数的零点u做为新的零点的估计值。线性分式插值法主要就是为了求有水平或竖直方向渐近线的函数，这类函数亲测采用二次反插会有格外的困难，经常莫名其妙地跑到无穷远处去；仔细分析时证明，因为这类函数在很宽的范围内较为平坦，二次反插每每会造成很是尖锐的抛物线，该抛物线零点甚至很容易朝远离原点方向移动；如此数次迭代则趋于无穷。而线性分式插值法则能够解决这个问题。

　　代码实现：

 

function [ zeropt, iteration ] = InveQuadra( func, start0, start1, start2, prec )
%   二次反插求零点。输入函数句柄func，三个起始点start0-start2，要求精度prec
%   返回两个值：零点zeropt，迭代步数iteration
a = start0; b = start1; c = start2;
fa = func(a); fb = func(b); fc = func(c);
diffab = fa - fb; diffbc = fb - fc; diffca = fc - fa;
next = a*fb*fc/(diffab*diffca) + b*fa*fc/(diffab*diffbc) + c*fa*fb/(diffbc*diffca);
next = - next;
iteration = 0;
while abs(next - c) > prec && iteration < 500
    a = b; b = c; c = next;
    fa = func(a); fb = func(b); fc = func(c);
    diffab = fa - fb; diffbc = fb - fc; diffca = fc - fa;
    next = a*fb*fc/(diffab*diffca) + b*fa*fc/(diffab*diffbc) + c*fa*fb/(diffbc*diffca);
    next = - next;
    iteration = iteration + 1;
end
if iteration >= 500
    error('Method fails to converge.');
end
zeropt = next;
end

 

　　试验以下：

% 这个函数刚好是分式的形式，然而用二次反插困难重重
func = @(x)1 - 3/x;
% 输入
[zero, iter] = LinFraction(func, 1, 5, 10,  0.0001)
% 输出
zero = 3, iter = 1

func = @(x)9 - 1/x^2;
% 输入
[zero, iter] = LinFraction(func, 1, 5, 10,  0.0001)
% 输出
zero = 0.3333, iter = 6

　　

　　【迭代步复杂度】 和二次反插基本相同，三次函数求值及四则运算若干。

　　【收敛率】惊了，竟然也和二次插值/反插相同，为 $r\approx 1.839$ 

　　该算法的特色已经如上所述。