使用“牛頓迭代法”求解方程


使用牛頓迭代法求解方程

盡管通過因式分解和利用求根公式可以很方便的得出多項式方程的根,但大多數時候這個多項式的次數都很高,計算將變得非常復雜,因此,我們必須轉向一些近似解法。

牛頓迭代法是其中最好的方法之一。從根本上說,牛頓迭代法通過一系列的迭代操作使得到的結果不斷逼近方程的實根。

首先,要選擇一個初始值x=x0,使得該初始值接近實根的值。然后,迭代計算如下的公式:

xi+1 = xi - f(xi) / f '(xi)

直到xi+1達到一個滿意的近似結果為止。在這個公式中,f(x)是要求解的多項式方程,而f '(x)是f(x)的導數。

多項式求導

多項式求導是微積分的基礎,現在讓我們來看看針對多項式求導的公式化描述。

要計算出多項式的求導結果,只需要對多項式的每一項套用如下兩個公式

d/dx * k = 0, d/dx *kxr = krx r-1

這里的k是為常數,r是有理數,x是未知數。符號d/dx表示求導,其中x是多項式中的變量

對於多項式中的每一常數項,套用第一個公式;否則,就用第二個公式。假設有如下函數:

f(x) = x3 + 5x2 +3x +4

要得到求導后的結果f '(x),對該多項式的前三項套用第二個公式,最后一項套用第1個公式,得到結果如下:

f '(x) = 1 * 3x(3-1) + 5 * 2x(2-1) + 3 * 1x(1-1) + 0 = 3x2 + 10x +3

有時候也有必要進行高階求導,即導數的導數。比如,f(x)的2階求導可記為f ''(x),它是對f '(x)的求導結果。同理,f(x)的3階求導可記為f '''(x),這是對f ''(x)的求導結果,以此類推。因此,在前面的例子中,如果要計算f(x)的2階導數的話,我們按照如下的方式對f '(x)求導即可:

f ''(x) = 3 * 2x(2-1) + 10 * 1x(1-1) + 0 =6x +10

理解1階和2階導數

理解1階和2階導數的意義,是正確使用牛頓迭代法非常重要的一點。

f(x)在點x=x0處的1階導數表示函數f(x)在點x0處的斜率。

1階導數決定函數f(x)是遞增的(從左到右上升)還是遞減的(從左右到下降)。如果求導結果為正,f(x)就是遞增的;而如果求導結果為負,f(x)就是遞減的;如果求導結果為0,則f(x)即不是遞增的也不是遞減的。求導結果的大小表示f(x)遞增或遞減的速率。如下圖所示,陰影部分表示函數的遞增區間,因此這些區域里函數的1階導數都為正值。而中間穿過橫軸的部分里函數是單調遞減的,因此這個區域內的1階導數所表示的斜率都為負值。

圖:函數f(x)的1階和2階導數的意義

f(x)在點x=x0處的2階導數表示函數在該點處的凹凸性,也就是說函數圖像是向上凸的還是向下凹的。

2階導數值的大小表示函數圖像的凹凸程度。在上圖的a和c中,虛線表示函數凹凸性發生改變的位置(曲線上凸部和下凹部的分界點,也稱為拐點),拐點必然是2階導函數與x軸的交點

另一種表示函數f(x)在某點x=c處的導數的方法是:按照點斜式表示出與f(x)相切於點c的直線方程。直線點斜式方程可寫作:

y - f(c) = f ' (c)(x-c)

因此,如果f(x) = x3 - x2 - 3x + 1.8如上圖a所示,那么與f(x)相切於c=1.5的直線方程就可以表示為:

y - ((1.5)3 - (1.5)2 - 3*1.5 + 1.8) = (3 * 1.52 -2 * 1.5 - 3)(x - 1.5)

y + 1.575 = 0.75(x - 1.5)

上圖d中畫出了這條切線。

為牛頓迭代法確定迭代初始值

牛頓迭代法中至關重要的一點是選出一個合適的初始迭代值x0為了使牛頓迭代法能夠收斂至我們要求的根,初始迭代值必須要足夠接近於所求的根。下面有兩條規則必須滿足:

1、為x0確定一個區間[a,b],使得該區間內有且只有一個根存在。為了實現這個目的,需要選擇兩個值a和b,使得f(a)和f(b)的符號不同,且f '(x)的符號不會改變。如果f(a)和f(b)有不同的符號,則表示該區間內至少存在一個根。如果f '(x)的符號在區間[a,b]上不變,則可確定該區間內只包含一個根,因為該函數在此區間內只可能單調的遞增或遞減。

2、使x0=a或x0=b,這樣f(x0)的符號就和f ''(x)在區間[a,b]上的符號相同這也說明了f ''(x)在區間[a,b]上的符號不會改變。回顧一下,f(x)的2次導數表示函數的凹凸性,如果f ''(x)的符號不變,則f(x0)就和f ''(x)擁有相同的符號。每經過一次牛頓迭代,就會在區間[a,b]上逐步逼近於根。如下圖。

圖:牛頓迭代法的收斂

 在上圖的4個部分中,f(x)以粗實線表示,a和b都表示為垂直的虛線。如果f(a)滿足前面給出的規定,則從a開始迭代,且切線會朝着根收斂的方向傾斜。如果f(b)滿足前面給定的規定,迭代就從b開始。且切線也會朝着根收斂的方向傾斜。

牛頓迭代法的工作過程

假設我們想找出方程f(x) = x3 - x2 - 3x + 1.8的根。根據函數圖像(如下圖),該方程會有三個根:一個在區間[-2,-1],另一個在區間[0,-1],第三個在區間[2,3]。一旦知道了方程的根的個數以及它們所在的區間,就根據上面的要求對每個區間做測試,以此挑出一個初始迭代值。要實現這一點,首先得知道如下信息:

f(x) = x3 - x2 - 3x + 1.8,f '(x) = 3x2 - 2x - 3,f ''(x) = 6x - 2

利用這些信息,我們發現[-1,-2]滿足第一個要求,因為f(-2) = -4.2,f(-1) = 2.8,且f '(x)在該區間內沒有改變符號(一直都是正的)。到這,我們知道在區間[-2,-1]內有且只有一個根。現在需要檢查是否滿足第二個要求,我們發現f ''(x)在區間上並沒有改變過符號(一直都是負的)。因為f(-2)=-4.2也是負值,所以選擇x0 = -2作為迭代的初始值。下圖展示了在該區間內迭代計算根的過程,其精度達到了0.0001。只經過5次迭代就得到了這個近似值。

圖:計算方程f(x) = x3 - x2 - 3x + 1.8 = 0的3個根的近似值,精確度0.0001

再來看看如何求出區間[0,1]上的根。我們能夠看到第一個要求在該區間上是能夠滿足的,但是在該區間的f''(x)的符號並不是不變的,因此這個區間不能滿足第二個要求。我們懷疑這個根更接近於1而不是0,因此我們接下來嘗試使用區間[0.5,1],這個可以解決不滿足第二個要求的問題了。f(0.5) = 0.175,f(1) = -1.2,f '(x)在這個區間內都是負值,因此第一個要求可以滿足。要完成第二個要求讓初始值x0 = 0.5,因為f(0.5)=0.175為正值,且同f ''(x)在區間[0.5,1]上的符號相同。上圖也展示了在區間[0.5,1]上迭代計算根的過程,其精度達到了0.0001。最后用了4次迭代就達到了根的近似值。計算第3個根的過程也與此類似。

方程求解的接口定義

root


int root(double (*f)(double x), double (*g)(double x), double *x, int *n, double delta)

返回值:如果找到了根返回0;否則返回-1。

描述:采用牛頓迭代法,根據給定的初始值來計算方程f的根。初始值由x[0]指定。函數f的導數由參數g指定。參數n代表迭代的最大次數。參數delta表示逐次逼近的差值,用該值來決定何時應該結束迭代。函數返回后,迭代過程中計算出的近似值都保存在數組x中,此時n代表數組x中的元素個數。由調用者負責管理同x相關聯的存儲空間。

復雜度: O(n),這里n表示調用者希望進行的最大迭代次數。

方程求解的實現與分析

求解形如f(x) = 0的方程,本質上就是要找出它的根。函數root采用牛頓迭代法以給定的初始迭代值開始逐次迭代逼近,從而找到實際根。

root函數包含單個循環,該循環采用牛頓迭代公式來連續計算根的近似值。在本節給出的實現中,f就代表需要求解的方程,g代表f的導函數。每一次迭代,我們要判斷當前得到的近似值是否滿足需求。如果當前得到的近似值與前一輪迭代得到的近似值之差小於delta所指定的值,則認為當前的近似值滿足要求。如果經過n次迭代后還是沒有找到滿足要求的近似根,則從函數中返回-1以表示沒有找到近似根。

root函數的時間復雜度為O(n),這里n表示調用者希望進行迭代的最大次數。最壞情況是當沒有找到所需要的近似根時,此時一共會進行n次迭代。

示例:方程求解的實現

#include <math.h>
#include "nummeths.h"

int root(double (*f)(double x), double (*g)(double x), double *x, int *n, double delta)
{
    int satisfied,i;
    
    i = 0;
    satisfied = 0;
    
    while(!satisfied && i+1 < *n)
    {
        /*x的下一個迭代值*/
        x[i+1] = x[i] - (f(x[i]) / g(x[i]));
        
        /*確定是否獲得期望的近似值*/
        if(fabs(x[i+1] - x[i]) < delta)
            satisfied = 1;
        
        /*為下一次迭代做准備*/
        i++;
    }
    /*即使沒有迭代,也表明一個值已經存儲在x中*/
    if (i==0)
        *n=1;
    else 
        *n=i+1;
    
    /*找到實根返回0,或者達到最大迭代次數仍沒有找到實根,返回-1*/
    if(satisfied)
        return 0;
    else 
        return -1;
}

 


免責聲明!

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



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