拉格朗日插值原理及實作（Python）
目錄- 拉格朗日插值原理及實作（Python）
- 一. 前言
- 二. 3種形式的Lagrange插值函式推導
- 1. 原始形態的Lagrange插值
- 2. 第一形式Lagrange插值
- 3. 第二形式的Lagrange插值（重心插值公式）
- 三. 利用Python編程實作這三種Lagrange插值
一. 前言
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插值計算的流程如下：
- 預處理：根據已知資料點，利用公式(4)計算\(w_0,w_1,\cdots,w_n\);
- 插值：計算\(g(x)=(x-x_0)(x-x_1)\cdots(x-x_n)\)，然后根據公式(5)，先計算中括號內的加式，然后乘以\(g(x)\)得到插值點的值，這樣計算一個新的點的時間復雜度就是\(O(n)\)；
- 補充資料點：如果資料集有更新，只需要更新\(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插值采用這種形式，
它的計算程序如下：
- 預處理：計算\(w_i\)；
- 插值：對需要插入的點(x,y)，計算（x-x_i)（可以存到一個串列中，計算時直接取用）；
- 補充資料點：補充新的點到資料集只需要更新\(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 = https://www.cnblogs.com/yang-ding/archive/2022/09/25/{'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
"""
轉載請註明出處，本文鏈接：https://www.uj5u.com/houduan/509492.html
標籤：其他
上一篇：爬蟲高性能相關
下一篇：創建一個Django專案總結