数值分析实验报告

1、.数值分析实验报告 目录第二章11.【第一题】1第三章32.【第二题】3第四章6第六章14第七章32第七章习题43第八章46第二章1. 【第一题】在区间-1,1上分别取n=10、20用两组等距节点对龙格函数做多项式插值及三次样条插值，对每个n值，分别画出插值函数及的图形。解：先编制Lagrange插值函数的m文件：function y1=lagrange(x0,y0,x1)n=length(x0);syms x;for k=1:n l(k)=x/x; for p=1:n if p=k l(k)=l(k)*(x-x0(p)/(x0(k)-x0(p); end endendz=0;for k=1:

2、n z=z+l(k)*y0(k);endy1=subs(z,x,x1);然后对n=10的情形做Lagrange插值并画图：x0=-1:.2:1;y0=1./(1+25*x0.2);x=-1:.001:1;y1=lagrange(x0,y0,x);y=1./(1+25*x.2);plot(x,y,x,y1,'-.',x0,y0,'p');legend('Runge Function','插值函数','插值节点');title('n=10时的Lagrange插值的龙格现象');xlabel('x

3、');ylabel('y');再对n=20的情形做Lagrange插值并画图：x0=-1:.1:1;y0=1./(1+25*x0.2);x=-1:.001:1;y1=lagrange(x0,y0,x);y=1./(1+25*x.2);plot(x,y,x,y1,'-.',x0,y0,'p');legend('Runge Function','插值函数','插值节点');title('n=20时的Lagrange插值的龙格现象');xlabel('x');yla

4、bel('y');第三章2. 【第二题】下列数据节点的插值：01491625364964012345678可得到平方根函数的近似。x=0 1 4 9 16 25 36 49 64;y=0 1 2 3 4 5 6 7 8;x0=0:.5:64;y0=lagrange(x,y,x0);y1=sqrt(x0);plot(x0,y1,x0,y0,'-.',x,y,'p')legend('开方函数','插值函数','插值节点'); x=0 1 4 9 16 25 36 49 64;y=0 1 2 3 4 5

5、6 7 8;x0=0:.5:64;y0=spline(x,y,x0);y1=sqrt(x0);plot(x0,y1,x0,y0,'-.',x,y,'p')legend('开方函数','三次条插值函数','插值节点');x=0.0 .1 .2 .3 .5 .8 1.0;y=1.0 .41 .50 .61 .91 2.02 2.46;>> p3=polyfit(x,y,3);>> p4=polyfit(x,y,4);>> x0=0:.001:1;y3=polyval(p3,x0);y

6、4=polyval(p4,x0);plot(x,y,'p',x0,y3,'-.',x0,y4,'-');legend('拟合节点','三次拟合','四次拟合');xlabel('x');ylabel('y');title('多项式拟合');grid on第四章1.【第一题】用不同的数值方法计算积分.在MATLAB中先求精确解：a=int('sqrt(x)*log(x)',0,1)a =-4/9(1)取不同的步长h，分别用复合梯形与复合辛

7、普森求积公式计算积分，给出误差中关于h的函数，并与精确积分值比较两个公式的精度，是否存在一个最小的h，使得精度不能再改善？采用复合梯形公式for k=50:50:8000x=linspace(0,1,k+1);x=x(2:k+1);y=x.(1/2).*log(x);b=trapz(x,y);a=a;k,b;endformat longa(:,2)+4/9ans = ；第 32 页 共 32 页略对于复合辛普生公式，先编制复合辛普生公式的M文件：function q=combinesimpson(F,a,b,n)%复合Simpson求积公式%F被积函数%a,b积分区间%n区间份数x=linsp

8、ace(a,b,n+1);q=0;for k=1:n q=q+F(x(k)+4*F(x(k)+x(k+1)/2)+F(x(k+1);endq=q*(b-a)/n/6;用辛普森公式有：F=inline('x.(1/2).*log(x)');a=;for k=100:100:5000b=combinesimpson(F,1/k,1,k);a=a;k,b;end>> a(:,2)+4/9ans = 略（2）用龙贝格求积计算完成问题（1）.function q,n,T=Romberg(f,x,a,b,tol)%Romberg求积公式%f被积函数%x积分变量%a,b积分区间%

9、tol积分精度%n迭代次数T=;k=0;h=b-a;T(0+1,0+1)=h/2*(subs(f,x,a)+subs(f,x,b);while k=0|abs(T(0+1,k+1)-T(0+1,k-1+1)>tol k=k+1;h=h/2;s=0; for jj=1:2:2k-1 s=s+subs(f,x,a+jj*h); end T(k+1,0+1)=1/2*T(k-1+1,0+1)+h*s; for jj=1:k T(k-jj+1,jj+1)=4jj/(4jj-1)*T(k-jj+1+1,jj-1+1)-1/(4jj-1)*T(k-jj+1,jj-1+1); endendn=k;q=

10、T(0+1,k+1);for jj=2:n+1 T(jj:n+1,jj)=T(1:n-jj+2,jj); T(1:jj-1,jj)=zeros(jj-1,1);End123（3）用自适应辛普生积分，使其精度达到.>> F=inline('sqrt(x).*log(x)');>> s=quad(F,0,1,1e-4)s = -0.44402.【第二题】计算二重积分.（1）若区域，试分别用复合辛普森公式及Gauss求积公式（取n=4）求积分。对于多重积分的复合辛普生公式，先编制多重积分复合辛普生公式的M文件：function q=combinesimpson

11、2(F,x0,a,b,n)%复合Simpson多元求积公式%F被积函数%x0被积函数自变量%a,b积分区间%n区间份数x=linspace(a,b,n+1);q=0;for k=1:n q=q+subs(F,x0,x(k)+4*subs(F,x0,(x(k)+x(k+1)/2)+subs(F,x0,x(k+1);endq=q*(b-a)/n/6;用辛普森公式计算多重积分有：>> syms x y;F=exp(-x.*y);>> s=combinesimpson2(combinesimpson2(F,'x',0,1,4),y,0,1,4)s = 0.796

12、6对于多重积分的Gauss公式，先编制Gauss公式的M文件，其中要用到函数GAUSSPOINTS，使用格式及作用如下：GAUSSPOINTS points for Gaussian quadrature nodes, weights = gausspoints(nquad) computes the nodes and weights for Gaussian quadrature with nquad (<=5) quadrature points on the standard interval -1,1编制的M文件为：function q=gaussquad(F,x0,a,b,n

13、)%Gauss求积公式%F被积函数%x0被积函数自变量%a,b积分区间%n节点个数syms t;F=subs(F,x0,(b-a)/2*t+(a+b)/2);x,A=gausspoints(n);q=(b-a)/2*sum(A.*subs(F,t,x);用Gauss求积公式计算多重积分有：>> clear>> syms x y;F=exp(-x.*y);>> s=gaussquad(gaussquad(F,x,0,1,4),y,0,1,4)s =0.7966（2）若区域，用复合辛普森公式（取n=4）求此积分.>> syms x y;F=exp(-

14、x.*y);>> s=combinesimpson2(combinesimpson2(F,y,0,sqrt(1-x2),4),x,0,1,4)s = 0.6701第五章1.function A,b,det=liezhuyuan(A,b)det=1;m=size(A);n=m(1);u=zeros(1,n);for k=1:n-1 ik=k; for jj=k+1:n if abs(A(jj,k)>abs(A(ik,k) ik=jj; end end if A(ik,k)=0 x=0;det=0;disp('无解');return; end if ik=k fo

15、r jj=k:n u(jj)=A(k,jj);A(k,jj)=A(ik,jj);A(ik,jj)=u(jj); end v=b(k);b(k)=b(ik);b(ik)=v;det=-det; end for ii=k+1:n A(ii,k)=A(ii,k)/A(k,k); for jj=k+1:n A(ii,jj)=A(ii,jj)-A(ii,k)*A(k,jj); end b(ii)=b(ii)-A(ii,k)*b(k); endendif A(n,n)=0 x=0;det=0;disp('无解');return;endb(n)=b(n)/A(n,n);for ii=n-1:

16、-1:1 s=0; for jj=ii+1:n s=s+A(ii,jj)*b(jj); end b(ii)=(b(ii)-s)/A(ii,ii);enddet=A(n,n)*det;function A,b,det=liezhuyuan(A,b)det=1;m=size(A);n=m(1);u=zeros(1,n);for k=1:n-1 ik=k; for jj=k+1:n if abs(A(jj,k)>abs(A(ik,k) ik=jj; end end if A(ik,k)=0 x=0;det=0;disp('无解');return; end if ik=k for

17、 jj=k:n u(jj)=A(k,jj);A(k,jj)=A(ik,jj);A(ik,jj)=u(jj); end v=b(k);b(k)=b(ik);b(ik)=v;det=-det; end for ii=k+1:n A(ii,k)=A(ii,k)/A(k,k); for jj=k+1:n A(ii,jj)=A(ii,jj)-A(ii,k)*A(k,jj); end b(ii)=b(ii)-A(ii,k)*b(k); endendif A(n,n)=0 x=0;det=0;disp('无解');return;endb(n)=b(n)/A(n,n);for ii=n-1:-

18、1:1 s=0; for jj=ii+1:n s=s+A(ii,jj)*b(jj); end b(ii)=(b(ii)-s)/A(ii,ii);enddet=A(n,n)*det;function A,b,Ip=liezhuyuanLU(A,b)m=size(A);n=m(1);s=zeros(1,n);tt=s;Ip=zeros(n,1);if det(A)=0 disp('无解');return;endfor r=1:n for ii=r:n A(ii,r)=A(ii,r)-A(ii,r:n)*A(r:n,r); end u=max(abs(A(ii,r) if ir=r

19、tt=A(r,:);A(r,:)=A(ir,:);A(ir,:)=tt; end A(r,r)=s(r); if r=n for ii=r+1:n A(ii,r)=A(ii,r)/A(r,r); u=0; for k=1:r-1 u=u+A(r,k)*A(k,ii); end A(r,ii)=A(r,ii)-u; end endendfor ii=1:n-1 u=Ip(ii); if ii=u v=b(u);b(u)=b(ii);b(ii)=v; endendfor ii=2:n u=0; for k=1:ii-1 u=u+A(ii,k)*b(k); end b(ii)=b(ii)-u;end

20、b(n)=b(n)/A(n,n);for ii=n-1:-1:1 u=0; for k=ii+1:n u=u+A(ii,k)*b(k); end b(ii)=(b(ii)-u)/A(ii,ii);end第六章1.【第一题】给出线性方程组，其中系数矩阵为Hilbert矩阵，假设，.若取n=6,8,10，分别用Jacobi迭代法及SOR迭代（=1,1.25,1.5）求解，比较计算结果。解：首先编制Jacobi迭代求线性方程解的m文件jacobiiteract.m：function x=jacobiiteract(A,b,x,n,flag)%Jacobi迭代求线性方程Ax=b的解%x-迭代初值%n-

21、迭代次数%flag-控制过程是否显示L=0 x' norm(A*x-b);h=size(A);p=h(2);for k=1:n for s=1:p w=0; for t=1:p if t=s w=w+A(s,t)*x(t); end end z(s)=(b(s)-w)/A(s,s); end x=z' L=L;k x' norm(A*x-b);endif flag Lend顺便编制Gauss Seidel迭代求线性方程解的m文件gauss_seidel_iteract.m：function x=gauss_seidel_iteract(A,b,x,n,flag)%Gau

22、ss Seidel迭代求线性方程Ax=b的解%x-迭代初值%n-迭代次数%flag-控制过程是否显示L=0 x' norm(A*x-b);h=size(A);p=h(2);for k=1:n for s=1:p w=0; for t=1:s-1 w=w+A(s,t)*x(t); end for t=s+1:p w=w+A(s,t)*x(t); end x(s)=(b(s)-w)/A(s,s); end L=L;k x' norm(A*x-b);endif flag Lend再编制超松弛迭代求线性方程解的m文件SOR.m：function x=SOR(A,b,x,omega,n,

23、flag)%超松弛迭代求线性方程Ax=b的解%x-迭代初值%omega-松弛因子%n-迭代次数%flag-控制过程是否显示L=0 x' norm(A*x-b);h=size(A);p=h(2);for k=1:n for s=1:p w=0; for t=1:p w=w+A(s,t)*x(t); end x(s)=x(s)+omega*(b(s)-w)/A(s,s); end L=L;k x' norm(A*x-b);endif flag Lend编制文件xiti6_1.m：for k=6:2:10xstar=ones(k,1);H=hilb(k);b=H*xstar;x=0.

24、8*ones(k,1);kdisp('Jacobi迭代结果');J=jacobiiteract(H,b,x,10,1)for v=1:.25:1.5disp('松弛因子'); vdisp('Gauss Seigel迭代结果');G=SOR(H,b,x,v,10,1)endend执行>> xiti6_1结果略第七章1.【第一题】求下列方程的实根：（1）；（2）.要求（1）设计一种不动点迭代法，要使迭代序列收敛，然后用Steffensen加速迭代，计算到为止。>> x=-10:.1:10;y1=x.2-3*x+2;y2=exp

25、(x);>> plot(x,y1,x,y2)>> x=-10:.1:6;y1=x.2-3*x+2;y2=exp(x);plot(x,y1,x,y2)>> x=-2:.05:2;y1=x.2-3*x+2;y2=exp(x);plot(x,y1,x,y2)>> x=0.255:.00005:0.26;y1=x.2-3*x+2;y2=exp(x);plot(x,y1,x,y2);grid on>> x0=0.257:0.00001:0.258;syms x;>> f1=sqrt(3*x-2+exp(x);>> f2=

26、(x2+2*x-exp(x)/3;>> f3=log(x2-3*x+2);>> f4=(exp(x)-2)/(x-3);>> f5=(exp(x)-2)/x+3;>> y1=subs(diff(f1),x,x0);>> y2=subs(diff(f2),x,x0);>> y3=subs(diff(f3),x,x0);>> y4=subs(diff(f4),x,x0);>> y5=subs(diff(f5),x,x0);>> subplot(3,1,1);plot(x0,y1)>&g

27、t; subplot(3,1,2);plot(x0,y2)>> subplot(3,1,3);plot(x0,y3)>> subplot(2,1,1);plot(x0,y4)>> subplot(2,1,2);plot(x0,y5)这表明f1 f3 f5迭代不收敛，而f4的收敛性要稍比f2好，所以采用f4迭代。根据作图，迭代初值设为0.2575。function xx,n,L=fixedpoint_iteract(f,x,x0,tol)%不动点迭代求解非线性方程%f-方程字符串%x-方程变量%x0-迭代初值%tol-迭代精度%xx-迭代结果%n-迭代次数%L

28、-迭代过程y=subs(f,x,x0);n=1;L=n,y,y-x0;while abs(y-x0)>tol x0=y; y=subs(f,x,x0);n=n+1; U=n,y,y-x0;L=L;U;endxx=y;>> xx n L=fixedpoint_iteract(f4,x,0.2575,1e-8);>> xxxx = 0.25753028364411>> nn = 10>> LL = 1.00000000000000 0.25754172799259 0.00004172799259 2.00000000000000 0.2575

29、2596199544 -0.00001576599715 3.00000000000000 0.25753191898097 0.00000595698553 4.00000000000000 0.25752966823006 -0.00000225075090 5.00000000000000 0.25753051864311 0.00000085041304 6.00000000000000 0.25753019732752 -0.00000032131559 7.00000000000000 0.25753031873177 0.00000012140425 8.000000000000

30、00 0.25753027286101 -0.00000004587076 9.00000000000000 0.25753029019259 0.00000001733157 10.00000000000000 0.25753028364411 -0.00000000654847采用Steffensen加速迭代：function xx,n,L=steffensen_iteract(f,x,x0,tol)%steffensen迭代求解非线性方程%f-方程字符串%x-方程变量%x0-迭代初值%tol-迭代精度%xx-迭代结果%n-迭代次数%L-迭代过程F=x-(f-x)2/(subs(f,x,f

31、)-2*f+x);y=subs(F,x,x0);n=1;L=n,y,y-x0;while abs(y-x0)>tol x0=y; y=subs(F,x,x0);n=n+1; U=n,y,y-x0;L=L;U;endxx=y;>> xx n L=steffensen_iteract(f4,x,0.2575,1e-8);>> xxxx = 0.25753028543986>> nn = 2>> LL = 1.00000000000000 0.25753028534589 0.00003028534589 2.00000000000000 0.2

32、5753028543986 0.00000000009397（2）用牛顿迭代，同样计算到。输出迭代初值及各迭代值和迭代次数k，比较方法的优劣。用Newton迭代法function xx,n,L=newton_iteract(f,x,x0,tol)%牛顿迭代求解非线性方程%f-方程字符串%x-方程变量%x0-迭代初值%tol-迭代精度%xx-迭代结果%n-迭代次数%L-迭代过程F=x-f/diff(f,x,1);y=subs(F,x,x0);n=1;L=n,y,y-x0;while abs(y-x0)>tol x0=y; y=subs(F,x,x0);n=n+1; U=n,y,y-x0;L

33、=L;U;endxx=y;>> f=x2-3*x+2-exp(x);>> xx n L=newton_iteract(f,x,0.2575,1e-8);>> xxxx = 0.25753028543986>> nn = 2>> LL = 1.00000000000000 0.25753028535414 0.00003028535414 2.00000000000000 0.25753028543986 0.00000000008572（1）>> x=-10:.1:10;y=x.3+2*x.2+10*x-20; plot(

34、x,y),grid on>> x=-4:.05:4;y=x.3+2*x.2+10*x-20; plot(x,y),grid on>> x=1:.005:2;y=x.3+2*x.2+10*x-20; plot(x,y),grid on>> x=1.3:.0001:1.4;y=x.3+2*x.2+10*x-20; plot(x,y),grid on>> x=1.36:.00005:1.38;y=x.3+2*x.2+10*x-20; plot(x,y),grid on>> x=1.365:.00001:1.370;y=x.3+2*x.2+1

35、0*x-20; plot(x,y),grid on>> x0=1.3685:0.000001:1.3690;syms x;>> f1=20/(x2+2*x+10);y1=subs(diff(f1),x,x0);>> plot(x0,y1)>> xx n L=fixedpoint_iteract(f1,x,1.3688,1e-8);>> xxxx = 1.36880810541671>> nn = 10>> LL =略>> f=x3+2*x2+10*x-20;>> xx n L=newto

36、n_iteract(f,x,1.3688,1e-8);>> xxxx = 1.36880810782137>> nn = 2>> LL = 1.00000000000000 1.36880810784040 0.00000810784040 2.00000000000000 1.36880810782137 -0.000000000019032.>> r=1 2 3 4 5 6 7 8 9 10;p=poly(r)p =略>> F=p;F(2)=F(2)+1e-6;>> x1=roots(F)x1 =略>> F

37、=p;F(2)=F(2)+1e-8;>> x2=roots(F)x2 =略>> F=p;F(2)=F(2)+1e-10;>> x3=roots(F)x3 =略第七章习题1. function xx,n,L=bicut(f,x,x1,x2,tol)%二分法迭代求解非线性方程%f-方程字符串%x-方程变量%x1,x2-迭代区间%tol-迭代精度%xx-迭代结果%n-迭代次数%L-迭代过程if subs(f,x,x1)*subs(f,x,x2)>=0 disp('迭代区间选取不合适，请重新选取');return;endn=0;L=;tol=m

38、in(abs(x2-x1)/2,tol);while abs(x2-x1)>tol xx=(x1+x2)/2;n=n+1; if subs(f,x,x1)*subs(f,x,xx)>0 x1=xx; elseif subs(f,x,xx)=0 return; else x2=xx; end U=n,xx,subs(f,x,xx);L=L;U;endsyms x;f=x2-x-1;xx n L=bicut(f,x,0,2,.05);xx = 1.59375000000000n = 6L = 1.00000000000000 1.00000000000000 -1.0000000000

39、0000 2.00000000000000 1.50000000000000 -0.25000000000000 3.00000000000000 1.75000000000000 0.31250000000000 4.00000000000000 1.62500000000000 0.01562500000000 5.00000000000000 1.56250000000000 -0.12109375000000 6.00000000000000 1.59375000000000 -0.053710937500002.（1）>> syms x;>> f=1+1/x2

40、;xx n L=fixedpoint_iteract(f,x,1.5,1e-5)略（2）>> f=(1+x2)(1/3);xx n L=fixedpoint_iteract(f,x,1.5,1e-5)略（3）function xx,n,L=fixedpoint_iteract(f,x,x0,tol)%不动点迭代求解非线性方程%f-方程字符串%x-方程变量%x0-迭代初值%tol-迭代精度%xx-迭代结果%n-迭代次数%L-迭代过程if abs(subs(diff(f,x,1),x,x0)>=1 disp('迭代发散'); xx=x0; n=0; L=; ret

41、urn;endy=subs(f,x,x0);n=1;L=n,y,y-x0;uu=x0;while abs(y-x0)>tol x0=y; y=subs(f,x,x0); if abs(y-uu)>5|n>1e4 disp('迭代发散');xx=y;return; end n=n+1; U=n,y,y-x0;L=L;U;endxx=y;>> f=1/sqrt(x-1);xx n L=fixedpoint_iteract(f,x,1.5,1e-5)迭代发散第八章2【第二题】function V,lambda,L=mifa(A,n)%幂法求矩阵A的主特征

42、值lambda与其对应的特征向量V%n迭代次数%L迭代过程k=size(A);v=ones(k(1),1);u=v;mu=max(v);L=0,u',mu;for k=1:n v=A*u;mu=max(v);u=v/mu;L=L;k,u',mu;endV=u;lambda=mu;>> A=1 1 .5;1 1 .25;.5 .25 2;>> V lambda L=mifa(A,20)略>> A=5 4 1 1;4 5 1 1;1 1 4 2;1 1 2 4;format long>> V lambda L=mifa(A,30)略f

43、unction V,lambda,L=mifa_rayleigh(A,n)%幂法求矩阵A的主特征值lambda与其对应的特征向量V,并用Raleigh商加速%n迭代次数%L迭代过程k=size(A);v=ones(k(1),1);u=v;mu=max(v);L=0,u',mu;for k=1:n v=A*u;mu=(u'*A*u)/(u'*u); u=v/mu; L=L;k,u',mu;endV=u;lambda=mu;>> V lambda L=mifa_rayleigh(A,30)略function V,lambda,L,n=mifa1(A,to

44、l)%幂法求矩阵A的主特征值lambda与其对应的特征向量V%n迭代次数%L迭代过程%tol特征值误差限k=size(A);v=ones(k(1),1);u=v;mu0=max(v);L=0,u',mu0;v=A*u;mu=max(v);u=v/mu;L=L;1,u',mu;n=1;while abs(mu-mu0)>tol n=n+1; mu0=mu; v=A*u;mu=max(v);u=v/mu;L=L;n,u',mu;endV=u;lambda=mu;>> V lambda L n=mifa1(A,1e-8);>> nn =28fun

45、ction V,lambda,L,n=mifa_rayleigh1(A,tol)%幂法求矩阵A的主特征值lambda与其对应的特征向量V,并用Raleigh商加速%n迭代次数%L迭代过程%tol特征值误差限k=size(A);v=ones(k(1),1);u=v;mu0=max(v);L=0,u',mu0;v=A*u;mu=(u'*A*u)/(u'*u);u=v/mu;L=L;1,u',mu;n=1;while abs(mu-mu0)>tol n=n+1; mu0=mu; v=A*u;mu=(u'*A*u)/(u'*u); u=v/mu; L=L;n,u',mu;endV=u;lambda=mu;>> V lambda L n=mifa_rayleigh1(A,1e-8);>> nn = 15第九章function T,X=ode_euler(f,t