11级数值方法计算实习题.doc

信计、数应专业11级数值方法计算实习题要求：1、用Matlab语言或你熟悉的其他语言编写程序，使之尽可能具有通用性；2、根据上机计算实践，对所使用的数值方法的特点、性质、有效性、误差和收敛性等方面进行必要的讨论和分析；3、完成计算后写出实验报告，内容包括：课题名称、解决的问题、采用的数值方法、算法程序、数值结果、对实验结果的讨论和分析等；4、特别说明：严禁抄袭，否则一经发现，所有雷同实验报告最多评为及格。一、在区间-1,1上分别取用两组等距节点对龙格函数作多项式插值及三次样条插值，对每个值，分别画出插值函数即的图形。解：作m文件，文件名：lagrange.mfunction lagrange(n)x0=linspace(-1,1,n+1);y0=1./(1+25*x0.2);x=linspace(-1,1,100);y=1./(1+25*x.2);y3=interp1(x0,y0,x,spline);for k=1:length(x) im=ones(1,n+1); for i=1:n+1 for j=1:n+1 if x0(i)=x0(j) im(i)=im(i)*(x(k)-x0(j)/(x0(i)-x0(j); end end end yp(k)=im*y0;endplot(x0,y0,*,x,y,b,x,y3,-,x,yp,-);legend(原插值点,原函数,三次样条插值,多项式插值)输入lagrange(10)调用函数，得到n=10时的插值函数与f（x）的函数图形输入lagrange(20)调用函数，得到n=20时的插值函数与f（x）的函数图形n=10时的插值函数与f（x）的函数图形n=20时的插值函数与f（x）的函数图形二、已知Wilson矩阵，且向量，则方程组有准确解。用Matlab 内部函数求，的所有特征值和；解： M=det(A) D=eig(A) N=cond(A)M = D = N = 1 0.0102 2.9841e+003 0.8431 3.8581 30.2887令，解方程组，并求出向量和，从理论结果和实际计算结果两方面分析方程组解的相对误差与的相对误差的关系；解：为了方便MATLAB的编辑，令B=A+A;M=A;X=x+x;t=x; B=10 7 8.1 7.2;7.08 5.04 6 5;8 5.98 9.89 9;6.99 4.99 9 9.98 ; b=32;23;33;31; X=Bb norm(t) X = ans = -81.0000 163.9695 % x的结果 137.0000 norm(t)/norm(x) -34.0000 ans = 22.0000 81.9848 %的结果 (3.7347e+003)/(-11.7943) m=B-A;norm(m)ans = ans = -316.6530 0.2308 %|A |的结果 X=Bb norm(m)/norm(A)X = ans = -81.0000 0.0076 %的结果 137.0000 d=inv(A);norm(d)*norm(m)*norm(x) -34.0000 ans = 22.0000 45.4829 % x=1;1;1;1;t=X-x; d=inv(A); t 1-norm(d*(m)t = ans = -82.0000 -11.7943 % 136.0000 (45.4829)/(-11.7943) -35.0000 ans = 21.0000 -3.8563 x=1;1;1;1;t=X-x norm(t)t = ans -82.0000 0.0074 136.0000 所以xAAx1AA -35.0000 21.0000 X=Bb X = -81.0000 137.0000 -34.0000 22.0000 x=1;1;1;1;t=X-x; tt = -82.0000 136.0000 -35.0000 21.0000 x=1;1;1;1;t=X-xt = -82.0000 136.0000 -35.0000 21.0000再改变扰动矩阵（其元素的绝对值不超过0.005），重复第2问。改变 ，取a2=0 0.002 0.001 0.003;0.001 0.004 0 0.001;0.003 -0.004 -0.001 0;-0.001 -0.002 0 -0.003；B2=a2+A; C2=B b; b=32;23;33;31;rank(B2)ans=4rank(C2)ans=4所以：rank(X2)=rank(C2)所以扰动矩阵有唯一解x2=B2bx2=1.06490.8893norm(x2)/norm(x)ans=6.3276norm(a2)ans=0.0071norm(a2)/norm(A)ans=1.3336e-004所以 =1.3336e-004norm(x2)ans=12.6552（norm(d)*norm(a2)*norm(x2)）/(1-norm(d*(a2)ans=13.0887因为12.655213.0887，所以即：x format rat a=0,-1,-1,-1,-1; b=2,2,2,2,2; c=-1,-1,-1,-1; f=1,0,0,0,0; n=length(a); m(1)=c(1)/b(1); y(1)=f(1)/b(1); for i=2:n-1 m(i)=c(i)/(b(i)-a(i)*m(i-1); y(i)=(f(i)-a(i)*y(i-1)/(b(i)-a(i)*m(i-1); end y(n)=(f(n)-a(n)*y(n-1)/(b(n)-a(n)*m(n-1); x(n)=y(n); for i=n-1:-1:1 x(i)=y(i)-m(i)*x(i+1); end i=1:n;fprintf(三角方程式的解为n);x(i)三角方程式的解为ans = 5/6 2/3 1/2 1/3 1/6 format short结果与和符合四、公元1225年，比萨的数学家Leonardo研究了方程，得到一个根，没有人知道他用什么方法得到这个值。对于这个方程，分别用下列方法：迭代法；迭代法；对的Steffensen加速方法；对的Steffensen加速方法；Newton法。求方程的根（可取），计算到Leonardo所得到的准确度。解：（1）迭代法编写m文件，文件名：diedai（）function x0,k,er,x=diedai(g,x0,wucha,max)%g是给定的迭代函数%x0是给定的初始值%wucha是规定的误差范围%max是迭代次数加1%x是不动点近似值%x（x1,x2.,xn）X(1)=x0;for k=2:max X(k)=feval(g,X(k-1); k；er=abs(X(k)-X(k-1); x=X(k) if(er diedai(g,1,10(-9),15)超出迭代次数k = 15x = 1.36880377314363ans = 1取解为：1.36880377314363（2）迭代法只需添加g1.m问价即可function y=g1(x)y=(20-2*x2-x3)/10;在窗口输入 diedai(g1,1,10(-9),15)超出迭代次数k = 15x = 1.36880377314363ans = 1取解为：1.36880377314363（3）对的Steffensen加速方法编写steffensen.m文件i=2;x0=1;f=inline(20/(x2+2*x+10);y=f(x0);z=f(y);x1=x0-(y-x0).2/(z-2*y+x0);s.result=x0;x1;while abs(x1-x0)=1e-9 x0=x1; y=f(x0); z=f(y); x1=x0-(y-x0).2/(z-2*y+x0); i=i+1; s.result(i)=x1;ends.step=(0:i-1);fprintf(迭代步数为：t%dn,i-1);for j=1:i fprintf(%10d,s.step(j);fprintf(t);end在命令窗口中输入： steffensen迭代步数为：4 01.0000000 11.3708139 21.3688082 31.3688081 41.36880815 分析结果知，steffensen迭代加速步数减少了，到第五步已经达到了精确精度要求（4）对的Steffensen加速方法把（3）中的迭代格式改为（20-2*x2-x*3）/10,保存，在命令窗口输入 steffensen迭代步数为：0 1.0000001 1.33349212 1.36841543 1.36880814 1.36880815 1.3688081分析结果知，steffensen迭代加速步数减少了，到第五步已经达到了精确精度要求（5）Newton法编写m文件，文件名：ndf.mfunction p1,er,k,y=ndf(f,df,p0,tol,max)p0;feval(ff,p0);for k=1:max p1=p0-feval(ff,p0)/feval(df,p0); er=abs(p1-p0); p0=p1; p1;er;k;y=feval(ff,p1); if(er ndf(ff,df,1,10(-9),10)k = 5ans = 1.36880810782137取解为：1.36880810782137五、用不同的数值方法计算积分的近似值，其中 取不同的步长，分别用复合梯形公式和复合辛普森公式计算积分，比较两个公式的计算效果，是否存在一个最小的，使得精度不能再被改善？ 用龙贝格求积公式，取，并打印出T-表。解：（1）用复合梯形公式，编写fhtx.m文件 用复化辛普森公式，编写fhxps.m文件function s=fhtx(f,a,b,n) function s=fhxps(f,a,b,n)%f是被积函数 %f是被积函数%a b分别为积分的上下限 %a b分别为积分的上下限%n是子空间的个数 %n是子空间的个数%s是梯形总面积 %s是梯形总面积h=(b-a)/n; h=(b-a)/(2*n);s=0; s1=0;s2=0;for k=1:n-1 for k=1:n x=a+k*h; x=a+(2*k-1)*h; s=s+feval(f,x); s1=s1+feval(f,x);end end s=h*(feval(f,a)+feval(f,b)/2+h*s; for k=1:n-1 x=a+2*k*h; s2=s2+feval(f,x); end编写被积函数文件 f.m s=h*(feval(f,a)+feval(f,b)+4*s1+2*s2)/3;function y=f(x);y=(10/x)2*sin(10/x);在命令窗口输入%n值取小的数如：10，15，20时 fhtx(f,1,3,10) fhxps(f,1,3,10)ans = ans = -4.7789 -1.3680 fhtx(f,1,3,15) fhxps(f,1,3,15)ans = ans = -2.8604 -1.4136 fhtx(f,1,3,20) fhxps(f,1,3,20)ans = ans = -2.2208 -1.4220 %n取较大的数如1000，2000时 format long fhtx(f,1,3,1000) fhxps(f,1,3,1000)ans = ans = -1.42633620719670 -1.42602475569236 fhtx(f,1,3,2000) fhxps(f,1,3,2000)ans = ans = -1.42610261856845 -1.42602475630535 format short分析结果知道两种方法有微小误差，当步长越大，精确度就越大，h取2000时两种方法误差精确到小数点第三位。（2）编写lbg.m文件function r,quad,tol,h=lbg(f,a,b,n,wucha)%f是被积函数%a b分别为积分的上下限%n+1是t数表的列数%wucha是容许的误差%quad是所求积分值m=1;h=b-a;tol=1;j=0;r=zeros(4,4);r(1,1)=h*(feval(f,a)+feval(f,b)/2;while(tolwucha)&(jn)|(j lbg(f,1,3,6,0.5*(10(-3)quad = -1.4260ans