数值分析实验报告二2汇总

1、数学与信息工程学院实验报告课程名称： 计算方法实验室：5206实验台号：24班 级：姓 名：实验日期： 2014 年 3 月26 日实验名称非线性方程求解实验目的和要求实验目的：1.计算机实现最小二乘法2. 最小二拟合与插值法比较3. 计算机实现数值积分实验要求：1.每种算法要求达到给定的精度，输岀数值积分结果及所需分半次数；2.利用逐次分半递推公式；实验内容和步骤：实验内容：1、分别用最小二乘法，拉格朗日插值，牛顿插值，求满足f(1)=-2,f(2)=3,f(3)=4,f(4)=-1,的三次拟合多项式，并在图中画出它们的图像。2、地球卫星轨道是一个椭圆，椭圆周长的计算方式是兀cs=4af屮-

2、(一)2sin2e d。，这里a是椭圆的半长轴，c是地球中心与轨道中va心(椭圆中心)的距离，记h为近地点距离，H为远地点距离，R=6371 (km为地球半径，则a=(2R + H +h)/2,c = (H -h)/2。我国第一颗人造地球卫星近地点距离h=439 ( km，远地点距离H=2384( k"，试求卫星轨道的周长。第一题:最小二乘法实验步骤：1) 实验编程 建立m文件：fun cti on c d l=leasts(x,y, n)% in put x y m维% output n次多项式m=le ngth(x)d=;A=zeros(m, n+1); b=zeros(m,1)

3、;for k=1:n+1A(:,k)=x.A(k-1);endb=A'*y'A=A'*Ac=Abfor i=1:md(i)=c(m+1-i);enddI=poly2sym(d)运行下列语句：x=1 2 3 4;y=-2 3 4 -1;n=3;leasts(x,y, n)2) 运行结果d =-0.3333 -0.0000 7.3333 -9.0000l =-3002399751579999/9007199254740992*xA3-311/1125899906842624*xA2+4128299658423301/562949953421312*x-25332747903

4、96013/281474976710656拉格朗日插值 实验步骤：1)实验编程建立M文件：fun cti onp,pO=Lagra nge1(x,y,x0)syms t;n=len gth(x);p=0;for i=1:nL=1;for j=1:i-1L=L*(t-x(j)/(x(i)-x(j);en d;for j=i+1:nL=L*(t-x(j)/(x(i)-x(j);en d;L=L*y(i);p=p+L;endp=simplify(p);p0=subs(p, 't',x0)p=vpa(p)运行下列语句：x=1 2 3 4;y=-2 3 4 -1;S=Lagra nge1

5、(x,y,1)运行结果：p0 =-2p =-.33333333333333333333333333333333*L3+7.3333333333333333333333333333333*t-9S =-.33333333333333333333333333333333*tA3+7.3333333333333333333333333333333*t-9牛顿插值法: 实验步骤：1)实验编程建立M文件：fun ctiony,R,A,C,L=newt on (X,Y,x,M)n=le ngth(X);m=le ngth(x);for t=1:mz=x(t);A=zeros (n,n);A(:,1)=Y&

6、#39;s=0.0;p=1.0;q1=1.0;c1=1.0;for j=2: nfor i=j:nA(i,j)=(A(i,j-1)- A(i-1,j-1)/(X(i)-X(i-j+1);endq1=abs(q1*(z-X(j-1);c1=c1*j;endC=A( n,n );q仁abs(q1*(z-X( n);for k=(n-1):-1:1C=co nv(C,poly(X(k);d=le ngth(C);C(d)=C(d)+A(k,k);end y(k)= polyval(C, z);endR=M*q1/c1;L(k,:)=poly2sym(C);L=vpa(L)运行下列语句：syms M,

7、X=1 2 3 4;Y =-2 3 4 -1; x=1;y,R,A,C,P=newt on (X,Y,x,M)2)运行结果：P=-2.6666666666666666666666666666667*xA3+14.*xA2-18.333333333333333333333333333333*x+5.画图：运行下列命令：>> x1=1 2 3 4;y仁-2 3 4 -1;x=0:0.001:3;t=0:0.001:3;y=-3002399751579999/9007199254740992*x.A3-311/1125899906842624*x42+412829965 8423301/

8、562949953421312*x-2533274790396013/281474976710656;Pl=-.33333333333333333333333333333333*t.A3+7.3333333333333333333333333333333*t -9.;P2=-.33333333333333333333333333333333*x.A3-.55511151231257827021181583404541e- 16*x.A2+7.3333333333333333333333333333333*x-9.;subplot(3,1,1)plot(x1,y1, 'r*',x

9、,y, 'b-')legend('数据点(xi,yi)','最小二乘法拟合曲线y=f(x)' );xlabel( 'x' );ylabel( 'y');title('数据点(xi,yi)和最小二乘法拟合曲线y=f(x)的图形')subplot(3,1,2)plot(x1,y1, 'r*',t,p1,'k-')legend('数据点(xi,yi)','拉格朗日插值曲线 y=f(x)' );xlabel( 'x' );yla

10、bel( 'y');title('数据点(xi,yi)和拉格朗日插值曲线y=f(x)的图形')subplot(3,1,3)plot(x1,y1, 'r*',x,P2, 'y-')legend('数据点(xi,yi)','牛顿插值曲线 y=f(x)' );xlabel( 'x' );ylabel( 'y'); title('数据点(xi,yi)和牛顿插值曲线y=f(x)的图形')运行结果：数据点脚.yi)和最小二乘法拟合曲线y=f(x)的團形数擄点(u再

11、和拉格朗日插倩曲线y詔凶的图形数据点浹观和牛顿插值曲线严地)的團形WlnrLi1111牛顿插值曲线y=f(x) >1J111 1 100.511 522.533.54实验结果分析：最小二乘法拟合的曲线误差最小。也可以得到三图合一的图像：在以上命令的基础上运行命令 plot(x1,y1,'r*',x,y,'b-',t,p1,'k-',x,P2,'y-')得到| IIIILLLL* 00.611.S22.533 5第二题：逐次分半梯形求积公式建立f.m的m文件fun cti on f=f(x)R=6371.0;H=2384.0;

12、h=439.0;a=(2*R+H+h)/2;c=(H-h)/2;f=4*a*sqrt(1-(c/a)A2*(si n(x)A2);再建立 compT.m 的 m 文件fun cti on n ,t2=compT(a,b,eps)%利用复化梯形递推公式求 fun在a,b的定积分; %输入积分下限a;输入积分上限b;精度eps； %输出分半次数n和积分值t2 ;具体问题1的值s;h=b-a;%节点步长t1=(f(a)+f(b)*h/2;%梯形值 T0 :h=h/2;t2=t1/2+h*f(a+h);%经1次分半复合梯形值err=t2-t1;淞差n=1;while abs(err)>epssu

13、m=0;h=h/2;t1=t2;n=n+1;k=1;while k<2Ansum=sum+f(a+k*h);k=k+2;endt2=t1/2+h*sumerr=t2-t1;end在程序窗口输入：n ,t2=compT(0,pi/2,0.00001)得到结果为t2 =4.8707e+0042t2 =4.8707e+004由此我们可以得到：S =4a o'2 .仁(c)2sinJdv - 4.8707 104复合辛普森(Simpson)数值积分建立f.m的m文件fun cti on f=f(x)R=6371.0;H=2384.0;h=439.0;a=(2*R+H+h)/2;c=(H-

14、h)/2;f=4*a*sqrt(1-(c/a)A2*(si n(x)A2);建立 simps on .m 的 m文件fun cti on simps on=simps on (x0,x n, eps)S1= (x n-x0)/6*(f(x0)+4*f(x0+x n)/2)+f(x n);S2=S1;d=0;n=1;m=1;while abs(S2-S1)>eps|d=0d=1;S仁 S2;for k=0:2*m-1S2=S2+2*(-x0+x n)/(2A n*3)*f(x0+(2*k+1)*(x n-x0)/2A( n+1);endfor j=0:m-1S2=S2-(-x0+x n)/

15、(2A n*3)*f(x0+(2*j+1)*(x n-x0)/2A( n);endS2=S2-1/2*S1;n=n+1;m=2*m;S2-S1;end cs=nS2在程序窗口输入:clearxO=O;xn=pi/2;eps=0.01;sin pus in( x0,x n, eps)运行程序得到:cs =S2 =4.8707e+004 while err>=eps由此我们可以得到:S =4a二/21 -(£)2 sin2 日d日4= 4.8707 10龙贝格积分求解 建立 Romberg.m 的 m 文件 fun cti on R,k,T=Romberg(a,b,eps)% f积

16、分函数% a/b :积分上下限% tol :积分误差% R : Romberg积分值% k :二分次数k=1;h=b-a;%第一步T(k,1)=h/2*(f(a)+f(b); err=1;T(k,k)= T(k,1);h=h/2;%第二步求梯形值T0temp=0;i=1;while i<2Ak temp=temp+f(a+i*h);i=i+2;endT(k+1,k)=T(k,k)/2+h*temp;%第三步计算加速值T(k+1,k+1)=(4Ak*T(k+1,k)-T(k,k)/(4Ak-1);%第四步精度控制err=abs(T(k+1,k+1)-T(k,k);k=k+1;endk=k-1;