插值与多项式逼近的数组计算方法实验讲解.doc

插值与多项式逼近的数组计算方法实验郑发进 2012042020022【摘要】计算机软件中经常要用到库函数，如，它们是用多项式逼近来计算的。虽然目前最先进的逼近方法是有理函数（即多项式的商），但多项式逼近理论更适于作为数值分析的入门课程。在已知数据具有高精度的情况下，通常用组合多项式来构造过给定数据点的多项式。构造组合多项式的方法有许多种，如线性方程求解、拉格朗日系数多项式以及构造牛顿多项式的方分和系数表。关 键 字 泰勒级数、拉格朗日插值法、牛顿插值法、帕德逼近一、 实验目的1. 通过具体实验，掌握泰勒级数、拉格朗日插值法、牛顿插值法、帕德逼近的编 程技巧。2. 比较各插值方法的优劣并掌握。二、 实验原理1. 泰勒级数在数学中，泰勒级数（英语：Taylor series）用无限项连加式级数来表示一个函数，这些相加的项由函数在某一点的导数求得。如果  在点x=x0具有任意阶导数，则幂级数称为  在点x0处的泰勒级数。在泰勒公式中，取x0=0，得到的级数称为麦克劳林级数。函数的麦克劳林级数是x的幂级数，那么这种展开是唯一的，且必然与的麦克劳林级数一致。2.拉格朗日插值法如对实践中的某个物理量进行观测，在若干个不同的地方得到相应的观测值，拉格朗日插值法可以找到一个多项式，其恰好在各个观测的点取到观测到的值。这样的多项式称为拉格朗日（插值）多项式。数学上来说，拉格朗日插值法可以给出一个恰好穿过二维平面上若干个已知点的多项式函数。在平面上有(x1,y1)(x2,y2).(xn,yn)共n个点，现作一条函数f（x）使其图像经过这n个点。作n个多项式pi(x),i=1,2,3.,n,使得最后可得 3. 牛顿插值法插值法利用函数f (x)在某区间中若干点的函数值，作出适当的特定函数，在这些点上取已知值，在区间的其他点上用这特定函数的值作为函数f (x)的近似值。如果这特定函数是多项式，就称它为插值多项式。利用插值基函数很容易得到拉格朗日插值多项式，公式结构紧凑，在理论分析中甚为方便，但当插值节点增减时全部插值基函数均要随之变化，整个公式也将发生变化， 这在实际计算中是很不方便的，为了克服这一缺点，提出了牛顿插值。牛顿插值通过求各阶差商，递推得到的一个公式：牛顿插值与拉格朗日插值具有唯一性。4. 帕德逼近它不仅与逼近论中其他许多方法有着密切的关系，而且在实际问题特别是许多物理问题中有着广泛的应用。设是在原点某邻域内收敛的、具有复系数的麦克劳林级数。欲确定一个有理函数，式中,使得前次方的系数为0，即使得 此处约定qk0（k>n）。虽然所求得的Pm(z)和Qn(z)不惟一，但是比式却总是惟一的。有理函数称为F(z)的(m,n)级帕德逼近,记为(m/n)。由(m/n)所形成的阵列称为帕德表。三、 实验内容1. P154.1:用plot命令，在同一幅图中绘制区间-1x1上的sin(x)，以及P5(x),P7(x)和P9(x)。其中： ·········································式（1） ····································式（2） ····························式（3）2. P171.2:下表（表一）给出了11月8号美国洛杉矶的一个郊区在5小时内的测量温度。（a）对表一中的数据构造一个拉格朗日插值多项式。（b）估计这5小时内的平均温度。（c）在同一坐标系中画出表中的数据和由（a）得到的多项式。讨论用（a）中的多项式计算平均温度可能产生的误差。 时间（下午） 华氏度 1 66 2 66 3 65 4 64 5 63 6 63 表13. P178.1:用牛顿插值多项式计算实验P171.2的内容。4. P194.1:比较对于函数的逼近： 泰勒多项式逼近： ····························式（4） 帕德逼近： ································式（5）（a）在同一坐标系中画出f(x)，T4(x)，R2,2(x)的曲线。（b）分别求出在区间-1，1上用T4(x)和R2,2(x)逼近f(x)的最大误差。5. P194.3:比较对于函数f(x)=tan(x)的逼近： 泰勒多项式逼近： ·······················式（6） 帕德逼近： ·····························式（7）（a）在同一坐标系中画出f(x)，T9(x)，R5,4(x)的曲线。（b）分别求出在区间-1，1上用T9(x)和R5,4(x)逼近f(x)的最大误差。四、 实验结果及分析1. P154.1： 实验描述： （1）plot绘图的原理为连续点绘图，只需输入一组等间距的坐标点即可完成； (2)坐标点的计算使用C+完成，计算完成后输入文件中； (3)绘图使用matlab的plot函数完成，具体方法为从文件中读取出坐标点，之后使用plot函数绘图。 实验结果： 表2 x及sin(x)及5,7,9阶泰勒展开公式计算结果xy=sin(x)y= P5(x)y= P7(x)y= P9(x)-100000000- 0.84147098-0.84166667-0.84146825-0.84147101-0.90000000-0.78332691-0.78342075-0.78332585-0.78332692-0.80000000-0.71735609-0.71739733-0.71735572-0.71735609-0.70000000-0.64421769-0.64423392-0.64421758-0.64421769-0.60000000-0.56464247-0.56464800-0.56464245-0.56464247-0.50000000-0.47942554-0.47942708-0.47942553-0.47942554-0.40000000-0.38941834-0.38941867-0.38941834-0.38941834-0.30000000-0.29552021-0.29552025-0.29552021-0.29552021-0.20000000-0.19866933-0.19866933-0.19866933-0.19866933-0.10000000-0.09983342-0.09983342-0.09983342-0.099833420.000000000.000000000.000000000.000000000.000000000.100000000.099833420.099833420.099833420.099833420.200000000.198669330.198669330.198669330.198669330.300000000.295520210.295520250.295520210.295520210.400000000.389418340.389418670.389418340.389418340.500000000.479425540.479427080.479425530.479425540.600000000.564642470.564648000.564642450.564642470.700000000.644217690.644233920.644217580.644217690.800000000.717356090.717397330.717355720.717356090.900000000.783326910.783420750.783325850.783326921.000000000.841470980.841666670.841468250.84147101 图1 y=sin(x)及其5,7,9阶泰勒展开函数 实验结论： （1）由表二可知，随着泰勒展开阶数的增加，Pn(x)越来越接近于原函数，当展开阶数n=9时，在误差delta=1e-7的精度要求下可以认为P9(x)与sin(x)完全拟合；（2）当y=sin(x)的泰勒展开函数Pn(x)的展开阶数为5,7,9时，其函数图像与原函数图像基本相符；（3）综上，当y=sin(x)泰勒展开到达5阶时，其五阶泰勒展开函数P5(x)便近似于原函数，可在范围内用于替代计算。2. P171.2: 实验描述：（1）拉格朗日插值多项式的公式为： ·········································式（8）以及 ·················式（9）（2）进行计算时，因为现有坐标点关于(3.5,y(3.5)点呈中心对称，故取y(3.5)的值近似替代五小时温度的平均值；（3）拉格朗日插值多项式及已知点的绘图使用matlab的plot函数完成，为完成绘图需将拉格朗日插值多项式及已知坐标点的坐标输入到文件中，以便于使用matlab读取。实验结果：(1)拟合后的6次拉格朗日插值多项式为： 式（10）(2)这五小时的平均温度为：64.5（华氏度）(3) 图2 温度的拉格朗日插值多项式拟合曲线图3. P178.1: 实验描述：（1）拉格朗日插值多项式的公式为： ··············式（11）其中为给出坐标点的各阶差商。（2）进行计算时，因为现有坐标点关于(3.5,y(3.5)点呈中心对称，故取y(3.5)的值近似替代五小时温度的平均值；（3）牛顿插值多项式及已知点的绘图使用matlab的plot函数完成，为完成绘图需将牛顿插值多项式及已知坐标点的坐标输入到文件中，以便于使用matlab读取。 实验结果：(1)拟合后的6次牛顿插值多项式为：·······································式(12)(2)这5小时的平均温度为：64.5华氏度(3)图三 温度的牛顿插值多项式拟合曲线图 实验结论：（1）当日下午1至6时的平均温度估计为64.5华氏度。（2）由于采样点的数据容量较小，平均温度可能存在误差。（3）对比拉格朗日插值与牛顿插值，牛顿插值的计算量更小，且当取5次点的时候，两者拟合精度都比较高。4.P194.1:实验描述：（1）实验画图使用matlab中的plot函数完成 (2）误差的最大值使用matlab中的max函数完成实验结果:（1） 图4 f(x)及其泰勒逼近T4(x)及帕德逼近R2,2(x)曲线图(2) 图5 泰勒逼近及帕德逼近误差曲线图泰勒逼近f(x)的最大误差：E4(x)=0.0099帕德逼近f(x)的最大误差：E2,2(x)=0.004实验结论：（1）在区间内，T4(x)与R2,2(x)与原函数y=ex大致完全拟合；（2）相较与T4(x)，R2,2(x)的误差在区间内更小，且随|x|的增大，误差也逐渐增大。5.P194.3:实验描述：（1）实验画图使用matlab中的plot函数完成（2）误差的最大值使用matlab中的max函数完成实验结果：（1） 图6 f(x)及其泰勒逼近T9(x)及帕德逼近R5,4(x)曲线图(2) 图7 泰勒逼近及帕德逼近误差曲线图泰勒逼近f(x)的最大误差：E4(x)=-0.0149帕德逼近f(x)的最大误差：E2,2(x)= -3.1725e-07实验结论：（1）在区间内，T4(x)与R2,2(x)与原函数y=tan(x)大致完全拟合；（2）相较与T9(x)，R5,4(x)的误差在区间内更小，且随|x|的增大，误差也逐渐增大;（3）与实验四相对比可得出，随着逼近次数的增加，计算精度也随之增加附件（代码）：1. P154.1:#include<stdlib.h>#include<iostream>#include<fstream>#include<iomanip>#include<cmath>using namespace std; struct X_FX /定义一个用于存放x,sin(x),P5(x),P7(x),P9(x)的数据类型double x;double y1;double y2;double y3;double y4;int main()int i;X_FX Fx21; /定义一个一维数组，此数组用于储存x,f(x),P5(x),P7(x),P9(x)double x=-1;i=0;double Fsinx(double);double Psinx(double,int);ofstream outfile;while(x<=1) /依次产生各x值,x=-1+0.1k,k=0,1,2.20Fxi.x=x;x=x+0.1;i+;for(i=0;i<=20;i+) Fxi.y1=sin(Fxi.x); /依次计算各x对应的sin(x)的值Fxi.y2=Psinx(Fxi.x,5); /依次计算各x对应的P5(x)的值Fxi.y3=Psinx(Fxi.x,7); /依次计算各x对应的P7(x)的值Fxi.y4=Psinx(Fxi.x,9); /依次计算各x对于的P9(x)的值outfile.open("154_1.txt");/输出计算结果，并将结果保存到文件“154_1.txt”中cout<<std:left<<setw(15)<<"x"<<std:left<<setw(15)<<"sin(x)"<<std:left<<setw(15)<<"P5(x)"<<std:left<<setw(15)<<"P7(x)"<<std:left<<setw(15)<<"P9(x)"<<endl;for(i=0;i<=20;i+)cout<<std:left<<setw(15)<<setiosflags(ios:fixed)<<setprecision(10)<<Fxi.x;cout<<std:left<<setw(15)<<setiosflags(ios:fixed)<<setprecision(10)<<Fxi.y1;cout<<std:left<<setw(15)<<setiosflags(ios:fixed)<<setprecision(10)<<Fxi.y2;cout<<std:left<<setw(15)<<setiosflags(ios:fixed)<<setprecision(10)<<Fxi.y3;cout<<std:left<<setw(15)<<setiosflags(ios:fixed)<<setprecision(10)<<Fxi.y4<<endl;outfile<<std:left<<setw(15)<<setiosflags(ios:fixed)<<setprecision(10)<<Fxi.x;outfile<<std:left<<setw(15)<<setiosflags(ios:fixed)<<setprecision(10)<<Fxi.y1;outfile<<std:left<<setw(15)<<setiosflags(ios:fixed)<<setprecision(10)<<Fxi.y2;outfile<<std:left<<setw(15)<<setiosflags(ios:fixed)<<setprecision(10)<<Fxi.y3;outfile<<std:left<<setw(15)<<setiosflags(ios:fixed)<<setprecision(10)<<Fxi.y4<<endl;outfile.close();system("pause");return 0;/此函数为用于计算y=sin(x)的解泰勒展开y=Pn(x)的函数，n为展开阶数/输x为y=Pn(x)的自变量x/输入n为y=Pn(x)的展开阶数n/输入y为y=Pn(x)的因变量ydouble Psinx(double x,int n)double y=0;int i,j=1;int inf(int);for(i=1;i<=n;i=i+2) y=y+j*pow(x,i)/inf(i);j=-j;return y;/此函数用于计算x的阶乘/输入x为要求阶乘x！的数x/输出y为对x做阶乘运算的结果int inf(int x)int i,y=1;for(i=1;i<=x;i+)y=y*i;return y;Matlab绘制曲线函数Val=load('154_1.txt'); %从文件中取出计算结果X=(Val(:,1)' %将x值保存到矩阵X中SinX=(Val(:,2)' %将y=sin(x)保存到矩阵SinX中P5X=(Val(:,3)' %将y=P5(x)保存到矩阵P5X中P7X=(Val(:,4)' %将y=P7(x)保存到矩阵P7X中P9X=(Val(:,5)' %将y=P9(x)保存到矩阵P9X中figure(1)subplot(2,2,1) %绘制y=sin(x)曲线plot(X,SinX,'g')legend('y=sin(x)')subplot(2,2,2); %绘制y=P5(x)曲线plot(X,P5X,'r')legend('y=P5(x)')subplot(2,2,3); %绘制y=P7(x)曲线plot(X,P7X,'y')legend('y=P7(x)')subplot(2,2,4); %绘制y=P9(x)曲线plot(X,P9X,'b')legend('y=P9(x)')2. P171.2:#include<iostream>#include<stdlib.h>#include<cmath>#include<iomanip>#include<fstream>using namespace std;struct Coordinate /构建一个用于保存坐标点的数据类型double x;double y;int main()ofstream outfile;int N,i;double adv;double *lagran_build(Coordinate*,int);void lagran_output(double*,Coordinate*,int);double lagran(double*,Coordinate*,int,double);cout<<"请输入坐标点的个数：" cin>>N;Coordinate *Data; /构建变量Data，用于保存坐标点double *lagran_coefficient; /构建变量lagran_coefficient，用于保存拉格朗日插值多项式的系数Data=new Coordinate N-1;cout<<"请依次输入各点坐标："<<endl; for(i=0;i<N;i+) /输入已知坐标点坐标cin>>Datai.x>>Datai.y;lagran_coefficient=lagran_build(Data,N); /构建拉格朗日插值多项式lagran_output(lagran_coefficient,Data,N); /输出N阶的拉格朗日插值多项式函数adv=lagran(lagran_coefficient,Data,N,3.5); /计算5小时的平均温度cout<<"这五小时的平均温度估计为:"<<adv<<"华氏度"<<endl;outfile.open("171_2.txt");for(i=0;i<N;i+) /将坐标点及拉格朗日插值多项式系数输出到文件"171_2.txt"outfile<<std:left<<setw(10)<<setiosflags(ios:fixed)<<setprecision(5)<<Datai.x;outfile<<std:left<<setw(10)<<setiosflags(ios:fixed)<<setprecision(5)<<Datai.y;outfile<<std:left<<setw(10)<<setiosflags(ios:fixed)<<setprecision(5)<<lagran_coefficienti<<endl;outfile.close();system("pause");return 0;/此函数为用于计算N次拉格朗日插值多项是系数的函数/输入Data为已知坐标点的指针/输入N为要求的拉格朗日插值多项式的次数/输出拉lagran为N次拉格朗日插值多项式的系数矩阵的地址double *lagran_build(Coordinate *Data,int N)int i,j;double *lagran;lagran=new double N-1; /生成一个一维矩阵，用于保存N次拉格朗日插值多项式的系数for(i=0;i<N;i+) /计算N次拉格朗日插值多项式的系数lagrani=1;for(j=0;j<i;j+)lagrani=lagrani*(Datai.x-Dataj.x);for(j=i+1;j<N;j+)lagrani=lagrani*(Datai.x-Dataj.x);lagrani=1/lagrani;return lagran;/此函数用于输出一个N次拉格朗日插值多项式/输入lagran_coefficient为保存N次拉格朗日插值多项式的一维矩阵的地址/输入Data为保存已知坐标点的地址/输入N为N次拉格朗日插值多项式的次数void lagran_output(double *lagran_coefficient,Coordinate *Data,int N)int i,j;cout<<"拟合后的"<<N<<"次拉格朗日插值函数为:"<<endl;cout<<"P"<<N<<"(x)="for(i=0;i<N;i+)cout<<Datai.y*lagran_coefficienti;for(j=0;j<i;j+)cout<<"(x-"<<Dataj.x<<")"for(j=i+1;j<N;j+)cout<<"(x-"<<Dataj.x<<")"if(i<N+1)cout<<"+"cout<<endl;return;/此函数为计算N次拉格朗日插值多项式在某一点的值的函数/输入lagran_coefficient为保存N次拉格朗日插值多项式系数的一维矩阵的地址/输入Data为用于保存已知坐标点的地址/输入N为N此拉格朗日插值多项式的次数/输入x为要求N次拉格朗日插值的拟合点的x坐标/输出y为要求N次拉格朗日插值的拟合点的y坐标double lagran(double *lagran_coefficient,Coordinate *Data,int N,double x)double y=0,yk;int i,j;for(i=0;i<N;i+)yk=Datai.y*lagran_coefficienti;for(j=0;j<i;j+)yk=yk*(x-Dataj.x);for(j=i+1;j<N;j+)yk=yk*(x-Dataj.x);y=y+yk;return y;Matlab绘图函数主函数Val=load('171_2.txt'); %从文件172_2.txt中取出已知点的坐标及拉格朗日插值多项式系数X=(Val(:,1)' %将已知坐标点x值保存到矩阵X中Y=(Val(:,2)' %将拉格朗日插值多项式系数保存到矩阵Lagran中Lagran=(Val(:,3)' %将拉格朗日插值多项式系数保存到矩阵Lagran中corf=Y.*Lagran; x=linspace(1,6);y=lagran(X,corf,x); %计算绘图点坐标plot(X,Y,'o') %绘制已知点及拉格朗日插值多项式函数hold all;plot(x,y);legend('采样点','温度拟合曲线')被调用的函数lagran%此函数用于计算拉格朗日插值多项式的拟合值%X为已知点的x坐标%coef为拉格朗日插值多项式系数%x为被求点x坐标%y为被求点y坐标function y = lagran( X,coef,x )n=length(X);y=0;for i=1:n yk=coef(i); for j=1:i-1 yk=yk.*(x-X(j); end for j=i+1:n yk=yk.*(x-X(j); end y=y+yk;endend3.P178.1:#include<iostream>#include<stdlib.h>#include<cmath>#include<iomanip>#include<fstream>using namespace std;struct Coordinate /构建一个用于保存坐标点的数据类型double x;double y;int main()ofstream outfile;int N,i;double adv;double *newton_build(Coordinate*,int);void newton_output(double*,Coordinate*,int);double newton(double*,Coordinate*,int,double);cout<<"请输入坐标点的个数："cin>>N;Coordinate *Data; /构建数组Data，用于保存坐标点double *newton_coefficient; /构建数组new_coefficient，用于保存牛顿插值多项式的系数Data=new Coordinate N-1;Data=new Coordinate N-1;cout<<"请依次输入各点坐标："<<endl;for(i=0;i<N;i+) /输入已知点坐标cin>>Datai.x>>Datai.y;newton_coefficient=newton_build(Data,N); /构建次牛顿插值多项式newton_output(newton_coefficient,Data,N); /输出N次牛顿插值多项式adv=newton(newton_coefficient,Data,N,3.5); /计算5小数平均温度cout<<"这五小时的平均温度估计为:"<<adv<"华氏度"<<endl;outfile.open("178_1.txt"); /将坐标点及牛顿插值多项式系数输出到文件"178_1.txt"for(i=0;i<N;i+)outfile<<std:left<<setw(10)<<setiosflags(ios:fixed)<<setprecision(5)<<Datai.x;outfile<<std:left<<setw(10)<<setiosflags(ios:fixed)<<setprecision(5)<<Datai.y;outfile<<std:left<<setw(10)<<setiosflags(ios:fixed)<<setprecision(5)<<newton_coefficienti<<endl;outfile.close();system("pause");return 0;/此函数为用于计算N次牛顿插值多项式系数的函数/输入Data为已知坐标点的指针/输入N为要求的牛顿插值多项式的次数/输出newton_coefficient为N次牛顿插值多项式的系数矩阵的地址double *newton_build(Coordinate *Data,int N)int i,j;double deltax;double *newton_coefficient;newton_coefficient=new double N-1; /生成一个一维数组，用于保存N次牛顿插值多项式的系数for(i=0;i<N;i+) /计算N次牛顿插值多项式的系数newton_coefficienti=Datai.y;for(i=0;i<N;i+)deltax=Datai+1.x-Data0.x;for(j=N-1;j>i;j-)newton_coefficientj=(newton_coefficientj-newton_coefficientj-1)/deltax;return newton_coefficient;/此函数用于输出一个N次牛顿插值多项式/输入newton_coefficient为保存N次牛顿插值多项式的一维矩阵的地址/输入Data为保存已知坐标点的数组的地址/输入N为N次牛顿插值多项式的次数void newton_output(double *newton_co