拉格朗日插值原理及實現(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
  • 前言 當別人做大數據用Java、Python的時候,我使用.NET做大數據、數據挖掘,這確實是值得一說的事。 寫的並不全面,但都是實際工作中的內容。 .NET在大數據項目中,可以做什麼? 寫腳本(使用控制台程式+頂級語句) 寫工具(使用Winform) 寫介面、寫服務 使用C#寫代碼的優點是什麼? ...
  • 前言 本文寫給想學C#的朋友,目的是以儘快的速度入門 C#好學嗎? 對於這個問題,我以前的回答是:好學!但仔細想想,不是這麼回事,對於新手來說,C#沒有那麼好學。 反而學Java還要容易一些,學Java Web就行了,就是SpringBoot那一套。 但是C#方向比較多,你是學控制台程式、WebAP ...
  • 某一日晚上上線,測試同學在回歸項目黃金流程時,有一個工單項目介面報JSF序列化錯誤,馬上升級對應的client包版本,編譯部署後錯誤消失。 線上問題是解決了,但是作為程式員要瞭解問題發生的原因和本質。但這都是為什麼呢? ...
  • 本文介紹基於Python語言中TensorFlow的Keras介面,實現深度神經網路回歸的方法。 1 寫在前面 前期一篇文章Python TensorFlow深度學習回歸代碼:DNNRegressor詳細介紹了基於TensorFlow tf.estimator介面的深度學習網路;而在TensorFl ...
  • 前段時間因業務需要完成了一個工作流組件的編碼工作。藉著這個機會跟大家分享一下整個創作過程,希望大家喜歡,組件暫且命名為"easyFlowable"。 接下來的文章我將從什麼是工作流、為什麼要自研這個工作流組件、架構設計三個維度跟大家來做個整體介紹。 ...
  • 1 簡介 我們之前使用了dapr的本地托管模式,但在生產中我們一般使用Kubernetes托管,本文介紹如何在GKE(GCP Kubernetes)安裝dapr。 相關文章: dapr本地托管的服務調用體驗與Java SDK的Spring Boot整合 dapr入門與本地托管模式嘗試 2 安裝GKE ...
  • 摘要:在jvm中有很多的參數可以進行設置,這樣可以讓jvm在各種環境中都能夠高效的運行。絕大部分的參數保持預設即可。 本文分享自華為雲社區《為什麼需要對jvm進行優化,jvm運行參數之標準參數》,作者:共飲一杯無。 我們為什麼要對jvm做優化? 在本地開發環境中我們很少會遇到需要對jvm進行優化的需 ...
  • 背景 我們的業務共使用11台(阿裡雲)伺服器,使用SpringcloudAlibaba構建微服務集群,共計60個微服務,全部註冊在同一個Nacos集群 流量轉發路徑: nginx->spring-gateway->業務微服務 使用的版本如下: spring-boot.version:2.2.5.RE ...
  • 基於php+webuploader的大文件分片上傳,帶進度條,支持斷點續傳(刷新、關閉頁面、重新上傳、網路中斷等情況)。文件上傳前先檢測該文件是否已上傳,如果已上傳提示“文件已存在”,如果未上傳則直接上傳。視頻上傳時會根據設定的參數(分片大小、分片數量)進行上傳,上傳過程中會在目標文件夾中生成一個臨 ...
  • 基於php大文件分片上傳至七牛雲,使用的是七牛雲js-sdk V2版本,引入js文件,配置簡單,可以暫停,暫停後支持斷點續傳(刷新、關閉頁面、重新上傳、網路中斷等情況),可以配置分片大小和分片數量,官方文檔https://developer.qiniu.com/kodo/6889/javascrip ...