非線性方程(組)：一維非線性方程（一）二分法、不動點迭代、牛頓法 [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$ 對結果進行插值造成完整的函數。