数值分析实验报告.docx

1、数值分析实验报告数值分析实验报告（第二章）实验题目：分别用二分法、牛顿迭代法、割线法、史蒂芬森迭代法求方程的根，观察不同初始值下的收敛性，并给出结论。问题分析：题目有以下几点要求：1. 不同的迭代法计算根，并比较收敛性。2. 选定不同的初始值，比较收敛性。实验原理：各个迭代法简述二分法：取有根区间的重点，确定新的有根区间的区间长度仅为区间长度的一版。对压缩了的有根区间重复以上过程，又得到新的有根区间，其区间长度为的一半，如此反复，可得一系列有根区间，区间收敛到一个点即为根。牛顿迭代法：不动点迭代法的一种特例，具有局部二次收敛的特性。迭代格式为割线法：是牛顿法的改进，具有超线性收敛的特性，收敛阶

2、为1.618. 迭代格式为史蒂芬森迭代法：采用不动点迭代进行预估校正。至少是平方收敛的。迭代格式为这里可采用牛顿迭代法的迭代函数。实验内容：1. 写出该问题的函数代码如下：function py= f(x)syms k;y=(k2+1)*(k-1)5;yy=diff(y,k);py(1)=subs(y,k,x);py(2)=subs(yy,k,x);end2. 分别写出各个迭代法的迭代函数代码如下：二分法：function y=dichotomie(a,b,e)i=2;m(1)=a;while abs(a-b)e t=(a+b)/2; s1=f(a); s2=f(b); s3=f(t); if

3、 s1(1)*s3(1)=e s=f(x); t=x-s(1)/s(2); en=t-x; x=t; m(i)=t; i=i+1;endy=x,i+1,m;end牛顿割线法：function y=Secant(x1,x2,e)i=3;m(1)=x1,m(2)=x2;while abs(x2-x1)=e s1=f(x1); s2=f(x2); t=x2-(x2-x1)*s2(1)/(s2(1)-s1(1); x1=x2; x2=t; m(i)=t; i=i+1;endy=x2,i+1,m;end史蒂芬森迭代法：Function p=StephensonIterative(x,e)i=2;m(2)

4、=x;en=2*e;while abs(en)=e y=fai(x); z=fai(y); t=x-(y-x)2/(z-2*y+x); en=t-x; x=t; m(i)=t; i=i+1;endp=x,i+1,m;end3. 因为经常被使用，故可以写一个函数。代码如下：function y=fai(x) s=f(x); y=x-s(1)/s(2);end4. 可以绘制不同的图形来比较不同迭代法的收敛性和不同初值下的收敛性。代码如下：clear all;%相同初始值，不同迭代法下的收敛x1=dichotomie(0,3,1e-10);x2=NewtonIterative(0,1e-10);x3

5、=Secant(0,2, 1e-10);x4=StephensonIterative(0,1e-10);x1(2),x2(2),x3(2),x4(2)figure,subplot(2,2,1),plot(x1(3:x1(2),title(二分法);subplot(2,2,2),plot(x2(3:x2(2),title(牛顿迭代法);subplot(2,2,3),plot(x3(3:x3(2),title(牛顿割线法);subplot(2,2,4),plot(x4(3:x4(2),title(史蒂芬森迭代法);figure,subplot(2,2,1),plot(x1(4:x1(2)-1)-x

6、1(1)./(x1(3:x1(2)-2)-x1(1),title(二分法);subplot(2,2,2),plot(x2(4:x2(2)-1)-x2(1)./(x2(3:x2(2)-2)-x2(1),title(牛顿迭代法);subplot(2,2,3),plot(x3(4:x3(2)-1)-x3(1)./(x3(3:x3(2)-2)-x3(1),title(牛顿割线法);subplot(2,2,4),plot(x4(4:x4(2)-1)-x4(1)./(x4(3:x4(2)-2)-x4(1),title(史蒂芬森迭代法);%不同初始值，相同迭代法下的收敛性x5=dichotomie(-1,1

7、,1e-10);x6=dichotomie(-2,3,1e-10);x7=dichotomie(0,4,1e-10);x8=dichotomie(-4,4,1e-10);x9=NewtonIterative(-2,1e-10);x10=NewtonIterative(-4,1e-10);x11=NewtonIterative(4,1e-10);x12=NewtonIterative(6,1e-10);figure,subplot(1,2,1),plot(1:x1(2)-2,x1(3:x1(2),1:x5(2)-2,x5(3:x5(2),1:x6(2)-2,x6(3:x6(2),1:x7(2)-

8、2,x7(3:x7(2),1:x8(2)-2,x8(3:x8(2),title(二分法);subplot(1,2,2),plot(1:x2(2)-2,x2(3:x2(2),1:x9(2)-2,x9(3:x9(2),1:x10(2)-2,x10(3:x10(2),1:x11(2)-2,x11(3:x11(2),1:x12(2)-2,x12(3:x12(2),title(牛顿迭代法);x13=Secant(-1,1, 1e-10);x14=Secant(-4,5, 1e-10);x15=Secant(0,7, 1e-10);x16=Secant(-8,2, 1e-10);x17=Stephenso

9、nIterative(-1,1e-10);x18=StephensonIterative(-4,1e-10);x19=StephensonIterative(4,1e-10);x20=StephensonIterative(6,1e-10);figure,subplot(1,2,1),plot(1:x3(2)-2,x3(3:x3(2),1:x13(2)-2,x13(3:x13(2),1:x14(2)-2,x14(3:x14(2),1:x15(2)-2,x15(3:x15(2),1:x16(2)-2,x16(3:x16(2),title(牛顿割线法);subplot(1,2,2),plot(1:

10、x4(2)-2,x4(3:x4(2),1:x17(2)-2,x17(3:x17(2),1:x18(2)-2,x18(3:x18(2),1:x19(2)-2,x19(3:x19(2),1:x20(2)-2,x20(3:x20(2),title(史蒂芬森迭代法);实验结果：1. 各个迭代值分布图 1.1 不同迭代法下的得到的迭代值迭代值的情况如下：二分法牛顿迭代法牛顿割线法史蒂芬森迭代法00001.50000000000.20000000002.00000000001.35555555560.75000000000.37049180320.33333333330.98161652831.12500

11、000000.50764420760.38071968010.99994600030.93750000000.61461894470.49828334190.99999999951.03125000000.69738690980.57049963330.98437500000.76155380910.63938062441.00781250000.81154111860.69427858790.99609375000.85067638570.74116926531.00195312500.88144821230.78027159970.99902343750.90572974000.81329

12、27871当二分法的初始区间选为，误差限为，牛顿迭代法初值选为，误差限为，牛顿割线法初始点为，误差限为，史蒂芬森迭代法初始点选为，误差限为，迭代情况如图所示。迭代次数分别为38次，100次，140次，9次。故而，史蒂芬森迭代法速度最快，效果最好。2. 收敛情况图 1.2 不同迭代法下迭代值得收敛情况二分法收敛效果较差，牛顿迭代法和牛顿割线法相近，史蒂芬森迭代法收敛次数高于1，效果最好3. 不同初值的收敛情况图 1.3 二分法，牛顿迭代法下不同初值的收敛情况图 1.4 牛顿割线法，史蒂芬森迭代法下不同初值的收敛情况1. 二分法的五个初始区间分别为；2. 牛顿迭代法的五个初始值分别为；3. 牛顿割

13、线法的五个初始区间分别为；4. 史蒂芬森迭代法的五个初始值分别为；由图可知，它们最终均达到收敛。收敛性分析及结论：1. 二分法收敛较慢且不能求解崇根，但算法简单；此处牛顿法具有了平方收敛；从迭代次数上看，牛顿割线法较牛顿法的多，所以收敛性较差，是超线性收敛；史蒂芬森迭代法收敛效果最好。2. 因为牛顿迭代法是局部的二次收敛，所以要注重初值的选取，本次实验中选择的初值均得到了收敛，效果比较好。牛顿割线法也应注意初值的选取。（第三章）实验题目：1. 区间作等距划分：以为结点对函数进行插值逼近。（1） 分别取用牛顿插值对进行逼近，并在同一坐标系下做出函数的图形，进行比较。写出插值函数对的逼近程度与节点

14、个数的关系，并分析原因；（2） 试用三次样条插值对进行逼近，在同一坐标下画出图形，观察样条插值函数对的逼近程度与节点个数的关系；（3） 整体插值有何局限性？如何避免？2. 已知一组数据如下，求其拟合曲线.表 2.1 数据表01234567891023478101114161819106.42108.2109.5110109.93110.49110.59110.6110.76111111.2（1） 求以上数据形如的拟合曲线及其平方误差；（2） 求以上数据形如的拟合曲线及其平方误差；（3） 通过观察（1）（2）的结果，写出你对数据拟合的认识.问题分析：题目除上述要求之外还有以下几点：1. 明确整体

15、插值和分段插值的不同。牛顿插值多项式属于整体插值，三次样条插值属于低次分段插值。2. 将结果在同一坐标下绘制出。但是为了方便分析节点个数对于插值效果的影响，也可以单独绘制。3. 第二题中为了确定各个参数的大小，可以进行适当变换，转化为线性，运用最小二乘法，得到拟合。实验原理：牛顿插值多项式：对于给定的插值节点构造次数不超过的插值多项式使其满足插值条件这样得到的插值多项式称为插值多项式。系数为差商，可以通过构造差商表得到。三次样条插值：三次样条插值函数在每个小区间上为三次多项式；在全进上存在二阶连续导数；其次符合插值条件。中存在内置的三次样条插值函数，命令为 实验内容：第一题：1. 牛顿插值函数

16、的构造代码如下：function f=Newton(x0,y0) %牛顿多项式插值函数syms x;SZ=size(x0,2);a(1)=y0(1);y(:,1)=y0;for j=2:SZ nx1=1; for i=1:SZ-j+1; nx2=nx1+j-1; y(i,j)=(y(i,j-1)-y(i+1,j-1)/(x0(nx1)-x0(nx2); nx1=nx1+1; endendf=y(1,1);for j=2:SZ ff=y(1,j); for i=1:j-1 ff=ff*(x-x0(i); end f=f+ff;endend2. 牛顿和三次样条插值情况及比较：代码如下：clear

17、all;clc;syms x;fx=1/(5+x2);%牛顿多项式插值x0=-1:0.01:1;y0=subs(fx,x0);n1=1,5,10,20,25;n2=5,55,60,62,67;h1=2./n1;h2=2./n2;x1=-1:h1(1):1; y1=subs(fx,x1); f1=Newton(x1,y1); y01=subs(f1,x0); x2=-1:h1(2):1; y2=subs(fx,x2); f2=Newton(x2,y2); y02=subs(f2,x0); x3=-1:h1(3):1; y3=subs(fx,x3); f3=Newton(x3,y3); y03=s

18、ubs(f3,x0); x4=-1:h1(4):1; y4=subs(fx,x4); f4=Newton(x4,y4); y04=subs(f4,x0); x5=-1:h1(5):1; y5=subs(fx,x5); f5=Newton(x5,y5); y05=subs(f5,x0); figure,plot(x0,y0,x0,y01,x0,y02,x0,y03,x0,y04,x0,y05),title(所有结果);x6=-1:h2(1):1; y6=subs(fx,x6); f6=Newton(x6,y6); y06=subs(f6,x0); x7=-1:h2(2):1; y7=subs(f

19、x,x7); f7=Newton(x7,y7); y07=subs(f7,x0); x8=-1:h2(3):1; y8=subs(fx,x8); f8=Newton(x8,y8); y08=subs(f8,x0); x9=-1:h2(4):1; y9=subs(fx,x9); f9=Newton(x9,y9); y09=subs(f9,x0); x10=-1:h2(5):1; y10=subs(fx,x10); f10=Newton(x10,y10); y010=subs(f10,x0); figure,plot(x0,y0,x0,y06,x0,y07,x0,y08,x0,y09,x0,y01

20、0),title(龙格现象);%三次样条插值spline自带命令x0=-5:0.01:5;y0=subs(fx,x0);figure,subplot(2,3,1),plot(x0,y0),title(精确图);y11=spline(x1,y1,x0);subplot(2,3,2),plot(x0,y11),title(n=1);y12=spline(x2,y2,x0);subplot(2,3,3),plot(x0,y12),title(n=5);y13=spline(x3,y3,x0);subplot(2,3,4),plot(x0,y13),title(n=10);y14=spline(x4,

21、y4,x0);subplot(2,3,5),plot(x0,y14),title(n=20);y15=spline(x5,y5,x0);subplot(2,3,6),plot(x0,y15),title(n=25);figure,plot(x0,y0,x0,y11,x0,y12,x0,y13,x0,y14,x0,y15),title(所有结果);第二题：代码如下：clear all;clc;x=2 3 4 7 8 10 11 14 16 18 19;y=106.42 108.2 109.5 110 109.93 110.49 110.59 110.6 110.76 111 111.2;%Age

22、nda 1A=(x.*x) x ones(11,1);A1=A*A;y1=A*y;c1=A1y1x1=2:0.1:19;y1=polyval(c1,x1);plot(x,y,k*,x1,y1,r),gtext(曲线一),hold ony01=polyval(c1,x);delta1_2=sum(y-y01).2)%Agenda 2xx=-1./x;yy=log(y);B=ones(11,1) xx;B1=B*B;yy1=B*yy;c2=B1yy1;syms s;a=exp(c2(1);b=c2(2);a bf=a*exp(-b/s);x2=x1;y2=subs(f,x2);plot(x2,y2

23、,b),gtext(曲线二);y02=subs(f,x);delta2_2=sum(y-y02).2)实验结果：第一题：1. 牛顿插值结果图 2.1 不同节点数的牛顿插值多项式综合图图 2.2 龙格现象由图可知，2. 三次样条插值结果图 2.3 不同节点下的三次样条插值结果图 2.4 不同节点数的三次样条的综合图由图可知，牛顿插值多项式在较小的时候如差值效果良好，当变大时如就出现了龙格现象，三次样条在各个子区间内为三次多项式，拟合效果好.第二题1. 第一问的多项式拟合得到的拟合曲线为平方误差为第二问的拟合曲线为平方误差为2. 拟合曲线如图所示图 2.5 拟合情况从图中可以看出曲线二大体符合黑点

24、的分布情况，拟合效果较好。结论：整体插值随着节点个数的增多，多项式的次数也在升高。高次多项式的插值会出现龙格现象，震荡剧烈，逼近效果不理想。但是当节点很多时，低次插值的精度又不够，所以为了避免这一局限性采用分段低次插值。其中三次样条插值有良好的收敛性和光滑性，效果较好。数据拟合时，可以选择趋势相似的函数形式，在求解相关参量时，将函数进行适当变换，从而运用最小二乘法得到拟合结果。（第四章）实验题目：设计区间分半求积算法、龙贝格求积算法和自适应辛普森求积算法的程序，观察时，积分的结果，并给出相应的评价. 问题分析：1. 分半求积算法即为复合梯度求积。2. 因为越大，在内震荡地越严重，自适应辛普森求

25、积算法很难适用，计算复杂度会很高。所以选取辛普森求积算法代替。实验原理：分半求积算法：将区间等分为个小区间，分点为利用定积分性质，有每个小区间上均用梯形求积公式，有令即为复合梯形公式。龙贝格求积算法：将步长为的复合梯形公式进行査尔逊外推，就得到龙贝格求积公式。自适应辛普森求积算法：按照被奇函数在区间上的变化来安排求积结点的辛普森算法。实验内容：1. 写得各个求积算法的函数。代码如下：（1）分半求积算法：function f=CompositeTrapezoida(n)syms x;y=x*x*cos(n*x);m=1:100;a=-1;b=1;for i=1:100 h=(b-a)/m(i);

26、 k=0:m(i); xk=a+k*h; yk=subs(y,xk); f(i)=0; for j=1:m(i) f(i)=f(i)+(yk(j)+yk(j+1); end f(i)=h/2*f(i);endend（2）辛普森算法function f=CompositeSimpson(n)syms x;y=x*x*cos(n*x);m=1:100;a=-1;b=1;for i=1:100 h=(b-a)/(2*m(i); k=0:(2*m(i); xk=a+k*h; yk=subs(y,xk); f(i)=0; for j=1:m(i) f(i)=f(i)+(yk(2*j-1)+4*yk(2*

27、j)+yk(2*j+1); end f(i)=h/3*f(i);endm;f;end（4）龙贝格求积算法function f=Romberg(n)epsilon=0.0001;syms x;y=x*x*cos(n*x);a=-1;b=1;h=b-a;fa=subs(y,a);fb=subs(y,b);T(1)=h/2*(fa+fb);m=1;while 1 h=h/2; k=1:(2m-1); xk=a+k*h; yk=subs(y,xk); su=sum(yk); S(1)=h/2*(fa+fb)+h*su; for i=1:m S(i+1)=S(i)+(S(i)-T(i)/(4i-1);

28、end if(abs(S(m+1)-T(m)epsilon) break; end T=S; m=m+1;endf=S(m+1);end（5）自适应辛普森求积算法function f=AdaptiveSimpson(n)epsilon=0.001;syms x;y=x*x*cos(n*x);a=-1;b=1;h=b-a;p=a b;p0=p;ep=epsilon;m=0;q=0;f=0;while 1 n1=length(ep); n=length(p0); if(n=1) break; end h=p0(2)-p0(1); k=0:4; yk=subs(y,p0(1)+k*h/4); s0=h/6*(yk(1)+4*yk(3)+yk(5);s1=h/12*(yk(1)+4*yk(2)+yk(3);s2=h/12*(yk(3)+4*yk(4)+yk(5); if(abs(s0-s1-s2)15*ep(1)