拉格朗日插值及python实现

1. 拉格朗日插值及python实现

1.1. 返回首页

首页地址https://kz16.top/na
本文来源https://kz16.top/na/interp

1.2. 拉格朗日插值理论

1.2.1. 问题提出

求经过点列 $\{(x_i, y_i)\}_{i=0}^n(x_k\neq x_j,k\neq j)$ 的 $n$ 次多项式 $L_n(x)$ .

1.2.2. 问题分析

设 $L_n(x)=a_0+a_1x+\cdots + a_n x^n$ ，根据 $L_n(x_i)=y_i$ 得线性方程组
$\begin{cases} a_0+a_1x_0 +\cdots +a_n x_0^n=y_1\\ a_0+a_1x_1 +\cdots + a_n x_1^n=y_1\\ \cdots\\ a_0+a_1x_n +\cdots + a_n x_n^n = y_n \end{cases}$

上式线性方程组的左端系数矩阵的行列式是一个Vandermonde行列式，由于 $(x_k\neq x_j,k\neq j)$ ，即该行列式不为零，上述线性方程组存在唯一解，从而 $L_n(x)$ 存在唯一。

1.2.3. 拉格朗日插值基函数

考虑点列 $P_k = \{(x_i, \delta _{ik})\}_{i=0}^n, k=0,1,\cdots,n$ ，其中 $\delta _{ik}$ 满足

$\delta _{ik} = \begin{cases} 1,\quad i=k\\ 0,\quad i\neq k \end{cases}.$

记经过点列 $P_k$ 的 $n$ 次多项式为 $L^k_n(x)$ ，同时记
$Y_n(x) = y_0L_n^0(x) + y_1L_n^1(x) +\cdots + y_nL_n^n(x) .$

则多项式 $Y_n(x)$ 满足 $Y_n(x_i)=y_i, i=0,1,\cdots,n$ . 因此
$Y_n(x) = L_n(x).$

于是问题转化为求经过点列 $P_k$ 的多项式 $L^k_n(x)$ ，我们称 $L^k_n(x)$ 为拉格朗日插值基函数，它满足 $L^k_n(x_i)=\delta _{ik}$ ，所以多项式 $L^k_n(x)$ 的零点分别为 $x_0,x_1,\cdots,x_{k-1},x_{k+1},\cdots x_n.$

所以可以写出拉格朗日插值基函数 $L^k_n(x)$ 的具体表达式

$L^k_n(x) = \frac{(x - x_0)(x - x_1)\cdots (x-x_{k-1})(x-x_{k+1})\cdots (x-x_{n})}{(x_k - x_0)(x_k - x_1)\cdots (x_k -x_{k-1})(x_k -x_{k+1})\cdots (x_k -x_{n})}$

1.2.4. 拉格朗日插值基函数算法分析

要计算出所有的 $L^k_n(x),k=0,1,\cdots,n$ ，也就是要分别需要考虑它的分子和分母。很明显，它的分子的因子都在集合 $\{x-x_i\}_{i=0}^n$ 里，它的分母的因子都在集合 $\{x_i-x_j\}_{i,j=0}^n$ 中，因此算且仅算一次 $\{x-x_i\}_{i=0}^n$ 和 $\{x_i-x_j\}_{i,j=0}^n$ ，由于公式中涉及到连乘，因此在算法实现中，只需要把集合 $\{x-x_i\}_{i=0}^n$ 中某个多余的项设定为 $1$ ，同时把集合 $\{x_i-x_j\}_{i,j=0}^n$ 中的 $0$ 元素也设定为 $1$ 。按照这个思路可以很高效用矩阵运算实现该算法。类似可以利用该思路快速计算插值函数的导函数。

1.3. 拉格朗日插值算例

1.3.1. 问题

绘制函数 $f(x)=\dfrac{1}{1+25x^2}$ 在区间 $x\in[-1,1]$ 上的等距节点的插值，多项式次数分别取 $4,6,8,10$ .

UEsDBBQACAgIALSheE8AAAAAAAAAAAAAAAAXAAAAZ2VvZ2VicmFfZGVmYXVsdHMyZC54bWztmltz4yYUgJ+7v4LhqX2ILWTL9mai7GR3ptPMZLOZJrPTVyxhmQaBKlAs+9cvAlmS13bqW2onTR6CDuKm7xwOB/DFpzxm4ImkkgruQ9RyICA8ECHlkQ8zNTobwE+XHy4iIiIyTDEYiTTGyodeUbKqp6WW5/aKPJwkPgwYlpIGECQMq6KKDycQgFzScy5ucUxkggNyH4xJjG9EgJVpZaxUct5uTyaT1ry/lkijdhSpVi5DCPRYufRh+XCum1uoNOmY4q7joPZfX29s82eUS4V5QCDQ3xGSEc6YkvqRMBITroCaJsSHiaBcQcDwkDAf3hUS+HWUEvIbBGUljceBlx9+uZBjMQFi+DcJdJ5KM1LVM0K7KKNffxFMpCD1Yb8PQWSToQ9dz9OUWDLGPnRsYYanJAVPmFU5OFMiMPVN7ggzSeZldU9fRUjsm25ZntPYUARSEa0ABIFMCAnNk/1CZLQxNYpttBcIkYYS5D68xbcQTMt0ZlNTxLC5p7OyS6+Zq6aMNEZ+0S6xbgY4JAnhoS60QBntRLk3MJSLZGiT1wy5+9KQe++Q10FG21P+xpts3Z3YItczcE367igadK/5nyTSY24y7rwzPijjRQvu7kTXMWydV0rWFLEMZfFfhzEiThjJDwieUV5DvDFCBd3dLbpoQneOgtzZGXmBw8JTYxo8ciJlQbZut3j4g4Z69TL9CR02UqVbQv2BbYH8wxdURrXGqC7zvBpGGQ+UcSgl2i9Z+tTURafrHEMbdZuHVsa+pNezlCQqpIrL/VyuDXu3gO7/bdgiU6zo+Zorvcsixlzl0qc9EpI86Ka+8YcUc1lstRYtab3eUjx9Tmfeu85OT2dzv3X7HaeVJjId3I/02MOm+nYLj9Yu4C3XO7YOt/DmK4nsH9KclEFvaa2HMavebl7BdbqrMbb6J2xWT/rzRM3jeynWMcJ7vLa9i1wRZONUEUkx/7ctC5tGjTl+N5crffStPvYf49abSq9jdOqhJftGjv1D3Y8OQj3kHlvNzwNe2J7cVRk1YnQkxCc6adbTDAQvzr/n2wsrVRy7b8x1HGAfRyPCrceVAOSOKTZ1TOWZU95M5MjIU2TezpDNNvX1wFOagytb48oWvHJt0rFJ1yZeBWi3zaNRbaK9ViN8/mlp6O6243lNjuRNKv0/COF5FpO04Rpu53JlPJ51Drq9jCyodgNXsM5O1luFZDTUJhRTraQzrb0Y50aLeCgFyxS5D1JCeH1JZ814QkM1LgI73feI5oW52DbBWKR0JriqaIBiFlwxc523cMKxynzc5wLYBWPdzz1jHrF6Nl5ZqdaAPcA3hX4+21ulmCZDp0TYa7mDDhp4HaeP+h+9QW9DpGhQI7UvNia64G5KdWywniBnYzPa391s5TTcVU4Dp0F9CttxDmwYS3vJ36uMehd0ikeDxmSWir7YqR8TQSbrk2wrVYQGbyy+wVlOGcXpdLmnFyOsSF5HGA9GaPwQ4QQBr/8UjT2qh3ZtpcZ9v/2YEdUUOY51BdsJ5Z9x8BilIuPh8rp1kE9Hx7at9dCGQjCCa0f0eS437pmXIoV1gDZfDV5s9gVjEjwORb6wuD3vY6isZ8CNERr3vytmwD5r3tnRTWGXM71NryXPlk4A242fQrXnP7e6/AFQSwcI64/wTNEEAAAQJgAAUEsDBBQACAgIALSheE8AAAAAAAAAAAAAAAAXAAAAZ2VvZ2VicmFfZGVmYXVsdHMzZC54bWztmM1u2zgQx8/tUxC81yZlSYmCKIXRHnaBtkiRy14ZamxzK5EqScdWXm3fYZ+pI1Jx5G5cNEYaoO364OHXDMnfnxqLPn+9bWpyA9Ypo0vKJ4wS0NJUSi9LuvaLV6f09cXL8yWYJVxbQRbGNsKXNOtH7vywNsmSvG8TbVtSWQvnlKSkrYXvXUq6oYRsnTrT5oNowLVCwpVcQSPeGSl8iLLyvj2bTjebzeRuvomxy+ly6SdbV1GCa9WupEPhDMPtOW1mYXjCGJ/+9f5dDP9KaeeFlkAJ7qOChVjX3mERamhAe+K7FnDFRis5wzlqcQ11Sf/UHjcHsl8ZkWt7g/6Dc0lnPGP04uWLc7cyG2Ku/8ZxJfV2DTv/UJn2Y7D7jamNJbakCacEuXKG9hptkSCwul2JkrIJZ/HD04JxnvMk+teiA0tuBAZlsUWsvZEhZGhdiNrB3Vic/L2pIPakw3itmoCYOA+oDk7uWoAqlOL2WZCqC6qP4ykNV76rgfiVkp80OMSfjZz6wh+qqqA/PNEH1BL0DRIx1qHkLMzSsTD8lg1nbMtDveOh95bH5uCPS7VqS+bRYx4HzpNoZtGk0WQ7JPBZx3W6/rukrbB4yjCQ7PvPp4PY/5FdbJUbqT7vq2/3lGazo5RmQWgWZGb3Iv+kkh6mS4Yy4J7//efbsMNjJIX14JTQI+xv+o6vuee/O/fDIDG+hhG/y1Df44dJ8Ch+RREAJrwICIPdZajsqTBKY2zlyDYmgZgawvdmF3Ih+t+gYZaDufEhqOxIqKbuVlBZo++5jpru0c4GtMc8SY+Vg2ezoEfGvz7Rk3RAkhU5S/P0ybQ59og/iuzcypVqoAKxjxaFfS60CY8/xulJQNubX4PtZYcZWVX7XJ/vyIaUgYsvItfklzmzl1a5Zp8qf0aqeUzMkWqR/5RUNfjdPj/05XFWzf7Pqo9h+XktqvAGNmz14119zJQfeU05nBrztOg/JznPTnma8KcC9COuGg9eNPrGeJvoorlNdgEfe/cg8zyak2hOoykO3ktU09ZKKv9tad3aLvB+/NCr8tC1r3J6nMro9+DL8uTke4/9feBneV3m3/tmNx3d86d3/yVcfAFQSwcICiLpwzwDAADtEAAAUEsDBBQACAgIALWheE8AAAAAAAAAAAAAAAAWAAAAZ2VvZ2VicmFfamF2YXNjcmlwdC5qc0srzUsuyczPU0hPT/LP88zLLNHQVKiuBQBQSwcI1je9uRkAAAAXAAAAUEsDBBQACAgIALWheE8AAAAAAAAAAAAAAAAMAAAAZ2VvZ2VicmEueG1svVhtb9s4Ev7c/RWE0DtkbxObpKgX9+wu0t4Wt0C7V1z2DovD4gBaom1eZMkr0ondbv/7zpCUbMdJmzTttZEpksMZzjPDmaHG32+WFblSrdFNPYnYgEZE1UVT6no+idZ2dpZH3z//ZjxXzVxNW0lmTbuUdhIlSNmvg94g4SmOydVqEhWVNEYXEVlV0uKSSXQdEV1Oold5loj0B352njNxJuJXfzsbnb98dcZGyYtU/DAavXh5HhGyMfpZ3fwkl8qsZKEuioVaytdNIa2Tt7B29Ww4vL6+HnQ7GzTtfDifTwcbU0YEtKrNJAovz4DdwaLr2JFzStnwlzevPfszXRsr60JFBDVe6+ffPBlf67psrsm1Lu1iEqWjPCILpecLhICCukMkWgEOK1VYfaUMLN3rOp3tchU5Mlnj/BP/RqpenYiU+kqXqp1EdMCTOMlZmrJ8xBjLRBaRptWqtoGYBaHDjt34SqtrzxffnMiE8YjYpqmmEpmS3wkjCYWHsBE5JWkGI5ywhAgYyWEkIzGOJUyQmCAJi4kQ0AocZinM4DT8JgkljMEM4ZRwTjgjPIZukpAEyDJcy4E2HTl+FB6khh3BE+NYHMPjxmIBD8c3YJR4NrCPJE7dW4LUwD/hqIEbjHMiRiAIB5KMkRj2AP2MEuAYI3vm9BCU4B8jAtnzjPCcAD9QHTlTQOdKGz2t1CSaycqAzXU9a8Hf+r6x20o5+MLAzl7sFP4DhX6nght4N4EZSk/xSeERnX/s2UUcWgWMQEG3U2yYb7gfpb5LY99w3wjfJJ5G+JXCk3pFqfA0In6shp1+8UP0y/f0Y6gE2AN375qY4L6Z2z82InRT33WORhkNo7kfHWE3faQy8WcpI+goe4hg267vBDHmx3IzeOA0Hclle0I9z4co24l0seq+qj7oPNzUs0eXZdn9RbJHGvVWkyaAKP6550hk/Cg1P0dienDgv4zCIr+3eMbz/7vMjN4a43zLQvtlDDG6vyEeGw97IJKPixwPu2w8DiAQs0DacJCtWhqEJYtJyvu0nGLWDLk54yRLSJbuZehTzNFpskvTmKTzgzSd5Ie5OsXBzCV+SI2YZn3S5qLL26chc/9+lLkh0YpdroUNIitGCNQGJMWwHJIu7IL3aZcnmHl5SiA1J5ykGPrvyMBQFjZG98AuVLXqIXcY6nq1tge4Fcuye7UNUMvKFX2BvmyKyxc90oGTksbus4WCaVeX+QLqoGx7Mq7kVFVQ+F6gHxByJSs8s07CrKktCS7AReTYuQpxrNZFpUst63+D3btq7Kf1cqpa4l4bVNIxweXk1lJSiMSTFE3TlhdbA25CNv9RLSzmOVTlYsSFSKGGgpQEnrwNU2k6iAVUUTSLORMxBBtTyMpFYj5I85xlac5illIGgXl7x1wmvGx1daGsBf0NkRtleuzmLR6gvc6P5kVT7YZWja7tS7my69bdGmB/Lap1Xs8r5bB0ZoYau7icNpsLD2Lqef28XWFE9TuYzl82VdMSOIE8SYAgtFPfOhrcWk9FHQ11FMFSyLSfZyPuKFw79a2jAjP7rQVVGev0pJ0YbYjvH7qVcxIs2de1tq+7jtXF5U5VXOBdwASvPeTJvhTP8fCG+43DweiccdmUyjty7OkP5seXqq1V5R2vBtOvm7Xx5H5nbttro95Kuzivy3+qORzatxKjpoWNeNKdgqUq9BIW+vEAtUQ3+Bco5kdLNW9VB4jfjDdE2CUxq1bJ0iyUsr01/LHYkVGvTLf9sZUQ112sX2oIKmdg7aXcuLIGjtIqnMCxKVq9Qg8nUwjtl2rnw6U2yKLcUxwhMaBbgYEKjGHREHC/XdtF07o7m7Q4gjFgA3s2eBXuTDmD+LPBqHay+ZZMIMoPyQkj35GTE4iTfyEwiD/fgpGd/8/WtRPj5KpKLeG6d3Nqx9rFEnAH0kz/B+GrT5BhPoS/K7xtOPSA7NYT486WrFYL2QNdyS3Grb3I57i+OXQjiId6o8qb5t8dKrsA34XbKbhq0ulI/cvfdVkqH/gbuNpru8WIlAc7QxRBFB3jp79Wsp2r93+qPzwNzu6ROcKoXi9Vi58bAgS1gwiYrMOuAxBmu5w2lS72Nu4TyN1o7rS7C0W6w3AA9a77x8SIMpYy/lFU/zGbGWUJ+OkZRO6ta27DvEsfFX4i8B4uOgdH6VPTVGurLgo4WPXuI4nXI6QbjgUKLGB57kSlgL8zYu8w4NX6HUSBHTSIyrk/op80M6N32jlEAYgMzrbdgYSDrpQPhv3CFSDlcshe+Lm3TxTNcinrktSuprtQv60VfszZFRSSAsQgTjInVXKABQ5mHey7th2R3njegeORv1Xa2N5DgPh292D0/g7ysJPjPzph2sVk6pFJ9kd3kfKMfcL1PxPeV9q+bartDXRP9OaUzOAXQttNUOcfx/Qozs0fFee6+uHIBF89zrHkiwa6o9TyNqSWfazJcd448NG3n8DyY1HuFjd+GIZf37kfgaZVG8sCon/+bd3Yvz7FnD35Fa6GxXv24T37jieb//IPT/3sLRkaWUSH/O6Btndc+tUytDav5c/ql5u+6y4yBrLlrI/ygPKbvRsuDRfhnr+xsoXDrvEGhLmKDpI0y3JO45Ql6SjlPqNAcIPrCRfpKOO5yFicwPi7rvTaM8Jwv7ByF6nwCf/5H1BLBwgGbNBezAcAAI4YAABQSwECFAAUAAgICAC0oXhP64/wTNEEAAAQJgAAFwAAAAAAAAAAAAAAAAAAAAAAZ2VvZ2VicmFfZGVmYXVsdHMyZC54bWxQSwECFAAUAAgICAC0oXhPCiLpwzwDAADtEAAAFwAAAAAAAAAAAAAAAAAWBQAAZ2VvZ2VicmFfZGVmYXVsdHMzZC54bWxQSwECFAAUAAgICAC1oXhP1je9uRkAAAAXAAAAFgAAAAAAAAAAAAAAAACXCAAAZ2VvZ2VicmFfamF2YXNjcmlwdC5qc1BLAQIUABQACAgIALWheE8GbNBezAcAAI4YAAAMAAAAAAAAAAAAAAAAAPQIAABnZW9nZWJyYS54bWxQSwUGAAAAAAQABAAIAQAA+hAAAAAA

温馨提示： 可以用鼠标拖动滑动条，查看不同次数的插值，也可以点击滑动条，然后按空格键启动自动动画。

1.3.2. GeoGebra代码

f(x)=1/(1+25*x*x)
n = slider(4,20,1)
ix = Sequence(-1, 1, 2/n)
P = (ix, f(ix))
fitpoly(P)

不妨试试上面的代码逐行指令复制粘贴到下面的命令行中，

不妨再试试非等距节点插值, GeoGebra代码代码如下

f(x)=1/(1+25*x*x)
n = slider(4,20,1)
iix = Sequence(-1, 1, 2/n)
ix = sin(pi/2*iix)
P = (ix, f(ix))
fitpoly(P)

1.3.3. GeoGebra简介

GeoGebra（音译：几何不难）
是数形结合的神器，个人使用完全免费，并且可跨平台使用

官网https://www.geogebra.org/
经典版本下载 http://download.geogebra.org/installers/5.0/
跨平台的网页版说明
https://wiki.geogebra.org/en/Reference:GeoGebra_Apps_Embedding
跨平台的网页版下载
https://download.geogebra.org/package/geogebra-math-apps-bundle

1.3.4. python代码1

使用numpy.polyfit函数

import numpy as np
import matplotlib.pyplot as plt

f = lambda x:1/(1+25*x**2)
xx = np.linspace(-1,1,200)
fig = plt.figure(figsize=(8,8), dpi=72,facecolor="w")
axes = plt.subplot(111)
plt.plot(xx, f(xx), '--', color="r")
for i in [5, 7, 9, 11]:
    x = np.linspace(-1,1,i)
    y = np.poly1d(np.polyfit(x, f(x), len(x)))(xx)
    plt.plot(xx, y)
plt.show()

1.3.5. python代码2

根据拉格朗日插值公式实现，这种方法很实用，稳定性相对前者较好，而且给出可求插值基函数的梯度的方法InterpBase.gradValues()，可用于有限元编程

import numpy as np

class InterpBase(object):
    def __init__(self, ip):
        self.ip = np.array(ip, float)
        a, b = np.meshgrid(ip, ip) 
        c  = b - a
        ran = range(len(ip))
        c[ran, ran] = 1
        self.c = self.abssort(c)
        
    def mapfx(self, f, x):
        return np.array(list(map(f, x)))

    def abssort(self, y):
        f = lambda x: x[np.abs(x).argsort()]
        return self.mapfx(f, y)

    def values(self, ix):
        z = ix - self.ip
        n = len(self.ip)
        ran = range(n)
        zs = np.zeros((n, n))
        zs[ran] = z
        zs[ran, ran] = 1
        return (self.abssort(zs)/self.c).prod(axis=1)

    def valuesm(self, ixm):
        return self.mapfx(self.values, ixm)
 
    def gradValues(self, ix):
        z = ix - self.ip
        n = len(self.ip)
        ran = range(n)
        zs = np.zeros((n, n, n))
        zs[ran] = z
        zs[:, ran, ran] = 1
        zs[ran, :, ran] = 1
        zs[ran, ran, ran] = 0
        ans = zs.prod(axis=2).sum(axis=0)
        return ans/self.c.prod(axis=1)
        
    def gradValuesm(self, ixm):
        return self.mapfx(self.gradValues, ixm)
        
f = lambda x:1/(1+25*x**2)
xx = np.linspace(-1,1,200)

import matplotlib.pyplot as plt
fig = plt.figure(figsize=(8,8), dpi=72,facecolor="w")
axes = plt.subplot(111)
plt.plot(xx, f(xx), '--', color="r")
for i in [5, 7, 9, 11]:
    x = np.linspace(-1,1,i)
    plt.plot(xx, InterpBase(x).valuesm(xx).dot(f(x)))
plt.show()

上面插值部分的代码也可以写得简单点

class InterpBase(object):
    def __init__(self, ip):
        self.ip = np.array(ip, float)
        self.c  = np.array([ip]).T - self.ip
        ran = range(len(ip))
        self.c[ran, ran] = 1

    def valuesm(self, ix):
        z = np.array([ix]).T - self.ip
        n = len(self.ip)
        m = np.array([ix]).shape[-1]
        ran = range(n)
        zs = np.zeros((n, m, n))
        zs[ran] = z
        zs[ran, :, ran] = 1
        tmp = np.einsum('ijk->jik', zs, optimize=True)
        return (tmp/self.c).prod(axis=2)

或者直接得到两个矩阵

import numpy as np
def interpBase(ip, iy, pureComputeH=True):
    c  = ip[:, None] - ip
    n = len(ip)
    ran = np.arange(n)
    c[ran, ran] = 1
    z = np.array([iy]).T - ip
    zs = np.einsum('ij,k->ikj', z, np.ones(n), optimize=True)
    zs[:, ran, ran] = 1
    matrixH = (zs/c).prod(axis=-1)
    if pureComputeH:
        return matrixH
    dzs = np.einsum('ij,k->ikj', c, np.ones(n), optimize=True)
    dzs[:, ran, ran] = 1
    dzs[ran, ran, :] = 1
    matrixD = (dzs/c).prod(axis=-1)
    matrixD[ran, ran] = 0
    matrixD[ran, ran] = -matrixD.sum(axis=-1)
    return matrixH, matrixH.dot(matrixD)

注意: 更快速的代码, 可参考这里

1.4. 返回首页

首页地址https://kz16.top/na
本文来源本文来源https://kz16.top/na/interp

仅供参考，部分漏洞在所难免
欢迎转载，转载请指明来源，
请勿用于商业
作者：Zhou Bingzhen
2019.11.24