实验二 插值与拟合

计算方法课程实验报告

实验名称 插值与拟合

班级 动创新13 姓名 封敏丽 学号 201302400104

序号 教师

赵美玲 地点

数学实验中心

评分

一、 实验目的

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

② 了解整体插值的局限性及分段插值的基本思想。 ③ 掌握最小二乘法拟合的基本原理和方法； ④ 培养运用计算机模拟解决问题的能力。

二、用文字或图表记录实验过程和结果

3.4.1 多项式插值

1． 给定sin110.190809,sin120.207912,sin130.224951,o o o ===构造插值多项式计算'

sin1130o

。

（1） 编程实现拉格朗日插值，并计算结果。 function f=Language(x,y,x0)

%求已知数据点的拉格朗日插值多项式 %已知数据点的x 坐标向量: x %已知数据点的y 坐标向量: y %插值点的x 坐标: x0

%求得的拉格朗日插值多项式或在x0处的插值: f x=[11 12 13];

y=[0.190809 0.207912 0.224951]; x0=11.5; 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);

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

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

运行结果如下：

所以'

sin1130o

=0.199368

通过对以上的程序修改可以得到插值函数：

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

查表的结果是sin11.5＝0.19936793441719718108604

由此可见，计算结果与查表的结果已经是什么相近了。

拉格朗日插值法的流程图如下：

Y

N

N

Y

开始

输入插值点（Xi ，Yi ），以

及所求值t ，i ：0到n-1

i=j ？ temp=temp*（t-Xi ） temp=temp*（Xi-Xj ） j=j+1

result=result+temp

初始化result=0

i=0 临时变量temp=Yi i=n 输出结果result 结束 i=i+1

2.区间[]5,5-作等距划分：10

1(0,1,,),k x kh k n h n

=-+==

　，以k x （0,1,,k n =）

为节点对函数2

1

()1f x x =

+进行插值逼近。（分别取5,10,20n =） （1） 用多项式插值对()f x 进行逼近，并在同一坐标系下作出函数的图形，进行比较。

写出插值函数对()f x 的逼近程度与节点个数的关系，并分析

本题重点在于画图，所以使用相对而言绘图很方便的Mathmatic 软件进行绘图： 先利用原函数取点，然后分别以步长为2,1,0.5，得到拟合函数,计算结果如下：

原函数图形：

N=5时：

N=10时：

N=20时：

通过Show语句进行图形组合得到：

结果分析：高次插值稳定性差，而低次插值对于较大区间逼近精度又不够，而且，随着节点的加密，采用高次插值，插值函数两端会发生激烈震荡。解决这一矛盾的有效方法就是采用分段低次代数插值。

f x进行逼近，在同一坐标下画出图形，观察分（2）试用分段插值（任意选取）对()

f x的逼近程度与节点个数的关系。

段插值函数对()

用Mathmatic进行分段插值表示不会，故而用Matlab对分段插值进行了运行结果如下：

结果分析：通过采用分段线性插值，发现随着节点数的增多，插值计算结果的误差越来越小，而且分段线性插值的计算简单，曲线连续和一致收敛，但是不具有光滑性。

3.4.2 数据拟合

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

x 1 2 3 4 5

i

y 4 4.5 6 8 8.5

i

（1）编程实现最小二乘算法，并画出其拟合曲线

（2）求出其平方误差。

本题第一小题和第二小题同时用Matlab运行时，我们可以先建立一个M文件，代码如下：N=input('please put in how many times the power will you overfit:');

M=input('how many couples of statistics are there in the table:');

%读入数据文件

S=[1 2 3 4 5;4 4.5 6 8 8.5],

%显示数据文件，确保正确输入

disp('S(x,y)=');

disp(S);

C=zeros(N+1,M);

for i=1:N+1

for j=1:M

if S(1,j)==0

C(i,j)=0;

else

C(i,j)=S(1,j).^(i-1);

end

end

end

%建立法方程组

A=C*C';

Y=zeros(M,1);

for i=1:M

Y(i,1)=S(2,i);

end

b=C*Y;

%用列主元高斯消元法解法方程组

A=[A,b];

for i=1:N+1

max=abs(A(i,i));

for j=j+1:N+1

if abs(A(j,i))>max

flag=j;

max=A(j,i);

end

end

for kh=i+1:N+1

m=-A(kh,i)/A(i,i);

A(kh,i)=0;

for kl=i+1:N+2

A(kh,kl)=A(kh,kl)+m*A(i,kl);

end

end

end

X=zeros(N+1,1);

for i=N+1:-1:1

for j=i-1:-1:1

m=-A(j,i)/A(i,i);

A(j,N+2)=A(j,N+2)+m*A(i,N+2);

end

X(i,1)=A(i,N+2)/A(i,i);

end

disp(X);

%根据系数求的待定曲线

syms x;

expr=0;

for i=1:N+1

expr=expr+X(i,1)*x.^(i-1);

end

%输出得到的曲线表达式

disp(expr);

%计算偏差

bias=zeros(M,1);

for j=1:M

for i=1:N+1

bias(j,1)=bias(j,1)+X(i,1)*S(1,j)^(i-1);

end

bias(j,1)=bias(j,1)-S(2,j);

end

%计算平方差

rms=0;

for i=1:M

rms=rms+bias(i,1)^2;

end

disp('the bias is:');disp(rms);

%制图

a=S(1,1):0.01:S(1,M);

y=subs(expr,x,a);

plot(a,y);

hold on;

grid on;

for i=1:M

x=S(1,i);

y=S(2,i);

plot(x,y,'*');

hold on;

end

调用M函数

这个时候运行结果如下：

同时得到曲线图形如下：

其中得到最小二乘解，a=1.25 b=2.45 所以拟合曲线是y=1.25x+2.45 同时所得的平方误差是0.6750

求最小二乘解的流程图：

False

False

开始 输入拟合的次数N ，数据点的个数M 0

C

T

==

(Y 为数据点中yi 组成的列向量)

输出表达式

∑=+=N

i i

x i a r 0)1(exp

对数据点S(1，1）到S （1，M ）区间按0.01步长作图

计算

∑=

M rms

i

1

2

δδ

；

δδ

i

i max 1max

<=<==。

结束

i+=1 j+=1

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

i

0 1 2 3 4 5 6 7 8 9 10 i

x

2 3 4 7 8

10

11

14

16

18

19

i

y 106.42 108.2 109.5 110 109.93 110.49 110.59 110.6 110.76 111 111.2

(1) 求以上数据形如2012()y x c c x c x =++的拟合曲线，及其平方误差。

分析：本题和上面一题是类似的，都是用多项式来拟合所给的数据点，所以将不再继续呈现

求最小二乘解的原理。在此直接用Matlab 自带函数来拟合图形并求得平方误差，减少工作量。 代码及运行结果如下：

所得到拟合曲线图形如下：

它的拟合函数是：x x x y 2

020454.062635.02927.106)(-+=

它的平方误差是：2.7796

(2) 求以上数据形如()b

x y x ae =的拟合曲线，及其平方误差。

运用Matlab 进行求解，代码如下： >> x=[2,3,4,7,8,10,11,14,16,18,19];

>> y=[106.42,108.2,109.5,110,109.93,110.49,110.59,110.6,110.76,111,111.2]; >> b0=[111.849816,-0.099527]; %初始参数值 >> fun=inline('b(1)*exp(b(2)./x)','b','x'); >> [b,r,j]=nlinfit(x,y,fun,b0); >> b %最佳参数 >> R=sum(r.^2)%求平方误差

>> x1=[2,3,4,7,8,10,11,14,16,18,19]; >> y1=111.4922*exp((-0.0902)./x);

>> plot(x,y,'*',x1,y1,'-b')%做出原数据离散的点和新拟合的曲线图形

其中要求得b0的初始参数值是任取图中的两组数据，带入到拟合函数中求得。可用Matlab 中的Solve 语句求得，如下所示：

如果得到的结果是复数则只取实数部分作为b0的初始参数。

最终运行结果如下图所示：

所画图形如下：

它的拟合函数是：

e

x

x y 0902.04922.111)(-

=

它的平方误差是：0.4719

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

答：从（1）（2）的图形可以看出，（2）的拟合函数比（1）更加的贴近数据点，其拟合度更高，同时从平方误差中也可以看出，（2）的平方误差更小，这也说明了（2）的拟合度更高。对于数据拟合不仅仅是可以用多项式来对其进行拟合，有时候选择指数函数，对数函数，幂函数或者三角函数进行非线性拟合，可达到更加高度的图形吻合，从而减小拟合所带来的误差。所以在对数据点进行拟合前，我们可以先画出它的离散图形，对图形进行定性的分析选择贴切的函数，然后再进行拟合，可以提高拟合精度。

三、练习与思考题分析解答

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

答：整体插值的过程中，若有无效数字则整体插值后插值后曲线的平方误差会比较大，即在该数据附近插值曲线的震动幅度较大。在插值处理前，应对原始数据进行一定的筛选，剔除无效数据。

2、简述数据拟合与插值的异同。 答：（1）不同点：拟合是指已知某函数的若干离散函数值)2,1(fn f f ???，通过调整该函数中若干待定系数),(

2

1

λλλn

f ???，使得该函数与已知点集的差别最小。而插值是指已知某函

数在若干离散点上的函数值或者导数信息，通过求解该函数中待定形式的插值函数以及待定系数，使得该函数在给定离散点上满足约束。

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

四、本次实验的重点难点分析

本次上机中，学习了多项式插值法的基本思路和步骤，同时了解整体插值的局限性及分段插值的基本思想，熟悉了最小二乘法拟合的基本原理和方法。

在本次上机中，我觉得对于插值逼近的理解比较吃力，在用matlab进行插值逼近时，最终失败了，没有用matlab画出图形，所以不得不用对于画图更方便的mathmatic软件。其中，对于分段插值不是非常的理解，这也是导致不能成功画出图形的原因。

对于最小二乘法的拟合，尤其是最后一题的非线性拟合，在此提出，可以对拟合的函数两边先取对数，从而转化成线性拟合函数。由于matlab有自带的非线性拟合函数，所以就直接用自带的函数拟合了，而没有先转化成线性函数，在此进行说明。