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

符號看得有點花…… 如果你想做的是「在 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 先修一點 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 要 f[x,y] 連續至少要 C^1，因為沿著 y=x 靠近 (a,a) 的時候要 f' 連續。 ※ 編輯: Vulpix (163.13.112.58 臺灣), 08/08/2023 17:31:27 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 找到比較一勞永逸的方式： 首先，這次既不是用 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 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。 那就把次方都改成 max(0, 原本的次方) 好了…… 或者乾脆改成原次方的絕對值。 不過反正負次方項的係數都是 0 啦，當成 notation 看就好，沒有什麼代值的功能。 首先，不管 Newton、Lagrange 插出來的多項式都跟 https://i.imgur.com/sPmLIUo.png 解出來的多項式是同一個。 所以雖然我算的是這個多項式的退化，但可以直接當成 Newton form 的極限，沒問題。 你把這個矩陣方程式展開，其實每一條都只是 p(x_i) = f(x_i) 而已。 因為我想要都用 divided difference 寫，所以會比 CV 多除以一些階乘。 也是因為這樣，我的係數才是組合數 C 而不是 CV 的排列數 P。 對，但我不要。這樣的行列式裡面會有一大堆階乘。 我的證明過程須要算行列式，不想讓他變醜。 數學歸納法。 不過我還在想到底有哪些東西要定義跟證明XD 以這個流程來說，至少就是 f[x_1, ..., x_n] 收斂到 f^{(n)}(x_0)/n! 吧。 是啊。 ※ 編輯: Vulpix (1.160.12.97 臺灣), 08/13/2023 22:58:39 既然反矩陣麻煩，那為什麼不用克拉瑪公式呢？ 作者: 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