精华内容
下载资源
问答
  • 拉格朗日插值和分段线性插值
    千次阅读
    2021-01-30 18:35:19

    《用MATLAB实现拉格朗日插值和分段线性插值》由会员分享,可在线阅读,更多相关《用MATLAB实现拉格朗日插值和分段线性插值(6页珍藏版)》请在人人文库网上搜索。

    1、用MATLAB实现拉格朗日插值和分段线性插值1、 实验内容:用MATLAB实现拉格朗日插值和分段线性插值。2、 实验目的:1) 学会使用MATLAB软件;2) 会使用MATLAB软件进行拉格朗日插值算法和分段线性差值算法;3、实验原理:利用拉格朗日插值方法进行多项式插值,并将图形显式出来。4、实验步骤及运行结果(1)实现lagrange插值1)定义函数: f = 1/(x2+1) 将其保存在f.m 文件中,具体程序如下:function y = f1(x)y = 1./(x.2+1);2) 定义拉格朗日插值函数:将其保存在lagrange.m 文件中,具体实现程序编程如下:function y。

    2、 = lagrange(x0,y0,x)m = length(x); /区间长度/n = length(x0);for i = 1:nl(i) = 1;endfor i = 1:mfor j = 1:nfor k = 1:nif j = kcontinue;endl(j) = ( x(i) -x0(k)/( x0(j) - x0(k) )*l(j);endendendy = 0;for i = 1:ny = y0(i) * l(i) + y;end3) 建立测试程序,保存在text.m文件中,实现画图:x=-5:0.001:5; y=(1+x.2).-1; p=polyfit(x,y,n); 。

    3、py=vpa(poly2sym(p),10) plot_x=-5:0.001:5; f1=polyval(p,plot_x); figureplot(x,y,r,plot_x,f1)输入n=6,出现下面的图形:通过上图可以看到当n=6是没有很好的模拟。于是重新运行text.M并选择n=11由此可见n=11时的图像是可以很好的实现模拟(2)分段线性插值:建立div_linear.m文件。具体编程如下/*分段线性插值函数:div_linear.m 文件*/function y = div_linear(x0,y0,x,n)%for j = 1:length(x)for i = 1:n-1if (x。

    4、 = x0(i) & (x = x0(i+1)y = (x - x0(i+1)/(x0(i) - x0(i+1)*y0(i) + ( x - x0(i)/(x0(i+1) - x0(i)*y0(i+1);elsecontinue;endend%end测试程序(text2.m):n = input(输入n =:);x0 = linspace( -5,5,n);for x = -5:0.01:5y = div_linear(x0,f(x0),x,n);hold on;plot(x,y,r);plot(x,f(x),b);end2)运行测试程序,这是会出现:输入n=:2)输入n=6,并按Enter键。

    5、,出现:4)关掉图形界面后,重新运行程序,输入n=11,并按enter键后出现:5)再次关掉图形界面,输入n=100,并按enter键,出现:此时。图形将于原函数图形基本吻合,说明分割区间越多,图像接近真实的图像。(3)用lagrange插值观察y = |sin(k*x)|的误差分析:1)编写 函数文件,保存在f2.m 中x=0:0.01:1;k= input(输入k:)n= input(输入n:);y=abs(sin(k*pi*x);p=polyfit(x,y,n-1);py=vpa(poly2sym(p),8);plot_x=0:0.01:1;f1=polyval(p,plot_x);pl。

    6、ot(x,y,plot_x,f1); 2)运行该程序:输入k=:1输入n=:2出现如下图形界面:关掉图形界面后重新运行f2.m,输入k=:1,n=:3出现如下界面:再次关掉图形界面,输入k=:1,n=:6 后出现:此时图形基本吻合。类推,输入k=2, n=3后出现:k =2, n =11,出现如下图形:k =2,n =15,这时出现:k =2,n =19,出现:当k=2,n=21时,图形如下:此时基本吻合。5、实验总结:通过本次课程设计,我初步掌握了MATLAB运用,加深了对于各种线性插值的理解;培养了独立工作能力和创造力;综合运用专业及基础知识,解决实际数学问题的能力;在本次课程设计中,在老师的精心指导下,收益匪浅。同时对数学的研究有了更深入的认识。

    更多相关内容
  • 用 MATLAB 实现拉格朗日插值和分段线性插值 1 实验内容 用 MATLAB 实现拉格朗日插值和分段线性插值 2 实验目的 1 学会使用 MATLAB 软件 2 会使用 MATLAB 软件进行拉格朗日插值算法和分段线性 差值算法 3实验原理 ...
  • 用MATLAB 实现拉格朗日插值和分段线性插值 1 实验内容 用MATLAB 实现拉格朗日插值和分段线性插值 2 实验目的 1 学会使用MATLAB 软件 2 会使用MATLAB 软件进行拉格朗日插值算法和分段线 性差值算法 3实验原理 利用...
  • 用MATLAB实现拉格朗日插值和分段线性插值.pdf
  • 用MATLAB实现拉格朗日插值和分段线性插值.doc
  • 实验四 用 MATLAB实现拉格朗日插值分段线性插值 一实验目的 1学会使用 MATLAB软件 2 会使用 MATLAB软件进行拉格朗日插值算法分段线性差值算法 二实验内容 1 用 MATLAB实现 y = 1./(x^2+1; -1的拉格朗日插值分段...
  • 实验四 用 MATLAB 实现拉格朗日插值分段线性插值 一实验目的 1学会使用MATLAB 软件 2会使用MATLAB 软件进行拉格朗日插值算法分段线性差值算法 二实验内容 1 用 MATLAB 实现 y = 1./(x^2+1; -1的拉格朗日插值分段...
  • 实验四用MATLAB实现拉格朗日插值分段线性插值.pptx
  • 计算方法的三种基础插值方法,拉格朗日插值线性插值;三样条插值
  • 这个是用C语言编的关于插值的代码,主要是三种插值方式,为拉格朗日插值法,分段线性插值法和三次样条插值法,三次样条采用追赶法。
  • 该程序由c++编写,主要用于实现基于函数y=e^(-2x)在区间[0,6]的插值函数,开发工具为VS2015,请用此IDE或者更高版本的IDE打开工程文件~~~~~
  • 实验四用MATLAE实现拉格朗日插值分段线性插值 一实验目的 )学会使用MATLAB软件 )会使用MATLAB软件进行拉格朗日插值算法分段线性差值算法 二实验内容 1用MATLAB实现y二l./(x.A2+l; (-1二xUl )的拉格朗日插值分段...
  • 实验四 用MATLAB实现拉格朗日插值分段线性插值
  • 实验四 用 MATLAB 实现拉格朗日插值分段线性插值 一实验目的 学会使用 MATLAB 软件 会使用 MATLAB 软件进行拉格朗日插值算法分段线性差值算法 二实验内容 1 用 MATLAB 实现 y = 1./(x^2+1;-1的拉格朗日插值分段...
  • 实验四用MATLAB实现拉格朗日插值分段线性插值.doc
  • 拉格朗日插值分段线性插值、三次样条插值

    万次阅读 多人点赞 2018-08-30 21:20:05
    本篇主要介绍在三种插值方法:拉格朗日插值分段线性插值、三次样条插值,以及这三种方法在matlab中如何实现。 1.拉格朗日插值:  1.1基本原理:先构造一组基函数:   是次多项式,满足 令 上式称为...

    本篇主要介绍在三种插值方法:拉格朗日插值、分段线性插值、三次样条插值,以及这三种方法在matlab中如何实现。

    1.拉格朗日插值:

      1.1基本原理:先构造一组基函数:

    l_{i}(x)=\frac{(x-x_{0})...(x-x_{i-1})(x-x_{i+1})...(x-x_{n})}{(x_{i}-x_{0})...(x_{i}-x_{i-1})(x_{i}-x_{i+1})...(x_{i}-x_{n})}

             =\prod_{j=0j\neq i}^{n}\frac{x-x_{j}}{x_{i}-x_{j}}             \left ( i=0,1,...,n \right )

    l_{i}(x)n次多项式,满足

    l_{i}(x_{j})=\left\{\begin{matrix} 0 & j\neq i & \\ 1& j=i & \end{matrix}\right.

    L_{n}(x)=\sum_{i=0}^{n}y_{i}l_{i}(x)=\sum_{i=0}^{n}y_{i}\left \{ \prod_{j=0j\neq i}^{n}\frac{x-x_{j}}{x_{i}-x_{j}} \right \}

    上式称为n次Lagrange插值多项式。

    1.2用Matlab作Lagrange插值:

    matlab没有现成的lagrange函数,需要手动写,如下:

    x0,y0为原始坐标点,维度必须相同。

    x为待插值的点。

    y是返回值,是最终插值结果。

    function y=lagrange(x0,y0,x)   %x0,y0为初始坐标,x为需要插值的点,返回的y为插值结果
    n=length(x0);m=length(x);
    for i=1:m
        z=x(i);
        s=0;
        for j=1:n
            p=1;
            for k=1:n
                if k~=j
                    p=p*((z-x0(j))/(x0(k)-x0(j)));
                end
            end
            s=p*y0(k)+s;
        end
        y(i)=s;
    end

     2.分段线性插值:

    2.1基本原理:

    将每两个相邻的节点用直线连起来,如此形成的一条折线就是分段线性插值函数。计算x点的插值时,只用到x左右的两个节点,计算量与节点个数n(初始值x0,y0的长度,n=length(x0))无关,而拉格朗日插值与n值有关。分段线性插值中n越大,分段越多,插值误差越小。

    2.2Matlab实现分段线性插值:

    用matlab实现分段线性插值不需要自己手动编写函数,matlab有现成的一维插值函数interp1

    y=interp1(x0,y0,x,'method')

    method指定插值方法,其值可为:

    linear:线性插值(默认)

    nearest:最近项插值

    spline:逐次3次样条插值

    cubic:保凹凸性 3 次插值

    所有插值方法都要求x0单调。

    3.三次样条插值:

    使用三次样条插值有两种方法:其中一种就是第二种插值方式(分段线性插值),只需要将method修改为spline即可实现。

    还有一种实现方式如下:

    pp=csape(x0,y0);

    y=ppval(pp,x);

    其中x0,y0,x与前面含义相同,返回值y即插值结果。

    4.彩蛋(例题):

    表1给出的 x, y 数据位于机翼断面的下轮廓线上,假设需要得到 x 坐标每改变0.1 时的 y 坐标。试完成加工所需数据,画出曲线。要求用 Lagrange、分段线性和三次样条三种插值方法计算。

    表1

    x   0 3 5 7 9 11 12 13 14 15
    y   0 1.2 1.7 2.0 2.1 2.0 1.8 1.2 1.0 1.6

    解:编写代码如下:

    clear,clc
    x0=[0,3,5,7,9,11,12,13,14,15];
    y0=[0,1.2,1.7,2.0,2.1,2.0,1.8,1.2,1.0,1.6];
    x1=[0:0.1:15];
    
    %拉格朗日插值
    y1=lagrange(x0,y0,x1);
    figure
    plot(x0,y0,x1,y1,'.')
    title('拉格朗日插值')
    
    %分段线性插值
    y2=interp1(x0,y0,x1);
    figure
    plot(x0,y0,x1,y2,'.')
    title('分段线性插值')
    
    %三次样条插值
    y3=interp1(x0,y0,x1,'spline');
    figure
    plot(x0,y0,x1,y3,'.')
    title('三次样条插值')
    
    

    结果图为:

    由图可知:拉格朗日插值误差较大,插值大量点的时候不建议使用。分段线性插值结果最好,使用的是method是linear(默认)。三次样条插值有点误差,但不是很大,得到的曲线比较光滑。

     

    学习插值函数时的一点小心得,希望能帮到大家!

    展开全文
  • 线性插值 原理 流程图 代码 x0=0.2; y0=21; x1=0.4; y1=25; x=0.7; L0=(x-x1)/(x0-x1); L1=(x-x0)/(x1-x0); y=y0*L0+y1*L1; 抛物插值 原理 流程图 代码 x0=0.2; x1=0.4; x2=0.6; y0=21; y1...

    目录

     

    线性插值 

    原理

    流程图

     代码

     抛物插值

    原理

    流程图

     代码

     拉格朗日插值

    代码

    牛顿插值

    原理

    代码

    分段线性插值

    代码


    线性插值 

    原理

    流程图

     单个点的线性插值代码

    X=[0.2 0.4];
    Y=[21 25];
    x=0.7;
    x0=X(1)
    y0=Y(1);
    x1=X(2);
    y1=Y(2);
    
    
    L0=(x-x1)/(x0-x1);
    L1=(x-x0)/(x1-x0);
    y=y0*L0+y1*L1;
    

    多个点的线性插值代码

    time = [1 3 8 12 15 20 24];
    tem = [8 9 16 23 22 18 10];
    time_i = 1:0.01:24;
    tem_i = self_interp1(time,tem,time_i,'linear');
    plot(time,tem,'o',time_i,tem_i);
    
    %自己写一个interp1类似功能的接口
    %在这个接口中参数x需要从大到小排列,y随意
    function [yi]=self_interp1(x,y,xi,method)
    % 初始化yi,给它xi对应的列
    col_xi = size(xi,2);
    yi = zeros(1,col_xi);
    % 检测使用的插值方法 这里期望的是'linear'
    if strcmp(method,'linear')
        % 找到每个xi在x序列中的位置
        col_x = size(x,2);
        for i = 1:col_xi,
            for j = 1:col_x-1, 
                % 假如需要计算插值公式
                if x(j+1) > xi(i), 
                    yi(i) = y(j)+(y(j+1)-y(j))/(x(j+1)-x(j))*(xi(i)-x(j));
                    break;
                end
                % 假如插值处的数据已经测得了,就直接把值给它,节约计算资源
                if x(j) == xi(i),
                    yi(i) = y(j);
                    break;
                end
            end
            % 以上没有把最后一个数据点考虑进去,需要加上
            yi(col_xi) = y(col_x);
        end
    else
        error('插值方法请选择(linear)\n');
    end
    end

     抛物插值

    原理

    流程图

     代码

    X=[0.2 0.4 0.6];
    Y=[21 25 23];
    x0=X(1);
    x1=X(2);
    x2=X(3);
    y0=Y(1);
    y1=Y(2);
    y2=Y(3);
    x=0.7;
    
    L0=(x-x1)*(x-x2)/(x0-x1)/(x0-x2);
    L1=(x-x0)*(x-x2)/(x1-x0)/(x1-x2);
    L2=(x-x0)*(x-x1)/(x2-x0)/(x2-x1);
    
    y=y0*L0+y1*L1+y2*L2

     拉格朗日插值

    代码

    x=[0.2 0.4 0.6 0.8];
    y=[21 25 23 20];
    yh=lagrange(x,y,0.7)
    
    function yh=lagrange (x,y,xh)
        n = length(x);
        m = length(xh);
        yh = zeros(1,m); 
        c1 = ones(n-1,1);
        c2 = ones(1,m);
        for i=1:n
          xp = x([1:i-1 i+1:n]);
          yh = yh + y(i)*prod((c1*xh-xp'*c2)./(x(i)-xp'*c2));
        end
    end

    牛顿插值

    原理

     

    代码

    xi=[1 4 9];
    yi=[1 2 3];
    x=7;
    p= Newton_fun(x,xi,yi)
    
    
    function p= Newton_fun(x,xi,yi)
    n=length(xi);
    f=zeros(n,n);
    
    % 对差商表第一列赋值
    for k=1:n      
        f(k)=yi(k);
    end
    % 求差商表
    for i=2:n       % 差商表从0阶开始;但是矩阵是从1维开始存储!!!!!!
        for k=i:n
            f(k,i)=(f(k,i-1)-f(k-1,i-1))/(xi(k)-xi(k+1-i));  
        end
    end
    disp('差商表如下:');
    disp(f);
    
    %求插值多项式
    p=0;          
    for k=2:n
        t=1;
        for j=1:k-1
            t=t*(x-xi(j));
            disp(t)
        end
        p=f(k,k)*t+p;
        disp(p)
    end
    p=f(1,1)+p;
    
    end
    

    分段线性插值

    原理

     

    代码

    x = [1 3 8 12 15 20 24];
    y = [8 9 16 23 22 18 10];
    
    yy=fdxx(x,y,7)
    
    function yy=fdxx(x,y,xx)
        n=size(x,2);
        for i=1:n-1
            if x(i)<xx&&xx<x(i+1)
                L1=(xx-x(i+1))/(x(i)-x(i+1));
                L2=(xx-x(i))/(x(i+1)-x(i));
                yy=L1*y(i)+L2*y(i+1);
                break;
            elseif x(i)==xx
                 yy=y(i);      
            end
        end
        
    end

    分段抛物插值

    原理

    代码

    x = [1 3 8 12 15 20 24];
    y = [8 9 16 23 22 18 10];
    y=fenduanpaowu(x,y,7)
    
    function y=fenduanpaowu(xi,yi,x)
        n=size(xi,2);
        if x<xi(2)
            L1=(x-xi(1))*(x-xi(3))/(xi(1)-xi(2))/(xi(1)-xi(3));
            L2=(x-xi(1))*(x-xi(3))/(xi(2)-xi(1))/(xi(2)-xi(3));
            L3=(x-xi(1))*(x-xi(2))/(xi(3)-xi(1))/(xi(3)-xi(2));
            y=L1*yi(1)+L2*yi(2)+L3*yi(3);
        elseif x>xi(end-1)
            L1=(x-xi(end-1))*(x-xi(end))/(xi(end-2)-xi(end-1))/(xi(end-2)-xi(end));
            L2=(x-xi(end-2))*(x-xi(end))/(xi(end-1)-xi(end-2))/(xi(end-1)-xi(end));
            L3=(x-xi(end-2))*(x-xi(end-1))/(xi(end)-xi(end-2))/(xi(end)-xi(end-1));
            y=L1*yi(1)+L2*yi(2)+L3*yi(3);
       else
            for k=2:n-1
                if xi(k+1)>x  
                    if abs(x-xi(k))<abs(x-xi(k+1))
                        i=k-1;
                        
                    else
                        i=k;
                    end
                    L1=(x-xi(i+1))*(x-xi(i+2))/(xi(i)-xi(i+1))/(xi(i)-xi(i+2));
                    L2=(x-xi(i))*(x-xi(i+2))/(xi(i+1)-xi(i))/(xi(i+1)-xi(i+2));
                    L3=(x-xi(i))*(x-xi(i+1))/(xi(i+2)-xi(i))/(xi(i+2)-xi(i+1));
                    y=L1*yi(i)+L2*yi(i+1)+L3*yi(i+2);
                end
            end
       end
        
    end

    展开全文
  • 文章目录一、基础概念插值是什么拟合是什么插值拟合的相同点插值拟合的不同点二、常用的基本插值方法高次多项式插值拉格朗日多项式插值牛顿插值差商矩阵低次多项式插值(不易震荡)分段线性插值Hermite插值三次...

    摘要(必看)

    本文特别长(长到需要些摘要·····),这对可读性造成很大负面影响,但我觉得都是插值,放在一起比较方便对比,所以目录逻辑做的比较清晰,实际上又臭又长的本文仅仅讲了两大类插值方法,一类是只考虑最基本插值条件的多项式插值,一类是还要考虑导数的插值,第一类又主要细分为拉格朗日多项式插值和牛顿多项式插值,且分别讲了基于他俩的分段低次插值,第二类主要细分为只考虑一阶导数相等的hermite插值和考虑一二阶导数相等的样条插值,并讲了基于他俩的三次分段插值。所有插值方法都给出了示例,并有图和代码。文末给出了精华总结。

    另外,本文参考了很多清风老师的数学建模教程中有关插值的资料(部分公式、代码和文字描述),引用的文字部分都统一用灰框框起来的,清风老师的这套资料和视频讲解特别详细有用,可以在B站找到,虽然要花点小钱,但亲测真的非常值得。

    文章太长,滚轮滚起来太慢,可以用crtl + home或者ctrl + end快速到达顶部和底部哦

    0 基础概念

    什么是插值

    求出已知个有限数据点的近似函数

    插值法是数值分析中最基本的方法之一。 在实际问题中遇到的函数是许许多多的,有的甚至给不出表达式,只供给了一些离散数据,例如,在查对数表时,需要查的数值在表中却找不到,所以只能先找到它相邻的数,再从旁边找出它的更正值,按一定的关系把相邻的数加以更正,从而找出要找的数,这种更正关系事实上就是一种插值。采用不同的插值函数,逼近的效果也不同。

    当然啦,到后面你就会知道,已知数据点上取值相同只是插值的最基本条件,是必须满足的基础条件,这个都不满足就不是插值了。

    在这里插入图片描述在这里插入图片描述在这里插入图片描述

    插值用途

    补全缺失数据
    基于已知数据进行预测

    什么是拟合

    求出一个不要求过已知数据点的近似函数,不要求过那些数据点,只要求在这些点上的总偏差最小。
    在这里插入图片描述

    插值和拟合的相同点

    都是根据一组数据构造一个近似函数

    插值和拟合的不同点

    近似的要求不同;数学方法完全不同

    一个问题到底应该插值还是拟合,有时容易确定有时不能明显看出。
    我个人觉得数据点给的很多的时候就不再需要插值了,只需要拟合,看数据的走向趋势即可,因为本身插值就是在数据点太少的情况下看不出数据趋势才需要去生成新的可靠数据点的。另一方面是因为,数据点本身很多的话,多项式插值则次数很高,龙格现象会造成不准确插值。当然这一点可以通过分段低次插值克服,下面会讲


    1 常用的基本插值方法

    1.1 多项式插值法

    不同的多项式插值方法都是去构造一个n次插值多项式

    多项式插值是函数插值中最常用的一种形式。在一般的插值问题中,插值条件可以唯一地确定一个次数不超过 的插值多项式。从几何上可以解释为:可以从多项式曲线中找出一些不超过 次的点通过平面上 个不同的点。插值多项式有两种常用的表达式形式,一种是拉格朗日插值多项式,另一种是牛顿插值多项式,此外拉格朗日插值公式与牛顿插值公式永远相等。

    在这里插入图片描述

    只要n+1个节点互不相同,就一定存在唯一的n次多项式使得这n+1个点过这个多项式,厉害!!啊哈哈

    若不限制多项式次数为n,就不唯一了哦。也是从证明里可以看出来的
    在这里插入图片描述在这里插入图片描述

    1.1.1 拉格朗日多项式插值法

    拉格朗日插值多项式是一种特殊的多项式形式,刚好过所有插值节点,是一种很凑巧的多项式在这里插入图片描述在这里插入图片描述
    例子
    在这里插入图片描述在这里插入图片描述在这里插入图片描述

    所以拉格朗日多项式一共n项,n是已知数据点(插值节点)的数目

    多项式插值并不是次数越大越好(龙格现象)

    在这里插入图片描述在这里插入图片描述

    分段低次线性插值以提高精度,减小插值误差

    分段二次插值

    只考虑最近的三个点,进行构造分段的拉格朗日插值多项式,以避免高次多项式插值带来的龙格效应,实现更高的精度
    在这里插入图片描述
    这张图的公式有误,应该改为
    f ( x ) ≈ L 2 ( x ) = ∑ k = i − 1 i + 1 [ y k ∏ j = i − 1 j = k̸ i + 1 ( x − x j ) ( x k − x j ) ] f(x)\approx L_2(x)=\sum_{k=i-1}^{i+1}[y_k \prod_{j=i-1 \atop j=\not k}^{i+1}\frac{(x-x_j)}{(x_k-x_j)}] f(x)L2(x)=k=i1i+1[ykj=kj=i1i+1(xkxj)(xxj)]

    实际就是
    在这里插入图片描述

    分段线性插值

    分段线性插值本质上是多项式插值,而且是加权式的插值方法。

    在这里插入图片描述在这里插入图片描述
    第一个图中的 F i F_i Fi通分后:
    x i f ( x i + 1 ) − x i + 1 f ( x i ) x i − x i + 1 + f ( x i ) − f ( x i + 1 ) x i − x i + 1 x \frac{x_if(x_{i+1})-x_{i+1}f(x_i)}{x_i-x_{i+1}}+\frac{f(x_i)-f(x_{i+1})}{x_i-x_{i+1}}x xixi+1xif(xi+1)xi+1f(xi)+xixi+1f(xi)f(xi+1)x
    这其实就是 ( x i , f ( x i ) ) (x_i,f(x_i)) (xi,f(xi)) ( x i + 1 , f ( x i + 1 ) ) (x_{i+1},f(x_{i+1})) (xi+1,f(xi+1))两点的连线的一次函数

    所以分段线性插值的子插值函数没什么高深的,就是每段两个端点的连线。

    l i ( x ) l_i(x) li(x)就像是一个权重函数,我们用一个配巨丑手画图的例子来具体看看它为什么是权重

    假设有5点已知数据点, x 1 , x 2 , ⋯   , x 5 x_1,x_2,\cdots,x_5 x1,x2,,x5,那么子插值函数 F 1 , F 2 , F 3 , F 4 F_1,F_2,F_3,F_4 F1,F2,F3,F4就分别是相邻两点的连线函数

    现在对一个未知函数值的取值在 x 1 , x 2 x_1,x_2 x1,x2之间的点 x c x_c xc进行插值,由于
    l 1 ( x ) = { x − x 0 x 1 − x 0 , x ∈ [ x 0 , x 1 ] x − x 2 x 1 − x 2 , x ∈ [ x 1 , x 2 ] 0 , x ∉ [ x 0 , x 2 ] l 2 ( x ) = { x − x 1 x 2 − x 1 , x ∈ [ x 1 , x 2 ] x − x 3 x 2 − x 3 , x ∈ [ x 2 , x 3 ] 0 , x ∉ [ x 1 , x 3 ] l_1(x)=\left\{ \begin{aligned} &\frac{x-x_0}{x_1-x_0}&,&x \in[x_0,x_1]\\ &\frac{x-x_2}{x_1-x_2}&,&x \in[x_1,x_2]\\ &0&,& x \notin [x_0,x_2] \end{aligned} \right.\quad l_2(x)=\left\{ \begin{aligned} &\frac{x-x_1}{x_2-x_1}&,&x \in[x_1,x_2]\\ &\frac{x-x_3}{x_2-x_3}&,&x \in[x_2,x_3]\\ &0&,& x \notin [x_1,x_3] \end{aligned} \right. l1(x)=x1x0xx0x1x2xx20,,,x[x0,x1]x[x1,x2]x/[x0,x2]l2(x)=x2x1xx1x2x3xx30,,,x[x1,x2]x[x2,x3]x/[x1,x3]

    所以 F ( x c ) = F 1 ( x c ) x c − x 2 x 1 − x 2 + F 2 ( x c ) x c − x 1 x 2 − x 1 F(x_c)=F_1(x_c)\frac{x_c-x_2}{x_1-x_2}+F_2(x_c)\frac{x_c-x_1}{x_2-x_1} F(xc)=F1(xc)x1x2xcx2+F2(xc)x2x1xcx1

    其中 x c − x 2 x 1 − x 2 和 x c − x 1 x 2 − x 1 \frac{x_c-x_2}{x_1-x_2}和\frac{x_c-x_1}{x_2-x_1} x1x2xcx2x2x1xcx1分别是 x c x_c xc x 2 x_2 x2的距离在 x 1 , x 2 x_1,x_2 x1,x2之间距离的比例, x c x_c xc x 1 x_1 x1的距离在 x 1 , x 2 x_1,x_2 x1,x2之间距离的比例

    F 1 ( x c ) F_1(x_c) F1(xc)则是 x 1 , x 2 x_1,x_2 x1,x2连线函数上 x c x_c xc对应函数值,即下图图示中的a

    F 2 ( x c ) F_2(x_c) F2(xc)则是 x 2 , x 3 x_2,x_3 x2,x3连线函数上 x c x_c xc对应函数值,即下图图示中的b
    在这里插入图片描述

    n越大,分段越多,插值越好,误差越小

    缺点:分段线性插值不能保证在节点处的插值函数的导数的连续性,即不光滑(看后面的示例)

    拉格朗日插值的局限

    从定义式可看到,每个基函数 l i ( x ) l_i(x) li(x)都需要利用所有插值节点来计算,所以一旦新增一个数据点,所有的基函数全部需要重新计算,所以计算量上没有优势,下面介绍的牛顿插值法在这一点上就没有这个缺点
    在这里插入图片描述

    示例

    在这里插入图片描述

    人口数据给了30个已知数据点,属于是高次插值,可以看出效果并不理想
    代码:

    clear all 
    clc
    % 人口数据
    T=1:30;  % +1970=年份
    P=[33815 33981 34004 34165 34212 34327 34344 34458 34498 34476 34483 34488 34513 34497 34511 34520 34507 34509 34521 34513 34515 34517 34519 34519 34521 34521 34523 34525 34525 34527];%人口
    %plot(T,P)% 折线
    plot(T,P,'r o')% 只画数据点
    hold on
    x1=1.5:30.5; 
    y1=Lagrange(T,P,x1);
    plot(x1,y1,'b +')
    axis([0 30 33600 34800])
    title('拉格朗日插值&人口预测模型')
    
    function y1=Lagrange(x,y,x1);
    % x   插值节点
    % y   被插值函数在节点x处的函数值
    % x1  待插值节点
    % y1  待插值节点的插值结果
    n=length(x); % 数据点个数
    m=length(x1); % 待插值点数
    for k=1:m
        z=x1(k);
        sum=0.0;
        for i=1:n
            prod = 1.0;
            for j=1:n
                if j~=i
                    prod = prod * (z - x(j)) / (x(i) - x(j));
                end
            end
            sum = sum + prod * y(i);
        end
        y1(k) = sum;
    end              
    

    1.1.2 牛顿多项式插值法(差商)

    插值的思想和目的还是一样的

    牛顿插值法的公式是另一种n次插值多项式的构造形式,然而它却克服了拉格朗日插值多项式的缺陷,它的一个显著优势就是每当增加一个插值节点,只要在原牛顿插值公式中增加一项就可形成高一次的插值公式。此外,如果在实际应用中遇到等距分布的插值节点,牛顿插值公式就能得到进一步的简化,从而得到等距节点的插值公式,这样为缩短实际运算时间做出了很大的贡献。
    在这里插入图片描述

    差商

    一个挺有意思的东西
    在这里插入图片描述在这里插入图片描述##### 差商矩阵
    理解差商矩阵是理解牛顿插值代码的灵魂之门

    (以5个插值数据点为例,则差商矩阵为5*5)

    matlab中索引从1开始,我就保持一致,不从0开始了

    y(1)
    y(2) y ( 2 ) − y ( 1 ) x ( 2 ) − x ( 1 ) \frac{y(2)-y(1)}{x(2)-x(1)} x(2)x(1)y(2)y(1), i.e. f [ x 1 , x 2 ] f[x_1,x_2] f[x1,x2]
    y(3) y ( 3 ) − y ( 2 ) x ( 3 ) − x ( 2 ) \frac{y(3)-y(2)}{x(3)-x(2)} x(3)x(2)y(3)y(2), i.e. f [ x 2 , x 3 ] f[x_2,x_3] f[x2,x3] f [ x 2 , x 3 ] − f [ x 1 , x 2 ] x ( 3 ) − x ( 1 ) , i . e . f [ x 1 , x 2 , x 3 ] \frac{f[x_2,x_3]-f[x_1,x_2]}{x(3)-x(1)},i.e. f[x_1,x_2,x_3] x(3)x(1)f[x2,x3]f[x1,x2]i.e.f[x1,x2,x3]
    y(4) y ( 4 ) − y ( 3 ) x ( 4 ) − x ( 3 ) \frac{y(4)-y(3)}{x(4)-x(3)} x(4)x(3)y(4)y(3), i.e. f [ x 3 , x 4 ] f[x_3,x_4] f[x3,x4] f [ x 2 , x 3 , x 4 ] f[x_2,x_3,x_4] f[x2,x3,x4] f [ x 1 , x 2 , x 3 , x 4 ] f[x_1,x_2,x_3,x_4] f[x1,x2,x3,x4]
    y(5) y ( 5 ) − y ( 4 ) x ( 5 ) − x ( 4 ) \frac{y(5)-y(4)}{x(5)-x(4)} x(5)x(4)y(5)y(4), i.e. f [ x 4 , x 5 ] f[x_4,x_5] f[x4,x5] f [ x 3 , x 4 , x 5 f[x_3,x_4,x_5 f[x3,x4,x5] f [ x 2 , x 3 , x 4 , x 5 ] f[x_2,x_3,x_4,x_5] f[x2,x3,x4,x5] f [ x 1 , x 2 , x 3 , x 4 , x 5 ] f[x_1,x_2,x_3,x_4,x_5] f[x1,x2,x3,x4,x5]

    计算关键代码:

     A(i,j) = (A(i,j-1) - A(i-1,j-1))/(X(i)-X(i-j+1));
    

    在这里插入图片描述在这里插入图片描述

    示例

    此代码被撸了仨小时,参考的一位博主的非常好的代码,把差商用一个矩阵存放,还计算了插值余项,链接
    牛顿插值的理论可以在链接里看,讲的很详细

    %newton.m
    %求牛顿插值多项式、差商、插值及其误差估计的MATLAB主程序
    %输入的量:X是n+1个节点(x_i,y_i)(i = 1,2, ... , n+1)横坐标向量,Y是纵坐标向量,
    %x是以向量形式输入的m个插值点,M在[a,b]上满足|f~(n+1)(x)|≤M
    %注:f~(n+1)(x)表示f(x)的n+1阶导数
    %输出的量:向量y是向量x处的插值,误差限R,n次牛顿插值多项式L及其系数向量C,
    %差商的矩阵A
    function[y,R,A,C,L] = newton(X,Y,x,M)
    n = length(X);
    m = length(x);
    for t = 1 : m
        z = x(t);
        A = zeros(n,n);
        A(:,1) = Y';
        s = 0.0; p = 1.0; q1 = 1.0; c1 = 1.0;
            for j = 2 : n
                for i = j : n
                    A(i,j) = (A(i,j-1) - A(i-1,j-1))/(X(i)-X(i-j+1));
                end
                q1 = abs(q1*(z-X(j-1))); % 余项计算中的差值连乘项
                c1 = c1 * j; % 余项计算中的分母,阶乘
            end
            C = A(n, n); q1 = abs(q1*(z-X(n)));
            for k = (n-1):-1:1
                C = conv(C, poly(X(k))); 
                % poly函数把根转换为求得此根的多项式,X(k)是一个数,所以转换为一次多项式,得到其系数向量
                % conv函数的输入是多项式的系数向量时,相当于直接相乘
                d = length(C);
                C(d) = C(d) + A(k,k);%在最后一维,也就是常数项加上新的差商
            end
            y(t) = polyval(C,z); % 插值结果
            R(t) = M * q1 / c1; % 用M代替f~(n+1)(x),即f(x)的n+1阶导数,这样求得的余项比真实值略大
    end
    L = vpa(poly2sym(C),3);% vpa让输出结果系数为小数而非分数,3表示用3位小数
    

    测试脚本:

    clear all 
    clc
    X = [0 pi/6 pi/4 pi/3 pi/2];  
    Y = [0 0.5 0.7071 0.8660 1];  
    x = linspace(0,pi,50);  
    M = 1;  
    [y,R,A,C,L] = newton(X, Y, x, M); 
    y1 = sin(x);  
    errorbar(x,y,R,'.g')  
    hold on  
    plot(X, Y, 'or', x, y, '.k', x, y1, '-b');  
    legend('误差','样本点','牛顿插值估算','sin(x)'); 
    

    这里已知数据点少,多项式次数低,插值效果就比较好

    插值效果:
    在这里插入图片描述
    计算出的牛顿插值多项式,多项式系数,差商矩阵:

    >> C
    C =
        0.0290   -0.2048    0.0217    0.9955         0
    >> L
    L =
    0.029*x^4 - 0.205*x^3 + 0.0217*x^2 + 0.996*x
    >> A
    A =
             0         0         0         0         0
        0.5000    0.9549         0         0         0
        0.7071    0.7911   -0.2086         0         0
        0.8660    0.6070   -0.3516   -0.1365         0
        1.0000    0.2559   -0.4469   -0.0910    0.0290
    

    另一套数据的插值结果:
    在这里插入图片描述
    代码:

    X=[0.4 0.55 0.65 0.80 0.90 1.05];
    Y=[0.41075 0.57815 0.69675 0.88811 1.02652 1.25386];
    x=0.4:0.01:1.05;
    M = 1;  
    [y,R,A,C,L] = newton(X, Y, x, M);   
    errorbar(x,y,R,'.g')  
    hold on  
    plot(X, Y, 'or', x, y, '.k');  
    legend('误差','样本点','牛顿插值估算'); 
    

    用于人口数据插值:

    惨不忍睹的结果······

    所以并没有一种插值可以通用,完美拟合任何数据,要选择适合数据的插值方式

    这个人口数据呈明显的对数函数趋势,要用对数型函数去拟合(见这篇博客)。而**不能插值!!!**除非用于预测,某年数据还说得过去

    这里的数据点有30个,很多,牛顿插值得到了29次多项式, 共30个系数,高次多项式的函数图像很容易震荡,下图也可以看出插值数据的震荡性,这使得得到的插值函数并不光滑,而人口数据本身是很光滑的对数形状的曲线,这无法保证插值函数处处收敛于被插值函数(就是原数据点真正遵循的那个函数),这是高次多项式插值的通病。
    在这里插入图片描述

    代码:

    clear all 
    clc
    % 人口数据
    T=1:30;  % +1970=年份
    P=[33815 33981 34004 34165 34212 34327 34344 34458 34498 34476 34483 34488 34513 34497 34511 34520 34507 34509 34521 34513 34515 34517 34519 34519 34521 34521 34523 34525 34525 34527];%人口
    P=P/10000;
    M = 1;  
    x=1:.1:30;
    [y,R,A,C,L] = newton(T, P, x, M);   
    errorbar(x,y,R,'.g')  
    hold on  
    axis([1 30 3.35 3.50])
    plot(T, P, 'or', x, y, '.k');  
    legend('误差','样本点','牛顿插值估算'); 
    

    低次多项式插值(不易震荡)的其他示例


    保凹凸性即hermite插值
    在这里插入图片描述
    最差无疑了
    在这里插入图片描述
    分段线性插值的子插值函数都是一次的,从图像上也可以看得出来,一段一段的,很不光滑
    在这里插入图片描述
    三次样条的结果明显要光滑很多,子插值函数都是三次的,还考虑了一二阶导数在这里插入图片描述

    代码:

    clear all 
    clc
    % 人口数据
    T=1:30;  % +1970=年份
    P=[33815 33981 34004 34165 34212 34327 34344 34458 34498 34476 34483 34488 34513 34497 34511 34520 34507 34509 34521 34513 34515 34517 34519 34519 34521 34521 34523 34525 34525 34527];%人口
    P=P/10000;  
    x=1:.1:30;
    %y=interp1(T,P,x,'*nearest') ; 
    %y=interp1(T,P,x,'*linear') ; 
    %y=interp1(T,P,x,'*spline') ; 
    y=interp1(T,P,x,'*cubic') ; 
    axis([1 30 3.35 3.50])
    plot(T, P, 'or', x, y, '.k');  
    legend('样本点','插值估算','Location','southeast'); 
    %title('最近项分段线性插值')
    %title('线性  分段线性插值')
    % title('逐段三次样条  分段线性插值')
    title('保凹凸性3次插值 分段线性插值')
    
    
    % 'north' 坐标轴中的顶部
    % 'south' 坐标轴中的底部
    % 'east' 坐标轴中的右侧区域
    % 'west' 坐标轴中的左侧区域
    % 'northeast' 坐标轴中的右上角(二维坐标轴的默认值)
    % 'northwest' 坐标轴中的左上角
    % 'southeast' 坐标轴中的右下角
    % 'southwest' 坐标轴中的左下角
    % 'northoutside' 坐标轴的上方
    % 'southoutside' 坐标轴的下方
    % 'eastoutside' 到坐标轴的右侧
    % 'westoutside' 到坐标轴的左侧
    % 'northeastoutside' 坐标轴外的右上角(三维坐标轴的默认值)
    % 'northwestoutside' 坐标轴外的左上角
    % 'southeastoutside' 坐标轴外的右下角
    % 'southwestoutside' 坐标轴外的左下角
    % 'best' 坐标轴内与绘图数据冲突最少的地方
    % 'bestoutside' 到坐标轴的右侧
    % 'none' 由 Position 属性决定。可使用 Position 属性在自定义位置显示图例。
    

    尝试用在样本数据点少的情况:
    对sin函数插值(前面用牛顿插值取得了不错的效果),可以看到和真实sin函数在后面差距越来越大。
    在这里插入图片描述
    在这里插入图片描述
    这时候最近项线性插值是完全不能用的:

    可以看出最近项分段线性插值只能用于数据点非常非常多的情况
    在这里插入图片描述
    linear分段线性插值也不适合数据点少:
    在这里插入图片描述

    1.1.3 拉格朗日插值和牛顿插值的对比

    不同

    前者不具有继承性导致计算量大
    后者有继承性,数据点增多只需要在原来的多项式基础上增加几项,计算量要小一些
    在这里插入图片描述

    相同

    都存在龙格效应
    都没有考虑导数也要相同,所以可能不能模拟出被插值函数的形态,所以就引入了后面的考虑导数的插值法
    在这里插入图片描述

    1.2 考虑导数的插值方法

    1.2.1Hermite插值

    在这里插入图片描述

    分段三次hermite插值

    有前途,就用这个了,比分段线性插值好呀
    在这里插入图片描述

    优点:关于插值函数和被插函数的贴合程度,埃尔米特插值比多项式的好。
    缺点:埃尔米特插值只有在被插值函数在插值节点处的函数值和导数值已知时才可以使用,而这在实际问题中是无法实现的,因为在一般情况下我们是不可能也没必要知道函数在插值节点处的导数值。因此成为能否运用埃尔米特插值的一个重要因素就是:我们知不知道插值函数在节点处的导数值。

    1.2.2 三次样条插值

    定义

    每个子区间的多项式函数是三次的,所以名字里有“三次”,就像前面拉格朗日插值的分段二次插值的每个子区间的插值函数是二次的一样

    整个插值函数曲线由分段三次曲线并接而成,并且在连接点也就是样点上必须要二阶连续可导
    在这里插入图片描述

    n+1个数据点,共有n个分段多项式
    相邻两个子多项式之间必须在函数值,一阶导数,二阶导数上都相等
    在这里插入图片描述
    对比发现,分段三次hermite插值只考虑了一阶导数,分段三次样条考虑了一阶导数和二阶导数,所以后者得到的结果更加精确

    示例

    在这里插入图片描述
    代码

    x = -pi:pi;
    y = sin(x);
    x_ = -pi:.1:pi;
    p1 = pchip(x,y,x_); % 分段三次Hermite插值
    p2 = spline(x,y,x_); % 分段三次样条插值
    plot(x,y,'ko',x_,p1,'r-',x_,p2,'b-')
    legend('插值节点','分段三次Hermite插值','分段三次样条插值','location','southeast')
    

    可见,三次样条的结果更加光滑,也更接近sin函数本身的图像

    用插值进行数据预测

    在这里插入图片描述
    代码

    % 用插值进行预测
    % 2009-2018 10年的人口数据
    population=[133126,133770,134413,135069,135738,136427,137122,137866,138639, 139538];
    year = 2009:2018; 
    % 预测 2019-2021三年的人口数据
    p1 = pchip(year, population, 2019:2021);
    p2 = spline(year, population, 2019:2021);
    plot(year, population,'bo',2019:2021,p1,'r*-',2019:2021,p2,'bx-')
    legend('原始数据','三次hermite插值预测','三次样条插值预测','location','southeast')
    

    总结(必看+1)

    1. 拉格朗日插值、牛顿插值、分段线性插值、分段三次Hermite插值和分段三次样条插值具有不同的优势和适用范围,对于方法的选择至关重要,我们需要对它们进行差异化的了解与认知。

    2. 分段插值函数比如分段线性插值、3次分段hermite和三次分段样条插值函数在每个单元区间上收敛性强,数值稳定性好且易于计算机编程实现

    3. 一般来说,不要直接使用多项式插值,尤其是数据点较多(n>7)时,高次多项式震荡很厉害;
      分段多项式插值可以避免震荡带来的误差;
      分段多项式插值又不如分段hermite,因为后者考虑了一阶导数;
      分段hermite又不如分段三次样条,因为后者除了一阶导数,还考虑了二阶导数,更好地保持数据的光滑性和连续性,减少信息量的损失

    4. 回顾上面的所有示例图,一般来说,使用逐段三次样条可以获得最优插值结果,由于我们不知道数据本身的分布,所以分段三次hermite也可以尝试。其他的当然也是可以使用的,根据自己的问题自己判断吧

    一般来说,使用逐段三次样条可以获得最优插值结果
    一般来说,使用逐段三次样条可以获得最优插值结果
    一般来说,使用逐段三次样条可以获得最优插值结果

    还是要分情况的,最近发现,数据量很大时,比如一条数据向量三千多维,给一千个索引,得到的插值结果在两端会出现强烈的龙格现象,很大的甩尾swings,也是一个棘手的山芋

    如果给的索引不是从1开始,linear 和nearest竟然不给两端插值!!!!!我也是服气,整出NAN玩儿

    K>> interp1([3 5 7],[4 3 6],1:9,'linear')
    
    ans =
    
           NaN       NaN    4.0000    3.5000    3.0000    4.5000    6.0000       NaN       NaN
    
    K>> interp1([3 5 7],[4 3 6],1:9,'nearest')
    
    ans =
    
       NaN   NaN     4     3     3     6     6   NaN   NaN
    
    K>> interp1([3 5 7],[4 3 6],1:9,'cubic')
    警告: 在以后的版本中将会更改 INTERP1(...,'CUBIC')。请改用 INTERP1(...,'PCHIP')。 
    > In interp1>sanitycheckmethod (line 265)
      In interp1>parseinputs (line 401)
      In interp1 (line 80)
      In feature_extract (line 393) 
    
    ans =
    
       11.0000    6.3750    4.0000    3.1250    3.0000    3.8750    6.0000    8.6250   11.0000
    
    K>> interp1([3 5 7],[4 3 6],1:9,'spline')
    
    ans =
    
         9     6     4     3     3     4     6     9    13
    
    展开全文
  • 计算方法中三种基本插值方法,拉格朗日插值分段线性插值;样条插值
  • 文章目录创建项目拉格朗日插值 创建项目 使用VS2019创建C++控制台应用 C++菜鸟教程 拉格朗日插值 https://blog.csdn.net/u013871100/article/details/41249285/ ...拉格朗日插值 牛顿插值 分段线性插值
  • 已知 sin(0.32)=0.314567,sin(0.34)=0.333487,sin(0.36)=0.352274, sin(0.38)=0.370920。 请采用线性插值、二次插值、三次插值分别计算 sin(0.35)的值。
  • 1、数值分析作业姓名:虞驰程题目:函数:fx=11+x2在-5,5上,取n=10,对其进行分段线性插值和拉格朗日插值,在Matlab中实现且绘图。Matlab实现:首先定义函数f,在Matlab中用function.m文件编...
  • 1、拉格朗日插值法考虑全局信息的比较经典的插值方法,编程简单,计算量大。#coding=utf-8 from matplotlib import pyplot as pltdef Lg(data,testdata): predict=0 data_x=[data[i][0] for i in range(len(data))...

空空如也

空空如也

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

拉格朗日插值和分段线性插值