Re: [分析] Hermite內插演算法的證明

[Vulpix (Sebastian)]

符號看得有點花……


如果你想做的是「在 x_1 和 x_2 分別趨近 x_0 後所得的極限 = Taylor 多項式」，

那你需要的就是 MVT of divided differences。

https://en.wikipedia.org/wiki/Mean_value_theorem_(divided_differences)

直接套上去就馬上做完了。

也不必去算新多項式的導數。

但是如果要一步一步來就沒那麼好算了。

  (x_n - x_0)lim_{x_1→x_0} f[x_0,...,x_n]

= f[x_0,x_2,...,x_n] - lim_{x_1→x_0} f[x_0,...,x_{n-1}]

上面這條遞迴式是用來算極限的。

本來的插值多項式是 f[x_0] + f[x_0,x_1](x-x_0) + f[x_0,x_1,x_2](x-x_0)^2。

在 x_1→x_0 之後，變成

    f(x_0) + f'(x_0)(x-x_0) + {f[x_0,x_2]-f'(x_0)}/(x_2-x_0) * (x-x_0)^2。

然後 {f[x_0,x_2] - f'(x_0)}/(x_2 - x_0) 在 x_2→x_0 下的極限 = f"(x_0)/2，

所以多項式的極限就變成 f(x_0) + f'(x_0)(x-x_0) + f"(x_0)/2 * (x-x_0)^2。


不過我本來在想的是用 Lagrange 觀點。

e_0(x) := Π_{i=1}^n (x-x_i)/(x_0-x_i)，其他 e_j 類推。

還是先用 n = 2 來觀察，

插值多項式 = f(x_0)e_0(x) + f(x_1)e_1(x) + f(x_2)e_2(x)

然後也先讓 x_1→x_0，多項式變成

    f(x_0){e_0(x)+e_1(x)} + {f(x_1)-f(x_0)}e_1(x) + f(x_2)e_2(x)。

所以我們分成三項來看：

1. e_0(x)+e_1(x) = (x-x_1)(x-x_2)/(x_0-x_1)(x_0-x_2)

                      + (x-x_0)(x-x_2)/(x_1-x_0)(x_1-x_2)    公因式

                 = (x-x_2)/(x_1-x_0)

                      * (x_1-x_0)(x_0-x_2+x_1-x)/(x_0-x_2)(x_1-x_2)

                 = (x-x_2)(x_0-x_2+x_1-x)/(x_0-x_2)(x_1-x_2)

                 → (x-x_2)(2x_0-x_2-x)/(x_0-x_2)^2

                 = 1 - (x-x_0)^2/(x_0-x_2)^2

   最後這個多項式，他代 x_0 得 1、導數得 0，而代 x_2 得 0。

2. {f(x_1)-f(x_0)}e_1(x) = {f(x_1)-f(x_0)}(x-x_0)(x-x_2)/(x_1-x_0)(x_1-x_2)

                         → f'(x_0)(x-x_0)(x-x_2)/(x_0-x_2)

   (x-x_0)(x-x_2)/(x_0-x_2) 代 x_0 得 0、導數得 1，而代 x_2 得 0。

3. e_2(x) = (x-x_0)(x-x_1)/(x_2-x_0)(x_2-x_1) → (x-x_0)^2/(x_2-x_0)^2

   最後這個多項式也是代 x_0 得 0、導數得 0，而代 x_2 得 1。

我們得到三個可以各自突顯 f(x_0), f'(x_0), f(x_2) 的多項式，

剛好跟 Lagrange 觀點有謀而合。

最後再讓 x_2→x_0，

  f(x_0){1-(x-x_0)^2/(x_0-x_2)^2}

      + f'(x_0)(x-x_0)(x-x_2)/(x_0-x_2) + f(x_2)(x-x_0)^2/(x_2-x_0)^2

= f(x_0)+f'(x_0)(x-x_0)

      + { f(x_2) - f(x_0) - f'(x_0)(x_2-x_0) }(x-x_0)^2/(x_2-x_0)^2

→ f(x_0) + f'(x_0)(x-x_0) + f"(x_0)(x-x_0)^2/2

其實仔細看，1, x-x_0, (x-x_0)^2/2 也是

在函數值、一階導數、二階導數之中各自突顯一項，而消滅其他兩項的多項式函數，

同樣符合 Lagrange 觀點的插值概念。


真正麻煩的還是 general case：

    有資料的點是 x_0, ..., x_n，每個點的高階導數已知階數不盡相同。

像是已知 f(-1), f(0), f'(0), f"(0), f(5), f(100), f'(100) 這樣。

然後先用 -1, 0, a, b, 5, 100, c 插值，再讓 a,b→0 和 c→100，

之後要檢查在 x = 0 的一階二階導數和在 x = 100 的一階導數。

不過我想，應該也是這樣一步步算極限就好。

但是那個 general form 就真的很難看，所以平常都是給 algorithm。

--
※ 發信站: 批踢踢實業坊(ptt.cc), 來自: 1.162.224.247 (臺灣)
※ 文章網址: https://www.ptt.cc/bbs/Math/M.1691321340.A.CC3.html
推 znmkhxrw    : 嗨V大, 我想證的是「在 x_1 和 x_2 分別趨近 x_0 後
→ znmkhxrw    : 所得的極限L(x)」會滿足L(x_0)=f(x_0),
→ znmkhxrw    : L'(x_0)=f'(x_0), L''(x_0)=f''(x_0)
→ znmkhxrw    : 不過今天我舉的特例剛好是泰勒多項式, 因此我想證的
→ znmkhxrw    : 可以直接去對泰勒多項式做微分檢查得證
→ znmkhxrw    : 但是general case得到的L(x)就不知道怎麼證會符合
→ znmkhxrw    : 微分條件
→ Vulpix      : 你那句話跟二階泰勒一樣啊。
推 znmkhxrw    : 嗨V大我回了wiki的例子一篇 推文不好排版 謝啦
先修一點 cap 的錯漏。
※ 編輯: Vulpix (163.13.112.58 臺灣), 08/07/2023 04:56:47

把 divided differences 算一算，感覺很有趣。

在 f is smooth enough 的前提下，divided differences 其實都可以用極限延拓。

說的是類似 f[a,a] 這種東西可以自然定義成 f'(a)。

不過這樣一來，如果 f 是 C^1，就保證 f[x,y] 有 C^0。

而要讓 f[x,y,z] 能處處連續，f 至少也得是 C^2。


總之，f[x_0] + f[x_0,x_1](x-x_0) + f[x_0,x_1,x_2](x-x_0)^2 的極限，

就是 f[x_0] + f[x_0,x_0](x-x_0) + f[x_0,x_0,x_0](x-x_0)^2。

f[x_0,x_0] 就是 f'(x_0) 沒問題，而 f[x_0,x_0,x_0] 則是 f"(x_0)/2。

有公式為證：

    d/dx f[z_1,z_2,...,z_n,x] = f[z_1,z_2,...,z_n,x,x]

用數學歸納法甚至可以得到：

    (D_x^3/3!)(D_y^2/2!)(D_z^2/2!) f[x,y,z] = f[x,x,x,x,y,y,y,z,z,z]。

f[x_0,x_0,x_0] 的情況比較簡單，就是 f[x_0] 在 x_0 的二階導數再除以 2!。

如果是已知 f, f', f" 在 -1, 0, 1 上的值，

 x      f(x)   f'(x)  f"(x)
-1       2      -8      56
 0       1       0       0
 1       2       8      56

那這種 p(x) 的 Newton form 就是

    f[-1] + f[-1,-1](x+1) + f[-1,-1,-1](x+1)^2 + f[-1,-1,-1,0](x+1)^3

  + f[-1,-1,-1,0,0]x(x+1)^3 + f[-1,-1,-1,0,0,0]x^2(x+1)^3

  + f[-1,-1,-1,0,0,0,1]x^3(x+1)^3 + f[-1,-1,-1,0,0,0,1,1]x^3(x+1)^3(x-1)

  + f[-1,-1,-1,0,0,0,1,1,1]x^3(x+1)^3(x-1)^2

然後 wiki 上那張表，就只是很理所當然地計算這些係數而已。

所以你應該是卡在：

1. f[1,2,3,4,5] 我會，但 f[-1,-1,-1,0,0] 到底是什麼？

2. 為什麼 f[-1,b,c,0,e] 在 b,c→-1 且 e→0 的時候會收斂到 f[-1,-1,-1,0,0]？

根據前面的脈絡，兩個問題的回答是一起的：

    f[a,b,c,d,e] 在 Ｒ^5 上無法直接用 divided difference 寫下的位置，

    以其極限取代之。則 f[a,b,c,d,e] 在 Ｒ^5 上連續。

總之就是要確認 a divided difference with repeated arguments is well defined。

    f[a,a] = lim_{x→a} { f(x)-f(a) }/{ x-a } = f'(a)

而一般一點的情況，建議從 expanded form 下手。

以 f[0,0,0,1,1] 為例，

  f[0,b,c,1,e]

= f(0)/(-b)(-c)(-1)(-e) + f(b)/b(b-c)(b-1)(b-e) + f(c)/c(c-b)(c-1)(c-e)

    + f(1)/1(1-b)(1-c)(1-e) + f(e)/e(e-b)(e-c)(e-1)

後二項在 b,c→0 的時候，會趨近於 -f(1)/(e-1) + f(e)/e^3(e-1)。

然後在 e→1 的時候會再趨近於 f'(1)-3f(1)。

前三項在 e→1 的時候，會→ f(0)/bc - f(b)/b(c-b)(b-1)^2 + f(c)/c(c-b)(c-1)^2。

最後在 b,c→0 下會趨近於 3f(0) + 2f'(0) + f"(0)/2。


關於連續性，適合的參考資料應該是 https://ftp.cs.wisc.edu/Approx/deboor2.pdf。

f[x,y] 在 Ｒ^2 上連續，這直接做就好，沒有很難做。

更多變數的情況下就要一些技巧了。
※ 編輯: Vulpix (163.13.112.58 臺灣), 08/08/2023 13:45:45
推 znmkhxrw    : 謝謝V大的分享! 關於連續性我有兩個看法:
→ znmkhxrw    : (1) 我自己對於f[x,y]跟f[x,y,z]都是一直羅畢達XD
→ znmkhxrw    : 但是general case我就羅不下去了, 太醜了, 你給的
→ znmkhxrw    : pdf應該就是general解決這件事吧
→ znmkhxrw    : (2) 對於微分條件, 用mean value thm for divided
→ znmkhxrw    : difference來看的話, 要處理f[x,y,z]確實需要f€C^2
→ znmkhxrw    : 才能讓3點退化成1點的(f''(ε)取極限把極限搬入)
→ znmkhxrw    : 但是我總覺得有辦法只要"f€C^1, f'€diff"就可以
→ znmkhxrw    : 以普通MVT舉例, (f(x)-f(a))/(x-a)=f'(ε)
→ znmkhxrw    : 如果f€C^1, 當然可以x→a讓f'(ε)趨近於f'(a)
→ znmkhxrw    : 但是其實f€diff即可, 因為根本不用透過MVT
要 f[x,y] 連續至少要 C^1，因為沿著 y=x 靠近 (a,a) 的時候要 f' 連續。
※ 編輯: Vulpix (163.13.112.58 臺灣), 08/08/2023 17:31:27
推 znmkhxrw    : V大我回應如下連結, 有數學式跟說明, 謝謝
→ znmkhxrw    : https://www.desmos.com/calculator/8eh2iuqzdf
https://i.imgur.com/kDaTMXM.png

對，我知道因為 MVT 的要求是 "f conti. on [a,b], f is diff. on (a,b)."，

所以如果只是要 f[0,b,c,1,e] 會收斂到 f[0,0,0,1,1] 上，可能可以放寬條件，

甚至上面這個應該不用到四階導數存在，只要二階可導就可以了。

可是大家應該也都知道……

(partially) derivable, differentiable, continuously differentiable 很煩人。

但我的確是想先建構函數再考慮，畢竟，f[x,y,z,u,v] 那種函數很美嘛。
※ 編輯: Vulpix (163.13.112.58 臺灣), 08/09/2023 17:34:42
推 znmkhxrw    : 同意你說的, 謝謝這串分享!

找到比較一勞永逸的方式：


首先，這次既不是用 Newton form，也不是用 Lagrange form，

畢竟 p(x) = Σ_{i=1}^n (c_i x^i) 這個 standard form 還是最容易求導數的長相。

已知 f(x_i) for i = 0, 1, ... , n，而且 x_i 各不相同。

那麼，https://i.imgur.com/sPmLIUo.png。

因為 Vandermonde determinant 不是 0，所以這組 c_i 有唯一解。


下一步是讓 x_1, x_2, ... , x_{m_0} 都趨近於 x_0，

不過在那之前，要先處理一下前 1 + m_0 列。

把上圖的等式拿來做以下列運算：

for(int i=1; i<=1+m_0; i++)

    for(int j=1+m_0; j>=i; j--)

        (第 j 列 -= 第 j-1 列) /= (x_j - x_{j-i});

總之，經過這串列運算以後，等式會變成 https://i.imgur.com/7UsrbCz.png。

        其中的「１」代表 1 函數，各個「ｘ^k」則各自代表 k 次方函數。

        其實１[x_0,x_1] = 0，上部矩陣的下三角都是 0。

        然後ｘ[x_0,x_1] = 1，上部矩陣的主對角線上都是 1。

        下部矩陣沒有動到，照抄。

接下來要確認一下他的行列式值。

雖然長相很兇惡，但是因為我們之前有紀錄列運算的過程，

所以實際上是 Vand. det. / sub-Vand. det. of {x_0, ... , x_{m_0}}，

所以這個行列式 = Π_{j>i>m_0} (x_j - x_i) * Π_{j>m_0≧i} (x_j - x_i)，

在 x_1, x_2, ... , x_{m_0} 都趨近於 x_0 的時候，

收斂到 Π_{j>i>m_0} (x_j - x_i) * Π_{j>m_0} (x_j - x_0)^{1+m_0} ≠ 0。

這表示如果那個方陣的極限存在的話，行列式值非零。

終於要算極限了，前面那個方程式左側的方陣和右側的行矩陣各自都是收斂的，

而且方陣極限的行列式值非零，所以 c 那一個行矩陣也收斂。

整個方程式的極限是 https://i.imgur.com/AoSFHbG.png。

        跟剛剛一樣，其實上部矩陣的下三角都是 0，而且主對角線上都是 1。

        只是為了能有個通式的長相就拉他們下水。

        下部矩陣沒有動到，照抄。

然後反覆把想拿來簡併的 x_k 併在一起，我們就得到了退化多項式 p(x) 的係數 c_i。

前 1 + m_0 列已經不會再被動到了。

p(x_0) = f(x_0) 可直接參照矩陣方程式的第一列。

觀察第二列可得 p'(x_0) = 1 + 2x_0 + 3x_0^2 + ... + nx_0^{n-1} = f'(x_0)，

同理可得其他高階導數的等式。


這個作法就是從最初的插值多項式直接退化成 p(x)，

並且證明了直到第「簡併數」階之前，p 在簡併點上的導數 = f 在簡併點上的導數。


-------------------------
至於 f[a,b,c] 收斂到 f"(a)/2!，似乎真的可以用 f" 存在來證。

總之先對 b,c 做 MVT：

如果 [ f'(ξ)(ξ-a)-f(ξ)+f(a) ]/(ξ-a)^2 收斂到 f"(a)/2，

那 f[a,b,c] 就收斂到 f"(a)/2!。

可是羅下去會卡在 f" 的連續性上，所以雖然我愛羅但是不能羅，

因為羅後不收斂不代表羅前也不收斂。

  [ f'(ξ)(ξ-a)-f(ξ)+f(a) ]/(ξ-a)^2

= [ f'(ξ)-f'(a) ]/(ξ-a) - [ f(ξ)-f(a)-f'(a)(ξ-a) ]/(ξ-a)^2

→　f"(a) - f"(a)/2 = f"(a)/2

前項直接算極限，後項則是先羅一次。

前面做 MVT 也很辛苦，f[a,x] 的連續性是顯然的，

而 f[a,x] 的可微性則要考慮是否在 a 微分。

不在 a 的時候，很簡單，商法則可搞定。

在 a 的話，[ f[a,x]-f'(a) ]/(x-a) 把分母處理好以後就是先羅一次，

其實跟上面那個先羅一次的後項是一樣的東西，所以也是收斂到 f"(a)/2。

f[a,b,c] = d/dx f[a,x] |x=ξ for some ξ between b and c，

這個 ξ 有可能是 a，

所以 f[a,b,c] 可能是 [ f'(ξ)(ξ-a)-f(ξ)+f(a) ]/(ξ-a)^2 或 f"(a)/2。

最後當 b,c→a 的時候，也讓 ξ→a，然後就回到前頭的計算了。

可是再多一個 d 的時候……

感覺上應該是可以做 MIT 的，但我覺得累。

f[a,b,c,d] = (f[a,b,c]-f[a,b,d])/(c-d) 且 c≠d，進行 MVT：

1. 檢查 f[a,b,x] 的連續性。

   雖然 b≠a，但是 x 可能是兩者之一。

   所以對於 x 的連續性是依賴於 f \in C^1。

   f 三階可微，所以 f 當然 \in C^1。

2. 檢查 f[a,b,x] 的可微性。

   當 x 不是 a 也不是 b 的時候，直接用公式算導函數。

   然後在 a 上計算 (f[a,b,x]-f[a,b,a])/(x-a) 的極限，

   也就是 f[a,a,x,b] 的極限，這個應該可以由 C^2 保證。

   （f[b,b,x,a] 在 b 的情形是類似的，因為此時 a,b 有對稱性。）

到此，應用 MVT 得 f[a,b,c,d] = d/dx f[a,b,x] |x=ξ

                             = lim_{x→ξ} (f[a,b,x]-f[a,b,ξ])/(x-ξ)

                             = lim_{x→ξ} f[a,b,ξ,x]

而這個 ξ 只是介於 c,d 之間而已，可能是兩數之間的任何數。

接下來分成 ξ 剛好是 a：此時 f[a,b,c,d] = f[a,b,a,a] （∵C^2）

           ξ 剛好是 b：此時 f[a,b,c,d] = f[a,b,b,b] （∵C^2）

           ξ 兩者皆非：此時 f[a,b,c,d] = f[a,b,ξ,ξ] （∵C^1）

所以綜合來說，f[a,b,c,d] = f[a,b,ξ,ξ] 都是對的。

然後計算 b,ξ→a 的時候 f[a,b,c,d] 的極限，

基於 b,ξ 的異同情況，f[a,b,c,d]=f[a,b,b,b] 或 f[a,η,η,η]，

統一寫作 f[a,b,c,d] = f[a,η,η,η]，其中 η 介於 b,ξ 之間或者就是 b，

最後利用三階可微能算出他收斂到 f[a,a,a,a] = f"'(a)/3!。

感覺上應該還是要善用 divided difference 的符號，

不過扯上簡併的連續性跟可微性都得驗證。

或者可以針對 b,c,d 直接做某種推廣版的 MVT：

    f[a,b,c,d] = f[a,η,η,η], where min(b,c,d)<η<max(b,c,d).
※ 編輯: Vulpix (1.160.12.97 臺灣), 08/13/2023 07:26:34
推 znmkhxrw    : 嗨V大辛苦了~好多資訊, 幾個問題跟回應:
→ znmkhxrw    : (1) "「１」代表 1 函數，各個「ｘ^k」則各自代表
→ znmkhxrw    : k 次方函數。"這句話看不太懂, 所以也不知道為什麼
→ znmkhxrw    : １[x_0,x_1] = 0以及ｘ[x_0,x_1] = 1
f[x_0,x_1] = [ f(x_1)-f(x_0) ]/(x_1-x_0)

f 代入 1 函數，得 (1-1)/(x_1-x_0) = 0。

f 代入 k 次方函數，得

        (x_1^k-x_0^k)/(x_1-x_0) = Σ_{i=0}^{k-1} x_0^i x_1^{k-1-i}。

例如 k = 2 的 2 次方函數，

ｘ^2[x_0,x_1] = (x_1^2-x_0^2)/(x_1-x_0) = x_1+x_0，

ｘ^2[x_1,x_2] = (x_2^2-x_1^2)/(x_2-x_1) = x_2+x_1，

ｘ^2[x_0,x_1,x_2] = [ (x_2+x_1) - (x_1+x_0) ]/(x_2-x_0) = 1 這樣。

所以其實（上部矩陣的）下三角都是 0、（上部矩陣的）主對角線都是 1。

→ znmkhxrw    : (2) 這個方程式 https://i.imgur.com/AoSFHbG.png
→ znmkhxrw    : 的x_0有負次方是要假設x_0不為零?
那就把次方都改成 max(0, 原本的次方) 好了……

或者乾脆改成原次方的絕對值。

不過反正負次方項的係數都是 0 啦，當成 notation 看就好，沒有什麼代值的功能。

→ znmkhxrw    : (3) 承上的方程式, 跟Confluent Vandermonde很像耶!
→ znmkhxrw    : V大幫我修連結XD https://i.imgur.com/rk7ijn2.jpg
→ znmkhxrw    : 以上連結的唯一解c_0~c_n所決定出來的多項式p(x)
→ znmkhxrw    : 如果可以證明出來他就是Newton/Lagrange的退化函數
首先，不管 Newton、Lagrange 插出來的多項式都跟

https://i.imgur.com/sPmLIUo.png 解出來的多項式是同一個。

所以雖然我算的是這個多項式的退化，但可以直接當成 Newton form 的極限，沒問題。

你把這個矩陣方程式展開，其實每一條都只是 p(x_i) = f(x_i) 而已。

→ znmkhxrw    : 或是證明他是difference table演算法的那個函數
→ znmkhxrw    : 那就解決了. 只是我沒有勇氣去寫出general case的
→ znmkhxrw    : Confluent Vandermonde的反矩陣然後再乘開...
→ znmkhxrw    : 然後剛剛看到V大你連結的矩陣, 發現跟CV很相似欸
→ znmkhxrw    : (簡稱CV為Confluent Vandermonde)
→ znmkhxrw    : 只是那個Ac=f的A矩陣跟f向量不太一樣, V大連結的f
→ znmkhxrw    : 的微分項帶有factorial
因為我想要都用 divided difference 寫，所以會比 CV 多除以一些階乘。

也是因為這樣，我的係數才是組合數 C 而不是 CV 的排列數 P。

→ znmkhxrw    : 我猜"(第 j 列 -= 第 j-1 列) /= (x_j - x_{j-i});"
→ znmkhxrw    : 這個操作如果加料一下應該就可以推到rk7ijn2.jpg了
對，但我不要。這樣的行列式裡面會有一大堆階乘。

我的證明過程須要算行列式，不想讓他變醜。

→ znmkhxrw    : (4) 最後有關微分條件放寬的思路謝謝分享
→ znmkhxrw    : 其中有"MIT"是指哪個定理呢?
數學歸納法。

不過我還在想到底有哪些東西要定義跟證明XD

以這個流程來說，至少就是 f[x_1, ..., x_n] 收斂到 f^{(n)}(x_0)/n! 吧。

→ PPguest     : @z大 difference table演算法不就是在算Newton form
→ PPguest     : 的退化函數
是啊。
※ 編輯: Vulpix (1.160.12.97 臺灣), 08/13/2023 22:58:39
推 znmkhxrw    : 08:54的回應理解了 謝謝
→ znmkhxrw    : 09:56 是的! 所以我自己當初研究退化問題就是在想說
→ znmkhxrw    : 那樣的表示法最好研究退化函數的長相, 然後Lagrange
→ znmkhxrw    : 跟矩陣就被我捨棄了, 當下下意識覺得矩陣最麻煩
→ znmkhxrw    : 因為還要寫出Vand反矩陣的通式還要取極限, 然後再跟
→ znmkhxrw    : f向量相乘, 覺得很費工, 想不到V大這樣整理蠻乾淨的
既然反矩陣麻煩，那為什麼不用克拉瑪公式呢？

→ znmkhxrw    : 09:56->09:06 打錯
→ znmkhxrw    : 08:56, 09:11, 09:16一起回: 原來V大是把x_0^(-N)
→ znmkhxrw    : 當作0的notation, 這樣確實就跟Confluent Vand差一
→ znmkhxrw    : 個左乘的對角矩陣D了, D的左上區是階層, 右下區是
→ znmkhxrw    : 單位矩陣
作者: Vulpix (Sebastian) 看板: Math
標題: Re: [分析] Hermite內插演算法的證明
時間: Sun Aug  6 19:28:57 2023

符號看得有點花……


如果你想做的是「在 x_1 和 x_2 分別趨近 x_0 後所得的極限 = Taylor 多項式」，

那你需要的就是 MVT of divided differences。

https://en.wikipedia.org/wiki/Mean_value_theorem_(divided_differences)

直接套上去就馬上做完了。

也不必去算新多項式的導數。

但是如果要一步一步來就沒那麼好算了。

  (x_n - x_0)lim_{x_1→x_0} f[x_0,...,x_n]

= f[x_0,x_2,...,x_n] - lim_{x_1→x_0} f[x_0,...,x_{n-1}]

上面這條遞迴式是用來算極限的。

本來的插值多項式是 f[x_0] + f[x_0,x_1](x-x_0) + f[x_0,x_1,x_2](x-x_0)^2。

在 x_1→x_0 之後，變成

    f(x_0) + f'(x_0)(x-x_0) + {f[x_0,x_2]-f'(x_0)}/(x_2-x_0) * (x-x_0)^2。

然後 {f[x_0,x_2] - f'(x_0)}/(x_2 - x_0) 在 x_2→x_0 下的極限 = f"(x_0)/2，

所以多項式的極限就變成 f(x_0) + f'(x_0)(x-x_0) + f"(x_0)/2 * (x-x_0)^2。


不過我本來在想的是用 Lagrange 觀點。

e_0(x) := Π_{i=1}^n (x-x_i)/(x_0-x_i)，其他 e_j 類推。

還是先用 n = 2 來觀察，

插值多項式 = f(x_0)e_0(x) + f(x_1)e_1(x) + f(x_2)e_2(x)

然後也先讓 x_1→x_0，多項式變成

    f(x_0){e_0(x)+e_1(x)} + {f(x_1)-f(x_0)}e_1(x) + f(x_2)e_2(x)。

所以我們分成三項來看：

1. e_0(x)+e_1(x) = (x-x_1)(x-x_2)/(x_0-x_1)(x_0-x_2)

                      + (x-x_0)(x-x_2)/(x_1-x_0)(x_1-x_2)    公因式

                 = (x-x_2)/(x_1-x_0)

                      * (x_1-x_0)(x_0-x_2+x_1-x)/(x_0-x_2)(x_1-x_2)

                 = (x-x_2)(x_0-x_2+x_1-x)/(x_0-x_2)(x_1-x_2)

                 → (x-x_2)(2x_0-x_2-x)/(x_0-x_2)^2

                 = 1 - (x-x_0)^2/(x_0-x_2)^2

   最後這個多項式，他代 x_0 得 1、導數得 0，而代 x_2 得 0。

2. {f(x_1)-f(x_0)}e_1(x) = {f(x_1)-f(x_0)}(x-x_0)(x-x_2)/(x_1-x_0)(x_1-x_2)

                         → f'(x_0)(x-x_0)(x-x_2)/(x_0-x_2)

   (x-x_0)(x-x_2)/(x_0-x_2) 代 x_0 得 0、導數得 1，而代 x_2 得 0。

3. e_2(x) = (x-x_0)(x-x_1)/(x_2-x_0)(x_2-x_1) → (x-x_0)^2/(x_2-x_0)^2

   最後這個多項式也是代 x_0 得 0、導數得 0，而代 x_2 得 1。

我們得到三個可以各自突顯 f(x_0), f'(x_0), f(x_2) 的多項式，

剛好跟 Lagrange 觀點有謀而合。

最後再讓 x_2→x_0，

  f(x_0){1-(x-x_0)^2/(x_0-x_2)^2}

      + f'(x_0)(x-x_0)(x-x_2)/(x_0-x_2) + f(x_2)(x-x_0)^2/(x_2-x_0)^2

= f(x_0)+f'(x_0)(x-x_0)

      + { f(x_2) - f(x_0) - f'(x_0)(x_2-x_0) }(x-x_0)^2/(x_2-x_0)^2

→ f(x_0) + f'(x_0)(x-x_0) + f"(x_0)(x-x_0)^2/2

其實仔細看，1, x-x_0, (x-x_0)^2/2 也是

在函數值、一階導數、二階導數之中各自突顯一項，而消滅其他兩項的多項式函數，

同樣符合 Lagrange 觀點的插值概念。


真正麻煩的還是 general case：

    有資料的點是 x_0, ..., x_n，每個點的高階導數已知階數不盡相同。

像是已知 f(-1), f(0), f'(0), f"(0), f(5), f(100), f'(100) 這樣。

然後先用 -1, 0, a, b, 5, 100, c 插值，再讓 a,b→0 和 c→100，

之後要檢查在 x = 0 的一階二階導數和在 x = 100 的一階導數。

不過我想，應該也是這樣一步步算極限就好。

但是那個 general form 就真的很難看，所以平常都是給 algorithm。

--
※ 發信站: 批踢踢實業坊(ptt.cc), 來自: 1.162.224.247 (臺灣)
※ 文章網址: https://www.ptt.cc/bbs/Math/M.1691321340.A.CC3.html
推 znmkhxrw    : 嗨V大, 我想證的是「在 x_1 和 x_2 分別趨近 x_0 後
→ znmkhxrw    : 所得的極限L(x)」會滿足L(x_0)=f(x_0),
→ znmkhxrw    : L'(x_0)=f'(x_0), L''(x_0)=f''(x_0)
→ znmkhxrw    : 不過今天我舉的特例剛好是泰勒多項式, 因此我想證的
→ znmkhxrw    : 可以直接去對泰勒多項式做微分檢查得證
→ znmkhxrw    : 但是general case得到的L(x)就不知道怎麼證會符合
→ znmkhxrw    : 微分條件
→ Vulpix      : 你那句話跟二階泰勒一樣啊。
推 znmkhxrw    : 嗨V大我回了wiki的例子一篇 推文不好排版 謝啦
先修一點 cap 的錯漏。
※ 編輯: Vulpix (163.13.112.58 臺灣), 08/07/2023 04:56:47

把 divided differences 算一算，感覺很有趣。

在 f is smooth enough 的前提下，divided differences 其實都可以用極限延拓。

說的是類似 f[a,a] 這種東西可以自然定義成 f'(a)。

不過這樣一來，如果 f 是 C^1，就保證 f[x,y] 有 C^0。

而要讓 f[x,y,z] 能處處連續，f 至少也得是 C^2。


總之，f[x_0] + f[x_0,x_1](x-x_0) + f[x_0,x_1,x_2](x-x_0)^2 的極限，

就是 f[x_0] + f[x_0,x_0](x-x_0) + f[x_0,x_0,x_0](x-x_0)^2。

f[x_0,x_0] 就是 f'(x_0) 沒問題，而 f[x_0,x_0,x_0] 則是 f"(x_0)/2。

有公式為證：

    d/dx f[z_1,z_2,...,z_n,x] = f[z_1,z_2,...,z_n,x,x]

用數學歸納法甚至可以得到：

    (D_x^3/3!)(D_y^2/2!)(D_z^2/2!) f[x,y,z] = f[x,x,x,x,y,y,y,z,z,z]。

f[x_0,x_0,x_0] 的情況比較簡單，就是 f[x_0] 在 x_0 的二階導數再除以 2!。

如果是已知 f, f', f" 在 -1, 0, 1 上的值，

 x      f(x)   f'(x)  f"(x)
-1       2      -8      56
 0       1       0       0
 1       2       8      56

那這種 p(x) 的 Newton form 就是

    f[-1] + f[-1,-1](x+1) + f[-1,-1,-1](x+1)^2 + f[-1,-1,-1,0](x+1)^3

  + f[-1,-1,-1,0,0]x(x+1)^3 + f[-1,-1,-1,0,0,0]x^2(x+1)^3

  + f[-1,-1,-1,0,0,0,1]x^3(x+1)^3 + f[-1,-1,-1,0,0,0,1,1]x^3(x+1)^3(x-1)

  + f[-1,-1,-1,0,0,0,1,1,1]x^3(x+1)^3(x-1)^2

然後 wiki 上那張表，就只是很理所當然地計算這些係數而已。

所以你應該是卡在：

1. f[1,2,3,4,5] 我會，但 f[-1,-1,-1,0,0] 到底是什麼？

2. 為什麼 f[-1,b,c,0,e] 在 b,c→-1 且 e→0 的時候會收斂到 f[-1,-1,-1,0,0]？

根據前面的脈絡，兩個問題的回答是一起的：

    f[a,b,c,d,e] 在 Ｒ^5 上無法直接用 divided difference 寫下的位置，

    以其極限取代之。則 f[a,b,c,d,e] 在 Ｒ^5 上連續。

總之就是要確認 a divided difference with repeated arguments is well defined。

    f[a,a] = lim_{x→a} { f(x)-f(a) }/{ x-a } = f'(a)

而一般一點的情況，建議從 expanded form 下手。

以 f[0,0,0,1,1] 為例，

  f[0,b,c,1,e]

= f(0)/(-b)(-c)(-1)(-e) + f(b)/b(b-c)(b-1)(b-e) + f(c)/c(c-b)(c-1)(c-e)

    + f(1)/1(1-b)(1-c)(1-e) + f(e)/e(e-b)(e-c)(e-1)

後二項在 b,c→0 的時候，會趨近於 -f(1)/(e-1) + f(e)/e^3(e-1)。

然後在 e→1 的時候會再趨近於 f'(1)-3f(1)。

前三項在 e→1 的時候，會→ f(0)/bc - f(b)/b(c-b)(b-1)^2 + f(c)/c(c-b)(c-1)^2。

最後在 b,c→0 下會趨近於 3f(0) + 2f'(0) + f"(0)/2。


關於連續性，適合的參考資料應該是 https://ftp.cs.wisc.edu/Approx/deboor2.pdf。

f[x,y] 在 Ｒ^2 上連續，這直接做就好，沒有很難做。

更多變數的情況下就要一些技巧了。
※ 編輯: Vulpix (163.13.112.58 臺灣), 08/08/2023 13:45:45
推 znmkhxrw    : 謝謝V大的分享! 關於連續性我有兩個看法:
→ znmkhxrw    : (1) 我自己對於f[x,y]跟f[x,y,z]都是一直羅畢達XD
→ znmkhxrw    : 但是general case我就羅不下去了, 太醜了, 你給的
→ znmkhxrw    : pdf應該就是general解決這件事吧
→ znmkhxrw    : (2) 對於微分條件, 用mean value thm for divided
→ znmkhxrw    : difference來看的話, 要處理f[x,y,z]確實需要f€C^2
→ znmkhxrw    : 才能讓3點退化成1點的(f''(ε)取極限把極限搬入)
→ znmkhxrw    : 但是我總覺得有辦法只要"f€C^1, f'€diff"就可以
→ znmkhxrw    : 以普通MVT舉例, (f(x)-f(a))/(x-a)=f'(ε)
→ znmkhxrw    : 如果f€C^1, 當然可以x→a讓f'(ε)趨近於f'(a)
→ znmkhxrw    : 但是其實f€diff即可, 因為根本不用透過MVT
要 f[x,y] 連續至少要 C^1，因為沿著 y=x 靠近 (a,a) 的時候要 f' 連續。
※ 編輯: Vulpix (163.13.112.58 臺灣), 08/08/2023 17:31:27
推 znmkhxrw    : V大我回應如下連結, 有數學式跟說明, 謝謝
→ znmkhxrw    : https://www.desmos.com/calculator/8eh2iuqzdf
https://i.imgur.com/kDaTMXM.png

對，我知道因為 MVT 的要求是 "f conti. on [a,b], f is diff. on (a,b)."，

所以如果只是要 f[0,b,c,1,e] 會收斂到 f[0,0,0,1,1] 上，可能可以放寬條件，

甚至上面這個應該不用到四階導數存在，只要二階可導就可以了。

可是大家應該也都知道……

(partially) derivable, differentiable, continuously differentiable 很煩人。

但我的確是想先建構函數再考慮，畢竟，f[x,y,z,u,v] 那種函數很美嘛。
※ 編輯: Vulpix (163.13.112.58 臺灣), 08/09/2023 17:34:42
推 znmkhxrw    : 同意你說的, 謝謝這串分享!

找到比較一勞永逸的方式：


首先，這次既不是用 Newton form，也不是用 Lagrange form，

畢竟 p(x) = Σ_{i=1}^n (c_i x^i) 這個 standard form 還是最容易求導數的長相。

已知 f(x_i) for i = 0, 1, ... , n，而且 x_i 各不相同。

那麼，https://i.imgur.com/sPmLIUo.png。

因為 Vandermonde determinant 不是 0，所以這組 c_i 有唯一解。


下一步是讓 x_1, x_2, ... , x_{m_0} 都趨近於 x_0，

不過在那之前，要先處理一下前 1 + m_0 列。

把上圖的等式拿來做以下列運算：

for(int i=1; i<=1+m_0; i++)

    for(int j=1+m_0; j>=i; j--)

        (第 j 列 -= 第 j-1 列) /= (x_j - x_{j-i});

總之，經過這串列運算以後，等式會變成 https://i.imgur.com/7UsrbCz.png。

        其中的「１」代表 1 函數，各個「ｘ^k」則各自代表 k 次方函數。

        其實１[x_0,x_1] = 0，上部矩陣的下三角都是 0。

        然後ｘ[x_0,x_1] = 1，上部矩陣的主對角線上都是 1。

        下部矩陣沒有動到，照抄。

接下來要確認一下他的行列式值。

雖然長相很兇惡，但是因為我們之前有紀錄列運算的過程，

所以實際上是 Vand. det. / sub-Vand. det. of {x_0, ... , x_{m_0}}，

所以這個行列式 = Π_{j>i>m_0} (x_j - x_i) * Π_{j>m_0≧i} (x_j - x_i)，

在 x_1, x_2, ... , x_{m_0} 都趨近於 x_0 的時候，

收斂到 Π_{j>i>m_0} (x_j - x_i) * Π_{j>m_0} (x_j - x_0)^{1+m_0} ≠ 0。

這表示如果那個方陣的極限存在的話，行列式值非零。

終於要算極限了，前面那個方程式左側的方陣和右側的行矩陣各自都是收斂的，

而且方陣極限的行列式值非零，所以 c 那一個行矩陣也收斂。

整個方程式的極限是 https://i.imgur.com/AoSFHbG.png。

        跟剛剛一樣，其實上部矩陣的下三角都是 0，而且主對角線上都是 1。

        只是為了能有個通式的長相就拉他們下水。

        下部矩陣沒有動到，照抄。

然後反覆把想拿來簡併的 x_k 併在一起，我們就得到了退化多項式 p(x) 的係數 c_i。

前 1 + m_0 列已經不會再被動到了。

p(x_0) = f(x_0) 可直接參照矩陣方程式的第一列。

觀察第二列可得 p'(x_0) = 1 + 2x_0 + 3x_0^2 + ... + nx_0^{n-1} = f'(x_0)，

同理可得其他高階導數的等式。


這個作法就是從最初的插值多項式直接退化成 p(x)，

並且證明了直到第「簡併數」階之前，p 在簡併點上的導數 = f 在簡併點上的導數。


-------------------------
至於 f[a,b,c] 收斂到 f"(a)/2!，似乎真的可以用 f" 存在來證。

總之先對 b,c 做 MVT：

如果 [ f'(ξ)(ξ-a)-f(ξ)+f(a) ]/(ξ-a)^2 收斂到 f"(a)/2，

那 f[a,b,c] 就收斂到 f"(a)/2!。

可是羅下去會卡在 f" 的連續性上，所以雖然我愛羅但是不能羅，

因為羅後不收斂不代表羅前也不收斂。

  [ f'(ξ)(ξ-a)-f(ξ)+f(a) ]/(ξ-a)^2

= [ f'(ξ)-f'(a) ]/(ξ-a) - [ f(ξ)-f(a)-f'(a)(ξ-a) ]/(ξ-a)^2

→　f"(a) - f"(a)/2 = f"(a)/2

前項直接算極限，後項則是先羅一次。

前面做 MVT 也很辛苦，f[a,x] 的連續性是顯然的，

而 f[a,x] 的可微性則要考慮是否在 a 微分。

不在 a 的時候，很簡單，商法則可搞定。

在 a 的話，[ f[a,x]-f'(a) ]/(x-a) 把分母處理好以後就是先羅一次，

其實跟上面那個先羅一次的後項是一樣的東西，所以也是收斂到 f"(a)/2。

f[a,b,c] = d/dx f[a,x] |x=ξ for some ξ between b and c，

這個 ξ 有可能是 a，

所以 f[a,b,c] 可能是 [ f'(ξ)(ξ-a)-f(ξ)+f(a) ]/(ξ-a)^2 或 f"(a)/2。

最後當 b,c→a 的時候，也讓 ξ→a，然後就回到前頭的計算了。

可是再多一個 d 的時候……

感覺上應該是可以做 MIT 的，但我覺得累。

f[a,b,c,d] = (f[a,b,c]-f[a,b,d])/(c-d) 且 c≠d，進行 MVT：

1. 檢查 f[a,b,x] 的連續性。

   雖然 b≠a，但是 x 可能是兩者之一。

   所以對於 x 的連續性是依賴於 f \in C^1。

   f 三階可微，所以 f 當然 \in C^1。

2. 檢查 f[a,b,x] 的可微性。

   當 x 不是 a 也不是 b 的時候，直接用公式算導函數。

   然後在 a 上計算 (f[a,b,x]-f[a,b,a])/(x-a) 的極限，

   也就是 f[a,a,x,b] 的極限，這個應該可以由 C^2 保證。

   （f[b,b,x,a] 在 b 的情形是類似的，因為此時 a,b 有對稱性。）

到此，應用 MVT 得 f[a,b,c,d] = d/dx f[a,b,x] |x=ξ

                             = lim_{x→ξ} (f[a,b,x]-f[a,b,ξ])/(x-ξ)

                             = lim_{x→ξ} f[a,b,ξ,x]

而這個 ξ 只是介於 c,d 之間而已，可能是兩數之間的任何數。

接下來分成 ξ 剛好是 a：此時 f[a,b,c,d] = f[a,b,a,a] （∵C^2）

           ξ 剛好是 b：此時 f[a,b,c,d] = f[a,b,b,b] （∵C^2）

           ξ 兩者皆非：此時 f[a,b,c,d] = f[a,b,ξ,ξ] （∵C^1）

所以綜合來說，f[a,b,c,d] = f[a,b,ξ,ξ] 都是對的。

然後計算 b,ξ→a 的時候 f[a,b,c,d] 的極限，

基於 b,ξ 的異同情況，f[a,b,c,d]=f[a,b,b,b] 或 f[a,η,η,η]，

統一寫作 f[a,b,c,d] = f[a,η,η,η]，其中 η 介於 b,ξ 之間或者就是 b，

最後利用三階可微能算出他收斂到 f[a,a,a,a] = f"'(a)/3!。

感覺上應該還是要善用 divided difference 的符號，

不過扯上簡併的連續性跟可微性都得驗證。

或者可以針對 b,c,d 直接做某種推廣版的 MVT：

    f[a,b,c,d] = f[a,η,η,η], where min(b,c,d)<η<max(b,c,d).
※ 編輯: Vulpix (1.160.12.97 臺灣), 08/13/2023 07:26:34
推 znmkhxrw    : 嗨V大辛苦了~好多資訊, 幾個問題跟回應:
→ znmkhxrw    : (1) "「１」代表 1 函數，各個「ｘ^k」則各自代表
→ znmkhxrw    : k 次方函數。"這句話看不太懂, 所以也不知道為什麼
→ znmkhxrw    : １[x_0,x_1] = 0以及ｘ[x_0,x_1] = 1
f[x_0,x_1] = [ f(x_1)-f(x_0) ]/(x_1-x_0)

f 代入 1 函數，得 (1-1)/(x_1-x_0) = 0。

f 代入 k 次方函數，得

        (x_1^k-x_0^k)/(x_1-x_0) = Σ_{i=0}^{k-1} x_0^i x_1^{k-1-i}。

例如 k = 2 的 2 次方函數，

ｘ^2[x_0,x_1] = (x_1^2-x_0^2)/(x_1-x_0) = x_1+x_0，

ｘ^2[x_1,x_2] = (x_2^2-x_1^2)/(x_2-x_1) = x_2+x_1，

ｘ^2[x_0,x_1,x_2] = [ (x_2+x_1) - (x_1+x_0) ]/(x_2-x_0) = 1 這樣。

所以其實（上部矩陣的）下三角都是 0、（上部矩陣的）主對角線都是 1。

→ znmkhxrw    : (2) 這個方程式 https://i.imgur.com/AoSFHbG.png
→ znmkhxrw    : 的x_0有負次方是要假設x_0不為零?
那就把次方都改成 max(0, 原本的次方) 好了……

或者乾脆改成原次方的絕對值。

不過反正負次方項的係數都是 0 啦，當成 notation 看就好，沒有什麼代值的功能。

→ znmkhxrw    : (3) 承上的方程式, 跟Confluent Vandermonde很像耶!
→ znmkhxrw    : V大幫我修連結XD https://i.imgur.com/rk7ijn2.jpg
→ znmkhxrw    : 以上連結的唯一解c_0~c_n所決定出來的多項式p(x)
→ znmkhxrw    : 如果可以證明出來他就是Newton/Lagrange的退化函數
首先，不管 Newton、Lagrange 插出來的多項式都跟

https://i.imgur.com/sPmLIUo.png 解出來的多項式是同一個。

所以雖然我算的是這個多項式的退化，但可以直接當成 Newton form 的極限，沒問題。

你把這個矩陣方程式展開，其實每一條都只是 p(x_i) = f(x_i) 而已。

→ znmkhxrw    : 或是證明他是difference table演算法的那個函數
→ znmkhxrw    : 那就解決了. 只是我沒有勇氣去寫出general case的
→ znmkhxrw    : Confluent Vandermonde的反矩陣然後再乘開...
→ znmkhxrw    : 然後剛剛看到V大你連結的矩陣, 發現跟CV很相似欸
→ znmkhxrw    : (簡稱CV為Confluent Vandermonde)
→ znmkhxrw    : 只是那個Ac=f的A矩陣跟f向量不太一樣, V大連結的f
→ znmkhxrw    : 的微分項帶有factorial
因為我想要都用 divided difference 寫，所以會比 CV 多除以一些階乘。

也是因為這樣，我的係數才是組合數 C 而不是 CV 的排列數 P。

→ znmkhxrw    : 我猜"(第 j 列 -= 第 j-1 列) /= (x_j - x_{j-i});"
→ znmkhxrw    : 這個操作如果加料一下應該就可以推到rk7ijn2.jpg了
對，但我不要。這樣的行列式裡面會有一大堆階乘。

我的證明過程須要算行列式，不想讓他變醜。

→ znmkhxrw    : (4) 最後有關微分條件放寬的思路謝謝分享
→ znmkhxrw    : 其中有"MIT"是指哪個定理呢?
數學歸納法。

不過我還在想到底有哪些東西要定義跟證明XD

以這個流程來說，至少就是 f[x_1, ..., x_n] 收斂到 f^{(n)}(x_0)/n! 吧。

→ PPguest     : @z大 difference table演算法不就是在算Newton form
→ PPguest     : 的退化函數
是啊。
※ 編輯: Vulpix (1.160.12.97 臺灣), 08/13/2023 22:58:39
推 znmkhxrw    : 08:54的回應理解了 謝謝
→ znmkhxrw    : 09:56 是的! 所以我自己當初研究退化問題就是在想說
→ znmkhxrw    : 那樣的表示法最好研究退化函數的長相, 然後Lagrange
→ znmkhxrw    : 跟矩陣就被我捨棄了, 當下下意識覺得矩陣最麻煩
→ znmkhxrw    : 因為還要寫出Vand反矩陣的通式還要取極限, 然後再跟
→ znmkhxrw    : f向量相乘, 覺得很費工, 想不到V大這樣整理蠻乾淨的
既然反矩陣麻煩，那為什麼不用克拉瑪公式呢？

→ znmkhxrw    : 09:56->09:06 打錯
→ znmkhxrw    : 08:56, 09:11, 09:16一起回: 原來V大是把x_0^(-N)
→ znmkhxrw    : 當作0的notation, 這樣確實就跟Confluent Vand差一
→ znmkhxrw    : 個左乘的對角矩陣D了, D的左上區是階層, 右下區是
→ znmkhxrw    : 單位矩陣
你誤會了，x_0^{-N} 之類的只是為了讓通式漂亮硬塞的而已。

→ znmkhxrw    : 不過我看到V大你對於組合數C有寫了C(n,m), n<m這種
→ znmkhxrw    : 形式, 比如C(1,2), 這種有定義參考資料嗎? 雖然跟他
→ znmkhxrw    : 相乘的都是x_0^(-N), 是都可以當他是0, 只是第一次
→ znmkhxrw    : 看到這種記號詢問一下
「C(n,m), n<m」才是 0。

組合解釋：從少數物品中，任取多數物品，方法數 = 0。

或者也可以把 (n-m)! 想像成 Γ(n-m+1) 這些發散到無限大的東西。

→ znmkhxrw    : 09:18 了解~多個點的f[a,b,c,d]這種要放弱微分條件
→ znmkhxrw    : 看起來要細寫組合就很多...用MIT也要找到一個好的
→ znmkhxrw    : induction方式來處理所有case, 例如f[a,b,c]就包含
→ znmkhxrw    : f[a,a,c], f[a,c,c], f[a,b,a], f[a,a,a]這幾種
→ znmkhxrw    : 廣義的f[x_0~x_n]就更多種了...
所以才麻煩，因為一直有長相不合群的在搗蛋。

不過畢竟要假設 3 階可微了，那 f[x,y] 就是可微分的函數應該沒有錯，

至少 f[x,y] 在所有位置上的兩個偏導數都存在。

那要從 f[a,a] 做出 f[a,a,c] 或者從 f[a,b] 做出 f[a,b,a] 都不是難事。

Ex. 1) f[a,a,c]

        f[a]
            〉f[a,a]
        f[a]        〉f[a,a,c]
            〉f[a,c]
        f[c]

       一開始我們的表只有白字，其中 f[a,a]=f'(a)。

       為了算 f[a,a,c] 加了一條紅字，因為 a 與 c 不同，

       所以用普通的均差遞迴關係就好。

Ex. 2) f[a,b,a]

        f[a]
            〉f[a,b]
        f[b]        〉f[a,b,a]
            〉f[b,a]
        f[a]

       這次因為 a 與 a 相同，所以要算 f[b,x] 的導數了。

       f[a,b,a] = lim_{x→a} (f[b,x]-f[a,b])/(x-a) = lim_{x→a} f[a,b,x]。

       雖然看起來很美好，但因為我們的假設很弱，所以這個極限突然下不了手。

       幸好均差是對稱函數，所以這個極限還可以不用爆開來算。

       f[a,b,a] = lim_{x→a} f[a,b,x] = lim_{x→a} f[a,x,b]

                = (f[a,b]-f[a,a])/(b-a)

                = f[a,a,b]

       回到了跟前一個例子相同的情況。

我是這樣想的：f[x,a_1,...,a_n,y] 都由那個均差表來定義。

如果 x 與 y 不同，那就是 (f[a_1,...,a_n,y]-f[x,a_1,...,a_n])/(y-x)。

但要是他們一樣，那就改成上式在 y→x 下的極限。


當然綜合來說都是 lim_{z→y} (f[a_1,...,a_n,z]-f[x,a_1,...,a_n])/(z-x)。

也因為是這樣定義的，所以至少 f[x,a_1,...,a_n,z] 對 z 是連續的。

然後靠對稱性（我懷疑這對於有簡併情形的均差是要證明一下的）移動 z 的位置，

就能知道均差對所有自變數都連續，只是要注意一點：都是視為單變數函數。

所以如果不要太貪心，一口氣讓 x_i 都趨近於 x_0，

而是一個一個來，那證明應該會輕鬆許多。


照上面的定義來看第三個例子。

Ex. 3) f[a,a,a]

        f[a,a,a] = lim_{x→a} (f[a,x]-f[a,a])/(x-a)

                 = lim_{x→a} [ (f(x)-f(a))/(x-a) - f'(a) ]/(x-a)

                 = lim_{x→a} (f(x)-f(a)-f'(a)(x-a))/(x-a)^2

                 = lim_{x→a} (f'(x)-f'(a))/2(x-a)

                 = f"(a)/2

       最後一步不能羅，要直接用導數定義。

反正怎麼搞這一步都跑不掉，索性就不要把 f[a,a] 定義成 f'(a) 了。

看著會覺得有趣的地方在分子那邊，其實會是 f(x)-泰勒多項式，

所以羅起來自然非常流暢。

→ znmkhxrw    : @P大 對, 令滿足微分條件的存在唯一多項式叫作p(x)
→ znmkhxrw    : 想證明:(1) Lagrange=Newton的退化函數就是p(x)
→ znmkhxrw    : 其中Newton型可以用不同點的差分表演算法得出
→ znmkhxrw    : (2) 重複點的差分表演算法所造的函數就是p(x)
→ znmkhxrw    : 目前來說, c大以及你當初給的google book給了第(2)
→ znmkhxrw    : 的general case證明. 而V大這裡給了(1)的special
→ znmkhxrw    : case的證明(n個不同點有m_0個點退化到同一個)
m_1 個點退化到另一個點上、m_2 個點退……

完成這個 modified CV 就可以了。

→ znmkhxrw    : 至於(1)的general case花點耐心應該可以仿造V大的
→ znmkhxrw    : 方式去寫出來
→ znmkhxrw    : 因為V大的推論結果就是讓不同點的Vand取極限(退化)
→ znmkhxrw    : 後變成Confluent Vand, 而CV可以直接跟存在唯一p(x)
→ znmkhxrw    : 連結上
→ znmkhxrw    : 因此"lim sol of V = sol of CV"就證得(1)
※ 編輯: Vulpix (163.13.112.58 臺灣), 08/14/2023 10:06:09
推 PPguest     : 剛想到可以用 Neville's algorithm 多項式的遞迴關
→ PPguest     : 係來看. https://imgur.com/arw1e82
→ PPguest     : P_1,2,...,n(x)是過那些點 x_1,x_2,...,x_n 函數值
→ PPguest     : 最高n-1次的多項式. 如果 x_i 有重複,那就是要符合
→ PPguest     : 高階導數的條件.
→ PPguest     : 在 x_i 有重複的情況,要和 Hermite 的 table 一樣先
→ PPguest     : 把重複的 x_i 合併,這樣子才容易檢查新生成的多項式
→ PPguest     : 會符合條件. 用這個看 Newton form 的 table 應會比
→ PPguest     : 較清楚.
→ PPguest     : 例如算 P_1,1,2(x),不要用 P_1,2(x)和 P_2,1(x)去算
→ PPguest     : 跟 divided difference 的情況類似要算極限. 改用
→ PPguest     : P_1,1(x)和 P_1,2(x)來算. 這樣子 general case 在
→ PPguest     : 圖中遞迴關係裡的 x_1 和 x_n 應該就會不相等
這個是直接插多項式出來，所以會在比較早的位置看到各階導數。

不過，就算先為了方便起見而把 x_i 依序排好，

他也跟均差一樣，遇到高階項就要背公式。

均差表的公式是 f[a,a,a] = f"(a)/2! 這些，

而 Neville's algorithm 的是 P_{0,2}(x) = f(a)+f'(a)(x-a)+f"(a)/2! (x-a)^2。

這都不是單純的「把前一層的結果代入遞迴公式並讓 x_1,x_2→a」。

要算這類極限，在算到這一層之前都不可以先算極限。

不過的確只要把這點先處理好，就可以比較輕鬆地檢查各點上的導數。

畢竟 Newton form 最後有「選」一條從零階均差直到 n 階均差的路徑：top row，

所以只有 x_0 那一點的導數好算，其他點的導數要找「其他路徑」去看。

所以使用 Newton form 的時候，要隨時記得「終點一樣的路就代表同一個多項式」。

也就是雖然我們只寫出一個 Newton form，

但是在把均差表完成的當下就已經有了同一個多項式的 2^n 個 Newton form。

※ 編輯: Vulpix (163.13.112.58 臺灣), 08/14/2023 16:49:24
→ PPguest     : 驗證 P_1,1,...,1(x)符合條件,會需要微分和極限互換
→ PPguest     : 註：圖片中微分的部分應該要 m<=n-2.
→ PPguest     : m個1,在微m-1次 之前 大概沒有什麼問題
→ PPguest     : 微m-1次時,用首項係數 f[x1,...,x1,x]應可以得到
→ PPguest     : f^(m-1). 但我圖中微分的式子好像怪怪的,少了階乘?
→ PPguest     : 啊沒事,階乘會在多項式首項微多次時掉出來
→ PPguest     : V大你說"不是單純的「把前一層..."是指 P_1,2(x)和
→ PPguest     : P_2,1(x)對吧,我那樣寫不太好,實際寫是用 P_1,2(x),
→ PPguest     : P_2,3(x),以及其他遞迴式裡的東西,再讓 x_3→x_1
→ PPguest     : P_2,3(x)確實沒有用到 P_1(x)
直接上實例：f(x) = x^3，選相異四數 a,b,c,d 算插值多項式。

P_aa
    〉P_ab
P_bb      〉P_ac
    〉P_bc      〉P_ad = P
P_cc      〉P_bd
    〉P_cd
P_dd

a^3
   〉(a^2+ab+b^2)x-ab(a+b)
b^3                       〉(a+b+c)x^2-(ab+bc+ca)x+abc
   〉(b^2+bc+c^2)x-bc(b+c)                            〉x^3
c^3                       〉(b+c+d)x^2-(bc+cd+db)x+bcd
   〉(c^2+cd+d^2)x-cd(c+d)
d^3

然後讓 b,c→a：

a^3
   〉   {3a^2}x-2a^3
a^3                       〉    {3a}x^2-{3a^2}x+a^3
   〉   {3a^2}x-2a^3                                   〉x^3
a^3                       〉(2a+d)x^2-(a^2+2ad)x+{a^2}d
   〉(a^2+ad+d^2)x-ad(a+d)
d^3

其中   {3a}x^2-{3a^2}x+a^3

    ≠ lim_{t→a} [(x-a)({3t^2}x-2t^3)-(x-t)({3a^2}x-2a^3)]/(t-a)

    ＝ {6a}x^2 - {9a^2}x + 4a^3

三項係數分別是正確值的 2、3、4 倍。

不像 {3a^2}x-2a^3 = lim_{t→a} [(x-a)(t^3)-(x-t)(a^3)]/(t-a) 那麼單純。
※ 編輯: Vulpix (163.13.112.58 臺灣), 08/14/2023 21:27:43
→ PPguest     : 感謝實例,我覺得我的認知和V大的一樣. 我是算下面
→ PPguest     : 這個：lim_{c→a} [(x-c)({3a^2}x-2a^3)-(x-a)
→ PPguest     : ({a^2+ac+c^2}x-ac{a+c})]/(a-c)
→ PPguest     : = {3a}x^2-{3a^2}x+a^3
推 znmkhxrw    : 嗨V大, 兩個回覆請教:
→ znmkhxrw    : (1) 02:02您提說用克拉瑪, 是另外一個方法嗎?
→ znmkhxrw    : 你這篇的矩陣方式應該不是用克拉瑪?
→ znmkhxrw    : 我把"用克拉瑪"翻譯成直接寫出各項係數的公式解
→ znmkhxrw    : 每個c_i的分子跟分母都是determinant, 然後取極限?
→ znmkhxrw    : 你的矩陣方式是Ac=f然後對A跟f取極限
→ znmkhxrw    : 如果用克拉瑪就是c=(A^-1)*f=cramer, 然後取極限
→ znmkhxrw    : 你提克拉瑪是指這樣算起來過程差不多?(結果論一樣)
→ znmkhxrw    : (2) 02:12您的回覆最後那邊寫了 "索性就不要把
→ znmkhxrw    :  f[a,a] 定義成 f'(a)", 這句話的"不要"是多打得?
→ znmkhxrw    : 您那段的意思我翻譯成"MVT要考慮很多, 索幸把f[a,a]
→ znmkhxrw    : 定義成f'(a)就好多了", 如果是這樣的話那我就沒問題
不是，就是乾脆不要先定義成這種形式。

因為 https://i.imgur.com/lylZcfd.png 會更方便。

至少在 f 有一定平滑性的時候這樣定義反而比較方便。

（對於不連續的 f，應該是有辦法把他也包容進來才對。）

1. 存在即連續，所以不用檢查太多次連續性。

2. 要是 x_n = x_0，還能利用 n 個變數的均差的對稱性快速得到

   https://i.imgur.com/88E2qLe.png 這條我們想要的微分公式。

這樣一來，每次均差的變數變多，首先就要證明 x_i 都可以對調。

pf): 1) x_n 可與 x_0 對調。

        如果 x_n 與 x_0 相同，那兩者對調確實不改變均差。

        即便兩者不同，定義中的極限變回普通的差分商，

        https://i.imgur.com/Mn5WvDa.png，兩者一樣可以對調。

     2) x_0 可與 x_1 對調。

        如果 x_0 與 x_1 相同，那兩者對調確實不改變均差。

        如果兩者不同，https://i.imgur.com/qt1duqL.png 說明還是可對調。

        其中有先假設變數不到 n+1 個的時候可以任意排列變數。

        以下補述少變數的情況，剩下的部份由數學歸納法遞迴可得。

        i) 雙變數的情況

           如果 a≠b：f[a,b] = (f(b)-f(a))/(b-a) = (f(a)-f(b))/(a-b) = f[b,a]

           至於 a=b：trivial case。

        ii) 三變數的情況

           如果 a,b,c 相異，可以直接寫出對稱式：

                f[a,b,c] = f(a)/(a-b)(a-c) + f(b)/(b-a)(b-c) + f(c)/(c-a)(c-b)

           如果 a＝b≠c：

                f[a,a,c] = f[c,a,a] = [(f(c)-f(a))/(c-a) -f'(a)]/(c-a)

                         = ( f(c)-f(a)-f'(a)(c-a) )/(c-a)^2

                f[a,c,a] = lim_{b→a} (f[c,b]-f[a,c])/(b-a)

                         = d/da f[a,c]

                         = ( f(c)-f(a)-f'(a)(c-a) )/(c-a)^2

           至於 a = b = c：trivial case。

     3) x_1,...,x_{n-1} 可任意排列。

        任取一個 1,...,n-1 之間的排列 σ，https://i.imgur.com/7SJV8cT.png。

     這些排列可以生成任意一個在 0,1,...,n 之間的排列，done。

→ znmkhxrw    : 最後就是您跟P大討論的涉及neville's algorithm
→ znmkhxrw    : , neville's表/差分表的top/not top row...之間的關
→ znmkhxrw    : 係, 目前對我資訊量有點大XD 我再慢慢看你們的討論
→ znmkhxrw    : 謝謝您這一系列的討論跟分享, 1000p聊表謝意
→ znmkhxrw    : 也謝謝P大的idea跟c大的代數證明, 666p奉上
沒事，那個你沒空可以暫時不去管，因為我覺得那個也有那個的麻煩。

不過解決 Newton form 的 top row 過於偏心的方案可能要了解一下。

※ 編輯: Vulpix (163.13.112.58 臺灣), 08/15/2023 02:55:44
推 znmkhxrw    : 嗨V大: (1) 用你https://i.imgur.com/lylZcfd.png
→ znmkhxrw    :  確實在一個一個點跑時方便很多. 不過我還是對於你
→ znmkhxrw    : 之前說的［反正怎麼搞這一步都跑不掉，索性就不要
→ znmkhxrw    : 把f[a,a]定義成f’(a) ］這句有點不理解, 這句話的
→ znmkhxrw    : 意思是［不管有沒有把f[a,a]定義成f’(a) 都還是避
→ znmkhxrw    : 不開最後一步的不使用羅必達直接求導數］嗎? 因為
→ znmkhxrw    : f[a,a]用你的定義算出來也是f’(a), 所以我不太清
→ znmkhxrw    : 楚這裡比較的對象是什麼…
是指在後面要計算 f[a,...,a] 的時候，該有的麻煩還是一樣在。

→ znmkhxrw    : (2) 總結來說是否如果要直接考慮多個點可能合併成
→ znmkhxrw    : 一個點的case很麻煩, 即便給很強的微分條件讓MVT所
→ znmkhxrw    : 出來的各階導數都是連續的, 還是要考慮各種形式的
→ znmkhxrw    : MVT(n個點有m個一樣blabla) 但是今天如果只考慮一
→ znmkhxrw    : 個一個點退化, 採用您lylZcfd.png的定義後就能有條
→ znmkhxrw    : 理的一個一個退化, 然後基於差分式的對稱性還能WL
→ znmkhxrw    : OG把退化順序排成退化的都在前面
推 znmkhxrw    : (3) Newton form偏心問題我沒follow到, 我一直都只
→ znmkhxrw    : 看top row, 但是發現top row對於［退化點是從f[a,
→ znmkhxrw    : a,b]這種型的就可以羅到爽］, 但是f[a,a,b]理應等
→ znmkhxrw    : 於f[b,a,a], 但是對f[b,a,a]用top row就很麻煩, 羅
→ znmkhxrw    : 下去有一堆東西要整理. 你說的top row偏心問題是指
→ znmkhxrw    : 這件事嗎? 如果是的話我大概從你跟P大討論不同row
→ znmkhxrw    : 猜測, 遇到f[b,a,a]又要用差分表的話, 就不要用to
→ znmkhxrw    : p row, 其他路徑都能得到一樣的結果?
top row 的問題是你另一篇回文一開始的問題。

因為 top row 就真的不能讓你直接看到 p(0)、p'(1) 那些東西。

這個只要看那篇的推文就好。

推 znmkhxrw    :  另外我剛剛回去看V大你最一開始提MVT那段, 你說f
→ znmkhxrw    : [a,b,c]對b,c做MVT可以寫成如下連結是基於哪種MVT
→ znmkhxrw    : 啊？我直接很誠實的對b,c做f[b,c]=f’(epsilon)然
→ znmkhxrw    : 後就沒了XDD
→ znmkhxrw    : https://i.imgur.com/oAtEbdl.jpg

對 b,c 做，那當然要先換成 f[b,a,c] 再做。

令 g(x) = f[a,x] = f[x,a]

因為 f'(a) 存在，所以 g 在 a 連續。

又 g[a,x] = f[a,a,x] = (f(x)-f(a)-f'(a)(x-a))/(x-a)^2 → f"(a)/2

其他 x 就不用說了，g(x) 肯定連續，用商法則就能算 g'(x) = f[a,x,x]。

總之 f 只要二階可微，那 g 就符合 MVT 的使用資格。

所以 f[b,a,c] = (f[a,c]-f[b,a])/(c-b) = g[b,c] = g'(ξ)

g'(ξ) 有兩種可能長相：

1. ξ = a，此時 g'(ξ) = g'(a) = f"(a)/2。

2. ξ≠a，此時 g'(ξ) = f[a,ξ,ξ] = (f'(ξ)(ξ-a)-f(ξ)+f(a))/(ξ-a)^2

在 b,c→a 的過程中，這兩個情況也可能交錯發生，

但只要 f[a,ξ,ξ] 也收斂到 f[a,a,a] = f"(a)/2 就好。

證明方式就是羅，上面兩串紫字都要羅一次。

不過 f[a,a,x] 可以直接羅，但 f[a,ξ,ξ] 要先換成 f'[a,ξ] - f[a,a,ξ] 才有羅。

這些該羅的，一個都跑不掉。


再開題外：f'[a,ξ] = f[a,ξ,ξ] + f[a,a,ξ] 這種形狀的公式，

          好像真的應該先做出來放著備用。
※ 編輯: Vulpix (163.13.112.58 臺灣), 08/15/2023 15:15:14

這應該是 f'[x_0,...,x_n] = Σ_{r=0}^n f[x_0,...,x_n, x_r]。

=> f"[x_0,...,x_n] = 2! Σ_{0≦r≦s≦n} f[x_0,...,x_n, x_r,x_s]

=> f"'[x_0,...,x_n] = 3! Σ_{0≦r≦s≦t≦n} f[x_0,...,x_n, x_r,x_s,x_t]

and so on.
※ 編輯: Vulpix (163.13.112.58 臺灣), 08/15/2023 16:26:02
推 znmkhxrw    : 了解~我花時間整理吸收一下, 再次謝謝V大
→ Vulpix      : 對稱性的證明怪怪的，我再想想。
推 znmkhxrw    : "用商法則就能算 g'(x) = f" 的 "=f"是typo?
對，大概本來是要打 f[a,x,x] 吧？

→ znmkhxrw    : by the way, 這些東西看起來好複交錯, 但是最終都能
→ znmkhxrw    : 拆成MVT跟L'hospital, 有種好複雜卻又有主幹的感覺
→ znmkhxrw    : 有那個Fu卻又還沒有完全理解的有趣感覺XDDD
直接 C^n 就沒有那麼多事了。只要求 Diff^n 是會釐清一些事情啦……

Def): The maximal degeneracy of (x_0,...,x_n)

      is the frequency of the mode of the sequence {x_i}_{i=0}^n.

總之就是，最大簡併數 := 眾數的出現次數。

f \in Diff^k => f[x_0,...,x_n] is well-defined at the points

                where the maximal degeneracy ≦ k+1.

                In particular, f[x_0,...,x_k] is well-defined everywhere.

f \in C^k => f[x_0,...,x_n] is continuous at the points

             where the maximal degeneracy ≦ k+1.

             In particular, f[x_0,...,x_k] is continuous everywhere.
※ 編輯: Vulpix (163.13.112.58 臺灣), 08/16/2023 04:38:19

其實放棄 Lagrange form 也是挺可惜的。

令 u_{ik}(x) 為 Π_{j≠i} {(x-x_j)/(x_i-x_j)}^{-1-m_j}

在 x = x_i 展開至 m_i-k 次項的多項式。

那麼 u_{ik}(x) * (x-x_i)^k/k! * Π_{j≠i} {(x-x_j)/(x_i-x_j)}^{1+m_j} 們

就是我們 Lagrange form 的基底……沒有很漂亮，但也沒有很不漂亮啦。

具體來說，https://i.imgur.com/ItFFMG9.png。

整個計算現在變成要算 https://i.imgur.com/LXZTetm.png，

而且因不看超過 m_i 次的項，最後一式也可變成 https://i.imgur.com/q85ceOt.png。

其實 u_{ik}(x) 可以用輾轉相除、長除法、泰勒展開等方式去做，方便就好。
※ 編輯: Vulpix (1.160.14.56 臺灣), 08/19/2023 08:45:44