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


免責聲明!

本站轉載的文章為個人學習借鑒使用,本站對版權不負任何法律責任。如果侵犯了您的隱私權益,請聯系本站郵箱yoyou2525@163.com刪除。



 
粵ICP備18138465號   © 2018-2025 CODEPRJ.COM