數值分析方法

webliu6發表於2024-08-04

數值分析方法

數值方法可以用近似方法來理解

對於一些得不到解析解否則解析解計算難度太大的問題,如何用較少的計算得到相對較好(滿足精度要求)的數值解、近似解。
img
數值分析最基本的立足點是容許誤差

例子:
img
img

近似計算基本方法:

  1. 離散化(e.g. 數值積分)
  2. 遞推法(將一個複雜的計算過程轉化為簡單過程的多次重複 e.g. 上面的第二個例子)
  3. 近似替代法(將無限過程的數學問題用計算機的有限次計算近似替代 e.g. 上面的第一個例子Taylor公式不可能計算無限次,必須要近似替代)

不是解析解買不起,而是數值解更有價效比

接下來會介紹最常用的數學模型的最基本的數值計算方法

img

有關誤差的一些基本概念

誤差的來源與分類

img
img
前兩種誤差和數值方法本身無關,數值方法應重點關注截斷誤差和舍入誤差

演算法穩定性的概念: 初始的小擾動迅速積累,誤差呈遞增走勢。——不穩定的演算法!
誤差逐步遞減, 這樣的演算法稱為穩定的演算法。

誤差和誤差限

絕對誤差: \(\epsilon(x) = x^* - x ,其中x^*為近似值\)
可正可負可為0

準確值\(x\)一般無法確切知道,所以我們一般估計一個絕對誤差的上界,\(|\epsilon(x)| \le \eta\)👉絕對誤差限\(\eta\)
顯然是正值

相對誤差\(\epsilon_r(x) = \frac{x^* - x}{x}\)
可正可負可為0

img
同樣的前提下
img

相對誤差、相對誤差限是無量綱的量
相對誤差限比絕對誤差限能更好地刻畫近值數的近似程度

有效數字

img

img
準確到小數點後第 j 位 🟰 保留k + j 位有效數字 🟰 絕對誤差限不超過\(0.5*10 ^{−j}\)

\(x^*\)可表示為\(\pm0.x_1x_2....x_n * 10^m, 其中x_1 \ne 0\),
\(x^*\)有效數字的位數為 n,則有\(x^*\)的相對誤差限\(\eta_r\)\(\frac{1}{2x_1}*10^{1-n}\)
\(x^*\)的相對誤差限\(\eta_r\)\(\frac{1}{2(x_1 + 1)}*10^{1-n}\),則\(x^*\)有效數字的位數為 n

證明:
img
img

算術運算的誤差與誤差限

加減法:

img
img

乘法:

img

除法:

img
當除數絕對值較小時,商的絕對誤差的絕對值將非常大

冪、對數:

img
img

近似計算應遵循的原則

img
為什麼?
img
img
img

若干值得注意的問題

  • 防止兩個相近的數相減 👉 會使有效數字的位數嚴重損失
    img
  • 避免除數的絕對值 << 被除數的絕對值 👉 除數稍有一點誤差,則計算結果會出現較大變化
    img
  • 防止大數“吃掉”小數
    計算🐔計算浮點數時對階後可能出現較小數的尾碼變成全0
  • 簡化計算步驟,減少運算次數
    img

插值方法

img
已知一組離散資料,人們希望建立一個簡單的而便於計算的函式 g(x)使其近似的代替 f(x)

插值多項式的存在性與唯一性

img
img
該行列式可以用範德蒙行列式的公式計算,因為\(x_i\)是互異的節點,所以行列式的值不為0

可以用求解線性方程組的方法可以確定插值多項式 pn(x),
但是當 n 較大時,這種方法的計算量非常大。

拉格朗日(Lagrange)插值

n + 1個插值節點 👉 n + 1個限制條件 👉 n + 1個插值基函式 👉 n次拉格朗日差值多項式

插值基函式
img
插值基函式很好寫,該函式是\(n\)次多項式,且已知它的\(n\)個零點,常量係數只需要滿足當代入\(x_i\)時基函式值為1即可

拉格朗日插值多項式是插值基函式的線性組合

這樣寫出的n次拉格朗日插值多項式一定滿足所有限制條件,且插值多項式存在且唯一,故拉格朗日插值多項式就是我們想要的插值多項式

誤差分析:

img

設在\(x\)處的誤差\(f(x) - L_n(x)\)可以表示為 \(c\omega(x)\ 也即\ c\prod_{i =0}^{n}(x-x_i), 這裡c並不是一個常量,它是x的函式\)

接下來構造一個具有n + 2個零點的函式,且該函式包含誤差的表示式\(f(x) - L_n(x)\)

在插值區間\([a, b]\)中找到一個不是插值節點的\(x\),在該點的誤差為\(c_x\omega(x)\)

構造的函式為\(F(t) = f(t) - L(t) - c_x\omega(t)\)
這樣的函式在所有的插值節點和x處函式值為0,一共n + 2個零點

img

對錶達式\(F(t) = f(t) - L(t) - c_x\omega(t)\)左右兩邊對t求n+1階導數 (t與x沒有任何函式關係)
img
代入點\(\xi\)
img
注意這裡的\(\xi\)與我們選取的x有關,存在某種未知的函式關係

img
img
img

牛頓(Newton)插值

Lagrange 插值雖然簡單,但每增加一個節點,原有的插值基函式\(l_i(x)\) 必須重新計算,從而不具有承襲性 👉 增加一個節點原本的所有插值基函式都多了一個新的零點,且係數k也隨之改變

因此引入Newton插值
img
多一個節點,只要基於前一個 \(\phi(x)\) 構造一個新的 \(\phi(x)\) 加進多項式中即可,之前的係數也不需要改變(新加入的 \(\phi\) 函式在前面的插值節點處的函式值為0)

特徵: 在 \(x_i\) 處的函式值僅由0~i個 \(c_i * \phi(x)\) 決定

引入差商的概念給出Newton插值的插值公式
img

計算差商時畫出差商表方便計算
img

差商的一些性質:

差商 \(f[x_0,x_1,...x_k]\)\(f(x)\)\(x_0,x_1,...x_n\) 處函式值的一個組合,組合的係數是一個連乘形式的倒數
img
該性質的證明可以用數學歸納法,首先0階差商就是函式值顯然成立,之後由k = n-1 時成立推出 k = n時成立

由差商是函式值的線性組合可以看出,差商的值與各點的排列順序無關

差商 \(f[x,x_0,x_1,...x_k]\) 同樣可以寫成組合形式,只不過其中一個節點要替換成x,則它應該是x的函式

根據定義
img
可以確定,如果 \(f[x,x_0,x_1,...x_k]\) 是x的n次多項式,那麼 \(f[x,x_0,x_1,...x_k,x_{k+1}]\) 是x的n-1次多項式

這個可以從組合形式中定性地看出,對於x來說,f(x)的次數不變,分母的 \((x-x_0)(x-x_1)...(x-x_k)\) 多了一項,那麼次數應該減少一次
也可以從定義式中推斷,因為 \(f[x,x_0,x_1,...x_k,x_{k+1}]\)可以寫成\(\frac{f[x_0,x_1,...x_k,x_{k+1}]-f[x,x_0,x_1,...x_k]}{x_{k+1} - x}\),且已知\(f[x,x_0,x_1,...x_k]\)是n次多項式,\(f[x_0,x_1,...x_k,x_{k+1}]\)是常數,且分子在\(x=x_{k+1}\)時為0,故n次多項式含有\((x-x_k+1)\)因子,可以和分母約掉,故整體是n-1次多項式

進一步地
img

Newton插值用差商表示可以寫成
img
該形式在\(x_0\)點顯然成立,在其他任意點由差商的定義\(f[x,x_0] = \frac{f(x) - f(x_0)}{x - x_0}\)也顯然成立
同時 \(f[x,x_0,x_1] = \frac{f[x_0,x_1] - f[x,x_0]}{x_1 - x}\)推出\(f[x,x_0] = f[x,x_0,x_1]*(x - x_1) + f[x_0,x_1]\)
所以 \(f(x) = f(x_0) + f[x_0,x_1](x-x_0)+f[x,x_0,x_1](x-x_0)(x-x_1)\)
推廣到n+1個點時,就是
img
可以發現前面的\(c_n\)就是Newton插值的形式,後面的\(R_n\)是牛頓插值的餘項

由插值多項式的唯一性可知,同為n次多項式的Newton插值和Lagrange,同為n次多項式的Newton插值和Lagrange插值是完全相等的,故餘項也是相等的
img

赫密特(Hermite)插值

就是限制條件變多的Lagrange插值
現在要求在插值節點處插值多項式的導數值要和f(x)一致

同樣的思路寫出插值基函式

n + 1個插值節點 👉 2n + 2個限制條件 👉 2n + 2個插值基函式
img
插值基函式同樣很好寫,在 \(\alpha_i\)\(x_i\) 以外的點值為0,導數值也為0,故包含 \((x - x_j)^2\) 因子,因為是2n+1次多項式,待定係數法 \((ax+b)\) ,使得在 \(x_i\) 處插值基函式的值為1,導數值為0

\(\beta_i\)同樣含有2n次的因子,且\(\beta_i\)\(x_i\) 處函式值為0,故還含有 \((x-x_i)\) 因子,故待定係數求一常數即可,該常數使 \(\beta_i\)\(x_i\) 處導數值為1

誤差分析:

和Lagrange插值多項式的證明類似,得到插值餘項為
img

注意 \(R_{2n+1}(x)\) 在插值節點處的值為0,且 \({R_{2n+1}}'(x)\) 在插值節點處也為0 👈 \(\omega(x)\) 的作用

分段插值

由於Rouge現象的出現,我們將大區間分段,每段用由較少插值節點得到的低階插值多項式進行擬合
img
缺點:在分段處不夠光滑

誤差分析:

在每段用和Lagrange插值同樣的誤差分析進行分析

img
其中 \(h_i\) 是第i個分段的長度
其他分段插值的誤差分析與之類似,都是找到所有段中Lagrange誤差餘項最大的那個作為整體的誤差上界

曲線擬合

曲線擬合是給定函式類H,按照某種準則,找到一條曲線,既能反映給定資料的總體分佈形式,又不致於出現區域性較大的波動

和插值不同,插值只要在插值節點處滿足插值條件即可,而曲線擬合是為了反映給定資料的分佈
插值在區域性的插值節點是精確的,而曲線擬合在給定的點處允許誤差,期望在整體分佈上是準確的

曲線擬合的最小二乘法

最小二乘法的出發點:
img
找到的g(x)在函式類H中計算得到的殘差平方和一定是最小的

做法:待定係數法,根據最小二乘條件解方程

直線擬合
img
img

多項式擬合
img
img
img
img
img
img

數值積分

雖然數學題中的積分往往可以找到原函式用牛頓-萊布尼茲公式求出解析解,但實踐中常常找不到原函式或者被積函式都是以離散形式給出的(例如插值問題中,就不知道f(x)是什麼,但我們可能需要求其積分)

故根據積分的定義,採用計算機累加的方法進行計算
img

機械求積公式和代數精度

img
機械求積公式嘗試給出 \(f(\xi)\) 的近似值
img
img
img

誤差分析:
img
引入代數精度的概念
img
對m次多項式成立,只需要求積公式分別對 \(f(x)為x^0,x^1,..x^m\) 成立即可 👈 積分的線性性質,即 \(∫(f(x)+h(x))dx=∫f(x)dx+∫h(x)dx\)

img
定理的意義:給定若干節點,機械求積公式的積分精度是有保證的
img
同時也說明f(x)用n次Lagrange插值來近似後再進行積分計算,代數精度是n
機械求積公式中的 \(A_i\) 可以用 \(\int_a^bl_i(x)\) 計算得到

求積公式的構造方法

從前面的分析可知,存在特定的係數使k+1個節點組合得到的求積公式的代數精度至少為k

可以根據代數精度的判斷條件,解方程求得
img

也可以用n次Lagrange插值近似被積函式,同樣的代數精度為n
img
當被積函式本身就是不大於n次的代數多項式,那麼n次Lagrange插值多項式實際上等價於被積函式(Lagrange餘項為0),那麼插值型求積公式就是精確的

可以證明
img
如果代數精度至少為n次,那麼該公式對於 \(l_j(x)\) 是精確的,因為 \(l_j(x)\) 是n次多項式,那麼 \(\int_{a}^{b} l_j(x) \,dx = \sum_{i=1}^{n}A_il_j(x_i)\) 因為 \(l_j(x_i)\) 只有在 \(i=j\) 時為1,其他點都為0,則有 \(\int_{a}^{b} l_j(x) \,dx = A_j\) ,對於所有的n+1個插值基函式都有這一結論,根據積分的線性性質,可以知道求積公式是插值型的

另一方面,如果求積公式是插值型的,說明我們用n次Lagrange插值對被積函式進行了近似。那麼我們可以將被積函式 \(f(x)\) 分為兩部分,一部分是不大於n次的多項式,另一部分是高於n次的多項式,因為不大於n次的部分在近似後是可以精確求積的,故代數精度至少為n。

Newton-Cotes 求積公式

上面的機械求積公式對 \(x_i\) 的選取沒做限制,因為 \(A_i = \int_a^bl_i(x)\),則當節點改變或者求積區間改變時,求積公式中的係數 \(A_i\) 需要重新計算

但如果我們將求積區間獨立出來,且規定 \(x_i\) 是區間的n等分點,將 \(A_i\) 用Cotes係數 \(C_i\) 代替,可以發現係數既與求積區間 \([a,b]\) 無關,也與節點 \(x_i\) 的值無關,僅與總的節點個數n有關

  • 求積區間獨立出來
    img
  • 節點是區間的n等分點
    img

最後的積分形式只有一個n是可變的量,所以Cotes係數對於給定的節點滿足條件的積分都適用

記住常見的Cotes係數
img

誤差分析:

img
注意這裡的 \(f^{(n+1)}(ξ)\) 不能視為常數提到積分式之外,回憶Lagrange插值的誤差分析,\(ξ\)\(x\) 的函式,故也是 \(t\) 的函式
這也是求積公式誤差分析的難點所在

img
因為 \((t-0)(t-1)...(t-n)\) 中n為偶數,故該函式關於 \((\frac{n}{2},0)\) 中心對稱,該函式積分值為0。

穩定性分析

img
截斷誤差在上一節給出了通式
舍入誤差是由機器字長有限導致的

接下來分析Newton-Cotes公式的穩定性
img
這裡第一行的 \(≈\) 是截斷誤差(代數精度至少為n,則高於n次的項可能積分值有誤)

img

常見Newton-Cotes公式

梯形公式

img
img
img
這裡由第一行變到第二行用到了積分第二中值定理來將 \(f''(ξ)\)提出來,\(t^2-t\)\([0,1]\) 區間上是小於0不變號的
積分第二中值定理如下
img

辛普森公式

(三個節點,二階Lagrange,代數精度為三)
img
img
因為Simpson公式的代數精度為3,分析其誤差則需要將 \(f(x)\) 表示為某個三次插值多項式與其餘項之和,然後再求餘項的積分
如果使用Lagrange插值多項式,辛普森公式沒法像梯形公式那樣,在分析誤差時直接使用積分第二中值定理來將 \(f'''(ξ)\) 提出來,因為另外一部分 \((t-0)(t-1)(t-2)(t-3)\) 不滿足積分割槽間上不變號的條件

所以改用Hermite插值多項式來滿足不變號的條件
img
img
img

柯特斯公式

(五個節點,四階Lagrange,代數精度為五)
img
img

三種求積公式及其餘項

img

復化求積法

在一定範圍內,求積公式的代數精度隨插值節點的增加而提高,但高次插值會出現Rouge現象,導致近似效果變差,則對近似函式求積得到的結果也變差。
且n 較大時,Newton-Cotes 係數不容易求解,且出現負數項,高階Newton-Cotes 公式數值不穩定。
又從餘項的討論看到,積分割槽間越小,求積公式的截斷誤差也越小。

因此,我們經常把積分割槽間分成若干小區間,在每個小區間上採用次數不高的插值公式,構造出相應的求積公式,然後再把它們加起來得到整個區間上的求積公式,這就是復化求積公式的基本思想。

復化T公式

img
img

誤差分析:
復化求積公式總的誤差等於所有小區間的誤差之和
注意這裡的h是復化梯形公式的步長 \(\frac{b-a}{n}\),而不是\(b-a\)
img
最後一行使用了介值定理

復化S公式

區間 \([a,b]\) 實際上劃分為m塊,每塊用Simpson公式,Simpson公式的步長是 \(\frac{1}{2}\) 區間長度,故總共 \(n\) 等分,\(n = 2m\),步長為 \(\frac{b-a}{n}\)
Simpson公式的係數是 \(\frac{1}{6} ,\, \frac{4}{6}, \, \frac{1}{6}\)
img
img

誤差分析:
注意共有m個小區間,步長h是每一小步的長度
img

復化C公式

Cotes公式的係數是 \(\frac{7}{90} ,\, \frac{32}{90} ,\, \frac{12}{90} ,\, \frac{32}{90} ,\, \frac{7}{90}\)
img

誤差分析:
注意共有m個小區間,步長h是每一小步的長度
img

Romberg 求積公式

從上述誤差分析可以看到,復化積分的截斷誤差隨著分割因子 n 的增大而減小
如何根據精度要求選取合適的步長呢?

二分法:
以梯形公式為例
img
img

該過程收斂速度較慢

對復化公式進行分析
img
注意從第一行公式到第二行公式有一個近似替換,將大區間n等分後小區間上離散的函式值乘步長的累加近似為大區間上的積分
img
\(\overline{T}\) 展開合併同類項可以發現
img

定性地比較一下復化T,復化S和復化C公式的誤差
img
img
img
可以發現
復化T每二分一次,誤差近似變為原來的 \(\frac{1}{4}\)
復化S每二分一次,誤差近似變為原來的 \(\frac{1}{16}\)
復化C每二分一次,誤差近似變為原來的 \(\frac{1}{64}\)

經過上面的分析,得到Richardson(理查德森)外推表
img
img
在復化T公式二分的同時,利用Richardson 外推進行加速

數值微分

給定某未知函式上若干離散的點的值,期望得到該函式的導函式(近似導函式)

差商型數值微分

img
img
img
img
誤差為
img
img

插值型數值微分

利用插值多項式近似函式,再求導
img

誤差分析:
img
當x取插值節點時,即求插值節點處的導數值, 因為 \(\omega_{n+1}(x) = 0\), 誤差中的後一項為0

常見的插值型數值微分(一般寫出Lagrange插值再求導即可):

兩點公式
img
帶有餘項的兩點公式
img

三點公式
img
img

常微分方程數值解法

該類問題是給定一個常微分方程,求滿足該方程的函式
n階常微分方程形如
img
n階常微分方程一般會給定n個初值或者n對邊值

引言

img

存在唯一解的充分條件:
img

如何計算 \(y(x)\)
img

約定 \(y(x_n)\) 為待求函式 y(x) 在 \(x_n\) 處的精確函式值,\(y_n\) 為待求函式 y(x) 在 \(x_n\) 處的近似函式值

補充的數學知識:

img
img
img

多元函式f(x,u,t)對x求導,其中u是x的函式,t與x無關 👉 則 \(f'(x,u,t) = f_x(x,u,t) + f_u(x,u,t)u'(x)\), 同樣的求n階導可以視為對n-1階導函式求導
則上式中 \(\frac{\partial f(x_n,y_n)}{\partial x}\) 的計算結果就是 \(f_x(x_n,y_n)+f_y(x_n,y_n)y'(x_n) = f_x(x_n,y_n)+f_y(x_n,y_n)f(x_n,y_n)\)

尤拉方法

尤拉方法是常微分方程的一種數值解法
其基本思路是不求 \(y(x)\) 的解析解,而是對求解區域 \([a,b]\) 做等距剖分,求 \(y(x)\) 在剖分節點 \(x_n\) 上的近似值 \(y_n\)

因為我們已經知道了導函式的表達形式
這一點可以由初值加上在剖分割槽間上的數值積分做到
(當然也可以將導數轉化為剖分節點的差商)
img
h為剖分割槽間的長度

途徑一:化導數為差商的方法

img
img
img
img

單步尤拉法僅需要當前步的資訊來計算下一步,二步尤拉法需要當前步和前一步的資訊來計算下一步

途徑二:初值+數值積分

img
img
img
img
梯形公式尤拉法
img

常微分方程數值解是一組離散的函式值資料,它的精確表示式很難求解得到,但可以進行插值計算後用插值函式逼近 y(x)

區域性截斷誤差分析

img
注意這裡和代數精度的區別
在上一步 \(y_n\) 精確的前提下,截斷誤差就是剖分割槽間上數值積分的誤差。顯然剖分割槽間越大,那麼積分帶來的誤差也越大,故誤差可以用h的次方來衡量。
而代數精度是對於次數不大於n的代數多項式,求積公式準確成立。

區域性截斷誤差可以用泰勒展開進行分析
因為上一步準確,所以 \(y(x_n) = y_n\),且有常微分方程可知 \(y'(x_n) = f(x_n,y(x_n))\)

泰勒展開
\(y(x_{n + 1}) = y(x_n + h) = y(x_n) + y'(x_n)h + \frac{y''(ξ)}{2}h^2\)

  • 顯式一步尤拉
    \(y_{n+1} = y_n + hf(x_n,y_n) = y(x_n) + hf(x_n, y(x_n)) = y(x_n) + hy'(x_n)\)
    故誤差為 \(T_{n+1} = \frac{y''(ξ)}{2}h^2 = O(h^2)\)
    一階精度
  • 隱式一步尤拉
    $y_{n+1} = y_n + hf(x_{n+1},y_{n+1}) $
    這裡做一個不甚精確的假設:我們計算的 \(y_{n+1}\) 可以近似地當作 \(y(x_{n+1})\) 代入到 \(f\)
    \(y_{n+1} = y(x_n) + hf(x_{n+1}, y(x_{n+1})) = y(x_n) + hy'(x_{n+1})\)
    我們可以對 \(y'(x_{n+1})\) 進行泰勒展開
    得到 \(y_{n+1} = y(x_n) + h(y'(x_{n}) + O(h))\)
    \(y_{n+1} = y(x_n) + hy'(x_{n}) + O(h^2)\)
    故誤差 \(T_{n+1} =y(x_{n+1}) - y_{n+1} = \frac{y''(ξ)}{2}h^2 - O(h^2) = O(h^2) - O(h^2) = O(h^2)\)
    一階精度
  • 顯式二步尤拉法
    做法是類似地,將 \(y(x_{n-1})\)\(x_n\) 處進行泰勒展開,之後作差
    img
    二階精度
  • 梯形公式尤拉法(隱式單步)
    梯形公式就是顯式和隱式一步尤拉法的算數平均,證明不做贅述
    img
    二階精度

img

改進的尤拉方法

從上述例子可以看到,梯形法由於具有二階精度,其區域性截斷誤差比顯式尤拉法和隱式尤拉法小,但梯形
法是一種隱式演算法,需要迭代求解,計算量大。
顯式尤拉法是一個顯式演算法,雖然計算量較小,但是精度不高。

改進的尤拉方法則先用顯式尤拉法計算 \(y_{n+1}\),代入到梯形法的右式中求更精確的 \(y_{n+1}\)
改進的尤拉法也被稱為預報-校正法,第一個 \(y_{n+1}\) 是預報值,第二個 \(y_{n+1}\) 是校正值
改進的尤拉法具有二階精度
img
img
當然也可以多次進行預報-校正過程
img
幾何意義
根據\((x_n,y_n)\)的資訊估計一個值,根據估計值的資訊重新估計一個值,取兩者的平均
img

龍格-庫塔方法

根據上面的分析,改進的尤拉法在剖分割槽間中預報了一個值,精度較顯式尤拉法有所提高
龍格-庫塔方法的思想是在剖分割槽間內多預報幾個點 \(k_i\) 的值,並用其加權平均值作為 \(k^*\) 的近似值

img
p階R-K公式如何構造可以透過泰勒展開進行求解

以二階R-K公式構造過程為例
img
展開過程略
img
img
四個未知量,三個方程,有無窮多解,且滿足方程組的任意一組係數構成的R-K公式都是2階精度的

其中包括了改進的尤拉法和變形的尤拉法
img

三階四階的R-K公式也可以用一樣的方法求解

三階R-K公式,8個未知量,6個方程,無窮多組解
最常用的是庫塔公式
img

四階最常用的是標準四階龍格庫塔公式和吉爾公式
img
img

實際應用中一般選四階及以下的龍格-庫塔
四階以上龍格-庫塔公式的計算量太大,並且精度不一定提高,有時反而會降低

變步長的龍格-庫塔

我們知道步長小精度高,步長大計算量小,那麼如何選擇合適的步長呢?

變步長的龍格-庫塔就是為了解決這個問題

img
img
這一點很好理解,假設 \(y_{n+1}^{(h/2)}\) 的誤差是 \(a\) ,因為 \(y_{n+1}^{(h)}\) 的誤差差不多是 \(y_{n+1}^{(h/2)}\) 誤差的 \(16\) 倍,且這兩個點不一定在精確解在數軸上的點的同一側,故兩點的距離最小是 \(15a\),即 \(|y_{n+1}^{(h)} - y_{n+1}^{(h/2)}|>15a\),又 \(|y_{n+1}^{(h)} - y_{n+1}^{(h/2)}|<15\epsilon\),故 \(a<\epsilon\),也就滿足了精度要求

有了判斷所選步長是否滿足精度要求的方法,我們就可以根據解的性態來調整步長的大小:在變化平緩的部分,數值求解使用較大步長;而在變化劇烈的部分,則使用較小步長

收斂性與穩定性

收斂性和穩定性的分析較難
我們接下來僅對較簡單模型的收斂性和穩定性進行分析

收斂性

收斂性的定義
img
雖然我們之前的誤差分析中,誤差最小是 \(O(h^5)\),當 \(\lim\limits_{h \to 0}\) 時,誤差顯然為0,但我們之前分析的誤差都是區域性誤差,建立在上一步結果精確的前提下。雖然我們現在每一步的誤差都近似為0,但同樣地 \(\lim\limits_{h \to 0}\) 也意味著 \(\lim\limits_{n \to \infty}\),現在有無窮多步。
\(n\)\(h\) 應滿足 \(n*h\) 是一個定值 \((x_n-x_0)\)

接下來我們僅對單步顯式方法進行分析

判斷收斂的方法是收斂性定理
img

證明
img
首先可以將誤差拆分成區域性誤差和 \(|\overline y_{n+1} - y_{n+1}|\) 兩部分誤差的和
區域性誤差是 \(O(h^{p+1})\) ,後一部分誤差利用Lipschitz條件,可以變成常數與前一步誤差的積
最終得到一個遞迴式

img

img
利用 \(e^x>(1+x)\) 進行放縮,使 \(nh\) 湊到一起,變成常數
同時 \(|\epsilon_0| = 0\),使得最終的誤差為0

穩定性

穩定性的定義我們已經介紹過了,穩定性用來判斷演算法的舍入誤差(微小擾動)是否會積累(惡性增長)的概念。

img

這裡只研究 \(y'=\lambda y \, (\lambda < 0)\),也即 \(y = e^{\lambda x} \, (\lambda < 0)\) 的情況

顯式尤拉法:
img
注意這裡的假設,可以看出對該模型穩定性的分析只是一個定性的分析
img

隱式尤拉法:
img

方程求根

rt,本章研究方程求實根的數值方法(一般是次數大於4的多項式方程或超越方程,低次的多項式方程有求根公式)

三個步驟:
img

二分法

一種很經典的解法,原理是零點存在定理,前提是給定一個有根區間
img
img
在迭代的過程中維護一個性質,即左端點和右端點滿足 \(f(left)*f(right)<0\),每次迭代取區間的中點進行判斷,將有根區間二分
img
結果取有根區間的中點即可

img

迭代法

img
具體步驟:
img
可以證明收斂的迭代序列一定收斂於方程的解
img
注意中間的變換要求 \(φ(x)\) l連續

迭代法的幾何意義:
img
注意 \(y = φ(x)\)\(y = x\) 的交點就是根

該方法還存在兩個問題:

  1. 不是所有迭代序列都是收斂的,如何找到收斂的迭代序列(即如何選擇初值 \(x_0\) 和迭代方程 \(x=φ(x)\)
  2. 收斂序列的收斂是用極限來判斷的,但我們不可能無限算下去(即判斷什麼是否停止計算,也即給出迭代到第k步的誤差分析)

這兩個問題在下一節解決

迭代法的收斂性與誤差估計

本節分析迭代法的收斂性和迭代過程的誤差

迭代法收斂性

預測序列是否收斂 👉 收斂定理(充分條件)
收斂定理保證了每迭代一次,\(x\) 離精確解 \(x^*\) 都更近
img
自對映是為了證明根的存在性、導數絕對值小於1是為了證明根的唯一性和收斂性
另外,魏爾斯特拉斯定理說明閉區間上的連續函式一定在該區間上取得最大值和最小值,故只要 \(|φ'(x)|<1\),就一定存在 \(|φ'(x)| \le L < 1\)
img
使用微分第一中值定理證明
img
使用微分第一中值定理證明
img

當不滿足迭代定理時,可能需要改變迭代方程或者初值所在區間 \([a,b]\) ,如何判斷是迭代函式取得不好,還是區間 \([a,b]\) 取得過大?

引入區域性收斂性的概念
img
可以發現,如果滿足區域性收斂性,只要 \([a,b]\) 區間取得夠好,一定收斂,也就能確定不是迭代函式取得不好,而是區間 \([a,b]\) 取得過大

img

定理很好證明,首先區域性收斂性只需要證明收斂皆可,解的存在和唯一是顯然成立的。又 \(φ'(x)\) 是連續的,\(|φ'(x^*)| < 1\) 可以推出存在一鄰域使得鄰域內的點都滿足 \(|φ'(x)| \le L < 1\),則 \(|φ(x) - x^*| = |φ(x) - φ(x^*)| \le L(x - x^*)\),也就說明迭代一次之後新的 $x $ 更接近精確解了,故序列是收斂的

該定理可以簡單地判斷出迭代方程是否具有區域性收斂性,但要求提前知道準確解,這顯然是不可能的

故採用如下不太嚴格的準則來判斷迭代方程選得好不好
img

迭代過程誤差估計

在滿足收斂定理的前提下,誤差的分析很簡單
img

證明
img

意義
img

迭代法的收斂速度及加速處理

img
img
img

加速過程:
img
img
類似於預報-校正,又有點像理查德森外推的過程

img
從這個式子中求解出 \(x^*\)
或者直接參照結論
img

牛頓法

牛頓法是一種通用的構造迭代方程的方法,構造出來的迭代方程滿足區域性收斂且收斂速度為平方收斂(即 \(φ'(x^*) = 0\)

做法是:對於求解 \(f(x) = 0\) 的問題,對 \(f(x)\) 在近似解 \(x_0\) 處進行泰勒展開且展開到二階,則 \(f(x) ≈ f(x_0) + f'(x_0)(x - x_0) = 0\) , 這樣得到一個線性方程可以直接求解,若 \(f'(x_0) \ne 0\),則有新的解 \(x_1 = x_0 - \frac{f(x_0)}{f'(x_0)}\),我們接下來在 \(x_1\) 處進行泰勒展開,並重復這一過程。那麼我們的迭代公式就是 \(x_{n+1} = x_n - \frac{f(x_n)}{f'(x_n)}\),其中 \(f'(x) \ne 0\)

這樣構造出來的迭代公式滿足下面的性質:
img
img
注意性質應滿足的條件:單根+ \(f(x)\) 有連續二階導數

證明:
img
img

幾何意義
img

計算步驟
img
注意迭代終止的條件是 \(|x_{n+1} - x_n| < \epsilon\)

牛頓法初值的選取

牛頓法是一種區域性收斂的演算法,如果初值 \(x_0\) 選擇的不恰當,就有可能得不到收斂的迭代序列
為使牛頓法收斂,必須滿足:用迭代公式算出的 \(x_1\)\(x_0\) 更靠近準確根 \(x^*\)
根據牛頓法的幾何意義可以發現 \(f'(x)\) 為0或者很小,都不能很快收斂

img
img
img
條件二是為了讓條件一有意義

牛頓法的soft版本
img

the end ......

相關文章