[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^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;
matlab求导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$
该算法的特点已经如上所述。
发布评论