数值分析实验2(wangwei)

实验2.1（多项式插值的振荡现象）问题提出：考虑在一个固定的区间上用插值逼近一个函数。

显然拉格朗日插值中使用的节点越多，插值多项式的次数就越高。

我们自然关心插值多项式的次数增加时，()n L x 是否也更加靠近被逼近的函数，龙格给出的一个例子是极著名并富有启发性的，设区间[-1，1]上函数 21()125f x x =+实验内容:考虑区间[-1，1]的一个等距划分，分点为 21,0,1,2,...,i ix i n n=-+= 则拉格朗日插值多项式为 201()()125nn i i iL x l x x==+∑其中的(),0,1,2,...,i l x i n =是n 次拉格朗日插值基函数。

实验要求：（1）选择不断增大的分点数目2,3...,n =画出原函数()f x 及插值多项式函数()n L x 在 [-1，1]上的图像，比较并分析实验结果。

（2）选择其他的函数，例如定义在区间[-5，5]上的函数4(),()arctan 1xh x g x x x==+ 重复上述的实验看其结果如何。

（3）区间[a,b]上切比雪夫点的定义为(21)cos ,1,2,...,1222(1)k a b b ak x k n n π⎛⎫+--=+=+ ⎪+⎝⎭以121,,...,n x x x +为插值节点构造上述各函数的拉格朗日插值多项式，比较其结果。

程序清单：1． 被逼近函数的函数文件func1.mfunction y=func1(x,c)%直接用被逼近函数计算函数值，c 用来选择函数 if c==1;y=1./(1+25*x.^2); end; if c==2;y=x./(1+x.^4); end; if c==3;y=atan(x);end;2．拉格朗日插值函数文件lagr.mfunction yy=lagr(x,y,xx)%用拉格朗日插值多项式进行插值,x为插值点的自变量列阵，y为对应的函数值列阵，xx为%待插值的自变量列阵，yy为求得的对应于xx的函数值列阵yy=xx-xx;%初始化为0向量n=length(x)-1;for i=0:n;z=yy;k=0;for j=0:n;if j~=ik=k+1;if k==1;z=(xx-x(j+1))/(x(i+1)-x(j+1));elsez=z.*(xx-x(j+1))/(x(i+1)-x(j+1));end;end;end;yy=yy+z*y(i+1);end;2．主程序main1.mm=input('请输入节点数目:');chof=input('请选择函数:(1:1/(1+25*x^2)) 2:x/(1+x^4) 3:arctanx)');div=input('请选择节点方式:(1.均布节点2.切比雪夫点)');%生成用于画图的自变量序列if chof==1;t=-1:0.01:1;else;t=-5:0.01:5;end;y1=func1(t,chof);%直接求t对应的函数值%按要求构造节点的自变量序列n=m-1;for i=0:n;if div==1;if chof==1;nod(i+1)=-1+2*i/n;else;nod(i+1)=-5+10*i/n;end;end;if div==2;if chof==1;nod(i+1)=cos((2*i+1)*pi/(2*(n+1))); elsenod(i+1)=5*cos((2*i+1)*pi/(2*(n+1))); end; end; end;nodv=func1(nod,chof);%求节点处函数值 y=lagr(nod,nodv,t); plot(t,y1,t,y); t=0; y=0; y1=0; nod=0; nodv=0;实验结果及其分析：（1） 分点数分别为2，3，5，8，11，15时的函数21()125f x x =+和插值多项式函数()n L x 的图像如下：n=2: n=3n=5: n=8:n=11: n=15:分析：从图上可以看出随着分点数目的增加区间端部的插值函数的值与原函数的误差越来达，而区间中部的误差越来越小，这就是所谓的龙格现象，这是等距节点的高次插值多项式的典型病态现象。

（2） 对于函数4()1xh x x =+情况如下： n=4: n=8:n=11: n=15:对于函数()arctan g x x =情况如下：n=3: n=7:n=11: n=15:分析：从图上可以看出对于上面这两种函数也有龙格现象出现。

（3） 若采用切比雪夫点则结果如下：21()125f x x =+：n=11: n=15:4()1xh x x =+: n=11： n=20:()arctan g x x =:n=11:分析：采用切比雪夫点有效的抑制了龙格现象，由于切比雪夫点在区间的端部较为密集而在区间中部较为稀疏，因此它有利于减小端部的误差，防止发生龙格现象。

实验2.3 (曲线逼近方法的比较)问题提出：曲线的拟合和插值，是逼近函数的基本方法，每种方法具有各自的特点和特定的适用范围，实际工作中合理选择方法是重要的。

实验内容：任然考虑实验2.1中的著名问题，用Matlab程序给出了该函数的二次和三次拟合多项式。

实验要求：（1）将拟合的结果与拉格朗日插值及样条插值的结果比较。

（2）归纳总结数值实验的结果，试定性地说明函数逼近各种方法的适用范围，及实际应用中选择方法应注意的问题。

程序清单：x=-1:0.2:1;y=1./(1+25*x.*x);xx=-1:0.02:1;%多项式拟合p2=polyfit(x,y,3);yy=polyval(p2,xx);plot(x,y,'o',xx,yy);xlabel('x');ylabel('y');hold on;%拉格朗日插值yy=lagr(x,y,xx);plot(xx,yy,'b');%样条插值yy=spline(x,y,xx);plot(xx,yy,'r');hold off；当需要分别画各种方法的图时，可对上面的程序进行相应的修改。

实验结果及其分析：对于实验2.1中的著名问题用各种逼近方法得到的结果分别如下：二次多项式拟合：六次多项式拟合：拉格朗日插值：三次样条插值：分析：从结果来看三次多项式拟合结果不理想，六次多项式拟合在中部稍好而在端部不理想， 拉格朗日插值在端部出现龙格现象，三次样条曲线拟合的结果较为理想。

（3） 曲线拟合适用于曲线不要求通过每一个数据点而只需反映数据点的分布趋势的情况，对于曲线拟合重要的是选好模型，例如前面的情况用多项式模型就不太合适，因此效果不理想。

而拉格朗日插值适用于严格要求曲线过每一个数据点，且在整个逼近区间上要求有一个统一表达式的情形，它的缺点是当通过的点数增多且点的分布不能满足一定要求时会出现龙格现象，此时在区间中部逼近较理想，端部的情况则往往很坏。

样条曲线插值适合于严格要求曲线过每个数据点，点数较多的情形，也就是拉格朗日插值不适用的情形。

它的适用范围较广，经常采用。

实验2.5 (高维积分数值计算的蒙特卡罗方法)问题提出：高维空间中的积分，如果维数不很高且积分区域是规则的或者能等价地写成多重积分的形式，可以用一元函数积分的数值方法来计算高维空间的积分。

蒙特卡罗方法对计算复杂区域甚至不连通的区域上的积分并没有特殊的困难。

实验内容：对于一般的区域Ω，计算其测度（只要理解为平面上的面积或空间中的体积）的一般方法是：先找一个规则的区域A 包含Ω，且A 的测度是已知的。

生成区域A 中m 个均匀分布的随机点,1,2,...,,i p i m =如果其中有n 个落在区域Ω中，则区域Ω的测度()m Ω为n/m 倍A 的测度。

函数()f x 在区域Ω上的积分可以近似为：区域Ω的测度与函数()f x 在Ω中n 个随机点上平均值的乘积。

1()()()k k p f x dx m f p nΩ∈Ω≈Ω⨯∑⎰实验要求：假设冰激凌的下部为一锥体而上面为一半球，考虑冰激凌体积问题：计算锥面222z x y =+上方和球面222(1)1x y z ++-=内部区域的体积。

如果使用球面坐标，该区域可以表示为如下的积分：2cos /422sin d d d ϕππρϕϕρθ⎰⎰⎰用蒙特卡罗方法可以计算该积分。

另一方面，显然这样的冰激凌可以装在如下立方体的盒子里 11,11,02x y z -≤≤-≤≤≤≤而该立方体的体积为8。

只要生成这个盒子里均匀分布的随机点，落入冰激凌锥点的个数与总点数之比再乘以8就是冰激凌锥的体积。

比较两种方法所得到的结果。

类似的方法可以计算复杂区域的测度（面积或体积）。

试求由下列关系所界定区域的测度：01,12,13(1)sin()0xx y z e y z y ≤≤≤≤-≤≤⎧⎪≤⎨⎪≥⎩33213,14(2)29201,01,01(3)sin 1x y x y x y y e x y z x y z x z e ≤≤-≤≤⎧⎪+≤⎨⎪≥-⎩≤≤≤≤≤≤⎧⎪+≤⎨⎪-+≤⎩程序清单：1． 计算冰激凌的体积n=input('请输入生成的随机点的数目:'); m=0; for i=1:n;p=rand(1,3); x=-1+2*p(1); y=-1+2*p(2); z=2*p(3);if x^2+y^2+(z-1)^2<1&x^2+y^2<z^2; m=m+1; end; end; 8*m/n2．(1)区域的测度n=input('请输入生成的随机点的数目:'); m=0; for i=1:n;p=rand(1,3); x=p(1); y=1+p(2); z=-1+4*p(3);if exp(x)<=y&sin(z)*y>=0; m=m+1;end;end;8*m/n3．(2)区域的测度n=input('请输入生成的随机点的数目:'); m=0;for i=1:n;p=rand(1,2);x=1+2*p(1);y=-1+5*p(2);if x^3+y^3<=29&exp(x)-2<=y;m=m+1;end;end;10*m/n4．（3）区域的测度n=input('请输入生成的随机点的数目:'); m=0;for i=1:n;p=rand(1,3);x=p(1);y=p(2);z=p(3);if x^2+sin(y)<=z&x-z+exp(y)<=1;m=m+1;end;end;m/n实验结果及其分析：（1）直接积分算得冰激凌的体积为。