拉格朗日插值基函数导出矩阵的快速计算

1. 问题

给定一组插值节点的横坐标{xi}i=1n\{x_i\}_{i=1}^{n}, 和赋值节点{yi}i=1m\{y_i\}_{i=1}^{\red{m}}, 需要求出定义在节点{xi}i=1n\{x_i\}_{i=1}^{n}上的n1n-1次拉格朗日基函数hi(x)h_{i}(x)在节点{yi}i=1m\{y_i\}_{i=1}^{\red{m}}处的值, 以及拉格朗日基函数的导数Li(x)L'_{i}(x)在节点{yi}i=1m\{y_i\}_{i=1}^{\red{m}}处的值. 也就是要求出矩阵H=[hj(yi)]i=1,mj=1,n\mathbb{H}=[h_{j}(y_i)]_{i=1,\red{m}}^{j=1,n}以及矩阵D=[hj(yi)]i=1,mj=1,n\mathbb{D}=[h'_{j}(y_i)]_{i=1,\red{m}}^{j=1,n}

1.1. 基本性质

D~=[hj(xi)]i=1,nj=1,n\widetilde{\mathbb{D}}=[h'_{j}(x_i)]_{i=1,n}^{j=1,n}, S~=[hj(xi)]i=1,nj=1,n\widetilde{\mathbb{S}}=[h''_{j}(x_i)]_{i=1,n}^{j=1,n} 以及 S=[hj(yi)]i=1,nj=1,n{\mathbb{S}}=[h''_{j}(y_i)]_{i=1,n}^{j=1,n}

证明: 根据多项式插值原理, 如果f(x)f(x)[1,1][-1,1]上的n1n-1多项式, 那么有

f(x)=k=1nf(xk)hk(x),f(x)=k=1nf(xk)hk(x),f(x)=\sum_{k=1}^n f(x_k) h_{k}(x), \quad f'(x)=\sum_{k=1}^n f(x_k) h'_{k}(x),

根据该基本性质, 实际计算, 只需计算出矩阵H\mathbb{H}D~\mathbb{\widetilde{D}}, 就可得到其它矩阵

1.2. 求解

1.2.1. 计算矩阵H

定义矩阵C=[cjk]j=1,nk=1,n\mathbb{C}=[c_{jk}]_{j=1,n}^{k=1,n}, 以及三维矩阵Z=[zijk]i=1:mj,k=1:n\mathbb{Z}=[z_{ijk}]_{i=1:\red{m}}^{j,k=1:n}, 其中

cjk={xjxk,jk1,j=k,以及 zijk={yixk,jk1,j=k.c_{jk}= \begin{cases} x_j-x_k, & \quad j\neq k\\ 1, & \quad j=k \end{cases},\quad \text{以及 } z_{ijk} = \begin{cases} y_i-x_k, & \quad j\neq k\\ 1, & \quad j=k \end{cases} .

根据公式

hj(yi)=(yix1)(yix2)(yixj1)(yixj+1)(yixn)(xjx1)(xjx2)(xjxj1)(xjxj+1)(xjxn)h_j(y_i)=\dfrac{(y_i-x_1)(y_i-x_2)\cdots(y_i-x_{j-1})(y_i-x_{j+1})\cdots(y_i-x_{n})}{(x_j-x_1)(x_j-x_2)\cdots(x_j-x_{j-1})(x_j-x_{j+1})\cdots(x_j-x_{n})}

可得

hj(yi)=zij1zij2zij3zijncj1cj2cj3cjnh_j(y_i)=\dfrac{z_{ij1}z_{ij2}z_{ij3}\cdots z_{ijn}}{c_{j1}c_{j2}c_{j3}\cdots c_{jn}}

Python的三维矩阵除以二维矩阵, 运算法则是这样的 Aijk/BjkA_{ijk}/B_{jk}, 因此可以有很简洁的表达式

H=(ZC)最后一个轴累乘\mathbb{H}=\left(\dfrac{\mathbb{Z}}{\mathbb{C}}\right)_{最后一个轴累乘}.

全部Python代码如下

import numpy as np
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, iy):
        z = np.array([iy]).T - self.ip
        n = len(self.ip)
        m = np.array([iy]).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)

注意:

事实上H\mathbb{H}矩阵乘以向量, 也就是多项式多点赋值问题, 可以有更快速的办法.

但是都比较复杂, 对于小规模问题, 感觉没什么必要, 可点我查看参考文献

1.2.2. 计算矩阵D

根据求导公式

hj(yi)=hj(yi)[1yix1+1yix2++1yixj1+1yixj+1+1yixn]h'_j(y_i)=h_j(y_i)\left[\frac{1}{y_i-x_1}+\frac{1}{y_i-x_2}+\cdots +\frac{1}{y_i-x_{j-1}}+\frac{1}{y_i-x_{j+1}}+\cdots \frac{1}{y_i-x_{n}}\right]

可分两种情况进行求解

  1. yi=xi,i=1,2,ny_i=x_i,\forall i=1,2,\cdots n, 且 m=nm=n,
    此时,
    hj(yi)={[1yix1++1yixj1+1yixj+1+1yixn],i=j(yix1)(yix2)(yixi1)(yixi+1)(yixj1)(yixj+1)(yixn)(xjx1)(xjx2)(xjxj1)(xjxj+1)(xjxn)ijh'_j(y_i)=\begin{cases} \left[\frac{1}{y_i-x_1}+\cdots +\frac{1}{y_i-x_{j-1}}+\frac{1}{y_i-x_{j+1}}+\cdots \frac{1}{y_i-x_{n}}\right], & i=j\\[1em] \frac{(y_i-x_1)(y_i-x_2)\cdots(y_i-x_{i-1})(y_i-x_{i+1})\cdots(y_i-x_{j-1})(y_i-x_{j+1})\cdots(y_i-x_{n})}{(x_j-x_1)(x_j-x_2)\cdots(x_j-x_{j-1})(x_j-x_{j+1})\cdots(x_j-x_{n})} & i\neq j\\ \end{cases}
    我们记 D~=[hj(xi)]i=1:nj=1:n\widetilde{\mathbb{D}}=[h'_j(x_i)]_{i=1:n}^{j=1:n}, 则有表达式
    D~=(RC)最后一个轴累乘,并替换对角线元素\widetilde{\mathbb{D}}=\left(\dfrac{\mathbb{R}}{\mathbb{C}}\right)_{最后一个轴累乘, 并替换对角线元素}
    这里矩阵C=[cjk]j=1,nk=1,n\mathbb{C}=[c_{jk}]_{j=1,n}^{k=1,n}, 以及三维矩阵R=[rijk]i=1:mj,k=1:n\mathbb{R}=[r_{ijk}]_{i=1:\red{m}}^{j,k=1:n}, 其中
    cjk={xjxk,jk1,j=k,以及 rijk={xixk,jki1,j=ki.c_{jk}= \begin{cases} x_j-x_k, & \quad j\neq k\\ 1, & \quad j=k \end{cases},\quad \text{以及 } r_{ijk} = \begin{cases} x_i-x_k, & \quad j\neq k 或 i\\ 1, & \quad j=k 或 i \end{cases} .
    注意: jhj(yi)=0,i\sum_{j} h'_j(y_i)=0, \forall i, 因此, 矩阵D~\widetilde{\mathbb{D}}的每行元素之和为0, 利用该性质可以很容易计算出D~\widetilde{\mathbb{D}}的对角线上的元素.

  2. yi=xj,,i,jy_i=x_j,\forall , i,j. 这种情况如果直接计算, 计算量太大, 这里给出一种比较简洁快速的方案, 假设插值点的纵坐标构成向量为v\bm v, 于是D~v\widetilde{\mathbb{D}}\bm v所得到的向量就是插值函数的导函数在插值点x\bm x处的值, 因此再做一次插值计算, 就可以得到插值函数的导数在y\bm y处的值. 于是有公式Dv=HD~v\mathbb{D}\bm{v}=\mathbb{H}\widetilde{\mathbb{D}}\bm v , 也就是,
    D=HD~.\mathbb{D}=\mathbb{H}\widetilde{\mathbb{D}}.

快速简洁的代码

import numpy as np
def interpBase(ip, iy, pureComputeH=True):
    c  = np.array([ip]).T - ip
    n = len(ip)
    ran = np.arange(n)
    c[ran, ran] = 1
    z = np.array([iy]).T - ip
    zs = np.tile(z, (n, 1, 1))
    zs[ran, :, ran] = 1
    zs = zs.transpose((1,0,2))
    abs_sort = lambda x: np.sort_complex(abs(x)+ 1j*x).imag 
    # abs_sort = lambda x: x
    matrixH = (abs_sort(zs)/abs_sort(c)).prod(axis=-1)
    if pureComputeH:
        return matrixH
    dzs = np.tile(c, (n, 1, 1))
    dzs[ran, :, ran] = 1
    dzs[ran, ran, :] = 1
    dzs = dzs.transpose((1,0,2))
    matrixD = (abs_sort(dzs)/abs_sort(c)).prod(axis=-1)
    matrixD[ran, ran] = 0
    matrixD[ran, ran] = -matrixD.sum(axis=-1)
    return matrixH, matrixH.dot(matrixD)

abs_sort 是绝对值排序, 但是实际用起来并不能明显提高精度, 最后给出测试代码.

# Compute LGL points and weights
from scipy import special as spc
deg = 60
x, _ = spc.roots_jacobi(deg-1, 1, 1)
x = np.hstack(([-1], x, [1]))
Lv = spc.eval_legendre(deg, x)
w = 2.0/(deg*(deg+1)*Lv**2)

# Compute matrix D 
a = np.array([Lv]).T / Lv
b = np.array([x]).T - x
ran = np.arange(deg+1)
a[ran, ran] = 0
b[ran, ran] = 1
D = a/b
D[0, 0] = - deg*(deg+1)/4.0
D[-1, -1] = deg*(deg+1)/4.0

D3 = interpBase(x, x, False)
err = abs(D - D3[-1])
print(err.max())

1.2.3. 计算矩阵D不稳定快速算法

如果完全硬算, 计算量太大, 而且也比较复杂. 这里提供一个矩阵向量乘积的快速算法, 也就是计算出Du\mathbb{D}\bm u. 需要借助FFT快速计算. 但是实际计算中该算法会丢精度.

全部Python代码如下:

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

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

def DMatrixVectorMultiplication(ip, iy, iv):
    # ip 插值点横坐标, iv 插值点纵坐标, iy需要赋值的点的横坐标
    # 直接算出矩阵向量的乘积
    n = len(ip)
    ran = np.arange(n)
    t = 2*np.pi/n
    w  = np.cos(t) - 1j*np.sin(t)
    w = w**ran
    v = InterpBase(ip).valuesm(w, np.complex).dot(iv)
    b = ifft(v)*ran
    c = np.zeros_like(b)
    c[:-1] = b[1:]
    c = fft(c)
    return InterpBase(w, np.complex).valuesm(iy, np.complex).dot(c).real

def DMatrix(ip, iy):
    # ip 插值点横坐标, iv 插值点纵坐标, iy需要赋值的点的横坐标
    # 只算出矩阵
    n = len(ip)
    ran = np.arange(n)
    t = 2*np.pi/n
    w  = np.cos(t) - 1j*np.sin(t)
    w = w**ran
    v = InterpBase(ip).valuesm(w, np.complex)
    b = ran*ifft(v.T)
    c = np.zeros_like(b)
    c[:, :-1] = b[:, 1:]
    c = fft(c).T
    return InterpBase(w, np.complex).valuesm(iy, np.complex).dot(c).real

注意:

2. 其它笔记(快速计算多项式乘法)

同样的道理, 可以借助FFT实现多项式乘法的快速计算

例如已知
f(x)=a0+a1x+anxmf(x)=a_0+a_1x+\cdots a_n x^m
g(x)=b0+b1x+bnxng(x)=b_0+b_1x+\cdots b_n x^n

h(x)=f(x)g(x)=c0+c1x+cnxm+nh(x)=f(x)g(x)=c_0+c_1x+\cdots c_n x^{m+n}的系数{ck}k=0m+n\{c_k\}_{k=0}^{m+n}.

快速算法的核心思想为:

Python 代码例如下

import numpy as np
from scipy.fftpack import fft,ifft

# 多项式乘积快速算法

a = np.array([1, 2, 0])
b = np.array([3, 4, 0])

c = ifft(fft(a)*fft(b))
print(c)


a = [1, 2, 3, 4]
b = [5, 6, 7, 8]

tz = [0]*(len(a)-1)
a = np.array(a + tz)
b = np.array(b + tz)

c = ifft(fft(a)*fft(b))
print(c.real)

注意: