插值算法

【问题】在比赛中,常常需要根据已知的函数点进行数据,模型的处理和分析,而有时候现有的数据又是极少的,不足以支撑分析的进行,这时候就需要使用一些数学的方法,“模拟产生”一些新的但又比较靠谱的值来满足需求,这就是插值的作用。
【引出】一维插值问题
已经有n+1个节点(xi,yi)(i=0,1,⋯ ,n)(x_i,y_i)(i=0,1,\cdots,n)(xi,yi)(i=0,1,,n)其中xix_ixi互不相同,不妨假设a=x0<x1<⋯<xn=ba=x_0<x_1<\cdots <x_n=ba=x0<x1<<xn=b,求任一插值点x∗(≠xi)x^*(\neq x_i)x(=xi)处的插值。
在这里插入图片描述思路:构造函数y=f(x)y=f(x)y=f(x),使得f(x)f(x)f(x)经过所有已知的节点,求f(x∗)f(x^*)f(x)即可得到y∗y^*y


插值法的定义

设函数y=f(x)y=f(x)y=f(x)在区间[a,b][a,b][a,b]上有定义,且已知在点a≤x0<x1<⋯<xn≤ba\leq x_0<x_1<\cdots <x_n\leq bax0<x1<<xnb上的函数值分别为:y0,y1,⋯ ,yny_0,y_1,\cdots,y_ny0,y1,,yn若存在一个简单函数P(x)P(x)P(x),使得P(xi)=yi,(i=0,1,⋯ ,n)P(x_i)=y_i ,(i=0,1,\cdots,n)P(xi)=yi,(i=0,1,,n)则称P(x)为f(x)P(x)为f(x)P(x)f(x)的插值函数,点x0,x1,⋯ ,xnx_0,x_1,\cdots,x_nx0,x1,,xn称为插值节点,f(x)f(x)f(x)叫做被插值函数,包含插值节点的区间[a,b][a,b][a,b]称为插值区间,求插值函数P(x)P(x)P(x)的方法称为插值法。


一般的多项式插值

  1. 插值法的原理:
  • 定理:设有n+1个不同的节点(xi,yi),(i=0,1,⋯ ,n)(x_i,y_i),(i=0,1,\cdots,n)(xi,yi),(i=0,1,,n)则存在唯一的多项式Ln(x)=a0+a1x+a2x2+⋯+anxnL_n(x)=a_0+a_1x+a_2x^2+\cdots+a_nx^nLn(x)=a0+a1x+a2x2++anxn使得Ln(xj)=yj(j=0,1,⋯ ,n)L_n(x_j)=y_j(j=0,1,\cdots,n)Ln(xj)=yj(j=0,1,,n)
  • 证明:n+1个方程解n+1个未知量(a0,a1,⋯ ,an)(a_0,a_1,\cdots,a_n)(a0,a1,,an)构造方程组:{a0+a1x0+a2x02+⋯+anx0n=y0a0+a1x1+a2x12+⋯+anx1n=y1⋯a0+a1xn+a2xn2+⋯+anxnn=yn\left\{\begin{matrix} a_0+a_1x_0+a_2x_0^2+\cdots+a_nx_0^n=y_0\\ a_0+a_1x_1+a_2x_1^2+\cdots+a_nx_1^n=y_1\\ \cdots\\ a_0+a_1x_n+a_2x_n2+\cdots+a_nx_n^n=y_n \end{matrix}\right.a0+a1x0+a2x02++anx0n=y0a0+a1x1+a2x12++anx1n=y1a0+a1xn+a2xn2++anxnn=yn令:A=[1x0⋯x0n1x1⋯x1n⋮⋮⋱⋮1xn⋯xnn],X=[a0a1⋮an],Y=[y0y1⋮yn]A=\begin{bmatrix} 1 & x_0 & \cdots & x_0^n\\ 1 & x_1 & \cdots & x_1^n\\ \vdots & \vdots & \ddots & \vdots\\ 1 & x_n & \cdots & x_n^n \end{bmatrix},X=\begin{bmatrix} a_0\\ a_1\\ \vdots\\ a_n \end{bmatrix},Y=\begin{bmatrix} y_0\\ y_1\\ \vdots\\ y_n \end{bmatrix}A=111x0x1xnx0nx1nxnn,X=a0a1an,Y=y0y1yn(A为范德蒙德行列式)方程组的矩阵形式如下:AX=YAX=YAX=Y由于∣A∣=∏i=1n∏j=0n−1(xi−xj)≠0|A|=\prod_{i=1}^{n}\prod_{j=0}^{n-1}(x_i-x_j)\neq 0A=i=1nj=0n1(xixj)=0所以该方程组有唯一解。从而Ln(x)=a0+a1x+a2x2+⋯+anxnL_n(x)=a_0+a_1x+a_2x^2+\cdots+a_nx^nLn(x)=a0+a1x+a2x2++anxn唯一存在,证毕。
  • 只要n+1个节点互异,满足上面插值条件的多项式是唯一存在的。
  • 如果不限制多项式的次数,插值多项式并不唯一。
  1. 拉格朗日插值法——一种多项式插值法
  • 两个点:(x0,y0),(x1,1)(x_0,y_0),(x_1,_1)(x0,y0),(x1,1)插值函数为y=x−x1x0−x1y0+x−x0x1−x0y1y=\frac{x-x_1}{x_0-x_1}y_0+\frac{x-x_0}{x_1-x_0}y_1y=x0x1xx1y0+x1x0xx0y1
  • 三个点:(x0,y0),(x1,y1),(x2,y2)(x_0,y_0),(x_1,y_1),(x_2,y_2)(x0,y0),(x1,y1),(x2,y2)插值函数为y=(x−x1)(x−x2)(x0−x1)(x0−x2)y0+(x−x0)(x−x2)(x1−x0)(x1−x2)y1+(x−x0)(x−x1)(x2−x0)(x2−x1)y2y=\frac{(x-x_1)(x-x_2)}{(x_0-x_1)(x_0-x_2)}y_0+\frac{(x-x_0)(x-x_2)}{(x_1-x_0)(x_1-x_2)}y_1+\frac{(x-x_0)(x-x_1)}{(x_2-x_0)(x_2-x_1)}y_2y=(x0x1)(x0x2)(xx1)(xx2)y0+(x1x0)(x1x2)(xx0)(xx2)y1+(x2x0)(x2x1)(xx0)(xx1)y2
  • 四个点:……
  • 缺点:龙格现象(Range phenomenon)
    在这里插入图片描述高次插值会产生龙格现象,即在两端处波动极大,产生明显的震荡,在不熟悉曲线的运动趋势的前提下,最好不要使用高次插值。

【解决办法】分段低次插值
插值多项式次数高,精度未必提高,反而误差可能越大。

  1. 分段低次插值
  • 分段低次插值主要是指分段线性插值和分段二次插值。
  • 分段线性插值:选取跟节点x最近的两个已知节点,求出经过这两个节点的直线即为节点x的插值函数
  • 分段二次插值:选取跟节点x最近的三个节点xi−1,xi,xi+1x_{i-1},x_i,x_{i+1}xi1,xi,xi+1进行二次插值,即在每一个区间[xi−1,xi+1][x_{i-1},x_{i+1}][xi1,xi+1]上取:f(x)=L2(x)f(x)=L_2(x)f(x)=L2(x)。在几何上就是每段都用抛物线代替f(x)f(x)f(x),故又叫做分段抛物插值。
  1. 牛顿插值法
  • f(x)=f(x0)+f[x0,x1](x−x0)+f[x0,x1,x2](x−x0)(x−x1)+⋯+f[x0,x1,⋯ ,xn−2,xn−1](x−x0)(x−x1)⋯(x−xn−2)]+f[x0,x1,⋯ ,xn](x−x0)(x−x1)⋯(x−xn−1)f(x)=f(x_0)+f[x_0,x_1](x-x_0)+f[x_0,x_1,x_2](x-x_0)(x-x_1)+\cdots+f[x_0,x_1,\cdots,x_{n-2},x_{n-1}](x-x_0)(x-x_1)\cdots(x-x_{n-2})]+f[x_0,x_1,\cdots,x_n](x-x_0)(x-x_1)\cdots(x-x_{n-1})f(x)=f(x0)+f[x0,x1](xx0)+f[x0,x1,x2](xx0)(xx1)++f[x0,x1,,xn2,xn1](xx0)(xx1)(xxn2)]+f[x0,x1,,xn](xx0)(xx1)(xxn1)
  • 差商的定义:称f[x0,xk]=f(xk)−f(x0)xk−x0f[x_0,x_k]=\frac{f(x_k)-f(x_0)}{x_k-x_0}f[x0,xk]=xkx0f(xk)f(x0)为函数关于点x0,xkx_0,x_kx0,xk的一阶差商(亦称均差)。
  • 二阶差商:f[x0,x1,x2]=f[x1,x2]−f[x0,x1]x2−x0f[x_0,x_1,x_2]=\frac{f[x_1,x_2]-f[x_0,x_1]}{x_2-x_0}f[x0,x1,x2]=x2x0f[x1,x2]f[x0,x1]
  • k阶差商:f[x0,x1,⋯ ,xk]=f[x1,⋯ ,xk−1,xk]−f[x0,x1,⋯ ,xk−1]xk−x0f[x_0,x_1,\cdots,x_k]=\frac{f[x_1,\cdots,x_{k-1},x_k]-f[x_0,x_1,\cdots,x_{k-1}]}{x_k-x_0}f[x0,x1,,xk]=xkx0f[x1,,xk1,xk]f[x0,x1,,xk1] (xk−1,xk可以不相邻)(x_{k-1},x_k可以不相邻)xk1,xk
  • 与拉格朗日插值法相比,牛顿插值法的计算过程具有继承性。(牛顿插值法每次插值只和前n项的有关,这样每次只要在原来的函数上添加新的项,就能够产生新的函数)但是牛顿插值法也存在龙格现象。
  1. 拉格朗日插值法与牛顿插值法的另一个缺点

    这两种插值仅仅要求插值多项式在插值节点处与被插函数有相等的函数值,而这种插值多项式却不能全面反映被插值函数f(x)f(x)f(x)的性态。
    然而在许多实际问题中,不仅仅要求插值函数P(x)P(x)P(x)与被插值函数f(x)f(x)f(x)在所有的节点处具有相同的函数值,它也需要在一个或者全部节点上插值多项式与被插函数有相同的低阶甚至高阶的导数值。
    对于这些情况,拉格朗日插值法和牛顿插值法都不能满足。


埃尔米特(Hermite)插值(多项式插值)

  1. 基本思路:不但要求在节点上函数值相等,而且还要求对应的导数值也相等。
    φ(xi)=yi,(i=0,1,⋯ ,n)φ′(xi)=yi′,(i=0,1,⋯ ,n)\begin{aligned} \varphi (x_i)=y_i,(i=0,1,\cdots,n)\\ \varphi ^{'}(x_i)=y_i^{'},(i=0,1,\cdots,n) \end{aligned}φ(xi)=yi,(i=0,1,,n)φ(xi)=yi,(i=0,1,,n)保持插值曲线在节点处有切线(光滑),使插值函数与被插值函数的密和程度更好。

  2. Hermite插值法的原理:
    f(x)f(x)f(x)在区间[a,b][a,b][a,b]上有n+1个互异节点a=x0<x1<x2<⋯<xn=ba=x_0<x_1<x_2<\cdots<x_n=ba=x0<x1<x2<<xn=b,定义在[a,b][a,b][a,b]上函数f(x)f(x)f(x)在节点上满足:f(xi)=yi,f′(xi)=yi′,(i=0,1,2,⋯ ,n)f(x_i)=y_i,f^{'}(x_i)=y_i^{'},(i=0,1,2,\cdots,n)f(xi)=yi,f(xi)=yi,(i=0,1,2,,n)(2n+2个条件)可唯一确定一个次数不超过2n+1的多项式H2n+1(x)=H(x)H_{2n+1}(x)=H(x)H2n+1(x)=H(x)满足:H(xj)=yj,H′(xj)=yj′,(j=0,1,2,⋯ ,n)H(x_j)=y_j,H^{'}(x_j)=y^{'}_j,(j=0,1,2,\cdots,n)H(xj)=yj,H(xj)=yj,(j=0,1,2,,n)其余项为:R(x)=f(x)−H(x)=f(2n+2)(ξ)(2n+2)!ω2n+2(x)R(x)=f(x)-H(x)=\frac{f^{(2n+2)}(\xi)}{(2n+2)!}\omega_{2n+2}(x)R(x)=f(x)H(x)=(2n+2)!f(2n+2)(ξ)ω2n+2(x)

  3. 分段三次埃尔米特插值

  • 直接使用Hermite插值得到的多项式的次数较高,也存在着龙格现象,因此在实际应用中,往往使用分段三次Hermite插值多项式(PCHIP)
  • 使用方法:调用Matlab内置函数
    p=pchip(x,y,new_x)
    x是已知样本点的横坐标,y是已知样本点的纵坐标,new_x是要插入处对应的横坐标,p是插入处对应的纵坐标。
    【例子】
    x=-pi:pi;
    new_x=-pi:0.1:pi;
    y=sin(x);
    p=pchip(x,y,new_x);
    plot(x,y,'o',new_x,p,'r-.')
    
    在这里插入图片描述

三次样条插值

  • 基本原理:
    y=f(x)y=f(x)y=f(x)在点x0,x1,x2,⋯ ,xnx_0,x_1,x_2,\cdots,x_nx0,x1,x2,,xn处的1值为y0,y1,y2,⋯ ,yny_0,y_1,y_2,\cdots,y_ny0,y1,y2,,yn,若·函数S(x)S(x)S(x)满足下列条件:
    (1)S(xi)=f(xi)=yi,i=0,1,2,⋯ ,nS(x_i)=f(x_i)=y_i,i=0,1,2,\cdots,nS(xi)=f(xi)=yi,i=0,1,2,,n
    (2)在每个子区间[xi,xi+1](i=0,1,2,⋯ ,n−1)[x_i,x_{i+1}](i=0,1,2,\cdots,n-1)[xi,xi+1](i=0,1,2,,n1)S(x)S(x)S(x)是三次多项式
    (3)S(x)S(x)S(x)[a,b][a,b][a,b]上二阶连续可微
    则称S(x)S(x)S(x)为函数f(x)f(x)f(x)的三次样条插值函数。

  • 构造方法:S(x)S(x)S(x)除了满足基本的插值条件S(xi)=f(xi)S(x_i)=f(x_i)S(xi)=f(xi)之外还应具有以下形式:S(x)={S0(x),x∈[x0,x1]S1(x),x∈[x1,x2]⋮Sn−1(x),x∈[xn−1,xn],Si(x)∈C3([xi,xi+1])S(x)=\left\{\begin{matrix} S_0(x),x\in [x_0,x_1]\\ S_1(x),x\in [x_1,x_2]\\ \vdots\\ S_{n-1}(x),x\in [x_{n-1},x_{n}] \end{matrix}\right.,S_i(x)\in C^{3}([x_i,x_{i+1}])S(x)=S0(x)x[x0,x1]S1(x)x[x1,x2]Sn1(x)x[xn1,xn],Si(x)C3([xi,xi+1])并且满足条件:{Si−1(xi)=Si(xi),插值的连续性Si−1′(xi)=Si′(xi)一阶微分的连续性Si−1′′(xi)=Si′′(xi)二阶微分的连续性\left\{\begin{matrix} S_{i-1}(x_i)=S_i(x_i),插值的连续性\\ S_{i-1}^{'}(x_i)=S_i^{'}(x_i)一阶微分的连续性\\ S_{i-1}^{''}(x_i)=S_i^{''}(x_i)二阶微分的连续性 \end{matrix}\right.Si1(xi)=Si(xi),Si1(xi)=Si(xi)Si1(xi)=Si(xi)

  • 使用方法:matlab内置函数:
    p=spline(x,y,new_x)
    x是已知点的横坐标,y是已知点的纵坐标,new_x是要插入处点的横坐标,p是要插入处点的纵坐标。

    x=-pi:pi;
    new_x=-pi:0.1:pi;
    y=sin(x);
    p1=pchip(x,y,new_x);
    p2=spline(x,y,new_x);
    plot(x,y,'o',new_x,p1,'r-',new_x,p2,'b-')
    legend('样本点','三次埃尔米特插值','三次样条插值','Location','SouthEast')
    %三次样条插值更加光滑,更接近被插值函数sin(x)
    

在这里插入图片描述由图可以看出:三次样条插值生成的曲线更加光滑,在实际建模中,由于我们不知道题目中数据的生成过程,因此三次埃尔米特插值与三次样条插值都可以使用。


n维数据的插值

p=interpn(x1,x2,…,xn,y,new_x1,new_x2,…,new_xn,method)
x1,x2,…,xn是已知样本点的横坐标,y是对应的纵坐标,new_x1,new_x2,…,new_xn是要插入点的横坐标,p是对应的纵坐标。
method 是要插值的方法:

  • ‘linear’:线性插值(默认方法)
  • ‘cubic’:三次插值
  • ‘spline’:三次样条插值(最为精准)
  • ‘nearest’:最邻近插值算法

【拓展】
上面学的插值算法有时候可以用来短期预测。

评论 1
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包

打赏作者

奇迹luanluan

你的鼓励将是我创作的最大动力

¥1 ¥2 ¥4 ¥6 ¥10 ¥20
扫码支付:¥1
获取中
扫码支付

您的余额不足,请更换扫码支付或充值

打赏作者

实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值