非線性方程(組)：一維非線性方程（二）插值迭代方法 [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$ 

　　該算法的特色已經如上所述。