插值法的程序实现

插值法的程序实现
一实验目的
1.熟悉Matlab编程；
2.学习插值方法及程序设计算法
二实验题目
分别用拉格朗日插值、牛顿插值、自然样条函数对1910、1965、2002的人口进行估算。

三实验原理与理论基础
1.拉格朗日插值算法设计
①利用已知条件得到xi，yi，i=0，1，2，…
②由Lk(x)=((x-x0)*…*(x-x(k-1))*(x-x(k+1))…(x-xn))/ ((xk-x0)*…*(xk-x(k-1))*(xk-x(k+1))…(xk-xn))得出Li(x);
③由Y=y1* L1(x)+…+yn*Ln(x)得出Y关于x的表达式。

④带值计算即可。

2.牛顿插值算法设计
①利用已知条件得到xi，yi，i=0，1，2，…
②利用差商公式
f[x0,…xk]=(f[x0,…,x(k-2),xk]-f[x0,…,x(k-1)])/(xk-x(k-1))各阶差商。

③利用牛顿插值公式
f(x)=f(x0)-f[x0,x1]*(x-x0)+…f[x0,x1,…xn]*(x-x0)*…(x-x(n-1)).
④带值计算即可。

3.自然样条函数算法设计
①利用已知条件得到xi，yi，i=0，1，2，…
②利用公式求出h(i) u(i) k(i);di;
③利用h(i) u(i) k(i);di解出向量M;
④将相关变量带入自然样条表达式中即可。

⑤带值计算即可。

四实验内容
（一）问题重述：
下面给出美国从1920年到1970年的人口表：
2002年的人口。

在1910年的实际人口约为91772000，请判断插值计算得到的1965年和2002
年的人口数据准确性是多少？
2用牛顿插值估计：
（1）1965年的人口数；
（2）2002年的人口数。

3用自然样条函数估计在1910，1965和2002年的人口数。

请比较以上三种方法所求值的效果。

那一种方法最优？
（二）实验代码：
1.①用matlab编写的拉格朗日插值M文件如下：
function [m,y]=cz5(x)
y1=105711;y2=123203;y3=131669;y4=150697;y5=179323;y6=203212;
x1=1920; x2=1930; x3=1940; x4=1950; x5=1960; x6=1970;
a1=((x-x2)*(x-x3)*(x-x4)*(x-x5)*(x-x6))/((x1-x2)*(x1-x3)*(x1-x4)*(x1-x5)*(x1-x6)); a2=((x-x1)*(x-x3)*(x-x4)*(x-x5)*(x-x6))/((x2-x1)*(x2-x3)*(x2-x4)*(x2-x5)*(x2-x6)); a3=((x-x1)*(x-x2)*(x-x4)*(x-x5)*(x-x6))/((x3-x1)*(x3-x2)*(x3-x4)*(x3-x5)*(x3-x6)); a4=((x-x1)*(x-x2)*(x-x3)*(x-x5)*(x-x6))/((x4-x1)*(x4-x2)*(x4-x3)*(x4-x5)*(x4-x6)); a5=((x-x2)*(x-x3)*(x-x4)*(x-x1)*(x-x6))/((x5-x1)*(x5-x3)*(x5-x2)*(x5-x4)*(x5-x6));
a6=((x-x2)*(x-x3)*(x-x4)*(x-x5)*(x-x1))/((x6-x1)*(x6-x2)*(x6-x3)*(x6-x4)*(x6-x5));
'实际值大约为：'
y=y1*a1+y2*a2+y3*a3+y4*a4+y5*a5+y6*a6;
if x==1910
m=91772;
'估计值大约为：'
else
m=(91772-31872)/31872*y+y;
'估计值大约为：'
End
②用c++编写程序如下：
#include<iostream.h>
#include<stdlib.h>
#include<stdio.h>
#include<math.h>
double H(double x,double a,double p[6]) //求得插值基函数
{
double h=1.0,f=1.0;
for(int i=0;i<6;i++)
{
if(fabs(a-p[i])<0.000001)
continue;
h*=(x-p[i]); //利用累乘求得插值基函数的分子部分
f*=(a-p[i]); //利用累乘求得插值基函数的分母部分
}
return h/f;
}
double L(double x,double p[6],double q[6])//求得拉格朗日插值多项式代入x取值时的函数值
{
double L=0;
for(int i=0;i<6;i++)
{
if(fabs(x-p[i])<0.000001)
return p[i];
L+=q[i]*H(x,p[i],p);// 利用累加实现L=∑(yi*∏[(x-xj)/(xi-xj)])
}
return L;
}
void N(double q[6],double f[6][5])
{
double l=1.0,h=0.0;
for(int i=0;i<6;i++)
{
f[i][0]=q[i];
}
for( i=1;i<6;i++)
{
for(int j=1;j<6;j++)
for(int k=0;k<=j;k++)
;
}
}
double Lagrange(double x,double p[6],double q[6])
{
return L(x,p,q);
}
void main(){
int i;
double x[6],y[6],l=1.0,h=1.0;
double f[6][5];
cout<<"请输入x的取值:"<<endl;
for(i=0;i<6;i++)
cin>>x[i];
cout<<"请依次输入x对应的y的取值:"<<endl;
for(i=0;i<6;i++)
cin>>y[i];
cout<<"Lagrange多插值项式:"<<endl;
cout<<"L=∑(yi*∏[(x-xj)/(xi-xj)]),其中(i=0,1,2,3,4,5;j=0,1,2,3,4,5;j≠
i)"<<endl;
cout<<endl;
cout<<"Lagrange插值多项式在1910，1965，2002处的函数值:"<<endl;
cout<<"当x值为1910时"<<Lagrange(1910,x,y)<<endl;
cout<<"当x值为1965时"<<Lagrange(1965,x,y)<<endl;
cout<<"当x值为2002时"<<Lagrange(2002,x,y)<<endl;
cout<<endl;
}
2.用matlab编写的牛顿插值M文件如下：
function [y]=niuden5(x)
y1=[105711 123203 131669 150697 179323 203212];
x0=[1920 1930 1940 1950 1960 1970];
for i=2:6
y2(i)=(y1(1)-y1(i))/(x0(1)-x0(i));
end
for i=3:6
y3(i)=(y2(2)-y2(i))/(x0(2)-x0(i));
end
for i=4:6
y4(i)=(y3(3)-y3(i))/(x0(3)-x0(i));
end
for i=5:6
y5(i)=(y4(4)-y4(i))/(x0(4)-x0(i));
end
y6(6)=(y5(5)-y5(6))/(x0(5)-x0(6));
y=y1(1)+y2(2)*(x-x0(1))+y3(3)*(x-x0(1))*(x-x0(2))+y4(4)*(x-x0(1))*(x-x0(
2))*(x-x0(3))+y5(5)*(x-x0(1))*(x-x0(2))*(x-x0(3))*(x-x0(4))+y6(6)*(x-x0(
1))*(x-x0(2))*(x-x0(3))*(x-x0(4))*(x-x0(5));
3.用matlab编写的自然样条函数M文件如下：
function []=yangtiaochazhi(x)
y1=[105711 123203 131669 150697 179323 203212];
x0=[1920 1930 1940 1950 1960 1970];
for i=1:5
h(i)=x0(i+1)-x0(i);
end
for i=1:4
u(i)=h(i)/(h(i)+h(i+1));
end
u(5)=0;
k(1)=0;
for i=2:5
k(i)=h(i)/(h(i-1)+h(i));
end
d(1)=0;d(6)=0;
for i=2:3
y21(i)=(y1(1)-y1(i))/(x0(1)-x0(i));
end
for i=3:4
y22(i)=(y1(2)-y1(i))/(x0(2)-x0(i));
end
for i=4:5
y23(i)=(y1(3)-y1(i))/(x0(3)-x0(i));
end
for i=5:6
y24(i)=(y1(4)-y1(i))/(x0(4)-x0(i));
end
d(1)=0;d(6)=0;
d(2)=(y21(2)-y21(3))/(x0(2)-x0(3));
d(3)=(y22(3)-y22(4))/(x0(3)-x0(4));
d(4)=(y23(4)-y23(5))/(x0(4)-x0(5));
d(5)=(y24(5)-y24(6))/(x0(5)-x0(6));
A=[2 k(1) 0 0 0 0;u(1) 2 k(2) 0 0 0;0 u(2) 2 k(3) 0 0;0 0 u(3) 2 k(4) 0;0 0 0 u(4)
2 k(5);0 0 0 0 u(5) 2];
B=[d(1);d(2);d(3);d(4);d(5);d(6)];
M=[];
M=A\B;
echo
if x<=1930
y=M(1)*((x0(2)-x)*(x0(2)-x)*(x0(2)-x))/(6*h(1))+M(2)*((x-x0(1))*(x-x0(1))*(x-x0(1)) )/(6*h(1))+(y1(1)-(M(1)*h(1)*h(1))/6)*((x0(2)-x)/h(1))+(y1(2)-(M(2)*h(1)*h(1))/6)*( (x-x0(1))/h(1));
elseif x<=1940
y=M(2)*((x0(3)-x)*(x0(3)-x)*(x0(3)-x))/(6*h(2))+M(3)*((x-x0(2))*(x-x0(2))*(x-x0(2)) )/(6*h(2))+(y1(2)-(M(2)*h(2)*h(2))/6)*((x0(3)-x)/h(2))+(y1(3)-(M(3)*h(2)*h(2))/6)*( (x-x0(2))/h(2));
elseif x<=1950
y=M(3)*((x0(4)-x)*(x0(4)-x)*(x0(4)-x))/(6*h(3))+M(4)*((x-x0(3))*(x-x0(3))*(x-x0(3)) )/(6*h(3))+(y1(3)-(M(3)*h(3)*h(3))/6)*((x0(4)-x)/h(3))+(y1(4)-(M(4)*h(3)*h(3))/6)*( (x-x0(3))/h(3));
elseif x<=1960
y=M(4)*((x0(5)-x)*(x0(5)-x)*(x0(5)-x))/(6*h(4))+M(5)*((x-x0(4))*(x-x0(4))*(x-x0(4)) )/(6*h(4))+(y1(4)-(M(4)*h(4)*h(4))/6)*((x0(5)-x)/h(4))+(y1(5)-(M(5)*h(4)*h(4))/6)*( (x-x0(4))/h(4));
else
y=M(5)*((x0(6)-x)*(x0(6)-x)*(x0(6)-x))/(6*h(5))+M(6)*((x-x0(5))*(x-x0(5))*(x-x0(5)) )/(6*h(5))+(y1(5)-(M(5)*h(5)*h(5))/6)*((x0(6)-x)/h(5))+(y1(6)-(M(6)*h(5)*h(5))/6)*( (x-x0(5))/h(5));
end
y
五实验结果
1. ①5次拉格朗日插值在matlab命令窗口中输入如下代码：
>>[m,y]=cz5(1910)
>>[m,y]=cz5(1965)
>>[m,y]=cz5(2002)
结果为：
ans =
实际值大约为：
ans =
估计值大约为：
m =
91772
y =
31872
ans =
实际值大约为：
ans =
估计值大约为：
m =
5.559574703016167e+005
y =
1.930815117187500e+005
ans =
实际值大约为：
ans =
估计值大约为：
m =
7.526371798548073e+004
y =
2.613874841600098e+004
②c++中的运行：
请输入x的取值:
1920
1930
1940
1950
1960
1970
请依次输入x对应的y的取值:
105711
123203
131669
150697
179323
203212
Lagrange多插值项式:
L=∑(yi*∏[(x-xj)/(xi-xj)]),其中(i=0,1,2,3,4,5;j=0,1,2,3,4,5;j≠i)
Lagrange插值多项式在1910，1965，2002处的函数值:
当x值为1910时31872
当x值为1965时193082
当x值为2002时26138.7
2. 牛顿插值在matlab命令窗口中输入如下代码：
>>niuden5(1965),niuden5(2002)
----------------------------------------------------------------------------------------------------------------- 结果为：
ans =
1.930815117187500e+005
ans =
2.613874841600133e+004
3. 样条插值在matlab命令窗口中输入如下代码：
>> yangtiaochazhi(1965),yangtiaochazhi(2002)
结果为：
y =
1.913745478468899e+005
y =
2.880973086315789e+005
六实验结果分析与小结
1. 5次拉格朗日插值与牛顿插值在估计人口时比较接近，但与真实值比起来误差较大，而用样条插值在区间内部估值比较准确，但在区间外估值误差相对
较大哦。

2.在matlab中数组名不能再用来作变量了哦。

忽略此处..。