实验二_插值与拟合

实验名称：插值与拟合                      

实验类型: 验证性实验

学　　时：2

3.1 实验环境

1操作系统：WindowsXP/Win7 

2编程环境：自定

3.2 实验目的

1掌握多项式插值法的基本思路和步骤；

2了解整体插值的局限性及分段插值的基本思想。

3掌握最小二乘法拟合的基本原理和方法；

4培养运用计算机模拟解决问题的能力。

3.3 实验原理和方法

3.3.1 多项式插值

若n次多项式在n+1个插值节点上满足插值条件：

则称这n＋1个n次多项式为插值节点上的n次插值基函数。

由于时，，故为的零点，从而可以设

由可得 

所以可得

。

若记，则有，从而有

。

上式是次数不超过的多项式，且满足所有的插值条件，称之为Lagrange插值多项式。

3.3.2 数据拟合

最小二乘法基本原理

已知数据对，求多项式

使得为最小，这就是多项式拟合的最小二乘法。

    最小二乘法的算法描述

 线性函数为例，拟合给定数据。

算法描述：

步骤1：输入值，及。

步骤2：建立法方程组，并求解。

步骤3：输出。

3.4 实验内容和步骤

3.4.1 多项式插值

1． 给定构造插值多项式计算。

（1）编程实现拉格朗日插值，并计算结果。

（2）将计算结果和查表结果进行比较。

2．区间作等距划分：，以（）为节点对函数进行插值逼近。（分别取）

（1）用多项式插值对进行逼近，并在同一坐标系下作出函数的图形，进行比较。写出插值函数对的逼近程度与节点个数的关系，并分析原因。

（2）试用分段插值（任意选取）对进行逼近，在同一坐标下画出图形，观察分段插值函数对的逼近程度与节点个数的关系。

解：1.

（1）拉格朗日插值程序代码：

function f = Language(x,y,x0)

syms t l;

if(length(x) == length(y))

    n = length(x);

else

    disp('x和y的维数不相等！');

    return;   %检错

end

h=sym(0);

for (i=1:n)

    l=sym(y(i));

    for(j=1:i-1)

        l=l*(t-x(j))/(x(i)-x(j));

    end;

    for(j=i+1:n)

        l=l*(t-x(j))/(x(i)-x(j));

    end;

    h=h+l;

end

simplify(h);

if(nargin == 3)

    f = subs (h,'t',x0); %计算插值点的函数值

else

    f=collect(h);

    f = vpa(f,6); %将插值多项式的系数化成6位精度的小数

end

x=[11,12,13];

y=[0.190809,0.207912,0.224951];

f=Language(x,y,11.5);

vpa(f,8)

程序结果：

（2）查表所得的精确解：0.199367934，所以插值所得到的结果精确到小数点后五位。

解：2.

（1）程序代码：

function y=lagrange(x0,y0,x);

n=length(x0);

m=length(x);

for i=1:m

    z=x(i);

    s=0.0;

    for k=1:n

        p=1.0;

        for j=1:n

            if j~=k

                p=p*(z-x0(j))/(x0(k)-x0(j));

            end

        end

        s=p*y0(k)+s;

    end

    y(i)=s;

end

x1=[-5:2:5];

y=1./(1+x1.^2);

x=[-5:0.001:5];

y1=L(x1,y);

x2=[-5:1:5];

y=1./(1+x2.^2);

y2=L(x2,y);

x3=[-5:0.5:5];

y=1./(1+x3.^2);

y3=L(x3,y);y=1./(1+x3.^2);

plot(x,subs(y1,'t',x),'r',x,subs(y2,'t',x),'g',x,subs(y3,'t',x),'b') 

hold on

plot(x,1./(1+x.*x),'m')

axis([-5 5 -0.5 2])

xlabel('x')

ylabel('y')

legend('n=5插值后的图','n=10插值后的图','n=20插值后的图','原图',0)

结果与原因分析：通过上图发现，高次插值稳定性差，而低次插值对于较大区间逼近精度又不够，而且，随着插值节点数目的增多，区域中间的确越来越逼近原函数，然而区域两端发生了较大震荡，即发生了Runge现象，这是由叉子余项的表达式所决定的。解决这一矛盾的有效方法就是采用分段低次代数插值。

（2）分段插值函数程序代码：

function output=L1

syms f x lx;

f=1/(1+x^2);

N=input('请输入插值节点数N=');

xx=-5:10/N:5;

ff=zeros(1,length(xx));

for i=1:(N+1)

    x=xx(i);

    ff(i)=eval(f);

end

M = -5:0.01:5;

output=zeros(1,length(M));

n = 1;

for i=2:N+1

    for x=-5:0.01:5

        if x=xx(i-1)

            lx(1)=ff(i-1)*(x-xx(i))/(xx(i-1)-xx(i));

            lx(2)=ff(i)*(x-xx(i-1))/(xx(i)-xx(i-1));

            output(n) = lx(1)+lx(2);

            n = n+1;

        end

    end

end

A =-5:0.01:5;

End

output=L1;

plot(A,output,'r');

hold on

output=L1;

plot(A,output,'g');

hold on

output=L1;

plot(A,output,'b');

hold on

ezplot(f,[-5,5])

xlabel('x')

ylabel('y')

legend('n=5分段差值后的图','n=10分段差值后的图','n=20分段插值后的图','原图',0)

请输入插值节点数N=5

请输入插值节点数N=10

请输入插值节点数N=20

程序结果：

结果分析：由图可以发现，插值节点数目越多，插值函数越接近原函数，函数没有出现Runge现象，这是因为分段插值将区间分为了n份。

3.4.2 数据拟合

1．已知一组数据如下，求它的线性拟合曲线。

（2）求出其平方误差。

2．已知一组数据如下，求其拟合曲线。

(2)求以上数据形如的拟合曲线，及其平方误差。

通过画出（1）（2）的图形，观察结果并结合其平方误差，写出你对数据拟合的认识。

解：1.

(1)最小二乘法程序代码

x=[1 2 3 4 5];

y=[4,4.5,6,8,8.5];

x1=0;

x12=0;

y1=0;

xy=0;

n=5;

for i=1:n

    x1=x1+x(i);

    x12=x12+x(i)^2;

    y1=y1+y(i);

    xy=xy+x(i)*y(i);

end

a0=(y1*x12-x1*xy)/(n*x12-x1*x1);

a1=(n*xy-x1*y1)/(n*x12-x1*x1);

yy =a0+a1*x

plot(x,y,'*')

hold on

plot(yy,'r')

axis([0,6,0,9])

结果：

（2）平方误差：

2.

（1）MATLAB程序代码：

xx=0;yy=0;x2=0;x3=0;x4=0;xy=0;y2=0;n=11;

x1=[2,3,4,7,8,10,11,14,16,18,19];

y1=[106.42,108.2,109.5,110,109.93,110.49,110.59,110.6,110.76,111,111.2];

for i=1:n

    xx=xx+x1(i);

    x2=x2+x1(i)^2;

    x3=x3+x1(i)^3;

    x4=x4+x1(i)^4;

    yy=yy+y1(i);

    xy=xy+x1(i)*y1(i);

    y2=y2+x1(i)*x1(i)*y1(i);

end

A=[n xx x2;xx x2 x3;x2 x3 x4];

B=[yy;xy;y2];

C=A\\B;

y2=C(1)+C(2)*x1+C(3)*x1.*x1;

plot(x1,y1,'*')

hold on

plot(x1,y2,'r')

结果：

平方误差：

程序代码：

syms a b

x=[2,3,4,7,8,10,11,14,16,18,19];

fi=a+b./x;

y1=[106.42,108.2,109.5,110,109.93,110.49,110.59,110.6,110.76,111,111.2];

y=log(y1);

fy=fi-y;

fy2=fy.^2;

J=sum(fy.^2);

Ja=diff(J,a); Jb=diff(J,b); 

Ja1=simple(Ja);Jb1=simple(Jb);

[a,b]=solve(Ja1,Jb1);

f=exp(a).*exp(b./x);

plot(x,y1,'*')

hold on

plot(x,f,'r')

xlabel('x')

ylabel('y')

legend('数值点','拟合后的图',0)

>> s=0;

for i=1:1:11

s=s+(f(i)-y1(i)).^2;

end 

vpa(s)

结果

平方误差 

拟合误差平方和为：0.47193209,通过计算比较发现，通过对数变换之后得到的拟合结果比进行二次函数拟合的结果要好很多。所以对于一些特殊的函数拟合，可以通过一些变换，进行拟合。  

3.5 练习思考

1、整体插值有何局限性？如何避免？

2、简述数据拟合与插值的异同。

解：1.线性插值的局限性：整体插值的过程中，若有无效数据则整体插值后插值曲线的平方误差会比较大，即在该数据附近插值曲线的震动幅度较大。在插值处理前，应对原始数据进行一定的筛选，剔除无效数据。对于函数高次插值稳定性差，而低次插值对于较大区间逼近精度又不够，而且，随着插值节点数目的增多，区域中间的确越来越逼近原函数，然而区域两端发生了较大震荡，即发生了Runge现象，这是由叉子余项的表达式所决定的。解决这一矛盾的有效方法就是采用分段低次代数插值。

2.①不同点：用插值的方法对一函数进行近似，要求所得到的插值多项式经过一直插值节点；在n比较大的情况下，插值多项式往往是高次多项式，这也就是容易出现震荡现象，即虽然在插值节点上没有误差，但在插值节点外插值误差变大，从整体上看，插值逼近效果将变得很差。

数据拟合是求一个简单的函数，例如是以的低次多项式，不要求通过已知的这些点，而是要求整体上尽量好的逼近原函数。

    ②相同点：通过已知一些离散点集M上的约束，求取一个定义在连续集合S(M包含于S)的未知连续函数，从而达到获取整体规律目的 。