拉格朗日插值原理及實現（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
  • Dapr Outbox 執行流程
    Dapr Outbox 是1.12中的功能。 本文只介紹Dapr Outbox 執行流程，Dapr Outbox基本用法請閱讀官方文檔 。本文中appID=order-processor，topic=orders 本文前提知識：熟悉Dapr狀態管理、Dapr發佈訂閱和Outbox 模式。 Outbo ...
  • Advanced .Net Debugging 7：托管堆與垃圾收集
    引言 在前幾章我們深度講解了單元測試和集成測試的基礎知識，這一章我們來講解一下代碼覆蓋率,代碼覆蓋率是單元測試運行的度量值,覆蓋率通常以百分比表示，用於衡量代碼被測試覆蓋的程度，幫助開發人員評估測試用例的質量和代碼的健壯性。常見的覆蓋率包括語句覆蓋率（Line Coverage）、分支覆蓋率（Bra ...
  • C#S7.NET實現西門子PLCDB塊數據採集的完整步驟
    前言 本文介紹瞭如何使用S7.NET庫實現對西門子PLC DB塊數據的讀寫，記錄了使用電腦模擬，模擬PLC，自至完成測試的詳細流程，並重點介紹了在這個過程中的易錯點，供參考。 用到的軟體： 1.Windows環境下鏈路層網路訪問的行業標準工具（WinPcap_4_1_3.exe）下載鏈接：http ...
  • WPF控制項：密碼框綁定MVVM
    從依賴倒置原則（Dependency Inversion Principle, DIP）到控制反轉（Inversion of Control, IoC）再到依賴註入（Dependency Injection, DI）的演進過程，我們可以理解為一種逐步抽象和解耦的設計思想。這種思想在C#等面向對象的編 ...
  • Python中的私有屬性與私有方法
    關於Python中的私有屬性和私有方法 Python對於類的成員沒有嚴格的訪問控制限制，這與其他面相對對象語言有區別。關於私有屬性和私有方法，有如下要點： 1、通常我們約定，兩個下劃線開頭的屬性是私有的（private）。其他為公共的（public）； 2、類內部可以訪問私有屬性（方法）； 3、類外 ...
  • C++ 訪問說明符詳解：封裝數據，控制訪問，提升安全性
    C++ 訪問說明符 訪問說明符是 C++ 中控制類成員（屬性和方法）可訪問性的關鍵字。它們用於封裝類數據並保護其免受意外修改或濫用。 三種訪問說明符： public：允許從類外部的任何地方訪問成員。 private：僅允許在類內部訪問成員。 protected：允許在類內部及其派生類中訪問成員。 示 ...
  • static_cast與dynamic_cast到底是什麼？
    寫這個隨筆說一下C++的static_cast和dynamic_cast用在子類與父類的指針轉換時的一些事宜。首先，【static_cast，dynamic_cast】【父類指針，子類指針】，兩兩一組，共有4種組合：用 static_cast 父類轉子類、用 static_cast 子類轉父類、使用 ...
  • C語言數據結構：雙向鏈表的增刪操作
    /******************************************************************************************************** * * * 設計雙向鏈表的介面 * * * * Copyright (c) 2023-2 ...
  • @ComponentScan註解的實現，Spring掃描包的過程
    相信接觸過spring做開發的小伙伴們一定使用過@ComponentScan註解 @ComponentScan("com.wangm.lifecycle") public class AppConfig { } @ComponentScan指定basePackage，將包下的類按照一定規則註冊成Be ...
  • opensips開啟python支持
    操作系統 ：CentOS 7.6_x64 opensips版本： 2.4.9 python版本：2.7.5 python作為腳本語言，使用起來很方便，查了下opensips的文檔，支持使用python腳本寫邏輯代碼。今天整理下CentOS7環境下opensips2.4.9的python模塊筆記及使用 ...
所有分類