文章目录
- 基本原理
- 算法实现
无论是Newton插值还是Lagrange插值，都只能在数值本身上满足插值函数与数据节点的重合，Hermite插值则要求其导数值相等。
基本原理
设在节点 a ⩽ x 0 ⩽ x 1 ⩽ … ⩽ x n ⩽ b a\leqslant x_0\leqslant x_1 \leqslant\ldots\leqslant x_n\leqslant b a⩽x0⩽x1⩽…⩽xn⩽b上， y j = f ( x j ) , m j = f ′ ( x j ) , j = 0 , 1 , … , n y_j=f(x_j),m_j=f'(x_j),j=0,1,\ldots,n yj=f(xj),mj=f′(xj),j=0,1,…,n，则插值多项式的条件为
H ( x j ) = y j , H ′ ( x j ) = m j , j = 0 , 1 , … , n H(x_j)=y_j,H'(x_j)=m_j,j=0,1,\ldots,n H(xj)=yj,H′(xj)=mj,j=0,1,…,n
总共有 2 n + 2 2n+2 2n+2个等式，可唯一确定一个次数不大于 2 n + 1 2n+1 2n+1的多项式 H 2 n + 1 ( x ) H_{2n+1}(x) H2n+1(x)，设其表达式可通过基函数 α j , β j \alpha_j,\beta_j αj,βj写成
H 2 n + 1 ( x ) = ∑ j = 0 n [ y j α j + m j β j ] H_{2n+1}(x)=\displaystyle\sum_{j=0}^n[y_j\alpha_j+m_j\beta_j] H2n+1(x)=j=0∑n[yjαj+mjβj]
则 α j , β j \alpha_j,\beta_j αj,βj需要满足
{ α j ( x k ) = δ j k α j ′ ( x k ) = 0 β j ( x k ) = 0 β j ′ ( x k ) = δ j k \left\{\begin{aligned} \alpha_j(x_k) &= \delta_{jk}\\ \alpha_j'(x_k)&= 0\\ \beta_j(x_k) &= 0\\ \beta_j'(x_k) & = \delta_{jk} \end{aligned}\right. ⎩ ⎨ ⎧αj(xk)αj′(xk)βj(xk)βj′(xk)=δjk=0=0=δjk
两个基函数的条件与Lagrange插值中的基函数相似，由于目标函数为 2 n + 1 2n+1 2n+1次，所以假设
{ α j = ( a x + b ) l j 2 ( x ) β j = ( c x + d ) l j 2 ( x ) \left\{\begin{aligned} \alpha_j = (ax+b)l_j^2(x)\\ \beta_j = (cx+d)l_j^2(x) \end{aligned}\right. {αj=(ax+b)lj2(x)βj=(cx+d)lj2(x)
结合二者所满足的基函数条件，可得
{ α j ( x j ) = ( a x j + b ) l j 2 ( x j ) = 1 α j ′ ( x j ) = l j ( x j ) [ a l j ( x j ) + 2 ( a x j + b ) l j ′ ( x j ) ] = 0 β j ( x j ) = ( c x j + d ) l j 2 ( x j ) = 0 β j ′ ( x j ) = l j ( x j ) [ c l j ( x j ) + 2 ( c x j + d ) l j ′ ( x j ) ] = 1 \left\{\begin{aligned} \alpha_j(x_j)&=(ax_j+b)l_j^2(x_j)=1\\ \alpha_j'(x_j)&=l_j(x_j)[al_j(x_j)+2(ax_j+b)l_j'(x_j)]=0\\ \beta_j(x_j)&=(cx_j+d)l_j^2(x_j)=0\\ \beta_j'(x_j)&=l_j(x_j)[cl_j(x_j)+2(cx_j+d)l_j'(x_j)]=1 \end{aligned}\right. ⎩ ⎨ ⎧αj(xj)αj′(xj)βj(xj)βj′(xj)=(axj+b)lj2(xj)=1=lj(xj)[alj(xj)+2(axj+b)lj′(xj)]=0=(cxj+d)lj2(xj)=0=lj(xj)[clj(xj)+2(cxj+d)lj′(xj)]=1
解得
{ α j = [ 1 − 2 （ x − x j ) ∑ k = 0 , k ≠ j n 1 x j − x k ] l j 2 ( x ) β j = ( x − x j ) l j 2 ( x ) \left\{\begin{aligned} \alpha_j &= [1-2（x-x_j)\displaystyle\sum_{k=0,k\not=j}^n\frac{1}{x_j-x_k}]l_j^2(x)\\ \beta_j &= (x-x_j)l_j^2(x) \end{aligned}\right. ⎩ ⎨ ⎧αjβj=[1−2（x−xj)k=0,k=j∑nxj−xk1]lj2(x)=(x−xj)lj2(x)
令 s j = ∑ k = 0 , k ≠ j 1 x j − x k s_j=\displaystyle\sum_{k=0,k\not=j}\frac{1}{x_j-x_k} sj=k=0,k=j∑xj−xk1,则Hermite插值公式可写为
H 2 n + 1 ( x ) = ∑ j = 0 n { [ y j + 2 s j y j x j − x j m j ] ⋅ l j 2 ( x ) + [ − 2 s j y j + m j ] ⋅ x l j 2 ( x ) } H_{2n+1}(x)=\displaystyle\sum_{j=0}^n\{ [y_j+2s_jy_jx_j-x_jm_j]\cdot l_j^2(x) +[-2s_jy_j+m_j]\cdot xl_j^2(x) \} H2n+1(x)=j=0∑n{[yj+2sjyjxj−xjmj]⋅lj2(x)+[−2sjyj+mj]⋅xlj2(x)}
算法实现
其代码为
# Hermite插值函数
# x,y分别为自变量、因变量，m为导数
function Hermite(x,y,m)
n = length(x)
unitMat = Matrix{Float64}(I,n,n)
xMat = x .- x' + unitMat
xInverse = 1./xMat - unitMat
A = zeros(n*2+1)
for i = 1:n
quot = polyDiv(polyMat,[-x[i],1])[1]
lag = quot ./ sumMulti(xMat[i,:]) #Lagrange基函数
lag = polyMulti(lag,lag)
sumInverse = sum(xInverse[i,:])
a0 = (1+2*x[i]*sumInverse)*y[i]-x[i]*m[i]
a1 = -2*y[i]+m[i]
α = 1-2
A += polyMulti(lag,[a0,a1])
end
return A
end
验证
> include("interpolation.jl");
> x = [i for i = 1:9];
> y = x.^8 + x.^6 + x.^2 .+ 1;
> m = [rand() for i = 1:9]
> a = Hermite(x,y,m);
> y1 = y
> y2 = [sum([a[i]*x[j]^(i-1) for i=1:18]) for j = 1:9]
> plot(x,[y1,y2],st=[:scatter :line])
> savefig("hermite.png")
由于引入了一组随即变量代表导数，所以拟合得到的多项式函数自然与输入的多项式不同
> f(x) = sum([a[i]*x^(i-1) for i=1:18])
> plot(f,xlims=(1,9))
> savefig("realfig.png")
