设函数$y=f(x)$在区间$[a,b]$上有定义，且在点$a⩽x_0⩽x_1\dots ⩽x_n⩽b$上的值$y_0,y_1\dots y_n$之间存在一个函数$P(x)$,使得

$$P(x_i)=y_i(i=0,1,\dots ,n)$$

成立，则称$P(x)$为$f(x)$的插值函数，在$x_0,x_1,\dots ,x_n$称为插值节点，包含节点的区间称为插值区间。如果$P(x)$是次数不超过$n$的代数多项式，则称之为插值多项式。

数学原理

若有$n+1$组一一对应的$(x,y)$点对，假设存在一多项式$L_n(x)$，使之满足$L_n(x_j)=y_j,j=0,1,\dots ,n$，即

$$\{\begin{array}{rlr}{a}_0+a_1{x}_0+\dots +a_n{x}_0^n& =& y_0\\ a_0+a_1{x}_1+\dots +a_n{x}_1^n& =& y_1\\ & ⋮& \\ a_0+a_1{x}_n+\dots +a_n{x}_n^n& =& y_n\end{array}$$

则$L_n$即为该组数据的插值多项式。对于任意$j=0,1,\dots ,n$，有$L_n(x_j)=y_j$,则可令

$$L_n(x)=\sum _{k=0}{y}_k{l}_k(x)$$

其中，

$$l_k(x_j)={\delta }_{jk}=\{\begin{array}{r}1,k=j\\ 0,k\ne j\end{array}(j,k=0,1,\dots ,n)$$

则可满足$L_n(x_j)=\sum _{k=0}^n{y}_k{l}_k(x_j)=y_j$，其中$l_k$被称为基函数。根据其所满足的要求，易得基函数的表达式为

$$l_k(x)=\frac{(x-x_0)\dots (x-x_{k-1})(x-x_{k+1})\dots (x-x_n)}{(x_k-x_0)\dots (x_k-x_{k-1})(x_k-x_{k+1})\dots (x_k-x_n)}$$

自此，就得到了Lagrange插值多项式。该式简洁明了，其插值基函数规律明显，能使人过目不忘。若记

$${\omega }_{n+1}(x)=(x-x_0)(x-x_1)\dots (x-x_n)$$

则

$${\omega }_{n+1}^{\prime }(x_k)=(x_k-x_0)\dots (x_k-x_{k-1})(x_k-x_{k+1})\dots (x_k-x_n)$$

那么Lagrange插值公式可写为

$$L_n(x)=\sum _{k=0}^n{y}_k\frac{{\omega }_{n+1}(x)}{(x-x_k){\omega }_{n+1}^{\prime }(x_k)}$$

特别地，对于两点而言，Lagrange插值退化为线性插值，即

$$L_2(x)=y_1\frac{x-x_2}{x_1-x_2}+y_2\frac{x-x_1}{x_2-x_1}$$

算法化

在Lagrange插值公式中，${\omega }_{n+1}(x)$为多项式乘法。Julia中虽然对函数式编程有着良好的支持，但如果将函数本身作为输出，则函数实现细节将被隐藏，因此并不能得到真正的拟合参数。故需对多项式进行抽象，即将其简化为数组的形式，其表现形式为[a_0,a_1,...,a_n],代表$a_0+a_{1}x+\dots +a_n{x}^n$。

那么多项式乘法可以通过进位的方式实现。对于常数项而言，并不会改变多项式的幂数；对于一阶项来说，乘以多项式数组，则将该数组向高位移一位。其实现方法为

# 多项式乘法，形式为an*x^n，建议a[1]为最低位常数项
function polyMulti(a,b)
    na = length(a); nb = length(b)
    c = zeros(na+nb-1)
    for i = 1:na
        c[i:i+nb-1] += a[i]*b
    end
    return c
end

Lagrange插值公式亦涉及到多项式除法，其实现方法与乘法类似，但运算顺序要从高位开始，并且被除多项式的阶数要不小于除数多项式。

Julia对象采用地址引用，在传参之后如果直接改变数组内容，将穿透作用域，使得函数之外的数组也发生变化，故而需要深拷贝。

function polyDiv(a,b)
    a = copy(a)         #此为余数
    na = length(a); nb = length(b)
    nc = na-nb+1
    c = zeros(nc)
    for i = 0:nc-1
        c[nc-i] = a[end-i] ÷ b[end]
        a[nc-i:end-i] += -c[nc-i]*b[1:end]
    end
    return c, a
end

此外，Lagrange插值需要频繁地进行连乘计算，其中，分子中的连乘为多项式连乘，分母中的多项式为值连乘，所以可建立连乘函数。

# 连乘函数，输入为数组，输出为数组中所有数的乘积
function sumMulti(x)
    return length(x) == 1 ? x : x[1]*sumMulti(x[2:end])
end

Lagrange插值公式中的多项式连乘形式单一，相比之下更像是多项式的两种表示形式的转换，所以再建立针对Lagrange插值公式的多项式转化函数

#将(x-x1)(x-x2)...(x-xn)化成an*x^n形式,A[1]为最低位(常数项)
#输入为[x1,x2,...,xn]输出为[a1,a2,...,an]
function polySimplify(x)
    if length(x) == 1
        return [-x[1],1]
    else
        return polyMulti([-x[1],1],polySimplify(x[2:end]))
    end
end

至此，可以创建Lagrange插值函数。

# Language插值函数
# x y: 维度相同的列变量
# A: 多项式系数，A[1]为最低位
using LinearAlgebra         #建立单位矩阵需要使用线性代数包
function Lagrange(x,y)
    n = length(x)
    xMat = x .- x' + Matrix{Float64}(I,n,n)
    polyMat = polySimplify(x)
    A = zeros(n)
    for i = 1:n
        quot = polyDiv(polyMat,[-x[i],1])[1]
        A += y[i] * quot ./ sumMulti(xMat[i,:])
    end 
    return A
end

函数中的xMat为矩阵

$$\left[\begin{array}{cccc}1& x_1-x_0& \dots & x_n-x_0\\ x_0-x_1& 1& \dots & x_n-x_1\\ & & ⋮& \\ x_0-x_n& x_1-x_n& \dots & 0\end{array}\right]$$

polyMat为多项式${\omega }_{n+1}(x)$。

测试

在命令行中调用得

julia> include("interpolation.jl")
Lagrange (generic function with 1 method)

julia> x = [i for i = 1:9];

julia> y = x.^8 + x.^6 + x.^2 .+ 1;

julia> a = Lagrange(x,y)
9-element Array{Float64,1}:
  1.0
  4.470348358154297e-8
  0.9999999850988388
 -1.4901161193847656e-8
 -7.450580596923828e-9
  2.7939677238464355e-9
  1.0
  0.0
  1.0

结果表明，其8次项、6次项、2次项和0次项分别为1或者接近于1，其他项分别为0或者接近于0，可见拟合结果尚可。通过Plots包对原始数据和拟合结果进行比较

julia> using Plots
julia> y1 = y;
julia> y2 = [sum([a[i]*x[j]^(i-1) for i=1:9]) for j = 1:9];
julia> plot(x,[y1,y2],st=[:scatter :line])
julia> savefig("Lagrange.png")

可见插值结果能够符合原始数据的趋势。