拉格朗日插值多項式的原理介紹及其應用

  插值，不論在數學中的數值分析中，還是在我們實際生産生活中，都不難發現它的身影，比如造船業和飛機制造業中的三次樣條曲線。那麼，什麼是插值呢？我們可以先看一下插值的定義，如下：

  （定義）如果對于每個1≤i≤n,P(xi)=yi1≤i≤n,P(xi)=yi,則稱函數y=P(x)y=P(x)插值資料點(x1,y1),...,(xn,yn)(x1,y1),...,(xn,yn).

   插值的定義無疑是清楚明了的，而在衆多的數學函數中，多項式無疑是最簡單，最常見的函數，關于它的理論研究也最為透徹。是以，我們可以不妨先考慮利用多項式來進行插值。那麼，這樣的多項式是否總是存在呢？答案是肯定的，因為我們有如下定理：

  （多項式插值定理）令(x1,y1),...,(xn,yn)(x1,y1),...,(xn,yn)是平面中的nn個點，各xixi互不相同。則有且僅有一個n−1n−1次或者更低的多項式PP滿足P(xi)=yi,i=1,2,...,n.P(xi)=yi,i=1,2,...,n.

  證明:先用歸納法證明存在性，再證明唯一性。

  當n=1n=1時，常函數(0次)P1(x)=y1P1(x)=y1即符合要求。假設當n−1n−1時存在一個次數≤n−2≤n−2的多項式Pn−1Pn−1,使得Pn−1(xi)=yi,i=1,2,...,n−1.Pn−1(xi)=yi,i=1,2,...,n−1.則令Pn(x)=Pn−1(x)+c(x−x1)(x−x2)...(x−xn−1)(x−xn)Pn(x)=Pn−1(x)+c(x−x1)(x−x2)...(x−xn−1)(x−xn)，其中cc為待定系數，利用Pn(xn)=ynPn(xn)=yn即可求出待定系數cc.此時，Pn(xi)=yi,i=1,2,...,n,Pn(xi)=yi,i=1,2,...,n,且Pn(x)Pn(x)的次數≤n−1≤n−1.這樣就證明了存在性。

  其次證明唯一性。假設存在兩個這樣的多項式，設為P(x)P(x)和Q(x)Q(x)，它們次數≤n−1≤n−1且都插值經過nn個點，即P(xi)=Q(xi)=yi,i=1,2,...,n.P(xi)=Q(xi)=yi,i=1,2,...,n.令H(x)=P(x)−Q(x)H(x)=P(x)−Q(x),HH的次數也≤n−1≤n−1，且有nn個不同的根x1,x2,...,xnx1,x2,...,xn.是以，由多項式基本定理可知，H(x)H(x)為0多項式，即恒等于0，故有P(x)=Q(x)P(x)=Q(x).這樣就證明了存在性。

  證畢。

  有了以上定理，我們可以放心地使用多項式進行插值，同時，通過上述定理，我們可以用歸納法來構造此多項式，但是，這樣的方法難免複雜麻煩。于是，天才的法國數學家拉格朗日（Lagrange）創造性地發明了一種實用的插值多項式方法來解決這個問題，那麼，他的方法是怎麼樣的？

  一般來說，如果我們有nn個點(x1,y1),...,(xn,yn)(x1,y1),...,(xn,yn)，各xixi互不相同。對于1到n之間的每個kk，定義n−1n−1次多項式

Lk(x)=(x−x1)..(x−xk−1)(x−xk+1)...(x−xn)(xk−x1)..(xk−xk−1)(xk−xk+1)...(xk−xn)Lk(x)=(x−x1)..(x−xk−1)(x−xk+1)...(x−xn)(xk−x1)..(xk−xk−1)(xk−xk+1)...(xk−xn)

Lk(x)Lk(x)具有有趣的性質：Lk(xk)=1,Lk(xj)=0,j≠k.Lk(xk)=1,Lk(xj)=0,j≠k.然後定義一個n−1n−1次多項式

Pn−1(x)=y1L1(x)+...+ynLn(x).Pn−1(x)=y1L1(x)+...+ynLn(x).

這樣的多項式Pn−1(x)Pn−1(x)滿足Pn−1(xi)=yi,i=1,2,...,n.Pn−1(xi)=yi,i=1,2,...,n.這就是著名的拉格朗日插值多項式！

  以上就是拉格朗日插值多項式的理論介紹部分，接下來我們就要用Python中的Sympy子產品來實作拉格朗日插值多項式啦~~

  實作拉格朗日插值多項式的Python代碼如下：

from sympy import *

def Lagrange_interpolation(keys, values):
    x = symbols('x')
    t = len(keys)
    ploy = []
    for i in range(t):
        lst = ['((x-'+str(_)+')/('+str(keys[i])+'-'+str(_)+'))' for _ in keys if _ != keys[i]]
        item = '*'.join(lst)
        ploy.append(str(values[i])+'*'+item)
    ploy = '+'.join(ploy)

    return factor(expand(ploy))

def main():
    #example 1, interpolate a line 
    x_1 = [1,2]
    y_1 = [3,5]
    if len(x_1) != len(y_1):
        print('The lengths of two list are not equal!')
    else:
        print('Lagrange_interpolation polynomials is:')
        print(Lagrange_interpolation(x_1,y_1))

    #example 2, interpolate a parabola
    x_2 = [0,2,3]
    y_2 = [1,2,4]
    if len(x_2) != len(y_2):
        print('The lengths of two list are not equal!')
    else:
        print('Lagrange_interpolation polynomials is:')
        print(Lagrange_interpolation(x_2,y_2))

    #example 3
    x_3 = [0,1,2,3]
    y_3 = [2,1,0,-1]
    if len(x_3) != len(y_3):
        print('The lengths of two list are not equal!')
    else:
        print('Lagrange_interpolation polynomials is:')
        print(Lagrange_interpolation(x_3,y_3))

main()           

函數Lagrange_interpolation()具體實作了拉格朗日插值多項式，參數(keys, values)為list形式的點對，在main()函數中舉了三個Lagrange_interpolation()函數的應用執行個體，一個是插值兩個點，即直線，一個是插值三個點，即抛物線，一個是插值四個點，但結果卻是一次多項式。該程式的運作結果如下：

  接下來，我們将介紹一個拉格朗日插值多項式的應用，即求

1k+2k+...+xk1k+2k+...+xk

的求和公式，其中x,kx,k為正整數。分析如下:

  首先，該求和公式應當是一個至多為k+1次的關于xx的多項式。然後，我們可以通過取k+2個不同的點，利用拉格朗日插值多項式的辦法來求解，這k+2個不同的點的橫坐标可以取x=1,2,...,k+2x=1,2,...,k+2,在求出其對應的縱坐标的值。

  以下代碼分别求出k=1,2,...,50k=1,2,...,50的求和公式，并将其插入到Redis中。

from sympy import *
import redis

def Lagrange_interpolation(keys, values):
    x = symbols('x')
    t = len(keys)
    ploy = []
    for i in range(t):
        lst = ['((x-'+str(_)+')/('+str(keys[i])+'-'+str(_)+'))' for _ in keys if _ != keys[i]]
        item = '*'.join(lst)
        ploy.append(str(values[i])+'*'+item)
    ploy = '+'.join(ploy)

    return factor(expand(ploy))

def degree_of_sum(k):
    x_list, y_list = [], []
    degree = k    # degree=k in expression of  1^k+2^k+...+x^{k}
    cul_sum = 0
    for i in range(1,degree+3):
        x_list.append(i)
        cul_sum += i**degree
        y_list.append(cul_sum)
    return Lagrange_interpolation(x_list,y_list)

def main(): 
    r = redis.Redis(host='localhost', port=6379,db=0)
    for k in range(1,51):
        expression = str(degree_of_sum(k))
        r.hset('sum_%s'%k,'degree',str(k))
        r.hset('sum_%s'%k,'expression',expression)
        print('Degree of %d inserted!'%k)

main()           

運作以上程式，結果如下：

在Redis中的儲存結果如下：

我們可以具體檢視當k=2k=2時的求和公式，如下：

  這樣我們就介紹完了一個拉格朗日插值多項式的應用了。看了上面的介紹，聰明又機智的你是否能想到更多拉格朗日插值多項式的應用呢？歡迎大家交流哦~~

  新的一年，新的氣象，就從這一篇開始~~

注意：本人現已開通兩個微信公衆号： 用Python做數學（微信号為：python_math）以及輕松學會Python爬蟲（微信号為：easy_web_scrape）， 歡迎大家關注哦~~