拉格朗日插值原理及實現(Python)

来源:https://www.cnblogs.com/yang-ding/archive/2022/09/25/16728202.html
-Advertisement-
Play Games

拉格朗日插值原理及實現(Python) 一. 前言 Lagrange插值是利用n次多項式來擬合**(n+1)個數據點**從而得到插值函數的方法。(註意n次多項式的定義是未知數最高次冪為n,但是多項式繫數有n+1個,因為還有個常數項) **Lagrange插值和Newton插值本質上相同,都是用(n- ...


拉格朗日插值原理及實現(Python)

目錄

一. 前言

Lagrange插值是利用n次多項式來擬合(n+1)個數據點從而得到插值函數的方法。(註意n次多項式的定義是未知數最高次冪為n,但是多項式繫數有n+1個,因為還有個常數項)

Lagrange插值和Newton插值本質上相同,都是用(n-1)次多項式來擬合n個數據點。所以這兩種插值方法得到的插值函數相同,因為多項式擬合的基本定理:同時通過n個數據點,且最高次冪小於(n-1)的多項式函數唯一。下麵順手證明一下這個重要的定理。

如果已知n+1個數據點\((x_0,y_0),(x_1,y_1),(x_2,y_2),\cdots,(x_n,y_n)\),假設\(L_1=k_0+k_1x+k_2x^2+\cdots+k_nx^n\)\(L_2=k_0’+k_1’x+k_2'x^2+\cdots+k_n’x^n\)都是通過這n個數據點的插值函數。那麼應該有\(L_1 - L_2\)通過所有\((x_1,0),(x_2,0),\cdots,(x_n,0)\)

代入這些點得到齊次線性方程組:

\[\begin{bmatrix} 1 & x_0 & x_0^2 & \cdots & x_0^n\\ 1 & x_1 & x_1^2 & \cdots & x_1^n\\ \vdots & \vdots & \vdots & & \vdots \\ 1 & x_n & x_n^2 & \cdots & x_n^n\\ \end{bmatrix} \begin{bmatrix} k_0-k_0'\\ k_1-k_1'\\ k_2-k_2'\\ \vdots \\ k_n-k_n'\\ \end{bmatrix} = \begin{bmatrix} 0 \\ 0 \\ 0 \\ \vdots \\ 0 \end{bmatrix} \]

它的繫數行列式是Vandermonde行列式,所以:

\[\begin{vmatrix} 1 & x_0 & x_0^2 & \cdots & x_0^n\\ 1 & x_1 & x_1^2 & \cdots & x_1^n\\ \vdots & \vdots & \vdots & & \vdots \\ 1 & x_n & x_n^2 & \cdots & x_n^n\\ \end{vmatrix} = \prod_{n\ge i > j \ge 0 }(x_i-x_j) \]

因為每個點都是不同的,所以\(x_i \ne x_j,i\ne j\),所以齊次線性方程組的繫數行列式不等於0,故方程解唯一且為0解:

\[k_i - k_i' = 0\\ k_i = k_i' \\ L_1 = L_2 \]

證畢。

二. 3種形式的Lagrange插值函數推導

1. 原始形態的Lagrange插值


為了用n次多項式擬合n+1個數據點:\((x_0,y_0),(x_1,y_1),(x_2,y_2),\cdots,(x_n,y_n)\)

Lagrange插值函數採用的方法是構建一個這樣的函數:

\[L(x) = l_0(x)y_0 + l_1(x)y_1 + \cdots +l_n(x)y_n \tag{1} \]

也就是用一組基函數\(\{l_0(x),l_1{x}, \cdots ,l_n(x)\}\)去構建插值函數\(L(x)\),那麼不難想到這樣的基函數需要滿足這樣的條件:

\[l_i(x_j)=\begin{cases} 0 , & j\ne i \\ 1 , & j = i \end{cases} \]

這樣對於\(L(x)\)就會有:

\[\begin{split} L(x_i) &= l_0(x_i)y_0 + \cdots l_i(x_i)y_i + \cdots +l_n(x)y_n \\ &=y_i \end{split} \]

這樣就實現了\(L(x)\)通過所有的數據點\((x_i,y_i)\)。接下來就構建這樣的基函數\(l_i(x)\)

首先實現\(l_i(x_j) = 0,i\ne j\)

\[l_i(x) = (x-x_0)(x -x_1)\cdots (x-x_{i-1})(x-x_{i+1}) \cdots(x-x_n) \]

這個函數實現了\(l_i(x)\)在所有非\((x_i,y_i)\)的點處為0,但是在\(x=x_i\)處:

\[l_i(x_i) = (x_i-x_0)(x_i -x_1)\cdots (x_i-x_{i-1})(x_i-x_{i+1})\cdots (x_i-x_n) \]

為了讓\(l_i(x_i)=1\),我們可以將\(l_i(x)\)除以這個值進行修正:

\[\begin{split} l_i(x_i)&= \frac{(x-x_0)(x -x_1)\cdots (x-x_{i-1})(x-x_{i+1})\cdots(x-x_n)}{(x_i-x_0)(x_i -x_1)\cdots (x_i-x_{i-1})(x_i-x_{i+1})\cdots (x_i-x_n)}\\ &=\prod_{n\ge j \ge0,j\ne i} \frac{x-x_j}{x_i -x_j} \end{split} \tag{2} \]

將(2)代入(1)就得到了原始形態的Lagrange插值函數.

\[\begin{split} L(x) &= l_0(x)y_0 + l_1(x)y_1 + \cdots +l_n(x)y_n \\ &= \prod_{n\ge j \ge0,j\ne 0} \frac{x-x_j}{x_0 -x_j}y_0 +\prod_{n\ge j \ge0,j\ne 1} \frac{x-x_j}{x_1 -x_j}y_1 +\cdots +\prod_{n\ge j \ge0,j\ne n} \frac{x-x_j}{x_n -x_j}y_n \end{split} \tag{3} \]

例:已知點(1,1),(2,2),(3,3),用原始Lagrange插值計算插值函數。

首先計算三個插值基函數:

\[l_0(x) = \frac{(x-2)(x-3)}{(1-2)(1-3)}=\frac{(x-2)(x-3)}{2}\\ l_1(x) = \frac{(x-1)(x-3)}{(2-1)(2-3)}=-{(x-1)(x-3)}\\ l_2(x) = \frac{(x-1)(x-2)}{(3-1)(3-2)}=\frac{(x-1)(x-2)}{2}\\ \]

從而得到插值函數:

\[L(X) = l_0(x)+ 2l_1(x)+3l_2(x)=\frac{(x-2)(x-3)}{2}-2{(x-1)(x-3)}+ \frac{3(x-1)(x-2)}{2} \]

插入點x=4試一下:\(L(4)=4\)

原始模式的Lagrange插值函數,每次計算插值點時需要進行n(n-1)次乘法,時間複雜度為\(O(n^2)\)

2. 第一形式Lagrange插值


為了降低計算的時間複雜度,我們對原始的Lagrange插值函數進行改進。

為了書寫方便,我們令\(w_i\)

\[\begin{split} w_i &= \prod_{n\ge j \ge0,j\ne i} (x_i -x_j) \\ &= (x_i-x_0)(x_i -x_1)\cdots (x_i-x_{i-1})(x_i-x_{i+1})\cdots (x_i-x_n) \end{split} \tag{4} \]

觀察(3)式,我們可以提取一個公因數\((x-x_0)(x-x_1)\cdots (x-x_n) = g(x)\)

\[\begin{split} L(x) &= (x-x_0)(x-x_1)\cdots (x-x_n)[\prod_{n\ge j \ge0,j\ne 0} \frac{y_0}{(x_0 -x_j)(x-x_0)} +\prod_{n\ge j \ge0,j\ne 1} \frac{y_1}{(x_1 -x_j)(x-x_1)} +\cdots +\prod_{n\ge j \ge0,j\ne n} \frac{y_n}{(x_n -x_j)(x-x_n)}]\\ &=g(x)[\frac{y_0}{w_0(x-x_0)}+\frac{y_1}{w_1(x-x_1)}+\cdots+\frac{y_n}{w_n(x-x_n)}]\\ &= g(x)\sum_{i=0}^{n}\frac{y_i}{w_i(x-x_i)} \end{split} \tag{5} \]

這樣就得到了第一形式Lagrange插值。這一形式的Lagrange插值計算的流程如下:

  1. 預處理:根據已知數據點,利用公式(4)計算\(w_0,w_1,\cdots,w_n\);
  2. 插值:計算\(g(x)=(x-x_0)(x-x_1)\cdots(x-x_n)\),然後根據公式(5),先計算中括弧內的加式,然後乘以\(g(x)\)得到插值點的值,這樣計算一個新的點的時間複雜度就是\(O(n)\)
  3. 補充數據點:如果數據集有更新,只需要更新\(w_i\)即可,加入一個新的點到數據集,只需要將每個\(w_i\)乘以\((x_i-x_{n+1})\),此外再增加一個\(w_{n+1}\)

例:已知點(1,1),(2,2),(3,3),使用Lagrange第一形式計算插值函數。

首先計算\(w_i\)

\[w_0=(1-2)(1-3)=2\\ w_1=(2-1)(2-3)=-1\\ w_2=(3-1)(3-2)=2 \]

插值函數就是:

\[L(x)=(x-1)(x-2)(x-3)[\frac{1}{2(x-1)}-\frac{2}{(x-2)}+\frac{3}{2(x-3)}] \]

插入點x=4試一下:\(L(4)=4\)

3. 第二形式的Lagrange插值(重心插值公式)


第一形式的Lagrange插值還要計算\(g(x)\),可以再進行優化。

第一形式:

\[L(x) = g(x)\sum_{i=0}^{n}\frac{y_i}{w_i(x-x_i)} \tag{6} \]

為了消掉\(g(x)\),我們取\(y=1\)這條直線上的點進行插值,即取n+1個點:\((x_0,1),(x_1,1),\cdots,(x_n,1)\)那麼這n+1個點的插值函數就是:

\[L'(x)=g(x)\sum_{i=0}^{n}\frac{1}{w_i(x-x_i)} \tag{7} \]

\(L'(x)=1,x=x_0,x_1,\cdots,x_n\) ,所以我們可以用(6)除以(7)消去g(x):

\[L(x)=\frac{L(x)}{1}=\frac{L(x)}{L'(x)}=\frac{\sum_{i=0}^{n}\frac{y_i}{w_i(x-x_i)}}{\sum_{i=0}^{n}\frac{1}{w_i(x-x_i)}} \tag{8} \]

這樣就得到了第二形式的Lagrange插值,也稱為重心插值,通常Lagrange插值採用這種形式。

它的計算過程如下:

  1. 預處理:計算\(w_i\)
  2. 插值:對需要插入的點(x,y),計算(x-x_i)(可以存到一個列表中,計算時直接取用);
  3. 補充數據點:補充新的點到數據集只需要更新\(w_i\)

例:已知點(1,1),(2,2),(3,3),利用重心插值公式計算插值函數。

首先計算\(w_i\)

\[w_0=(1-2)(1-3)=2\\ w_1=(2-1)(2-3)=-1\\ w_2=(3-1)(3-2)=2 \]

得到\(L(x)\)

\[L(x)=\frac{\frac{1}{2(x-1)}-\frac{2}{(x-2)}+\frac{3}{2(x-3)}}{\frac{1}{2(x-1)}-\frac{1}{(x-2)}+\frac{1}{2(x-3)}} \]

插入點x=4試一下:\(L(4)=4\)

三. 利用Python編程實現這三種Lagrange插值

import numpy as np


class LagrangeInterpolation:
    """
    There are three modes of the LagrangeInterpolation.The input data is supposed
    to be two lists.
    ---------------------------------------------------------------------------------------
    self.data : Contain the points known before.
    self.dataLength : Indicate the number of points in the data list.
    self.weight : The self.weight[i] is (xi - x1)(xi - x2)...(xi - x(i-1))(xi - x(i+1))
                  ...(xi - xn)
    self.items : The self.items[i] is (x - x1)(x - x2)...(x - x(i-1))(x - x(i+1))...
                 (x - xn)
    ---------------------------------------------------------------------------------------
    """
    def __init__(self, x, y):
        self.data = {'x': list(x), 'y': list(y)}
        self.dataLength = len(self.data['x'])
        self.weight = []
        self.items = []
        # control is a flag indicating if there is anything wrong with the data.
        self.control = True
        if len(self.data['x']) != len(self.data['y']):
            print("The length of x isn't equal to the length of y!")
            self.control = False
        else:
            self.__preprocess(order=0)

    # Appending function is used to add points to the data list.
    def data_append(self, ap):
        if ap[0] in self.data['x']:
            print("The point already exist.")
            self.control = False
        else:
            self.control = True
            self.data['x'].append(ap[0])
            self.data['y'].append(ap[1])
            self.dataLength = self.dataLength + 1
        self.__preprocess(order=1)

    """
    Preprocessing is used to update the self.weight and self.items.
    order = 0 : Initialize the self.weight.
    order = 1 : Update the self.weight when new point is added to the data list.
    order = 2 : Calculate the self.items for each point waiting for interpolation.
    """
    def __preprocess(self, order=0, x=0):
        if order == 0:
            self.weight = list(np.zeros(self.dataLength))
            for i in range(self.dataLength):
                weight_temp = 1
                for j in range(self.dataLength):
                    if i == j:
                        pass
                    else:
                        weight_temp = weight_temp * (self.data['x'][i] - self.data['x'][j])
                self.weight[i] = weight_temp
        elif order == 1:
            self.weight.append(1)
            for i in range(self.dataLength - 1):
                self.weight[i] = self.weight[i] * (self.data['x'][i] - self.data['x'][-1])
                self.weight[-1] = self.weight[-1] * (self.data['x'][-1] - self.data['x'][i])
        elif order == 2:
            self.items = list(np.zeros(self.dataLength))
            for i in range(self.dataLength):
                self.items[i] = x - self.data['x'][i]

    # The mode1 is the initial mode of Lagrange interpolation.
    def mode1(self, px):
        self.__preprocess(order=2, x=px)
        if self.control:
            dataCheck = False
            for w in self.weight:
                if w != 0:
                    dataCheck = True
                else:
                    dataCheck = False
            if dataCheck:
                py = 0.0
                for i in range(self.dataLength):
                    py_temp = 1
                    for j in range(self.dataLength):
                        if i == j:
                            pass
                        else:
                            py_temp = py_temp * self.items[j]
                    py = py + py_temp * self.data['y'][i] / self.weight[i]
                return py
            else:
                print("There is a same x!")
                return None
        else:
            return None

    def mode2(self, px):
        self.__preprocess(order=2, x=px)
        itemsProd = np.prod(self.items)
        itemsSum = 0
        for i in range(self.dataLength):
            itemsSum = itemsSum + self.data['y'][i] / (self.weight[i] * self.items[i])
        py = itemsProd * itemsSum
        return py

    def mode3(self, px):
        self.__preprocess(order=2, x=px)
        denomSum = 0
        numeSum = 0
        for i in range(self.dataLength):
            dtemp = self.weight[i] * self.items[i]
            numeSum = numeSum + self.data['y'][i] / dtemp
            denomSum = denomSum + 1 / dtemp
        py = numeSum / denomSum
        return py

demo.py :

from lagrange_interpolation import LagrangeInterpolation as lag


x = [1, 2, 3, 4]
y = [1, 2, 3, 4]

inter1 = lag(x, y)
inter1.data_append((5, 5)) #往數據集中追加一個點
z1 = inter1.mode1(6)
z2 = inter1.mode2(6)
z3 = inter1.mode3(6)
print(z1)
print(z2)
print(z3)
"""
Result:
6.0
6.000000000000002
6.0000000000000036
"""

您的分享是我們最大的動力!

-Advertisement-
Play Games
更多相關文章
  • 簡述 類型:結構型 目的:降低對象創建時大量屬性也隨之被新建而帶來的性能上的消耗 話不多說,我們看一個案例。 優化案例 最初版v0 現在需要採購一批辦公用的電腦,以下是Computer類的定義。 class Computer { private String sn; // 序列號,電腦的唯一識別碼 ...
  • 探索密碼學的奇妙之旅。介紹分組密碼常用模式CFB密文反饋模式的相關理論。並基於AES標準,使用golang crypto包的cipher模塊實現了加密、解密字元串的過程。 ...
  • 七牛雲文件上傳 @RequestMapping("/upload") public Result upload(MultipartFile imgFile) { try { //獲取原始文件名 String originalFilename = imgFile.getOriginalFilename ...
  • Stream流呢,以前我也有所瞭解,像一些面試題中也出現過,Java8的新特性,有一塊就是這個Stream操作集合,而且在看一些項目中也使用的比較多。但總感覺自己學的一知半解,所以今天打算系統的過一下,再鞏固鞏固。 ###概念 Stream是JDK8 API中的新成員,它允許以聲明性方式處理集合。 ...
  • 一、開發背景 您好,我是 @馬哥python說 ,這是我獨立開發的Python可視化大屏,看下演示效果: 截圖: 視頻演示效果: https://www.zhihu.com/zvideo/1556218745923821568 這個大屏,是通過pyecharts可視化開發框架實現。 下麵詳細介紹,這 ...
  • 來源:https://www.linuxmi.com/50-million-pc-linux.html 開源社區的一大勝利! 繼德國之後,中國現在想在 5000 萬台 PC 上拋棄 Windows 並運行 Linux! 如果您一直密切關註 Linux 新聞,您可能聽說過德國去年在超過 25000 台 ...
  • ###一、Scrapy 介紹 Scrapy是一個Python編寫的開源和協作的框架。起初是用於網路頁面抓取所設計的,使用它可以快速、簡單、可擴展的方式從網站中提取所需的數據。 Scrapy也是通用的網路爬蟲框架,爬蟲界的django(設計原則很像),可用於數據挖掘、監測和自動化測試、也可以應用在獲取 ...
  • 2022-09-25 首先,要安裝好虛擬環境,之後要切換到虛擬環境中,使用的命令 workon 創建好的虛擬環境的名稱 之後,創建一個Django項目使用的命令: django-admin startproject 項目名稱 進入到該項目的目錄下,創建一個子應用,使用的命令: python mana ...
一周排行
    -Advertisement-
    Play Games
  • 前言 插件化的需求主要源於對軟體架構靈活性的追求,特別是在開發大型、複雜或需要不斷更新的軟體系統時,插件化可以提高軟體系統的可擴展性、可定製性、隔離性、安全性、可維護性、模塊化、易於升級和更新以及支持第三方開發等方面的能力,從而滿足不斷變化的業務需求和技術挑戰。 一、插件化探索 在WPF中我們想要開 ...
  • 歡迎ReaLTaiizor是一個用戶友好的、以設計為中心的.NET WinForms項目控制項庫,包含廣泛的組件。您可以使用不同的主題選項對項目進行個性化設置,並自定義用戶控制項,以使您的應用程式更加專業。 項目地址:https://github.com/Taiizor/ReaLTaiizor 步驟1: ...
  • EDP是一套集組織架構,許可權框架【功能許可權,操作許可權,數據訪問許可權,WebApi許可權】,自動化日誌,動態Interface,WebApi管理等基礎功能於一體的,基於.net的企業應用開發框架。通過友好的編碼方式實現數據行、列許可權的管控。 ...
  • Channel 是乾什麼的 The System.Threading.Channels namespace provides a set of synchronization data structures for passing data between producers and consume ...
  • efcore如何優雅的實現按年分庫按月分表 介紹 本文ShardinfCore版本 本期主角: ShardingCore 一款ef-core下高性能、輕量級針對分表分庫讀寫分離的解決方案,具有零依賴、零學習成本、零業務代碼入侵適配 距離上次發文.net相關的已經有很久了,期間一直在從事java相關的 ...
  • 前言 Spacesniffer 是一個免費的文件掃描工具,通過使用樹狀圖可視化佈局,可以立即瞭解大文件夾的位置,幫助用戶處理找到這些文件夾 當前系統C盤空間 清理後系統C盤空間 下載 Spacesniffer 下載地址:https://spacesniffer.en.softonic.com/dow ...
  • EDP是一套集組織架構,許可權框架【功能許可權,操作許可權,數據訪問許可權,WebApi許可權】,自動化日誌,動態Interface,WebApi管理等基礎功能於一體的,基於.net的企業應用開發框架。通過友好的編碼方式實現數據行、列許可權的管控。 ...
  • 一、ReZero簡介 ReZero是一款.NET中間件 : 全網唯一開源界面操作就能生成API , 可以集成到任何.NET6+ API項目,無破壞性,也可讓非.NET用戶使用exe文件 免費開源:MIT最寬鬆協議 , 一直從事開源事業十年,一直堅持開源 1.1 純ReZero開發 適合.Net Co ...
  • 一:背景 1. 講故事 停了一個月沒有更新文章了,主要是忙於寫 C#內功修煉系列的PPT,現在基本上接近尾聲,可以回頭繼續更新這段時間分析dump的一些事故報告,有朋友微信上找到我,說他們的系統出現了大量的http超時,程式不響應處理了,讓我幫忙看下怎麼回事,dump也抓到了。 二:WinDbg分析 ...
  • 開始做項目管理了(本人3年java,來到這邊之後真沒想到...),天天開會溝通整理需求,他們講話的時候忙裡偷閑整理一下常用的方法,其實語言還是有共通性的,基本上看到方法名就大概能猜出來用法。出去打水的時候看到外面太陽好好,真想在外面坐著曬太陽,回來的時候好兄弟三年前送給我的鍵盤D鍵不靈了,在打"等待 ...