[分析] 拉格朗日插值的退化理論與實作

[znmkhxrw (QQ)]

 
請問各位板友一個對函數做Lagrange interpolation的問題, 先看這個連結
 
https://www.desmos.com/calculator/cegj2qpzez
 
先定義名詞: (1) 強N次多項式: x^N的領導係數不為0
            (2) N次多項式: 多項式至多N次方
 
再來回到圖, P(x)是一個通過(1,0), (s,log_2(s)), (2,1)三個點的強二次多項式
即通過三點的拉格朗日插值多項式
 
接著滑動s時, 會發現s越趨近於1^+時, P(x)越接近P_L(x)
而s越趨近於2^-時, P(x)越接近P_R(x)
先描述以上這種現象叫作退化 
 
然後從這個例子觀察出兩個現象:
 
(1) 退化成功, 即係數的極限值都存在
(2) 退化後還是強二次多項式, 並不是變成通過(1,0),(2,1)兩個點的一次多項式
 
然後從其他例子發現:
 
(1) 對f(x)=x做(0,0),(s,s),(1,1)的二次多項式內插, 出來的多項式仍是P(x)=x
    即便退化後仍是P(x)=x
    也就是說, 原本三個點決定的多項式式是強一次多項式, 退化後仍是強一次多項式
 
(2) 對log_2(x)做(1,0), (s,log_2(s)), (p, log_2(p)), 的強二次多項式內插
    先把s逼近到1, 仍是強二次多項式, 再把p逼近到1, 也還是強二次多項式
    不論是單退化(三點剩二點)還是雙退化(三點剩一點) 
    退化後仍是極限存在的強二次多項式 
 
(2) 對log_2(x)做(1,0), (s,log_2(s)), (p, log_2(p)), (r,log_2(r))
    的強三次多項式內插, 不論單退化、雙退化、三退化
    退化後仍是極限存在的強三次多項式 
 
因此我想請問跟猜測的兩個面向是:
 
【理論】 
 
上述的log_2(x)現象是巧合還是有定理支持?
我猜測的定理是:
 
令f(x)為定義在[a,b]的連續函數
在[a,b]中任取N個相異點a_1,...,a_N, 對(a_1,f(a_1)),...,(a_N,f(a_N))
造出N-1次的多項式P(x), (此P(x)唯一並且至多N-1次)
則我們有: (1) 對任何n<=N-1, n點的退化是均勻收斂
          (2) 若P(x)是強M次多項式, 上述的退化亦是強M次多項式
 
這裡的任何n點退化很廣泛, 即任挑n點、任意順序的退化, 都會收斂並且還是強M次 
(退化就是依序把某個a_i當變數, 其他a_j固定, 對多項式P(x)取lim{a_i→a_j})
也就是說強M次的M已經被當初的a_1,...,a_N決定了, 不管怎麼退化仍會是強M次多項式
 
【實作】 
 
我在 #1ajhtDXJ 問極限的電腦實作問題就是源自於這篇
我想對(1,0), (s,log_2(s)), (p, log_2(p)), (2,1)這四個點做三次多項式
但是實務上我無法控制使用者的s跟p會不會很靠近到浮點數的相等
也無法控制s or p會不會很靠近1 or 2到浮點數的相等
而即便浮點數不相等, 但是很靠近的話, 會導致這個三次多項式的closed form的係數
出現0除以0的情形
 
這個問題有其他解法嗎? 還是一樣只能靠我上篇問的方式, 設threshold/大數運算來解
----------------------------------------------------------
 
以上兩個面向詢問, 謝謝幫忙~
 
--
※ 發信站: 批踢踢實業坊(ptt.cc), 來自: 59.102.225.191 (臺灣)
※ 文章網址: https://www.ptt.cc/bbs/Math/M.1689965257.A.465.html
推 LPH66       : 沒有很仔細想, 不過你的逼近退化法看起來等同於
→ LPH66       : 取重合點的微分值等於中間的函數在該重合點的微分
 
確實很吃函數的微分值, 令f(x)=cos(2*pi*x), a=0, b=1, s€(a,b)
則造P(x)為通過(0,f(0)), (s,f(s)), (1,f(s))的強二次多項式
會發現s不管退化成0還是1, 最後退化出的多項式都變成常數
 
https://www.desmos.com/calculator/aqnxzlmcgr
 
光三個點要整理成ax^2+bx+c然後做羅畢達就花不少時間算
四個點我好沒有勇氣XDDDD
 
推 LPH66       : 又仔細想了一下應該沒錯, 你把 (t, g(t)) "退化"和
→ LPH66       : (0, g(0)) 合併之後的函數在 0 的微分會等於 g'(0)
→ LPH66       : 因此你想證的東西強多項式性質應該都有辦法調微分值
→ LPH66       : 去造出反例來
→ LPH66       : 在 g 是對數的狀況下可能因為對數函數的凸性
→ LPH66       : 造成一些可能會"弱化"的情形在對數函數裡找不到
→ LPH66       : 02:14 那兩行的結論應該能用拉格朗日插值公式
→ LPH66       : 把要經過的點給拆開, 然後就能只考慮 (0, g(0)) 和
→ LPH66       : (t, g(t)) 兩點的狀況, 這時極限狀況的直線
→ LPH66       : 就是 g 在 0 的切線 (即是這 g'(0) 是直接由定義得)
→ LPH66       : 咦等等, 插值公式沒有保證微分值為 0..這我要再想想
→ LPH66       : 啊有了, 因為是合併後的極限, 兩點合一變重根後
→ LPH66       : 這重根的重覆因子就能讓微分值為 0
→ LPH66       : 所以全部加起來之後微分值確實會變成 g'(0)
 
嗨L大, 謝謝分享, 不過我不太理解你02:29~02:34的意思
看起來你說的"合併"是點的合併, 但是這樣不就降次了嗎
以log_2(x)這個例子來說, 通過(1,0),(s,log_2(s)), (2,1)三個點的強二次多項式P(x)
對s→0還是得到強二次多項式P_L(x), 然後會通過(1,0), (2,1)兩個點
然後當然也可以針對(1,0),(2,1)這兩個點直接造出強一次多項式L(x)
 
你那幾行的回應是在討論如何從P(x)退化到L(x)?

推 LPH66       : 我是在講你開頭寫的"退化"沒錯
→ LPH66       : 算是非正式地在證明那個退化就是取"合併"點上的微分
了解~

推 PPguest     : 用divided difference, Newton form的interpolation
→ PPguest     : polynomial來看 f(x)=cos(2*pi*x) 的例子，
→ PPguest     : 因為f(0)=f(1),在退化時造成1st divided difference
→ PPguest     : 兩個都變成0(一個是因f'(0)=0),因此2nd divided
→ PPguest     : difference也是0,因此多項式的二次以及一次項係數都
→ PPguest     : 是0，結果就變成常數
→ PPguest     : 能確定的是,若n個點都退化到同一個點,對一般的f來說
→ PPguest     : n-1次項的係數即為(n-1)th divided difference,
→ PPguest     : 其值等於f微n-1次在那一點取值，然後除以(n-1)!
沒用過divided diff, 剛剛去查發現跟Lagrange Interpolation有等式上的關係式
想問P大你推文的內容, 感覺可以寫出任何退化的closed form多項式!?
如果可以的話我之後去研究一下, 謝謝!

如我上面給的連結 https://www.desmos.com/calculator/aqnxzlmcgr
光是算出三個點中的一個點的左右退化就好辛苦
因為我需要把Lagrange內插整理成a_n*x^n+...+a_0後再逐一對係數做羅畢達
我不敢想像四個點的各種退化情況會多複雜...

→ PPguest     : Lagrange form 的 interpolation 和 Newton form 的
→ PPguest     : 都是同一個多項式,只是寫法不同
Soga!

推 PPguest     : 另外不知道你的最後目的是什麼？是要求某點代入多項
→ PPguest     : 示的值？或是有很多不同的x要代入多項式？又或者只
→ PPguest     : 要a_0......a_n的值,不需要代值？
→ PPguest     : 某點代入多項式的話可用Neville's algorithm
→ PPguest     : 不同的x要代入多項式的話,先求Newton form的係數,再
→ PPguest     : 用Newton form算x代入多項式
我是嵌入式系統上要實作某個非線性函數, 需要用到log_2(x)
其中我用多項式P(x)來逼近log_2(x), 1<=x<=2
而這個非線性函數一般會有兩個不連續點, 經過計算我只要讓P(x)通過
(1,0), (s,log_2(s)), (p,log_2(p)), (2,1)這四個點, 那就可以讓那兩個不連續點消失
而s跟p的決定是使用者會輸入一些參數t而轉換來的, 也就是說s=U_1(t), p=U_2(t)

雖然U_1,U_2這兩個轉換函數已知, 但是我自己跑參數確實會讓一些t導致:
s=p or s=1 or p=1 or s=p=1...
這也是為什麼我會問Lagrange退化的公式與實做的原因
因為如果知道退化公式, 我可以設定某個threshold讓內插點相等或極近時代入退化公式
同時我又擔心內插點極近但是沒有過我的threshold所造成的P(x)誤差問題
(因為P(x)這個多項式的係數充滿著內插點相似而產生0+/0+情況)

也就是說, 使用者決定t→程式轉換成s,p→程式轉換成多項式P(x)的四個係數a_3~a_0
然後開機運行時會一直進來x, 就會讓P(x)跑a_3*x^3+a_2*x^2+a_1*x+a_0

→ PPguest     : 剛想了一下,07/23 16:00 講的是一般情況,不知道是否
→ PPguest     : 適合退化的情形.另外我想原po應該知道算P(x)時,用
→ PPguest     : ((a_3*x+a_2)*x+a_1)*x+a_0會比較穩定之類的吧
我是這樣實作沒錯, 不過不是因為精度因素, 而是這樣的乘法數量比較少XDD
你說的"穩定"是這樣也有精度比較高的好處?

推 PPguest     : 感覺要先確定一件事,如果參數t確定後變成2~3個點,多
→ PPguest     : 項式是退化的最好嗎？用一次或二次多項式會比較差嗎
在log_2(x)這個case中, 係數極限退化是比較好的沒錯
點數減少的退化方式會降次方, 誤差整個就上去了QQ

→ PPguest     : 三次多項式有4個未知係數,已知例如過2點,還有兩個條
→ PPguest     : 件的自由度可以選擇
喔喔 你的意思是如果退化, 那就隨便補點來維持三次多項式
這樣就繞過手動計算極限係數了, 這approach解決我一半問題耶XDD 謝謝

→ PPguest     : 我想表達的是就算是三次多項式還有很多選擇,不過會
→ PPguest     : 不會符合退化的那種是最"好"的？如果是的話,搞清楚
→ PPguest     : 其另外加入的條件是什麼,然後在看計算上用哪種方式
→ PPguest     : 會比較適合
→ PPguest     :                                    ^係數
維持三次多項式的話, 確實可以細部處理以下case:
(1) s=p=1    :選擇額外兩點(p1,log_2(p1)), (p2,log_2(p2))使得誤差最小
(2) s=1, s!=p:選擇額外一點(p1,log_2(p1))使得誤差最小
(3) p=1, s!=p:選擇額外一點(p1,log_2(p1))使得誤差最小

不過如果s,p,1,2這四點皆相異的係數精度是夠的話, 我是打算採取最簡單的方式:
如果這四點有誰相等, 就加個float_epsilon, 強制讓這三個點是浮點不相等
然後照原四點的內插公式算係數

推 PPguest     : Hermite interpolation 看起來是除了各點函數值要一
→ PPguest     : 樣,各點例如微一次以及微兩次後的函數值也要一樣
→ PPguest     : 看起來前面提的點變少的退化,猜測其他的條件應該就
→ PPguest     : 是對該點微分一次(甚至微分兩次)的函數值也要一樣
這我也沒聽過QQ 查了是數值分析領域的方法
P大提這個內插方式是優化哪個部份呢?

→ PPguest     : 如果跟點數一樣的比,就是用更多的資訊,更高次的多項
→ PPguest     : 式來逼近.如果跟一樣次數的多項式比,大概就是誤差項
→ PPguest     : 不太一樣
了解~

推 Vulpix      : Hermite插值就是你的問題啊。內容就是用函數值與各
→ Vulpix      : 階導數把多項式找出來。如果插值用的點只剩一個a（
→ Vulpix      : 其他a_i都一個個收斂到a了），那插出來的就是泰勒
→ Vulpix      : 多項式。如果插值用的資料不含導數，那就會得到拉
→ Vulpix      : 格朗日插值公式。
→ Vulpix      : 例如已知f(1),f'(1),f(2)就會插出P_L。
原來Taylor跟Lagrange都是Hermite的特例! 之前完全沒研究過Hermite說
我之後再理論研究, 謝謝V大的idea
※ 編輯: znmkhxrw (59.102.225.191 臺灣), 07/24/2023 00:15:54