南京邮电大学数值计算实践报告

1、数值计算实践I、方程求根一、实验目的熟悉和掌握Newton法,割线法，抛物线法的方法思路，并能够在 matlab上编 程实现二、问题描述(1) .给定一个三次方程，分别用Newton法,割线法，抛物线法求解.方程的构造方法：(a)根:方程的根为学号的后三位乘以倒数第二位加1再除以1000.假设你的学号为 B06060141，则根为141*(4+1)/1000=0.564(b)方程:以你的学号的后三位数分别作为方程的三次项，二次项，一次项的系数，根 据所给的根以及三个系数确定常数项.例如：你的学号是B06060141，则你的方程是x3+4x2+x+a0=0的形式.方程的根为0.564,因止匕有0

2、.5643+4*0.5642+0.564+a0=0，于是 a0=-2.015790144你的方程为 x3+4x2+x-2.015790144=0.(2)假设方程是sinx+4x2+x+a0=0的形式(三个系数分别是学号中的数字)，重新 解决类似的问题(3)构造一个五次方程完成上面的工作.四次方程的构造：将三次多项式再乘以(x-p*)2得到对应的五次多项式(p*为已经确 定的方程的根，显然，得到的五次方程有重根).(4)将(2)中的方程同样乘以(x-p*)得到一个新的方程来求解注:(1)Newton法取0.5为初值，割线法以0,1为初值，抛物线法以0,0.5,1为初值, (2)计算精度尽量地取高

3、.终止准则：根据| pn - pnFW来终止可供研究的问题:（一）8的取值不同对收敛速度有多大的影响（二）将注（1）中的初值该为其它的初值，对收敛性以及收敛速度有无影响（二）能否求出方程的所有一的根（4）实验报告的撰写实验报告包含的内容：（一）实验目的（二）问题描述（三）算法介绍（包 括基本原理）（四）程序（五）计算结果（六）结果分析（七）心得体会三、算法介绍在本问题中，我们用到了 newton法，割线法，抛物线法。1 .Newton法迭代格式为：f （xk）一“一国当初值与真解足够靠近,newton迭代法收敛，对于单根,newton收敛速度 很快，对于重根，收敛较慢。2 .割线法：为了回避导

4、数值的计算，使用上的差商代替，得到割线法迭代 公式：f （Xk）（Xk -Xk j）Xk 1 = xk -f（Xk） -f （XkJ割线法的收敛阶虽然低于newton法，但迭代以此只需计算一次函数值，不 需计算其导数，所以效率高，实际问题中经常应用。3 .抛物线法：可以通过三点做一条抛物线，产生迭代序列的方法称为抛物 线法。其迭代公式为：p（X）= fXk fXk,Xk_i（X-Xk）fXk,Xk_i,Xk/（X-Xk）（X-Xk_i）其中fXk,Xk, fXk, Xk山Xkj是一阶均差和二阶均差。收敛速度比割线法更接近于newton法。对于本问题的解决就以上述理论为依据。终止准则为：| Xn

5、 -Xn|< 8本题中所有精度取1e-8。四、程序计算结果问题一根据所给的要求，可知待求的方程为：X3 +2x2 +x-0.3356 =0牛顿法建立newton_1.m的源程序，源程序代码为：function y=newton_1（a,n,x0,nn,eps1）x（1）=x0;b=1;i=1;while(abs(b)>eps1*x(i)i=i+1;x(i)=x(i-1)-n_f(a,n,x(i-1)/n_df(a,n,x(i-1);b=x(i)-x(i-1);if(i>nn)errorreturn;endendy=x(i);建立n_f.m的源程序来求待求根的实数代数方程的函数

6、，源程序代码为: function y=n_f(a,n,x)% 待求根的实数代数方程的函数y=0.0;for i=1:(n+1)y=y+a(i)*xA(n+1-i);end建立n_df.m的源程序来方程的一阶导数的函数，源程序代码为：function y=n_df(a,n,x)% 方程的一阶导数的函数y=0.0;for i=1:ny=y+a(i)*(n+1-i)*xA(n-i);end在matlata件中执行下列语句并得到的最终结果截图：» a=l, 2, 1, -0.3366:» *3;» x0=0. 5;» nn=10C0：>> epsl

7、=le-8；>> y=nevt on_ 1 fa, nr xO, eps 1)0. 2240割线法建立gexian.m的源程序，源程序代码为function x=gexian(f,x0,x1,e)if nargin<4,e=1e-4;endy=x0;x=x1;i=0;while abs(x-y)>ei=i+1;z=x-(feval(f,x)*(x-y)/(feval(f,x)-feval(f,y);y=x;x=z;endi在matlata件中执行下列语句并得到的最终结果截图:>> fun= inline ( x 3+2+x 2+x-0. 3356：>&

8、gt; gexian (fun, 0, 1, le-S)ans =0. 2240抛物线建立paowuxian.m的源程序，源程序代码为：function x=paowuxian(f,x0,x1,x2,e)if nargin<4,e=1e-4;endx=x2;y=x1;z=x0;i=0;while abs(x-y)>ei=i+1;h1=y-z;h2=x-y;c1=(feval(f,y)-feval(f,z)/h1;c2=(feval(f,x)-feval(f,y)/h2;d=(c1-c2)/(h1+h2);w=c2+h2*d;xi=x-(2*feval(f,x)/(w+(w/abs(

9、w)*sqrt(wA2-4*feval(f,x)*d);z=y;y=x;x=xi;end在matlata件中执行下列语句并得到的最终结果截图:>> funFinlineCx'3+2*x 2+x-O,3356h);>> paowuxiandurij 必 0, 53 L le-8)CL 2240研究一：只改变初值由上述结果可知，方程的解在0.2附近，所以将牛顿法为0.2;割线法的初值设为0, 0.4;抛物线法的初值设为0, 0.2, 0.4;牛顿法根据问题1中牛顿法的程序，在matlab软件中执行下列语句并得到的最终结果截 图：» H 2, 1,-0,33

10、56;» n=3;» x0=0.2;» epsl=le-8;» nn=lDOO;» y=nevrton_1 区 延 nn, epsl)0. 2240割线法根据问题1中割线程序，在matlab软件中执行下列语句并得到的最终结果截图:>> finline ( x" 3+2*k * 2+x-0. 3356')>> gestianlfun, Cl, . % 1e-8)i -ans =0. 2240抛物线法根据问题1中抛物线法程序，在 matlab软件中执行下列语句并得到的最终结果 截图：» fun=i

11、nline( x 2+x-O. 3356)：>> paowuxian(fun, 0* 2t 0* 4, 1 s-8)ans 二Q. 2240研究二只改变精度将精度由1e-8改为1e-50和1e-100观察迭代次数有何变化牛顿法：根据问题1中的牛顿法的程序，在 matlab软件中执行下列语句并得到的最终结果截图：精度为1e-50时» 卒1,2. L-0, 3356;» n=3;» xO=O. E:» nn=1000;>> epsl=1e_50;» y=nevrt on_ 1 (n, x0, nrij ep s 1)I 22

12、40精度为1e-100时» epsl=le-lC0;» y=nevrton_l (a, n3 kOj nn, epsl)0. 2240割线法根据问题1中的割线法的程序，在 matlab软件中执行下列语句并得到的最终结果截图：精度为1e-50时>> fun= inline(工"3+2株"2+x-0. 3356 );>> gexian (fun, 0, 1, le-BO)g ans -精度为1e-100时>> fun= inline工-3+2"2+x-04 3356?);>> gexianffurij

13、 0? lf ls_100)ans =0.2240抛物线法根据问题1中的抛物线法的程序，在 matlab软件中执行下列语句并得到的最终 结果截图：精度为1e-50时>> fun=inline ( x* 3+2*x"2+x-0.3356');» paowuxian(fun, 0, 0» 5j L le-505i =10ars =0. 2240精度为1e-100时>> paowuxian(funj 0 0* 53 1 j le-100)10ans =0. 2240研究结论在只改变初值时，当初值定得越靠近初值，迭代次数就越少 在只改变精度

14、时，当精度越来越大时，迭代次数并几乎不变。综上所述，初值对迭代次数的影响比较大，精度对迭代次数影响不大问题二问题描述根据所给的要求，可知待求的方程为：sin x + 2x2 +x-0.3356 = 0问题分析仍然利用（1）中方法求解这一问题，并利用图解法找到初值，通过观察图像， 将newton法初值设为：0.1,割线法初值设为：0,0.2。抛物线法初值设为：0,0.1,02 图像见下图：43.532.521.510.50-0.500.10.20.30.40.50.60.70.80.91Newton 法问题一的牛顿法的求解只适用于线性方程，所以在问题二中用其他方法来求解方程。建立newton1.

15、m源程序，源程序代码为：function x=newton1(fn,dfn,x0,e) if nargin<4,e=1e-4;endx=x0;x0=x+2*e;i=1 ;while abs(x0-x)>ei=i+1 ；x0=x;x=x0-feval(fn,x0)/feval(dfn,x0);end在matlab软件中执行下列语句并得到最终结果截图» fun= inline (" sin (x) +2*x* 2-hc-O. 3356,)；» dfun=inline(J cos;» newt on! (fun, dfun (L I, le-8)i

16、 -5ans -0.1466割线法利用问题一中的割线法程序，在matlab软件中执行下列语句» gexiandurij Qj Q, 2,1 e-3)ans =0.1466抛物线法利用问题一中的抛物线法程序,在 matlab软件中执行下列语句>> paowuxian (fun, 0, 0. 1, 0. 2t 1 e-8)0.1466问题三问题描述按照题目要求对五次方程进行构造为：(x3 +2x2 +x-0.3356)(x-0.224)2 =0问题分析仍然利用一中方法求解这一问题，并利用图解法找到初值，通过观察图像， 将牛顿法的两组初值为0以及0.5,割线法初值设为：0,0.

17、5以及0, 0.2。抛物线 法初值设为：两组初值为 0,0.5,1以及0,0.25,05牛顿法利用问题二中的程序，在 matlab软件中执行下列语句：初值为0时> > fun=inline C (x 3+2*k 2+x-0. 3356)224 2');> > dfun= inlineC 2* (i-3+2*xR2-hc-0, 3356) * (x-Q. 224) + (3*x*2+4*x+l)* (x-0. 224) *Z );> > neirtonl (fuiij dfun, 0, le-8137ans =0. 2240初值为0.5时>>

18、; newt on 1 (fdfuiij 0.1 e-S)33ans =0. 2240割线法利用问题一中割线法的程序，在 matlab软件中执行下列语句: 初值为0, 0.5时» funFuilineC24*L-0, 3356) *(x-0, 224) ,2,» g exian (fun., Oj 0. 5j le-8)50ans =0. 2240初值为0, 0.2时>> fun=inline(, (x"3+2*x*2tK-0, 3356) * (xO. 224) *2");» gezianCfurij 0, 0, 2, le-8)

19、i -44ans -0. 2240抛物线法利用问题一中抛物线法的程序，在matlab软件中执行下列语句:初值为0, 0.5, 1时» fun=inlineC (x"3+2*x 2+x-0,3355)221);» paowuxian(funj 0, 0. 5；, lj le-B)i =63ans -0. 2240初值为0, 0.25, 0.5时» fun=inlineC (x'3+2*x 2-tx-O. 33S6) * (x-0. 224) 2J);» paowwcianCfurij 0j 0. 25, 0. 5,1 e-8)i -51a

20、naQ. 224Q问题四问题描述根据题目要求对方程进行构造为:(sin x 2x2 x-0.3356)( x-0.224) =0问题分析仍然利用问题一中方法求解这一问题，并利用图解法找到初值，通过观察图像， newton法初值选取了两组初值为0以及0.5,割线法初值设为：0, 0.5和0, 0.3 抛物线法初值设为：两组初值为0,0.2,0.4以及0, 0.5, 1。0.080.070.060.050.040.030.020.010-0.0100.050.10.150.20.250.3牛顿法利用问题二中的程序，在 matlab软件中执行下列语句:初值为0时» fun=inline C

21、 (sin(x)+2xx 2+-Hc-0. 3356) * (x-0. 224)' )：» dfun=inlinef (cos (k)+4*x+1)*(s-0. 224) + (sin(x)+2*s 2+so, 335E)');>> nevtonl (fun, dfun, 0： le-8)ans =0. 1466初值为0.5时» fun=inline(n (sin(x)+2*x"2i la 0. 3356) * (xO. 224)1);> > (ifun=inline C (cos (x)+44x+l) *(x-0, 224

22、) + (sin (x)+2*x''2+x-U4 3356)：> > nevtuni (fun, dfun, 0. 5, le-E)ans -0.2240割线法利用问题一中割线法的程序，在matlab软件中执行下列语句:初值为0, 0.5时» funpijtiline C (sin(x)+2*x*2-Hi-O4 3356) * tf-O* 224)1);> > gesianCfun, 0, 0, 5j le-8'113ans =初值为0, 0.3时» fun=inline C (sin(x)+2*x*2+x-0. 3356)*

23、(x-0.224)r)：>> gexian(fun, . 0. 3, le-S)i =10独3 =0. 2240抛物线法利用问题一中抛物线法的程序，在matlab软件中执行下列语句:初值为0, 0.1, 0.2时>> fun=inline (? (sin(x)+2*x 2+-hc-04 3356) * tx-0. 224), )：>> paowuxianffun, 口，0. 1, 0. 2, 1 e-8)i =14ans -0,1465 + 0.00001初值为0, 0.5, 1时>> funFinline( (sin+2*x 2+x-0

24、7; 3355)*(x-04 224)：>> paowuxian(fun, 0, 0. 5, 1,1 e-8)14ans -五、计算结果及分析Newton 法割线法抛物线法问题一0.2240.2240.224问题二0.14660.14660.1466问题三0.2240.2240.224问题四0.224 和0.224 和0.224 和0.1465772585571010.14660.1466问题一迭代步数问题二迭代步数问题三迭代步数问题四迭代步数Newton 法653789割线法855013 10抛物线法956314 13结果分析将Newton法，割线法，抛物线法进行比较可以看到在本

25、文题中，三种方法计算得到的最终结果基本相同，但是迭代步数有较大差别，综合看来Newton法迭代步数最少，割线法次之，抛物线法最次。在各个问题的研究中，我通常都会 采用不同的初值，发现不同初值会对应不同的迭代次数， 另外针对问题一，我选 用了不同的精度，发现迭代次数并没有很大的变化，因而一个初值的选定可以对 迭代次数产生很大的影响，而精度的改变对迭代次数的影响很小。对每个算法单独来看，显然选择初值不同对于迭代步数影响较大，对于找到 根也会有影响。因此应该先通过画图确定根的大致位置，给出在其附近的初值。六、心得体会在实现这三个算法的过程中，本身编程较易实现，最重要的是对算法本身的 理解，只有真正理

26、解算法的含义才能更快更好的实现程序。II、离散型最小二乘和连续型最小二乘问题一、实验目的掌握并能够利用离散型最小二乘和连续性最小二乘求解问题、问题描述1:以函数错误！未找到引用源。生成的6个数据点(0.25,23.1) ,(1.0,1.68), (1.5,1.0), (2.0,0.84) , (2.4,0.826) , (5.0,1.2576)为例，运行程序得到不同阶对应的曲线拟合产生的多项式函数。2.例题：计算f(x)=exp(x)在卜1, 1上的二、三次最佳平方逼近多项式。并画图进行比较。三、问题分析问题1是离散最小二乘问题最小离散最小二乘就是根据一批有误的数据(错误！未找到引用源。i=1

27、,2,，n确定参数，并通过均方误差来比较曲线拟合的优劣，在本题中通过 画图来比较不同阶方程拟合效果的优劣。选择两种方法实现离散最小二乘。方法一，建立normal equation(法方程组),求解k次多项式系数。法方程组构 造方法：m0a。' xi =0mnan"xi=0m_ 、 0- yi xii =0mm-1-2a0xi .a xi -i =0i =0m_ n，1anxii旬m1二 xxi =0ma°v xini =0mm- a/' xin 1 ' . aj xi2ni =0i =0m=£i =0nYiX方法二：由于在matlab中存在

28、ployfit函数，可以即为方便的用k次多项式拟问题2是一个连续型最小二乘法的应用实例。对于给定的函数f(x) w Ca, b如果存在S (x) Span( 0(x), 1(x), IH, n(x)使得2b2(P(x) f (x) -S (x) I dx = min f P(x) If (x) -s(x) Tdxa<<b一.一一4八一，、则称$(刈是£(刈在集合Span(50(x), %(x), | * (x)中的最佳平万逼近函数。n显然，求最佳平万逼近函数S (x) =£ aj ,p(x)的问题可归结为求它的系数j 0a；,a；,，a；,使多元函数bnI (a

29、。，ai, " , a；)=dx(P(x) |f (x) -Z aj啊(x)a-1/一取得极小值，也即点(a0, a1*,，a：)是I (a。，an)的极点。由于I (a。，ai,an)是关于a。, ai,，an的二次函数，利用多元函数取得极值的必要条件,=0 (k = 0, 1,2, , n) a即Hbn =2 ：(x) f(x) -y aj j(x) L ：k(x)dx=0-aka |L «得方程组n b'、aj ?(x) (x) j(x)dx aj=0b=7(x)f (x) k(x)dx, (k=0,1,2, ,n) a如采用函数内积记号b促k,%)= P(x

30、)%(x)Q(x)dx,q(f, k) =:?(x)f(x) k(x)dx,a那么，方程组可以简写为(Dn、 (),。问=(f, 1) (k 0,1,2,|l|,n)j f这是一个包含n + 1个未知元a0, ai,，an的n + 1阶线性代数方程组，写成矩阵形式为1'(%wjWH)iHi,%)WiWi)llllllllllll ©,%)(中 n,*)IIIIIIIll促呼n)、 促 1,Q)(中 nWn)Ja0ai、(f,九),(f,Q)(2)此方程组叫做求aj (j = 0, 1,2,，n)的法方程组。显然，其系数行列式就是克莱姆行列式 Gn = Gn (<p0,物

31、，<Pn)o由于邛0,明，旺线性无关，故Gn0 0 ,于是上述方程组存在唯一解ak =ak(k =0,i，n)。从而月止了函数 f (x)在 SpaK%(x),%(x),| Wn(x)中如果存在最佳平方逼近函数，则必是.(3)n_ *一* 2S(x)c aj j(x)j =0四、实验程序问题一法一：离散最小二乘法function f=normalequation(n)syms tr=0;xx=0.25,i.0,i.5,2.0,2.4,5.0;yy=23.i,i.68,i.0,0.84,0.826,i.2576;x=zeros(n,n);y=zeros(n,i); for i=i:nfor

32、 k=i:6x(i,i)=xx(k(2*i-2)+x(i,i);y(i,i)=yy(k).*xx(k(i-i)+y(i,i);endendfor i=2:nfor j=i:i-ifor k=i:6x(i,j)=xx(k).A(i+j-2)+x(i,j);x(j,i)=x(i,j);endendendz=xy;for i=i:nend vpa(r,5)法二:用matlab中的ployfit函数对k次多项式进行拟合。建立number2.m文件，带入如下：x=0.25,1.0,1.5,2.0,2.4,5.0;y=23.1,1.68,1.0,0.84,0.826,1.2576; plot(x,y,&#

33、39;o') p1=polyfit(x,y,1) %利用线性拟合 xi=-0.25:0.01:6.0;y1=polyval(p1,xi);plot(x,y,'o',xi,y1,'r:');hold on;p2=polyfit(x,y,2) % 二次拟合y2=polyval(p2,xi);plot(x,y,'o',xi,y2,'m');hold on;p3=polyfit(x,y,3) % 三次拟合y3=polyval(p3,xi);plot(x,y,'o',xi,y3,'b:');hold

34、on;p4=polyfit(x,y,4) % 四次拟合y4=polyval(p4,xi);plot(x,y,'o',xi,y4,'c');hold on;p5=polyfit(x,y,5) % 五次拟合y5=polyval(p5,xi);plot(x,y,'o',xi,y5,'g'); hold off;legend('初始点,'y=-2.3840*x+9.8499','初始点,'y=2.1793*xA2-15.8845*x+22.3999','初始点','y

35、=-2.0143*xA3+18.3499*xA2-45.2167*x+32.7386','初始点','y=1.4828*xA4-16.0435*xA3+56.3154*xA2-79.4994*x+39.6688',' 初始点 ','y=-0.9961*xA5+12.1743*xA4-55.2367*xA3+116.6262*xA2-116.7266*x+45.8090');问题二：最佳平方逼近算法(1)输入被逼近函数f(x)和对应的逼近区间a,b并选择逼近函数系q(x)和权函数；(2)解方程组(1)或(2),其中方程组的系

36、数矩阵和右端的项由式(3)得至(3)由式(3)得到函数的最佳平方逼近。将上述算法编写成MATLAB程序共需三个程序:第一个程序(函数名 Sapproach.m)计算最佳逼近函数的系数：function S=Sapproach(a,b,n)% 定义逼近函数global i;global j;if nargin<3 n=1;end% 判断X=zeros(n+1,n+1);for i=0:nfor j=0:n;X(i+1,j+1)=quad(fan,a,b); % 求 fan 积分 end endY=zeros(n+1,1);for i=0:nY(i+1)=quad(yb,a,b);%求 yb

37、 积分end s=XY第二个程序(函数名：fan.m)function y=fan(x) global i;global j;y=(poly(x,i).*poly(x,j);；第三个程序(函数名：yb.m)function y=yb(x) global i;y=(poly(x,i).*exp(x);编写多项式函数function y=poly(x,k)% 多项式函数if k=0y=ones(size(x);elsey=x.Ak;end五、实验结果及图像问题一：(1) 最小二乘法 normalequation.m运行结果为: 线性拟合：>> norjiial equat i on(

38、2)ans 二10. 68 - 2. SllS*t二次拟合：» normalequat ion(3) ans 二2. 5532M、2 - I6.962*t + 22.931三次拟合：» normalequationans =-2.2955*t13 + 19. 506*t-2 - 46. 507*+ + 33. 037四次拟合：» normalequal ion(b)ans =E6303*t'4 - 17. 152*t 3 + 5S. 393*t2 - 80, 932rt + 39. 917五次拟合：» no rm ale quation(6) a

39、ns =-1. 0853M"5 + 13. 027#t4 - 5；. 505+t"3 + 119.36粒 *2 - 113. 14*t + 46. 024用方法二得到的结果如下：number2.m运行结果为:pl =-2.9L1910. 6805p2 =2.5532-16.962022. 930919.5064-46. 507433, 03721. 6803 -17. 152253.3927 -80. 932439. 9168p5 =-1. 085313.0258 -57.5051 119. 3595 -118. 139545.0236C因此得到线性拟合多项式为：y - -

40、2.3840X 9.8499二次拟合多项式为：y =2.1793x2 -15.8845x 22.3999三次拟合多项式为：y = -2.0143x3 18.3499x2 -45.2167x 32.7386四次拟合多项式为：y =1.4828x4 -16.0435x3 +56.3154x2 -79.4994x+39.6688五 次 拟 合 多 项 式 为y = -0.9661x5 12.1743x4 -55.2367x3 116.6262x2 -116.7266x 45.8090绘制的各多项式拟合图像为：50403020100-10-2000,511.522.533.544.55线性拟合 二次拟

41、合 三次拟合 四次拟合 五次拟合 原数据点问题2:当求二次逼近时，有以下结果:>> a approach (-1, lj 2)0. 9963I, 10360. 6367绘制两者的图形：> > fun='exp(x)'> > fplot(fun,-1,1)> > hold on> > xi=-1:0.1:1;> > yi=polyval(S,xi);>> plot(xi,yi,'r:')得到以下结果：当三次逼近时，有以下结果» sapproach(1, lj 3) s =

42、0.99630,99800,53670.1761绘制图形如下：> > fun='exp(x)'> > fplot(fun,-1,1)> > hold on> > xi=-1:0.1:1;> > yi=polyval(S,xi);>> plot(xi,yi,'r:')得到如下所得图形：0 S -C 6-0.40 2 Q Q 20.40.60 8六、结果分析显然离散最小二乘中次数较高拟合的效果较好，但由于本问题中初始点较 少，并不能完全反应函数本身的特点，同时在本问题中，选择了两种方式，结果 接

43、近，也可以互相验证正确性。在连续最小二乘法的实验中较顺利的达到了预期的结果。从试验结果看出三次逼近没有二次逼近效果理想，验证了最佳平方逼近理论。七、心得体会在离散最小二乘与连续最小二乘程序设计的过程中，要么通过均差来表示函 数拟合的优劣，要么通过图像来评价，本问题中选择了图像，更加清晰直观。III、数值积分一、实验目的熟悉并掌握数值积分的方法，重要训练复化梯形公式，复化 simpson公式以及 romberg 积分。二、问题描述22问题三数值积分椭圆周长的计算。考虑椭圆与+%=1，为计算其周长，只要计算a b其第一象限的长度即可.一 .、一一. x = a cost用参数方程可以表示为(0 &

44、lt; t M二/ 2),y = bsint计算公式为.a2 sin21 b2 cos2 tdt0为计算方便，我们可以令a = 1，即计算下面的积分V2 :222sin t b cos tdt=f0 /l +(b2 -1)cos2 tdt/ /2 -212-22 ,(I. ;a sin t b cos tdt0=a(sin21 +(b)2cos tdt可以归结为上面的形式 采用复化梯形公式，复化 Simpson公式以及Romberg积分的方法计算积分72：-I (b) = p , 1 (b - 1)cos tdt给出通用程序，该通用程序可以计算任何一个函数在任意一个区间在给定的精度下的数值积分

45、。程序输出为计算出的数值积分值以及计算函数值的次数。利用你的程序计算5个椭圆的周长。三、算法介绍首先利用给出的各迭代公式，设计程序。在 matlab对话框中输入要计算的函数， 给出区间和精度。bhn 1旻化梯形的迭代公式为：a f(X)dX=21f(a)2 f (Xj) . f(b)2j 1bhr ,、f)-1 ,、,、,、复化 simpson迭代公式为：f f (x)dx =- f (a)+2£ f(x2j)+4£ f (x2j)+ f (b)a3j 1j 1Romberg迭代公式为： RJ=R+ Rk,jRk"四、算法程序复化梯形公式和复化simpso心式我们

46、放在jifen.m中。(标记后的程序可用来把b看为变量时的算法实现)%复化梯形公式function y=jifen(f,n,a,b)(说明：f表示任一函数，n精度，a, b为区间)fi=f(a)+f(b);h=(b-a)/n;d=1;%function f=jifen(n,a,b,c)%syms t%y=sqrt(1+(cA2-1)*cos(t)A2);%ya=subs(y,t,a);%yb=subs(y,t,b);%fi=ya+yb;for i=1:n-1x=a+i*h;fi=fi+2*f(x);d=d+1;%yx=subs(y,t,x);%fi=fi+2*yx;endf4=h/2*fi,d

47、%复化simposo心式f1=0;f2=0;dd=1;for i=1:n-1dd=dd+1;if rem(i,2)=0; x1=a+i*h; f1=f1+f(x1);else rem(i,2)=0; x2=a+i*h; f2=f2+f(x2); endendf3=(h/3)*(f(a)+4*f1+2*f2+f(b),ddromberg 积分建立romberg.m文件function y=romberg(f,n,a,b)(说明：f表示任一函数，n精度，a, b为区问)z=zeros(n,n);h=b-a;z(1,1)=(h/2)*(f(a)+f(b);f1=0;for i=2:nfor k=1:

48、2A(i-2)f1=f1+f(a+(k-0.5)*h);endz(i,1)=0.5*z(i-1,1)+0.5*h*f1;h=h/2;f1=0;for j=2:iz(i,j)=z(i,j-1)+(z(i,j-1)-z(i-1,j-1)/(4A(j-1)-1);endendz,n五、实验结果(1)对于复化梯形公式和复化simpson式，我们运行下列语句并得到结果:>> fun=inline(? sqrt 11+' C. 5 2-1 j. *coe (t). 2 1 );>> j if ent fun, 8j Ojpi/2)f4 =1.2111d -f3 =L2U1d

49、d =8题中将椭圆中的未知量a取为1, b取为0.5。f4为复化梯形公式得到的1/4个椭 圆周长，f3为复化simpson公式得到的1/4个椭圆周长。从而得到椭圆的周长为 4*1.2111=4.8444。(2)对于romberg,运行下列语句并最终得到结果为：» fun=inline(' sqrt(1+10. S"21). *cos(t), -2)J):» romberg fun, 31 0t pi/0. 5)2 -3. 141600000003. 14163.1416000004. 71245.23606.3756000004.83984. 88234.

50、85874. 8505000084J24,84574.84324. 84304,84290004.84424, 84424.S4414. 84414.84424.6442004.84424. 84424. 94424. 84424. 84424.E4424. 844204.84424, 34424.84424. 84424,84424.64424.84424,8442n =S题中将椭圆中的未知量a取为1, b取为0.5。用romberg算法对椭圆的周长进行 计算的到椭圆的周长为4.8442。六、结果分析我们计算了当椭圆长轴为1,短轴为0.5时的周长。通过上述三种方 法的计算可以看到，结果相差不

51、大。根据椭圆周长的一个计算公式我 们可以得到L=4.283。因此三种方法都较好的接近真值。七、心得体会我们应该熟练掌握这三种方法，才能在编程时正确快速的写出迭代公 式。同时在一种思想的前提下可以寻找多种方法实现算法，互相验证,从而得到更准确的结果。IV、newton差值和hermite差值方法一、实验目的掌握并能够利用newton差值和herm让e差值方法解决问题。二、问题描述二/222I(b) = f0 也 +(b2 -1)cos2 tdt22 bcos t上述函数的导数为 I'(b)： dt0221 (b - 1)cos t采用三种方法中最好的方法计算这一积分(1)利用数值积分白方

52、法给出I(b)在b = 0,0.1,0.2,|, 1 (可以直接计算精确值的，用精确值)，用Newton插值方法得到5个椭圆的周长(2)利用数值积分的方法给出I(b),I '(b)在b = 0,0.1, 0.2川|,1(可以直接计算精确值的，用精确值)，用Hermle插值方法得到5个椭圆的周长(3)选做题：利用I(b)以及导数更多的值来进行插值，插值误差会有什么变化？(4)选做题：采用其它的插值方法改进插值的效果。三、算法介绍血定，对于给定的b值都对应着一个椭圆，在本问题中用newton®值法和hermae 得到的多项式代替椭圆周长公式中的，1+(b21)cos2t进行积分，

53、首先画出图像，选择初始点。图彳O勺实现代码见picture1.m。newtonj雨值法迭代公式：nPn (X)= f Xo 二 f Xo, XiXk ( X - Xo )( X - Xi )(X - Xk)k =1Hermite法迭代公式：2n 1H 2n 1 (X) = f Zo% f Zi Zk( X - Zo) (X - Zk)k 1作图时建立picture.m文件画出和其导数图像。(注：此图像为b=0.5时)0.80.60.40.20-0.2-0.40.20.40.60.81.21.41.61.8sqrt(1+(0.5 2-1).*cos(x). 2)0.750./(1-0.75.*c

54、os(x).2).(1/2).*cos(x).*sin(x)根据图像，我们选取我们选取0,0.3,0.6,0.9,1.2,1.5为初始点四、算法程序(1)建立 newtondedail.m文件。function z=newtondedai1(f,n)syms xia=zeros(n,n);x=0 0.3 0.6 0.9 1.2 1.5;y=feval(f,x);a(:,1)=y;for i=2:nfor j=2:ia(i,j)=(a(i,j-1)-a(i-1,j-1)/(x(1,i)-x(1,i-j+1);endendt=xi-x(1,1);p=a(1,1);for i=2:np=p+a(i,

55、i)*t;t=t*(xi-x(1,i);endp=collect(vpa(p)(2)建立 hermite3.m文件。function z=hermite3(fn,dfn,n)syms xia=zeros(2*n,2*n);x=0.3 0.6 0.9;xx=0.3 0.3 0.6 0.6 0.9 0.9;y=feval(fn,x);yy=feval(dfn,x);for i=1:nz(1,2*i)=x(1,i);z(1,2*i-1)=x(1,i);a(2*i,1)=y(1,i);a(2*i-1,1)=y(1,i);a(2*i,2)=yy(1,i);endfor i=1:nif 2*i<6a(2*i+1,2)=(a(2*i+1,1)-a(2*i,1)/(z(1,2*i+1)-z(1,2*i);endendfor i=3:(2*n)for j=3:ia(iJ)=(a(iJ-1)-a(i-1J-1)/(z(1,i)-z(1,i-j+1);endendp=a(1.1);t=xi-xx(1,1);for i=2:2*np=a(i,i)*t+p;t=t*(xi-xx(1,i);endp=vpa(collect(p) %p=vpa