上机1成中涛080902.docx
《上机1成中涛080902.docx》由会员分享，可在线阅读，更多相关《上机1成中涛080902.docx（15页珍藏版）》请在冰点文库上搜索。
上机1成中涛080902
插值
1.用系数矩阵求方程组的方法求得拉氏插值多项式的方法
由给定的变量与自变量函数值表
-1
1
2
5
-7
7
-4
35
求次数不超过3的拉氏插值多项式。
在MATLAB中，输入命令
formatrat%使数据格式采用有理数形式
x=[-1,1,2,5]’,y=[-7,7,-4,35]’%输入列向量
A=[x.^0,x.^1,x.^2,x.^3]%构造出线性方程组（2.2）的系数矩阵
A\y
ans=[10,5,-10,2]’
结果：
y=
-7
7
-4
35
A=
1-11-1
1111
1248
1525125
ans=
10
5
-10
2
（运行文件名统一取为ff1.m）,其结果表明
2．拉格朗日插值函数的计算机实现
编制计算拉格朗日插值的M文件：
以下是拉格朗日插值的名为y_lagrl的M文件：
functiony=y_lagr1（x0,y0,x）
n=length（x0）;m=length（x）;
fori=1:
m
z=x（i）;
s=0.0;
fork=1:
n
p=1.0;
forj=1:
n
ifj~=k
p=p*（z-x0（j））/（x0（k）-x0（j））;
end
end
s=p*y0（k）+s;
end
y（i）=s;
end
运用此文件到具体插值问题中的例子：
选择函数y=exp（-x2）（-2≤x≤2）,在n个节点上（n不要太大，如5~11）用拉格朗日插值方法，计算m个插值点的函数值（m要适中，如50~100），并画出插值图形。
n=7;m=61;x=-2:
4/（m-1）:
2;
y=exp（-x.^2）;
z=0*x;
x0=-2:
4/（n-1）:
2;
y0=exp（-x0.^2）;
y1=y_lagr1（x0,y0,x）;
plot（x,y,’r’,x,y1,'b:
'）；
（文件名取为ff2.m）运行后结果：
3．龙格现象
在MATLAB中，通过编程可以观察到这种“龙格现象”。
%步1定义名为f.m的函数文件
functiony=f（x）
y=1./（1+x.^2）;
%步2编写计算拉格朗日插值多项式值并作图程序
%作被插函数图象
x=-5:
1:
5;
y=f（x）;
plot（x,y）;
%固化图形屏幕
Holdon
%计算下列节点处插值多项式的值
z=-5:
0.005:
5;
m=max（size（z））;
n=10;
fork=1:
m
L（k）=0;
fori=1:
n+1
t（i）=1;
forj=1:
n+1
ifj~=i
t（i）=t（i）*（z（k）-x（j））/（x（i）-x（j））;
end
end
L（k）=L（k）+t（i）*y（i）;
end
end
plot（z,L,'b'）
holdoff
在matlab中实现，（文件名为ff3.m）结果为：
4．分段线性插值函数逼近被插函数的计算机摸拟
作函数
在节点
的分段线性插值图象，并将它们进行比较，观察可发现，随着n的增大，两者吻合得越来越好，龙格现象并未发生。
fprintf（'\n'）%为屏幕清晰输出一空行
disp（'给出插值节点的个数n（正整数）'）
fprintf（'\n'）
n=input（'n='）;
fori=1:
200+1%作y=1/（1+x^2）的图象（多取了一些点作图）
u（i）=-5+0.05*（i-1）;
v（i）=1/（1+u（i）^2）;
end
plot（u，v，'r'）
holdon%固化图形屏幕
%给出插值条件
fori=1:
n+1
x（i）=-5+（10/n）*（i-1）;
y（i）=1/（1+x（i）^2）;
end
z=-5:
10/（2*n）:
5;%为作分段线性插值的图象，多取了一些点z
m=max（size（z））;
%计算分段线性插值函数在z处的值，并作图
fork=1:
m
i=1;
whilei<=n
ifz（k）>=x（i）&z（k）<=x（i+1）
Ih（k）=y（i）*（z（k）-x（i+1））/（x（i）-x（i+1））+y（i+1）*（z（k）-x（i））/（x（i+1）-x（i））;
break
else
i=i+1;
end
end
end
plot（z，Ih，'b'）
holdoff%释放固化的图形屏幕
在matlab中实现，输入不同的n值，比较其插值情况（文件名为ff4.m）
N=10:
N=15
N=20
N=30
从上面额图形可见，随着插值点的增加，拟合曲线和原曲线越来越接近。
需要大家自行完成的题目：
程序1下面给出美国从1920年到1970年的人口表：
年份
1920
1930
1940
1950
1960
1970
人口（千人）
105711
123203
131669
150697
179323
203212
用表中数据构造一个5次拉格朗日插值多项式，并用此估计1910，1965和2002年的人口。
在1910年的实际人口约为91772000，请判断插值计算得到的1965年和2002年的人口数据准确性是多少？
1）程序代码：
x0=1920:
10:
1970;
x=1910:
1:
2005;
m=length（x）;
n=length（x0）;
y0=[105711123203131669150693179323203212];
fori=1:
m
z=x（i）;
s=0.0;
fork=1:
n
p=1.0;
forj=1:
n
ifj~=k
p=p*（z-x0（j））/（x0（k）-x0（j））;
end
end
s=p*y0（k）+s;
end
y（i）=s;
end
plot（x,y）
2）曲线图
3）估计1910，1965和2002年的人口：
在上面的程序中，矩阵y已经包含了上述要求的年份的人口数，只要输出相应的值即可。
他们分别是y
（1）=31932,y（1965-1910+1）=1.9308e+005,y（2002-1910+1）=2.4499e+004.单位都是千人。
4）误差估计：
可以采用事后估计的方法来估计误差。
由于在1910年的实际人口约为91772000，这样实际上得到了7个插值点。
记从1910~1970年分别为x0,x1…x6.共7个插值点。
上面用x1~x6共6个点做了一次插值，结果用L1（x）表示；
若再用x0~x5做一次插值，结果若记作L2（x）.那么，误差值可以表示为：
F（X）-L（X）=|L1（X）-L2（X）|*[X-X0]/[X6-X0].
所以，还需要再做一次插值，方法与上面完全一样。
此次得到1965年的人口数为17836，，得到2002年的人口数为-396210（是个负值，可见是很不准确的）.
故1965年估计误差为：
|（17836-19308）|*55/60=1349（千人）。
这个误差还是可以接受的。
同样，2002年的误差也可以由此估计，但显然，该误差相当大。
这主要是由估计点在拟合点数据之外太远造成的。
一般地，只有在插值点范围内的估计值才较准确。
拟合
1．matlab自带最小二乘拟合程序
clc;
clear;
x=[1:
12];
fx=[1223434-13425233494523];
xi=[1:
0.2:
12];
y=polyfit（x,fx,11）%对于所给数据x，fx对，计算11阶拟合多项式的系数y
z=polyval（y,xi）;%求多项式值
plot（x,fx,'O',xi,z）
（文件名为ff5.m）运行结果：
2．自编最小二乘拟合程序
functiona=ployfit1（x,y）%最小二乘拟合
m=length（x）;
G=zeros（m）;
D=zeros（m,1）;
fori=1:
m
forj=1:
m
G（i,j）=sum（x.^（i+j-2））;
end
D（i,1）=sum（y.*（x.^（i-1）））;
end
G=G（1:
2,1:
2）
D=D（1:
2）
a=G^（-1）*D
例题：
已知一组试验数据如下，求它的拟合曲线
xi
1
2
3
4
5
fi
4
4.5
6
8
8.5
解答：
1）画出散点图，确定拟合曲线；作图如下：
不难发现，可以用三次多项式拟合。
设y=a0*x^3+a1*x^2+a2*x+a3;
2）用最小二乘求解，代码如下：
functionp=polyfit1（）
x=1:
5;
G=[1111;8421;27931;641641;1252551];
y=[4,4.5,6,8,8.5]';
p=[（G'*G）^（-1）*G'*y]';
t=0.5:
0.01:
5.5;
u=polyval（p,t）;
plot（t,u,x,y,'*'）;
结果如下：
ans=
-0.20831.9107-3.88106.2000
这四个值分别代表上面的a0a1a2a3
●● 上机作业格式
●● 题目
●● 所用方法
●● 实现算法
●● 输入、输出说明
●● 程序
●● 运行结果
●● 结果分析（选择项
光信1成中涛
二〇一〇年十一月二十日
升级会员 