插值法的程序实现

更新时间：2023-12-18 11:33:01 阅读量：教育文库文档下载

说明：文章内容仅供预览，部分内容可能不全。下载后的文档，内容与下面显示的完全一致。下载之前请确认下面内容是否您想要的，是否完整无缺。

插值法的程序实现

一实验目的

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年的人口表： 1920 1930 1940 1950 1960 1970 年份 105711 123203 131669 150697 179323 203212 人口 (千人) 1.用表中数据构造一个5次拉格朗日插值多项式，并用此估计1910，1965和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 #include #include #include

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 Lagrange(double x,double p[6],double q[6]) { }

void main(){

int i;

double x[6],y[6],l=1.0,h=1.0; double f[6][5];

cout<<\请输入x的取值:\return L(x,p,q); double l=1.0,h=0.0; for(int i=0;i<6;i++) { }

for( i=1;i<6;i++) { }

for(int j=1;j<6;j++)

for(int k=0;k<=j;k++) ; f[i][0]=q[i];

}

for(i=0;i<6;i++)

cin>>x[i];

cout<<\请依次输入x对应的y的取值:\for(i=0;i<6;i++)

cin>>y[i];

cout<<\多插值项式:\

cout<<\∑(yi*∏[(x-xj)/(xi-xj)]),其中(i=0,1,2,3,4,5;j=0,1,2,3,4,5;j≠cout<

cout<<\插值多项式在1910，1965，2002处的函数值:\cout<<\当x值为1910时\cout<<\当x值为1965时\cout<<\当x值为2002时\cout<

i)\

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