任意方程求根

簡介

方程和函式是代數數學中最為重要的內容之一，從初中直到大學，我們都在研究著方程與函式，甚至我們將圖形代數化，從而發展出了代數幾何、解析幾何的內容。而在方程與函式中，我們研究其性質最多的，往往就是方程的根（零點），即使是研究方程的極值點、鞍點等，我們無非也只是研究其微商的零點。
我們在初等數學中已經學過許多簡單初等函式、線性方程的求解方法，在本文中，我們重點討論任意方程，尤其是計算困難的非線性方程的求根方法。

方程

分類和介紹

方程就是指含有未知數的等式。是表示兩個數學式（如兩個數、函式、量、運算）之間相等關係的一種等式，使等式成立的未知數的值稱為“解”或“根”。在這裡，根據一些性質的不同，我們將方程分成以下幾類：

• 單個方程
  • 線性方程：本質是等式兩邊乘以任何相同的非零數，方程的本質都不受影響。通常認為只含有一次項的方程。
  • 非線性方程：是因變數與自變數之間的關係不是線性的關係的方程。
    • 多項式方程
    • 超越方程：指含有未知量的超越式（指數、對數、三角函式、反三角函式等）的方程。換言之，超越方程中都有無法用含有未知數的多項式、分式或開方表示的式子。
• 多個方程
  • 線性方程組
  • 非線性方程組

方程的零點（根、解）

若有一個值或一些值能夠使得方程 \(f(x)=0\) 成立，那麼這個值就被成為方程的解，也常常被叫做零點和根。

若方程有且只有一個解\(x^*\)，那麼我們稱方程有單根\(x^*\)。

若對於方程\(f(x)=0\)

PS：若方程是簡單冪函式多項式組成，那麼方程的解的數量應和最高此項的數值一致，因為存在虛根。

求根方法

求根的方法基本上大同小異，都是通過區間去逼近方程的根的點。

首先我們說一個定理1：對於實函式方程\(f(x)=0\)，當\(x\in(a,b)\)，且\(f(x)\)在\(x\in(a,b)\)時單調且連續，若\(f(a)\cdot f(b)<0\)，則方程在\(x\in(a,b)\)有且只有一個根。

二分法

二分法和演算法中的二分搜尋法非常的類似，取定有根區間\([a,b]\)進行對分，求得\(mid = \frac{a+b}{2}\)進行判斷含根區間，如果\(f(a)\cdot f(mid)<0\),則令\(b=mid\)；反之若\(f(b)\cdot f(mid)<0\),則令\(a=mid\)。當\(|b_n-a_n|<\epsilon\)停止計算返回結果。

產生的截斷誤差為\(|e_{n-1}| = x_{n+1} - x^*\leq[b_n - a_n] = \frac{b_0 - a_0}{2^n}\)。

可以計算出最小迭代次數為\(n = \frac{lg(c_0-a_0)-lg\epsilon}{lg2}\)
程式碼實現：

private static double epsilon = 0.001;
// func為函式，寫法如x=>x*x+2*x-1，a，b必須為有效的含根開區間
public static double Binary(Func<double, double> func, double a, double b)
{
    var f1 = func.Invoke(a);
    var f2 = func.Invoke(b);
    if (f1 * f2 > 0)
        throw new ArgumentException("此區間無根或根不唯一");
    double mid = (a + b) / (double)2;
    var fm = func.Invoke(mid);
    if (fm == 0)
        return fm;
    if (f1 * fm < 0)
        b = mid;
    else if (f2 * fm < 0)
        a = mid;
    if (Math.Abs(b - a) <= epsilon)
        return (a + b) / (double)2;
    return Binary(func, a, b);
}

一般迭代法

迭代法是將方程求根問題轉換成了函式求交點問題再轉換成數列求極限的問題，這個過程中，對方程進行一個巧妙的變換之後，方程就可以在若干次迭代之後解出一個近似解。

操作方法如下：首先將方程\(f(x)=0\)改寫成\(x = \varphi(x)\)的形式，這樣就可以將方程的解看成是函式\(y=x\)和\(y=\varphi(x)\)的交點。給定初始值\(x_0\)後，則\(x_1 = \varphi(x_0)\)，不斷重複這個過程，若\(\displaystyle \lim_{k \to \infty}x_k\)存在，則迭代便可以達到使得\(x_k\)趨於交點。

通常，我們需要保證迭代函式在指定的含根區間內的導數值\(|\varphi^{'}(x)|<1\)，否則迭代函式將會不收斂。

迭代過程如下圖所示：

迭代法的收斂性證明

在這裡，我們將證明迭代法求根的合理性和可行性。

前提條件：設函式\(\varphi(x), x\in[a,b]\)有連續的一階導數，並且滿足以下條件：

• \(\forall x\in [a,b],\varphi(x)\in[a,b]\)

• \(\exist L \in (0,1),\forall x\in[a,b],|\varphi^{'}(x)|\leq L\)

要證明和解決以下命題和問題：

• 在\(x\in[a,b],\exist x^*,\varphi(x^*) = 0\)

• \(\forall x_0\in[a,b]\)，迭代過程\(x_{k+1} = \varphi(x_k)\)均收斂與\(x^*\)

• 求解誤差分析式

現在開始證明

1：證明在區間內有且只有一個根存在：

證：在\(x \in [a,b]\)時，\(\varphi^{'}(x)\)存在，所以有\(\varphi(x)\)連續，於是可以作\(g(x) = x - \varphi(x)\)，易知\(g(x)\)連續。

又因為\(\varphi(x) \in [a,b]\)，且\(g(a)*g(b)<0\),故存在實根，使得\(x=\varphi(x)\)。

利用反證法：若在[a,b]上還有一實根\(\bar{x}\)，那麼通過拉格朗日中值定理必定有：

\[x^*-\bar{x} = \varphi(x^*)-\varphi(\bar{x}) = \varphi^{'}(\xi)(x^*-\bar{x})\Longrightarrow \varphi^{'}(\xi) = 1 \]

顯然與假設不符合。

2：證明這個根收斂於\(x^*\)

根據拉格朗日中值定理，容易知道

\[x^*-x_{k+1} = \varphi(x^*)-\varphi(x_k) = \varphi^{'}(\xi)(x^*-x_k) \]

又因\(|\varphi^{'}(x)|\leq L\)，故易得：

\[|x^*-x_{k+1}|\leq|L(x^*-x_k)|=|L^2(x^*-x_{k-1})=\cdots=|L^k(x^*-x_0)| \]

因為\(\displaystyle \lim_{k \to \infty}L^k=0\)，故\(\displaystyle \lim_{k \to \infty}|x^*-x_{k+1}|=0(絕對值必定非負)\)，得\(x^*=x_{k+1}(k\to\infty)\)

3:迭代法的誤差式：

設某一次迭代後誤差為\(\epsilon\)，則可以推出：

\[|x_{k+1}-x_{k}| = |(x^*-x_k)-(x^*-x_{k+1})\geq|x^*-x_k|-|x^*-x_{k+1}|\geq|x^*-x_k|-L|x^*-x_k| \]

其中從\(|x^*-x_{k+1}| \leq L|x^*-x_k|\)則是因為每一次迭代後，誤差總是減少的。

故可以輕鬆的計算出誤差估計式為：

\[|x^*-x_k|\leq\frac{1}{1-L}|x_{k+1}-x_k| \]

通過中值定理，又可以推出：

\[|x_{k+1}-x_k| = |\varphi^{'}(\xi)(x_k-x_{k-1})| \]

因為\(|\varphi^{'}(\xi)|\leq L\)，將上式遞推後放縮成第二種誤差估計式：

\[|x^*-x_k|\leq\frac{L^k}{1-L}|x_1-x_0| \]

// func是迭代函式而不是實際函式
public static double Iterative(Func<double, double> func, double x)
{
    double temp = func.Invoke(x);
    if (Math.Abs(temp - x) <= epsilon)
    {
        return temp;
    }
    x = temp;
    return Iterative(func, x);
}

牛頓法

牛爵爺在整個微積分、數值分析中都有著舉足輕重的地位，這裡闡述的牛頓法，就是Taylor展開式的一部分。牛頓迭代法的核心思想就是：設法將一個非線性方程\(f(x)=0\)轉化為某種線性方程去求解。在數學意義上，我們知道泰勒公式可以將任意函式以簡單冪函式的形式表示出來，而在幾何意義上，我們是利用切線與X軸的交點去進行迭代，在根處，切線必過零點。若設\(f(x)=0\)的近似解為\(x_k\)，則方程可以用一階Taylor展開進行近似表示，牛頓迭代法的核心公式如下：

\[p_1(x) = f(x_k)+f^{'}(x_k)(x-x_k) \]

普通牛頓法

從上述公式中我們知道，其中\(p_1(x)\)就是泰勒多項式的表達，將其看成是方程\(f(x)=0\)，通過迭代的思想，從而得到了一個線性方程：

\[x_{k+1} = x_k - \frac{f(x_k)}{f^{'}(x_k)} \]

這就是普通牛頓法的迭代公式。在幾何意義上，他就是一個切線與x軸的交點去逼近零點，如圖：

牛頓下山法

所有的迭代法都有一個無法逃過的缺點，如果選取的初始值離根太遠，則會導致迭代過程次數過多，並且有可能導致迭代發散，因為迭代法只具有區域性收斂性。

為了避免迭代失敗或時間過長，我們加上這樣一個條件用於保證數值穩定下降

\[|f(x_{k+1})|<|f(x_k)| \]

將這個條件和牛頓法結合再一起，即再保證數值穩定下降的過程中，利用牛頓法加快迭代，這就是牛頓下山法。具體的操作如下：

將牛頓法的結果\(\bar{x}_{k+1}\)與前一步的迭代值\(x_k\)進行加權平均作為新的迭代值\(x_{k+1}\)，

則有：

\[x_{k+1} = \lambda\bar{x}_{k+1} + (1-\lambda)x_k \]

或

\[x_{k+1} = x_k - \lambda\frac{f(x_k)}{f^{'}(x_k)} \]

其中\(\lambda(0\leq\lambda<1)\)稱為下山因子，它的值是一個逐步試探的過程，可以從1開始取值，一旦滿足\(|f(x_{k+1})|<|f(x_k)|\)則稱為下山成功，否則需要另選初始值\(x_0\)進行試算。

簡單牛頓法

牛頓法可以說是一個非常有效的計算方法，準確度和迭代次數上都要比普通的迭代法要好得多，但是牛頓法最大的問題是我們需要求方程的導數，對於某些極其複雜的函式而言，導數是無法通過人工的方式計算，假如我們使用微積分的方式去求解導數，這對整個程式的效能會有比較大的影響，因此我們可以利用一個常數值\(\lambda\)來代替方程中的導數項。此時迭代公式為：

\[x_{k+1} =x_k - \frac{f(x_k)}{\lambda} \]

不過對於常數\(\lambda\)的取值是有限制的，因為我們需要保證迭代函式的收斂性，如果函式不收斂於\(x^*\)，那麼一切都沒有意義。於是有：

\[\varphi(x) = x - \frac{f(x)}{\lambda}\Longrightarrow\varphi^{'}(x) = 1-\frac{f^{'}(x)}{\lambda} \]

牛頓迭代法的收斂性遵循前文中普通迭代法的收斂性，於是可以得到：

\[|\varphi^{'}(x)| = |1-\frac{f^{'}(x)}{\lambda}|\Longrightarrow0<\frac{f^{'}(x)}{c}<2 \]

這樣我們就可以很輕鬆的確定下c的值了。

簡單牛頓法的幾何意義就簡單許多了，和我們之前討論的普通迭代法一致，只不過普通迭代法是將函式值和\(y=x\)進行處理變換，而簡單牛頓法則是和\(y= \lambda(x-x_k)+f(x_k)\)進行變換，本質是一致的。

///這裡只寫普通牛頓法，另外的函式由讀者自己思考
// 其中f1x為導數
public static double Newtown(Func<double, double> fx, Func<double, double> f1x, double x)
{
    var temp = f1x.Invoke(x);
    if (temp == 0)
    {
        throw new ArgumentException();
    }
    x = x - fx.Invoke(x) / temp;
    if (Math.Abs(fx.Invoke(x)) <= epsilon)
    {
        return x;
    }
    return Newtown(fx, f1x, x);
}

弦截法

弦截法是牛頓法的一種改進，保留了牛頓法中收斂速度快的優點，還克服了再牛頓法中需要計算函式導數值\(f^{'}\)的缺點。

弦截法中，我們用差商

\[\frac{f(x_k)-f(x_{k-1})}{x_k-x_{k-1}} \]

去代替牛頓法中的導數值，於是可以得到以下離散化的迭代

\[x_{k+1} = x_k - \frac{f(x_k)}{f(x_k)-f(x_{k-1})}(x_k-x_{k-1}) \]

這種方法叫做雙點弦截法，如圖所示：

從圖中可以知道，弦截法一直利用兩點之間的連線作為迭代的內容，那麼，他的合理性在哪裡呢？

整個迭代法都離不開中值定理，這裡也是這樣，事實上，這個差商之所以可以對導數值進行替代，是因為中值定理中說過，連續函式中兩函式上的點的連線的斜率必為兩點之間某一點的導數值，並且迭代過程中，這兩點的中值會越來越逼近函式零點，事實上這已經說明了弦截法是牛頓法的改進方法了。

不過，如果將函式中\(x_{k-1}\)用一個定點\(x_0\)代替，這種方法叫做單點弦截法，幾何意義如圖所示：

//這裡就對雙點弦截法進行程式設計
public static double StringCut(Func<double, double> func, double x1, double x2)
{
    var temp = x1 - (func.Invoke(x1) / (func.Invoke(x1) - func.Invoke(x2))) * (x1 - x2);
    x2 = x1;
    x1 = temp;
    if (Math.Abs(func.Invoke(x1)) <= epsilon)
    {
        return x1;
    }
    return StringCut(func, x1, x2);
}

收斂迭代的加速

除了普通迭代法意外，我們介紹的牛頓法、弦截法都是針對迭代法的改進，那麼對於普通迭代法應該如何去進行加速呢？首先我們應該對收斂速度進行一個研究。

收斂速度的計算

對於收斂的速度，我們很難直接從數值看出速度，因為收斂過程中，誤差的變化是越來越小的，因此我們再次進行一個比階。

\[\exists p\in N \And C>0,\displaystyle \lim_{k \to \infty}\frac{e_{k+1}}{e_k^p} = C \]

其中：

\[e_k = |x_k-x^*| \]

\[p= \begin{cases} p=1,線性收斂\\ p=2,平方收斂\\ p>2,高階收斂 \end{cases} \]

並且有一個定理：若\(\varphi(x)\)在所求的根\(x^*\)鄰域有連續的二階導數，並且有\(|\varphi^{'}(x)|<1\)，有以下特點：

1. 當\(\varphi^{'}(x)\neq0\)，那麼迭代過程線性收斂

2. 當\(\varphi^{'}(x)=0,\varphi^{''}(x)\neq0\)，那麼迭代是平方收斂的。

證明過程留給讀者，只需要利用泰勒展開是便可以證明該定理。

普通迭代加速

對於迭代過程，利用中值定理有

\[x^*-x_{k+1} = \varphi(x^*)-\varphi(x_k) = \varphi^{'}(\xi)(x^*-x_k) \]

當\(\Delta x \rightarrow0\)，我們將\(\varphi^{'}(\xi)\)看成定值\(\lambda(\lambda<1)\)

於是有

\[\lambda(x^*-x_k) = x^*-x_{k+1}\longrightarrow x^* = \frac{1}{1-\lambda}\bar{x}_{k+1}-\frac{\lambda}{1-\lambda}x_k \]

故可以推出最終的迭代公式為

\[\begin{cases} \bar{x}_{k+1} = \varphi(x_k)\\ x_{k+1} = \bar{x}_{k+1} + \frac{\lambda}{1-\lambda}(\bar{x}_{k+1}-x_k) \end{cases} \]

利用這種方式的好處就是再求得\(x_k\)時，利用加權的方式，使得迭代法變得有點類似像牛頓迭代法一樣變成切線性質的線而不是x=c或y=c的線，你經過畫圖和簡單的代數運算後，你會發現\(\bar{x}_{k+1}<x_{k+1}\)，也就是說達到了我們的加速目的。

Atiken加速法

Atiken加速法其實就是再普通加速迭代公式上在進行一次迭代，這裡直接寫出公式：

\[\begin{cases} \bar{x}_{k+1} = \varphi(x_k)\\ \bar{\bar{x}}_{k+1} = \varphi(\bar{x}_{k+1})\\ x_{k+1} = \bar{\bar{x}}_{k+1} - \frac{(\bar{\bar{x}}_{k+1}-\bar{x}_{k+1})}{\bar{\bar{x}}_{k+1} - 2\bar{x}_{k+1}+x_k} \end{cases} \]

Github

BiliBili主頁

WarrenRyan'sBlog

部落格園