拉格朗日插值原理及實現（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; // 序列號，電腦的唯一識別碼 ...
  • 密碼學奇妙之旅、01 CFB密文反饋模式、AES標準、Golang代碼
    探索密碼學的奇妙之旅。介紹分組密碼常用模式CFB密文反饋模式的相關理論。並基於AES標準，使用golang crypto包的cipher模塊實現了加密、解密字元串的過程。 ...
  • Python中class內置方法__init__與__new__作用與區別探究
    七牛雲文件上傳 @RequestMapping("/upload") public Result upload(MultipartFile imgFile) { try { //獲取原始文件名 String originalFilename = imgFile.getOriginalFilename ...
  • Java8Stream流
    Stream流呢，以前我也有所瞭解，像一些面試題中也出現過，Java8的新特性，有一塊就是這個Stream操作集合，而且在看一些項目中也使用的比較多。但總感覺自己學的一知半解，所以今天打算系統的過一下，再鞏固鞏固。 ###概念 Stream是JDK8 API中的新成員，它允許以聲明性方式處理集合。 ...
  • 【可視化大屏教程】用Python開發智慧城市數據分析大屏！
    一、開發背景 您好，我是 @馬哥python說 ，這是我獨立開發的Python可視化大屏，看下演示效果： 截圖： 視頻演示效果： https://www.zhihu.com/zvideo/1556218745923821568 這個大屏，是通過pyecharts可視化開發框架實現。 下麵詳細介紹，這 ...
  • 棄用 Windows，政府機構 5000 萬臺電腦將替換為國產 Linux！
    來源：https://www.linuxmi.com/50-million-pc-linux.html 開源社區的一大勝利！ 繼德國之後，中國現在想在 5000 萬台 PC 上拋棄 Windows 並運行 Linux！ 如果您一直密切關註 Linux 新聞，您可能聽說過德國去年在超過 25000 台 ...
  • Scrapy框架
    ###一、Scrapy 介紹 Scrapy是一個Python編寫的開源和協作的框架。起初是用於網路頁面抓取所設計的，使用它可以快速、簡單、可擴展的方式從網站中提取所需的數據。 Scrapy也是通用的網路爬蟲框架,爬蟲界的django(設計原則很像)，可用於數據挖掘、監測和自動化測試、也可以應用在獲取 ...
  • 創建一個Django項目總結
    2022-09-25 首先，要安裝好虛擬環境，之後要切換到虛擬環境中，使用的命令 workon 創建好的虛擬環境的名稱 之後，創建一個Django項目使用的命令： django-admin startproject 項目名稱 進入到該項目的目錄下，創建一個子應用，使用的命令： python mana ...
一周排行
    -Advertisement-
    Play Games
  • 移動開發（一）：使用.NET MAUI開發第一個安卓APP
    移動開發（一）：使用.NET MAUI開發第一個安卓APP 對於工作多年的C#程式員來說，近來想嘗試開發一款安卓APP，考慮了很久最終選擇使用.NET MAUI這個微軟官方的框架來嘗試體驗開發安卓APP，畢竟是使用Visual Studio開發工具，使用起來也比較的順手，結合微軟官方的教程進行了安卓 ...
  • wpf ToggleButton選中效果和一個登錄界面
    前言 QuestPDF 是一個開源 .NET 庫，用於生成 PDF 文檔。使用了C# Fluent API方式可簡化開發、減少錯誤並提高工作效率。利用它可以輕鬆生成 PDF 報告、發票、導出文件等。 項目介紹 QuestPDF 是一個革命性的開源 .NET 庫，它徹底改變了我們生成 PDF 文檔的方 ...
  • 彈幕樹洞項目功能新增篇
    項目地址 項目後端地址： https://github.com/ZyPLJ/ZYTteeHole 項目前端頁面地址： ZyPLJ/TreeHoleVue (github.com) https://github.com/ZyPLJ/TreeHoleVue 目前項目測試訪問地址： http://tree ...
  • 第27篇 sqlserver2022詳細安裝步驟
    話不多說，直接開乾 一.下載 1.官方鏈接下載： https://www.microsoft.com/zh-cn/sql-server/sql-server-downloads 2.在下載目錄中找到下麵這個小的安裝包 SQL2022-SSEI-Dev.exe，運行開始下載SQL server； 二. ...
  • .NET 開源高性能 MQTT 類庫
    前言 隨著物聯網（IoT）技術的迅猛發展，MQTT（消息隊列遙測傳輸）協議憑藉其輕量級和高效性，已成為眾多物聯網應用的首選通信標準。 MQTTnet 作為一個高性能的 .NET 開源庫，為 .NET 平臺上的 MQTT 客戶端與伺服器開發提供了強大的支持。 本文將全面介紹 MQTTnet 的核心功能 ...
  • Serilog文檔翻譯系列（六） - 可用的接收器、增強器、格式化輸出
    Serilog支持多種接收器用於日誌存儲，增強器用於添加屬性，LogContext管理動態屬性，支持多種輸出格式包括純文本、JSON及ExpressionTemplate。還提供了自定義格式化選項，適用於不同需求。 ...
  • 警惕 Visual Studio 屬性求值副作用導致邏輯不符合預期
    目錄簡介獲取 HTML 文檔解析 HTML 文檔測試參考文章 簡介 動態內容網站使用 JavaScript 腳本動態檢索和渲染數據，爬取信息時需要模擬瀏覽器行為，否則獲取到的源碼基本是空的。 本文使用的爬取步驟如下： 使用 Selenium 獲取渲染後的 HTML 文檔 使用 HtmlAgility ...
  • [使用目前最新版]HybridCLR6.9.0+YooAsset2.2.4實現純C# Unity熱更新方案 (一)
    1.前言 什麼是熱更新 游戲或者軟體更新時，無需重新下載客戶端進行安裝，而是在應用程式啟動的情況下，在內部進行資源或者代碼更新 Unity目前常用熱更新解決方案 HybridCLR，Xlua，ILRuntime等 Unity目前常用資源管理解決方案 AssetBundles，Addressable， ...
  • 在 ASP.NET Core Web API 中使用操作篩選器統一處理通用操作
    本文章主要是在C# ASP.NET Core Web API框架實現向手機發送驗證碼簡訊功能。這裡我選擇是一個互億無線簡訊驗證碼平臺，其實像阿裡雲，騰訊雲上面也可以。 首先我們先去 互億無線 https://www.ihuyi.com/api/sms.html 去註冊一個賬號 註冊完成賬號後，它會送 ...
  • 第28篇 如何.net中實現高效可靠數據同步api
    通過以下方式可以高效，並保證數據同步的可靠性 1.API設計 使用RESTful設計，確保API端點明確，並使用適當的HTTP方法（如POST用於創建，PUT用於更新）。 設計清晰的請求和響應模型，以確保客戶端能夠理解預期格式。 2.數據驗證 在伺服器端進行嚴格的數據驗證，確保接收到的數據符合預期格 ...
所有分類