三次样条插值 订阅
三次样条插值(Cubic Spline Interpolation)简称Spline插值,是通过一系列形值点的一条光滑曲线,数学上通过求解三弯矩方程组得出曲线函数组的过程。实际计算时还需要引入边界条件才能完成计算。一般的计算方法书上都没有说明非扭结边界的定义,但数值计算软件如Matlab都把非扭结边界条件作为默认的边界条件。 展开全文
三次样条插值(Cubic Spline Interpolation)简称Spline插值,是通过一系列形值点的一条光滑曲线,数学上通过求解三弯矩方程组得出曲线函数组的过程。实际计算时还需要引入边界条件才能完成计算。一般的计算方法书上都没有说明非扭结边界的定义,但数值计算软件如Matlab都把非扭结边界条件作为默认的边界条件。
信息
应    用
数学工程等
外文名
Cubic Spline Interpolation
简    称
Spline插值
中文名
三次样条插值
性    质
数学术语
基    础
样条插值
三次样条插值基本概念
早期工程师制图时,把富有弹性的细长木条(所谓样条)用压铁固定在样点上,在其他地方让它自由弯曲,然后沿木条画下曲线。成为样条曲线。
收起全文
精华内容
下载资源
问答
  • 三次样条插值

    2018-05-24 09:19:21
    三次样条插值 三次样条插值 三次样条插值 三次样条插值 三次样条插值
  • Spline(三次样条插值)

    万次阅读 多人点赞 2015-09-09 09:59:58
    三次样条插值

    关于三次样条插值,计算方法比较复杂,但是静下心来仔细研究也是可以理解的。

    本文借鉴文章来源:http://www.cnki.com.cn/Article/CJFDTotal-BGZD200611035.htm

    定义:

         简单来说就是给定了一些在区间[a,b]的数据点{x1,x2,x3.....xn},对应函数值{y1,y2,y3.....yn},函数在[xj,xj+1]  (j=1,2,...n-1此处根据你的编译器所定,matlab数组下标从1开始的)上有表达式S(x),且满足下面条件:

    1. S(x)是一个三次多项式,在这里设为

    2. S(xj)=yj                    (2-1)

    3. S(xj-0)=S(xj+0)  (j=2,3....,n-1)这就是保证了在非端点处的其它点连续               (2-2)

    4. S'(xj-0)=S'(xj+0) (j=2,3....,n-1)一阶导数连续                              (3-1)

    5. S''(xj-0)=S''(xj+0) (j=2,3....,n-1)二阶导数连续                              (3-2)

    【注】3 4 5的区间从2开始到n-1是因为两个端点不需要判断是否连续,端点处没连续这一说。

           有一个说法“n个未知数需要n个方程才能确定唯一解”,具体对不对,可以参考线性代数的知识。我们的最终目标是求出每个区间的式(1)或者函数值。 共有n-1个区间,每个区间四个参数aj,bj,cj,dj,那么就总共需要4(n-1)个求未知数。在(2-1)中给出了n个方程,(2-2)(3-1)(3-2)总共给出了3(n-2)个方程,所以依据唯一解方法可知还需要4(n-1)-3(n-2)-n=2个方程。

           对于剩下的两个方程,三次样条插值给出了两个边界约束方程,刚好凑齐两个,并且有三种,可以依照自己的兴趣选择一种便于实现的。

    1. 给定了端点处的一阶导数值:S'(x1)=y1'    S'(xn)=yn'

    2.给定了端点处的二阶导数值:S''(x1)=y1''    S''(xn)=yn''。自然边界条件:y1''=yn''=0

    3.给定了一个周期性条件:如果f(x)是以b-a为周期的函数,在端点处便满足:S'(x1+0)=S'(xn-0),S''(x1+0)=S''(xn-0)

    下面的推导是以边界方程1为例的:

    推导过程:

    (一) 利用二阶导得到一些式子:

          S(x)为每个区间的三次多项式,那么S''(x)就是一次多项式。假设S''(xj)=Mj,S''(xj+1)=Mj+1的值已知,那么:
           
    其中hj=xj+1-xj,然后利用S''(x)求S(x):
     

    (二)依据S(x)得到一些式子

    按照上面的证明可以得到:
    其中:hj=xj+1-xj
    依据S(xj)=yj和S(xj+1)=yj+1得到:
    ——————更新日志2020-6-13————————
    感谢评论区老哥@meng_zhi_xiang ,公式(5)的A下标不是j+1,应该是j
    解得(求解过程就不写了,简单的二元方程组)
                                         (6-1)
              (6-2)
    将Aj和Bj带入S(x)中可以得到S(x)的最终式子:
    原文此处应该是有错误,Aj和Bj都没有xj的式子,但是原文的结果中包含。
     
    这里就得到了S(x)的雏形了,xj、xj+1、hj都是已知的,但是Mj和Mj+1是假设已知的,下面就是求它们了。

    (三)利用一阶导数得到一些式子

    依据式(7)求出一阶导:
    然后为了使在xj处连续平滑,那么在xj处的一阶导必须相等。即要满足S'(xj-0)=s'(xj+0)
    等式坐标表示S(x)在[xj-1,xj]的xj的一阶导值,右边表示[xj,xj+1]的xj一阶导值
    其中利用S(x)的表达式可以得到等式左右两边的值:
     
    由上得到两个式子:
     
    由两式相等整理可得:
    令:
    则μjMj-1+2Mj+(1-μj)Mj+1=dj(j=2,3,...,n-1)

    (四)带入边界条件

      此处选择边界条件1,即
    计算:
    结果写为2M1+M2=β1
     
    结果写为Mn-1+2Mn=βn

    (五)结果

    结合(三)、(四)可以看出所有的n个式子已经齐全:
    包括(三)结尾算得μjMj-1+2Mj+(1-μj)Mj+1=dj(j=2,3,...,n-1)的n-2个式子加上(四)得到的两个式子,刚好n个式子.
    用矩阵表示出这n个式子即为:
    方程中
       
    并且hj=xj+1-xj
    这样便解出了矩阵M,然后带入S(x)的式子即上面算得:
    这样便求得了每一个区间上的S(x)了。
     
    matlab代码:
    SplineThree.m
    <pre name="code" class="plain">function f = SplineThree(x,y,dy1,dyn,x0)
    %x,y为输入的已知点
    %x0为待求插值点的横坐标
    %dy1,dyn为约束条件,是端点处的一阶导数值
    n=length(x);
    for j=1:n-1
        h(j)=x(j+1)-x(j);
    end
    
    %得到式子右边的结果矩阵
    d(1,1)=6*((y(2)-y(1))/h(1)-dy1)/h(1);   %等式(11)的结果矩阵的上端点值
    d(n,1)=6*(dyn-(y(n)-y(n-1))/h(n-1))/h(n-1); %等式(11)的结果矩阵的下端点值
    for i=2:n-1
        u(i)=h(i-1)/(h(i-1)+h(i));
        d(i,1)=6*((y(i+1)-y(i))/h(i)-(y(i)-y(i-1))/h(i-1))/(h(i-1)+h(i));
    end
    
    %得到系数矩阵
    A(1,1)=2;
    A(1,2)=1;
    A(n,n-1)=1;
    A(n,n)=2;
    for i=2:n-1
        A(i,i-1)=u(i);
        A(i,i)=2;
        A(i,i+1)=1-u(i);
    end
    
    %解方程
    M=A\d;
    
    format long
    syms t;
    %得到每个区间的式子
    for i=1:n-1
       a(i)=y(i+1)-M(i+1)*h(i)^2/6-((y(i+1)-y(i))/h(i)-(M(i+1)-M(i))*h(i)/6)*x(i+1);
       b(i)=((y(i+1)-y(i))/h(i)-(M(i+1)-M(i))*h(i)/6)*t;
       c(i)=(t-x(i))^3*M(i+1)/(6*h(i));
       e(i)=(x(i+1)-t)^3*M(i)/(6*h(i));
       f(i)=a(i)+b(i)+c(i)+e(i);
        %f(i)=M(i)*(x(i+1)-t)^3/(6*h(i))+M(i+1)*(t-x(i))^3/(6*h(i))+(y(i)-M(i)*h(i)^2/6)*(x(i+1)-t)/h(i)+(y(i+1)-x(i+1)*h(i)^2/6)*(t-x(i))/h(i);
        % f(i)=((x(j+1)-x)^3)*M(i)/(6*h(i))+((x-x(i))^3)*M(i+1)/(6*h(i))+(y(i)-M(i)*(h(i)^2)/6)*((x(i+1)-x)/h(i))+(y(i+1)-(M(i+1)*(h(i)^2)/6))*((x-x(i))/h(i));
    end
    
    %化简
     f=collect(f);
     f=vpa(f,6);
    
     
    if(nargin==5)
       nn=length(x0);
    for i=1:nn
        for j=1:n-1
            if(x0(i)>=x(j)&x0(i)<=x(j+1))
                 yynum(i)=subs(f(j),'t',x0(i));   %计算插值点的函数值.subs是替换函数,把x0用t替换
            end
        end
    end   
    f=yynum;
    else
        f=collect(f);          %将插值多项式展开
        f=vpa(f,6);            %将插值多项式的系数化成6位精度的小数
    end
    end
    
     
    
     
    SplineThree.m
    <pre name="code" class="plain">% x=[-1.5 0 1 2];
    % y=[0.125 -1 1 9];
    % dy1=0.75;
    % dyn=14;
    % x0=1.5;
    % answer=SplineThree(x,y,dy1,dyn);
    % 
    % X=[-1.5 0 1 2];
    % Y=[0.125 -1 1 9];
    % dY=[0.75 14]
    % m=5;
    % spline3(X,Y,dY,x0,m)
    
    X=0:2*pi;
    Y=sin(X);
    dY=[1,1];
    dy1=1;
    dyn=1;
    xx=0:0.5:6;
    m=5;
    % for i=1:length(xx)
    % spline3(X,Y,dY,xx,m);
    % end
    yy=SplineThree(X,Y,dy1,dyn,xx);
    plot(xx,yy,'r')

    本文已经同步到微信公众号中,公众号与本博客将持续同步更新运动捕捉、机器学习、深度学习、计算机视觉算法,敬请关注

     
    展开全文
  • 内容提要:三次样条插值的定义、三次样条插值函数的建立方法、不同边值条件下求解三次样条插值函数的方法.《递推算法与多元插值》简介本书详细介绍了多元差商与多元逆差商的递推算法及其在多元多项式与连分式插值中...

    内容提要:三次样条插值的定义、三次样条插值函数的建立方法、不同边值条件下求解三次样条插值函数的方法.


    d568b1fef8f8a881fe80f96730164564.png

    c4e1d9e0837cfcb6d4c1f7f4881fd67f.png


    aea8433466cc4a1471cffdd904cfee6e.png


    eff7deebff76e9da93a645db8cca4631.png


    7bda2ab65e19790e79b49cb06c0644fc.png


    eb877611ecb82b1a6fddac76686590e8.png


    1241a9f866d922e277b9d933535b8945.png


    d28754c005ffad4e697549d6b0b878bf.png


    《递推算法与多元插值》简介

    本书详细介绍了多元差商与多元逆差商的递推算法及其在多元多项式与连分式插值中的应用。内容包括常用的张量积型二元多项式与连分式插值方法概述、直角三点组上的二元多项式与连分式插值及其比较研究、直角三点组上二元多项式插值余项等的进一步研究、非矩形网格上的二元多项式插值、基于二元递推多项式的散乱数据插值、基于二元连分式的散乱数据插值递推格式、非张量积型二元连分式插值、金字塔型网格点上的三元分叉连分式插值等。

    本书可作为计算数学、应用数学等学科高年级本科生、硕士生、博士生数值分析、数值逼近、计算几何等相关课程的参考书或专业著作,还可供从事数值逼近与计算几何、计算机辅助几何设计、图像处理及相关领域的科技工作者参考。

    本书获国家科学技术学术著作出版基金资助,入选“信息与计算科学丛书”。

    展开全文
  • 所谓三次样条插值对于一个区间(a,b)将区间分成x0 = a < x1 ......xn-1 < b = xn 的n-1个区间,我们需要通过已知的n+1个点来模拟一个未知的函数,在三次样条插值中我们采用分段的方法来做这件事情。三次样条...

    所谓三次样条插值对于一个区间(a,b)将区间分成x0 = a < x1 ......xn-1 < b = xn 的n-1个区间,我们需要通过已知的n+1个点来模拟一个未知的函数,在三次样条插值中我们采用分段的方法来做这件事情。

    三次样条插值得到的分段函数保证一下条件成立,而这些条件也是用来求解每一段样条插值的条件:

    1 模拟出来的函数在已知点的函数值等于f的函数值

    2模拟出来的分段函数是二阶连续的也就是说导数和二阶导数在分段的交界点是相等的

    3需要知道在a和b点的二阶导数的情况,或者二阶导数在这n+1个点的变化规律

    下面直接转载其内在的规律

    已知:

    a. n+1个数据点[xi, yi], i = 0, 1, …, n

    b. 每一分段都是三次多项式函数曲线

    c. 节点达到二阶连续

    d. 左右两端点处特性(自然边界,固定边界,非节点边界)

    根据定点,求出每段样条曲线方程中的系数,即可得到每段曲线的具体表达式。

    插值和连续性:

    , 其中 i = 0, 1, …, n-1

    微分连续性:

     , 其中 i = 0, 1, …, n-2

    样条曲线的微分式:

    将步长

     带入样条曲线的条件:

    a. 由 (i = 0, 1, …, n-1)推出

    b. 由 (i = 0, 1, …, n-1)推出

    c. 由  (i = 0, 1, …, n-2)推出

    由此可得:

    d. 由  (i = 0, 1, …, n-2)推出

     ,则

    a.  可写为:

     ,推出

    b. 将ci, di带入  可得:

    c. 将bi, ci, di带入 (i = 0, 1, …, n-2)可得:

    这样我们可以构造一个以m为未知数的线性方程组,而且在端点条件已知的情况下我们是知道其中几个未知数的值的

    端点条件

    由i的取值范围可知,共有n-1个公式, 但却有n+1个未知量m 。要想求解该方程组,还需另外两个式子。所以需要对两端点x0和xn的微分加些限制。 选择不是唯一的,3种比较常用的限制如下。

    a. 自由边界(Natural)

    首尾两端没有受到任何让它们弯曲的力,即 。具体表示为 和 

    则要求解的方程组可写为:

    b. 固定边界(Clamped)

    首尾两端点的微分值是被指定的,这里分别定为A和B。则可以推出

    将上述两个公式带入方程组,新的方程组左侧为

    c. 非节点边界(Not-A-Knot)

    指定样条曲线的三次微分匹配,即

    根据 和 ,则上述条件变为

    新的方程组系数矩阵可写为:

    右下图可以看出不同的端点边界对样条曲线的影响:

    1.3 算法总结

    假定有n+1个数据节点

    a. 计算步长 (i = 0, 1, …, n-1)

    b. 将数据节点和指定的首位端点条件带入矩阵方程

    c. 解矩阵方程,求得二次微分值。该矩阵为三对角矩阵,具体求法参见我的上篇文章:三对角矩阵的求解。

    d. 计算样条曲线的系数:

    其中i = 0, 1, …, n-1

    e. 在每个子区间 中,创建方程

    原创作者博客地址:https://blog.csdn.net/flyingleo1981/article/details/53008931

    展开全文
  • 分段三次样条插值3. 三对角线性系统的求解4. 固定边界三次样条的求解5. 附录5.1 龙格库塔现象1. 拟合与插值一个情形有时我们获得了某条曲线上若干个离散点的信息,却不知道具体的曲线方程,又希望得到这些离散点...

    三次样条(cubic spline)插值

    1. 拟合与插值2. 分段三次样条插值3. 三对角线性系统的求解4. 固定边界三次样条的求解5. 附录5.1 龙格库塔现象

    1. 拟合与插值

    一个情形

    有时我们获得了某条曲线上若干个离散点的信息,却不知道具体的曲线方程,又希望得到这些离散点之外的信息。

    两种方法

    • 拟合:不要求方程通过所有已知点,讲究神似,曲线整体趋势一致即可。

      85f52b29853ad1a47505c32303d3607d.png

    • 插值:要求方程通过所有已知点,讲究形似。

    • 0c621fdad5f596cca38123f668178c4d.png

    样条插值

    在数值分析这个数学分支中,样条插值是使用一种名为样条的特殊分段多项式进行插值的形式。由于样条插值可以使用低阶多项式样条实现较小的插值误差,这样就避免了使用高阶多项式所出现的龙格现象(见附录1),所以样条插值得到了流行。常用的是三次样条插值

    2. 分段三次样条插值

    所谓分段三次样条插值,顾名思义,就是将区间 分成 个区间 bb715d38366adc7cb68678cdd1cf44fb.png ,共有 个点,其中端点值 bbe0f43d51e0974f6d8783f8b4d07491.png 。三次样条就是说每个分段区间上的曲线都是一个三次多项式方程 ,并满足以下条件。

    • 1. 在每个分段区间 c1affe2269094ce6c2663e5ec64d05a0.png 上,89ea3941d28e9199ff7b67f2d3a2bd9b.png 都是一个三次多项式方程。

    • 2. 在所有已知点上满足插值条件,即 6107a278ef03f231665a71448faa6f1c.png

    • 3. 曲线光滑,即 cff90e9f91ea694c0616d942f18a78eb.png 平滑, 连续。

    我们可将构造如下的三次多项式方程:6aae8e154bf19bd9f53f54be8d9596e4.png
    我们称这个方程为三次样条函数 。显然,每个分段区间上的三次样条函数均有 个未知系数,有 个分段区间,共有 个未知系数,需要 个方程来求解。

    分析

    我们需要寻找出 个方程去求解 个未知系数。

    • 1. 所有已知点满足插值条件,即 6107a278ef03f231665a71448faa6f1c.png 。除端点外的 内部点 7ff4ddee59268923784fddb81e1ffcc2.png 满足 7fe7910702adb444ae493ad056fd367f.png ,共有 41592786eedc278b49f8152a5fcce16b.png 个方程;两个端点满足 46a697aff00b92e0e471d3deb6ec2d10.png ,共有 个方程。总计 个方程。

    • 2. 个内部点一阶导数连续,即 46fd8a4402eee784f382037a88385b4f.png 。总计 个方程。

    • 3. 个内部点二阶导数连续,即 68cd69655b3e4d9ffc3009da5d551bca.png 。总计 个方程。

    这样,我们就找出了 个方程了,还差 个方程就可以解出所有未知系数了。这 个方程我们通过边界条件给出。

    • 自然边界(natural spline):指定端点的二阶导数为 ,即 b437fe1b2ef36ea6684177e0b31233e4.png

    • 固定边界(clamped spline):指定端点一阶导数为特定值,即 ba52147e170b0988cb57f7b20be944c9.png

    • 非扭结边界(not-a-knot spline):使最前面两个插值点的三阶导数值相等、最后面两个插值点的三阶导数值相等,即 1e0e9caba8575725a733eeac2703292a.png

    推导

    为了方便推导,我们令:99c537aafac08f2c08144c5f3dcb9473.png

    • dbda363cb9dd85c0bc10e8a8e6afe0a2.png 可得:9638256e8b7b6c7c19ac17c3b3dde9fc.png

    • 82ae5426aaf46d0c8bdd862a676e0236.png 并令 c65b9aad2c7d95c07cf6fd4151e83b0b.png 可得:805fc37c0d25625ddcf80bf2c5c7e690.png

    • 46fd8a4402eee784f382037a88385b4f.png 可得:

      e765631e90450a1111da997900827399.png

    • 68cd69655b3e4d9ffc3009da5d551bca.png 可得:

      a10a59f8636b1f0143ab8f19e538cecd.png

    我们令 c466925c55c281a3612c07212e352569.png ,则 ,由 可得 7c7d63d87f34f8a4238903de78e1271c.png

    7c7d63d87f34f8a4238903de78e1271c.png 代入 中,可得:6782e67e678e00e554c26c33b07e1859.png
    至此,35dbb604360a30ccb46e32da5a6c75c1.png 均表示成二阶导的关系式,将其代入 可得:3379ead37b0421025d39db8813a70e2b.png
    这样,我们就能构造出一个以二阶导 为未知数的线性方程组了。

    • 自然边界e6e5e3b3a732c61c6aaa9ee66af7d9df.png

    75b14d87164e1b387c0738a3e59477c9.png

    可以看出,左侧的系数矩阵为严格对角占优矩阵,故线性方程组有唯一解。再根据 35dbb604360a30ccb46e32da5a6c75c1.png 关于二阶导的表达式就能求出每个分段区间上的具体方程了。

    • 固定边界ba52147e170b0988cb57f7b20be944c9.png

      1d88c24cc3cb885e3fdb29b3ef14fcf7.png 有: ,而 9b424b4572e0e3e0972d78fded3d2f5b.png ,则有:79784445b11bc749ed216ce90199d021.png
      893faf48de7c965b7697c688e5ffabb8.png 有:f163e8edab698f6263dd5d6a64effe21.png ,将 ff7b5ea025002fa3ed33f29bdfccab78.png 的表达式带入可得:7a46c52e9f2fa6dea6f34be874afe451.png
      则线性方程组为:5a49007cd9f8a329911a9981a4628136.png

    • 非扭结边界fa4a86ddc519717619f052b5d9d5fd06.png

      55a16f3a3e9cb112d2f695482d72e93f.png 有:9d77cac54f9016c7a25e564ed8b589a9.png ,代入 的表达式可得:fb44b01de5910e56241c3760f0fae6a6.png
      则线性方程组为:c690f73acfb87d3780846e9e1144a1f2.png

    算法

    输入: 个数据节点 1d3d1627aa8ff79475ab8b8868045feb.png

    [1] 计算步长 c65b9aad2c7d95c07cf6fd4151e83b0b.png

    [2] 配置并求解线性方程组;

    [3] 计算样条曲线的系数;

    [4] 在每个分段区间 c1affe2269094ce6c2663e5ec64d05a0.png 上创建三次多项式方程 c5c6db0abf9eb6ff4acaecd3012fdada.png

    3. 三对角线性系统的求解

    对于系数矩阵为严格对角占优三对角矩阵的线性方程组,有唯一解。且其 分解、前代及回代过程,只需 d89dfba738342d6dd0f7cae936a42cb6.png 次操作。f47b8927bbdf412e74c441e9960a7632.png
    在编程方面,通常无需存储整个 系数矩阵,使用三个列向量来代替。

    void tridag(int n, double a[], double b[], double c[], double r[], double u[]){
        double bet;
        double *gma = vector(1, n - 1, (double)(NULL));

        bet = b[0];
        u[0] = r[0] / bet;

        for (int i = 1; i// 分解及前代
            gma[i] = c[i - 1] / bet;
            bet = b[i] - a[i] * gma[i];
            u[i] = (r[i] - a[i] * u[i - 1]) / bet;
        }

        for (int i = n - 2; i >= 0; i--) {    // 回代
            u[i] -= gma[i + 1] * u[i + 1];
        }
    }

    4. 固定边界三次样条的求解

    仿照前述三对角线性系统的求解方法,对于固定边界的三次样条,编写如下样条插值的功能块。

    /**
     * 在处理列表函数时,程序spline仅调用一次。一旦调用过后,对任意的xx,插值函数的值通过splint获得。
     **/


    /**
     * 给定数组x[0...n-1](严格单调递增)与y[0...n-1]构成一个列表函数,即y[i] = f(x[i])。
     * yp1与ypn表示插值函数分别在第1个与第n-1个点处的一阶导数。
     * 函数修改y2[0...n-1],表示插值函数在表中点x_i处的二阶导数。
     **/

    void spline(double x[], double y[], int n, double yp1, double ypn, double y2[]){
        double bet;
        double *gma = vector(1, n-1, (double)(NULL));

        bet = 2*(x[1]-x[0]);
        y2[0] = 6*((y[1]-y[0])/(x[1]-x[0]) - yp1) / bet;

        for(int i=1; i-1; i++) {    // 分解及前代
            gma[i] = (x[i]-x[i-1]) / bet;
            bet = 2*(x[i+1]-x[i-1]) - (x[i]-x[i-1]) * gma[i];
            y2[i] = (6*((y[i+1]-y[i])/(x[i+1]-x[i]) - (y[i]-y[i-1])/(x[i]-x[i-1])) - (x[i]-x[i-1]) * y2[i-1]) / bet;
        }
        gma[n-1] = (x[n-1]-x[n-2]) / bet;
        bet = 2*(x[n-1]-x[n-2]) * gma[n-1];
        y2[n-1] = (6*(ypn - (y[n-1]-y[n-2])/(x[n-1]-x[n-2])) - (x[n-1]-x[n-2]) * y2[n-2]) / bet;for(int i=n-2; i>=0; i--) {    // 回代
            y2[i] -= gma[i+1] * y2[i+1];
        }
    }/**
     * 给定数组x[0...n-1]和y[0...n-1]作为列表函数,数组y2a[0...n-1]由spline程序输出给定。
     * 给定要计算的xx值,程序通过yy返回其三次样条插值结果。
     **/
    void splint(double x[], double y[], double y2[], int n, double xx, double *yy){int k, kl=0, kh=n-1;double a, b, c, d;// find the corresponding interval with binary searchingwhile(kh-kl > 1) {
            k = (kh+kl)/2;if (x[k] > xx) {
                kh = k;
            } else {
                kl = k;
            }
        }
        a = y[kl];
        b = (y[kh]-y[kl])/(x[kh]-x[kl]) - (x[kh]-x[kl])*y2[kl]/2 - (x[kh]-x[kl])*(y2[kh]-y2[kl])/6;
        c = y2[kl]/2;
        d = (y2[kh]-y2[kl])/(x[kh]-x[kl])/6;
        *yy = a + b*(xx-x[kl]) + c*(xx-x[kl])*(xx-x[kl]) + d*(xx-x[kl])*(xx-x[kl])*(xx-x[kl]);
    }

    测试

    针对上述功能块 fc63ad26c794a87d3cf5c71501c96250.png398457a909ba664c1de27768192161ea.png ,进行如下测试:在 cec2efce62f6fafe70951bd41bf0f262.png 上选取若干个控制点,计算控制点之外的插值点。

    #include 
    #include 
    #include "nrutil.h"

    int main(){
        FILE *fp;
        int n = 21;

        double xx;
        double yy;
        double yp1 = -20;
        double ypn = 20;

        double *x = vector(0, n-1, (double)(NULL));
        double *y = vector(0, n-1, (double)(NULL));
        double *y2 = vector(0, n-1, (double)(NULL));

        // select n control points from f(x)=x^2
        fp = fopen("control.dat","w");
        for (int i=0; i        x[i] = i-10.0;
            y[i] = (i-10.0)*(i-10.0);
            fprintf(fp,"%lf\t%lf\n",x[i],y[i]);
        }
        fclose(fp);

        spline(x, y, n, yp1, ypn, y2);

        // calculate the interpolating points
        fp=fopen("interpolate.dat","w");
        xx = -10.0;
        while (xx <= 10.0) {
            splint(x, y, y2, n, xx, &yy);
            fprintf(fp,"%lf\t%lf\n",xx,yy);
            xx += 0.1;
        }
        fclose(fp);

        return 0;
    }

    画图脚本:

    set xlabel 'x'
    set ylabel 'y'

    set multiplot layout 1,1 title 'cubic spline'

    plot "control.dat" u 1:2 t 'control' w points pointtype 7 linecolor 'red' linewidth 1 pointsize 1,\
         "interpolate.dat" u 1:2 t 'interpolatation' w points pointtype 6 linecolor 'blue' linewidth 1 pointsize 1 

    unset multiplot

    效果图:

    50dd39743138d9dc5293f0fbd27e4839.png

    5. 附录

    5.1 龙格库塔现象

    龙格现象是在一组等间距插值点上使用高次多项式插值时出现的区间边缘处的振荡问题,它表明使用高次多项式插值并不总是能提高准确性。随着多项式阶次的增加,以这种方式产生的多项式 0bab07c57b782790b3da9660322b8e48.png 实际上可能偏离 ,通常发生在靠近插值点的端点。

    解决龙格现象的方法:

    使用分段多项式样条可以避免龙格库塔现象。如果要减小插值误差,那么可以增加构成样条的多项式的数目,而不必是增加多项式的阶次。

    7f56269fb06be1d59a46b671c23be45f.png

    展开全文
  • 三次样条插值分段线性插值的优点:计算简单、稳定性好、收敛性有保证且易在计算机上实现缺点:它只能保证各小段曲线在连接点的连续性,却无法保证整条曲线的光滑性,这就不能满足某些工程技术的要求。三次Hermit插值...
  • 三次样条插值函数高次插值函数的计算量大,有剧烈振荡,且数值稳定性差;在分段插值中,分段线性插值在分段点上仅连续而不可导,分段三次埃尔米特插值有连续的一阶导数,如此光滑程度常不能满足物理问题的需要,样条...
  • 本篇介绍力求用容易理解的方式,介绍一下三次样条插值的原理,并附C语言的实现代码。1. 三次样条曲线原理假设有以下节点1.1 定义样条曲线是一个分段定义的公式。给定n+1个数据点,共有n个区间,三次样条方程满足以下...
  • MATLAB程序分享三次样条插值法求信号的包络线源程序-MATLAB三次样条插值法 求信号的包络线源程序代码.rar
  • 1 三次样条插值 早期工程师制图时,把富有弹性的细长木条(所谓样条)用压铁固定在样点上,在其他地方让它自由弯曲,然后沿木条画下曲线,成为样条曲线。 设函数S(x)∈C2[a,b] ,且在每个小区间[xj, xj+1]上是三次...
  • 三次样条插值算法的数学原理和公式推导在这里就不详细展开,主要是其在轨迹生成中的应用,具体了参考以下博客:https://blog.csdn.net/adamshan/article/details/80696881​blog.csdn.net数学公式:(1)求出点与点...
  • PAGE PAGE #/ 4 第一边界条件源代码 function y=yt1(x0,y0,f_0,f_n,x) (1) %第一类边界条件下三次样条插值 %xi 所求点 %yi 所求点函数值 %x已知插值点 %y已知插值点函数值 %f_0 左端点一次导数值 %f_n 右端点一次导...
  • 前面我们提到,轨迹即包含时间这一维度的路径,...其中,轨迹生成是非常重要的一步,在本节我们介绍一种基于三次样条插值的路径生成方法,后面我们将结合Moritz Werling提出的基于Frenet坐标系的最佳轨迹生成来完成...

空空如也

空空如也

1 2 3 4 5 ... 20
收藏数 1,207
精华内容 482
关键字:

三次样条插值