《信计10级数值分析计算实习题》由会员分享，可在线阅读，更多相关《信计10级数值分析计算实习题（17页珍藏版）》请在金锄头文库上搜索。
1、信计 101 郭英贵 学号：201000901031信计、数应专业 10 级数值方法计算实习题要求：1、用 Matlab 语言或你熟悉的其他语言编写程序，使之尽可能具有通用性；2、根据上机计算实践，对所使用的数值方法的特点、性质、有效性、误差和收敛性等方面进行必要的讨论和分析；3、完成计算后写出实验报告，内容包括：课题名称、解决的问题、采用的数值方法、算法程序、数值结果、对实验结果的讨论和分析等；4、特别说明：严禁抄袭，否则一经发现，所有雷同实验报告最多评为及格。1、在区间-1,1上分别取 用两组等距节点对龙格函数 作多项式10,2n 21()5fx插值及三次样条插值，对每个 值，分别画出插值
2、函数即 的图形。解：（1）多项式插值先建立一个多项式插值的 M-file；输入如下的命令（如牛顿插值公式）：function C,D=newpoly(X,Y)n=length(X);D=zeros(n,n) D(:,1)=Y for j=2:nfor k=j:nD(k,j)=(D(k,j-1)- D(k-1,j-1)/(X(k)-X(k-j+1);endendC=D(n,n);for k=(n-1):-1:1C=conv(C,poly(X(k) m=length(C);C(m)= C(m)+D(k,k);end当 n=10 时，我们在命令窗口中输入以下的命令：clear,clf,hold on
3、;X=-1:0.2:1;Y=1./(1+25*X.2); C,D=newpoly(X,Y);x=-1:0.01:1;y=polyval(C,x);plot(x,y,X,Y,.);grid on;xp=-1:0.2:1;z=1./(1+25*xp.2);plot(xp,z,r)得到插值函数和 f（x）图形：当 n=20 时，我们在命令窗口中输入以下的命令：clear,clf,hold on;X=-1:0.1:1;Y=1./(1+25*X.2); C,D=newpoly(X,Y);x=-1:0.01:1;y=polyval(C,x);plot(x,y,X,Y,.); grid on;xp=-1:0
4、.1:1;z=1./(1+25*xp.2);plot(xp,z,r)得到插值函数和 f（x）图形：（2）三次样条插值先建立一个多项式插值的 M-file；输入如下的命令:function S=csfit(X,Y,dx0,dxn) N=length(X)-1;H=diff(X); D=diff(Y)./H;A=H(2:N-1);B=2*(H(1:N-1)+H(2:N); C=H(2:N); U=6*diff(D);B(1)=B(1)-H(1)/2;U(1)=U(1)-3*(D(1);B(N-1)=B(N-1)-H(N)/2;U(N-1)=U(N-1)-3*(-D(N); for k=2:N-1t
5、emp=A(k-1)/B(k-1);B(k)=B(k)-temp*C(k-1); U(k)=U(k)-temp*U(k-1); end M(N)=U(N-1)/B(N-1);for k=N-2:-1:1M(k+1)=(U(k)-C(k)*M(k+2)/B(k);endM(1)=3*(D(1)-dx0)/H(1)-M(2)/2;M(N+1)=3*(dxn-D(N)/H(N)-M(N)/2;for k=0:N-1S(k+1,1)=(M(k+2)-M(k+1)/(6*H(k+1);S(k+1,2)=M(k+1)/2;S(k+1,3)=D(k+1)-H(k+1)*(2*M(k+1)+M(k+2)/6;
6、S(k+1,4)=Y(k+1);end当 n=10 时，我们在命令窗口中输入以下的命令：clear,clcX=-1:0.2:1;Y=1./(25*X.2+1); dx0= 0.0739644970414201;dxn= -0.0739644970414201;S=csfit(X,Y,dx0,dxn) x1=-1:0.01:-0.5;y1=polyval(S(1,:),x1-X(1);x2=-0.5:0.01:0;y2=polyval(S(2,:),x2-X(2); x3=0:0.01:0.5; y3=polyval(S(3,:),x3-X(3);x4=0.5:0.01:1;y4=polyval
7、(S(4,:),x4-X(4); plot(x1,y1,x2,y2,x3,y3,x4,y4, X,Y,.)结果如图：当 n=20 时，我们在命令窗口中输入以下的命令：clear,clcX=-1:0.1:1;Y=1./(25*X.2+1); dx0= 0.0739644970414201;dxn= -0.0739644970414201;S=csfit(X,Y,dx0,dxn) x1=-1:0.01:-0.5;y1=polyval(S(1,:),x1-X(1);x2=-0.5:0.01:0;y2=polyval(S(2,:),x2-X(2); x3=0:0.01:0.5; y3=polyval(
8、S(3,:),x3-X(3);x4=0.5:0.01:1;y4=polyval(S(4,:),x4-X(4); plot(x1,y1,x2,y2,x3,y3,x4,y4, X,Y,.)结果如图：二、已知 Wilson 矩阵 ，且向量 ，则方程组1078569A321b有准确解 。Axb1Tx用 Matlab 内部函数求 ， 的所有特征值和 ；A2()condA令 ，解方程组 ，并求出078.2.5465.96.A()xb向量 和 ，从理论结果和实际计算结果两方面分析方程组 解的相x2 A对误差 与 的相对误差 的关系；A2A再改变扰动矩阵 （其元素的绝对值不超过 0.005） ，重复第 2 问
9、。解：（1） A=10 7 8 7;7 5 6 5;8 6 10 9;7 5 9 10;b=32;23;33;31; C=det(A) %求A的行列式C =1 X=eig(A) %求A的所有特征值X =0.01020.84313.858130.2887 D=cond(A) %求A的2-范数条件数D =2.9841e+003所以 A 的行列式为 1，cond(A) 2=2.9695e+003(2) 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; rank(B),rank(B,b)ans =4
10、4 x=Bb x = -5.476111.4934 -1.42922.4838求 假设 X= xx x=Bb;x1=1;1;1;1;X=x-x1 X =-6.476110.4934-2.42921.4838%求 2xnorm(X) ans =12.655212.655 就是 。2x%求 norm(X)/norm(x1)ans =6.32766.3276 即为 2x norm(a)ans =0.2244 norm(A)ans =30.2887 norm(a)/norm(A) %求 2Aans =0.0074所以 =0.00742Ainv(A) %求 A 的逆矩阵，下求 |)(|1Ax d=inv
11、(A); norm(d)*norm(a)*norm(x)ans =288.4858 1-norm(d*(a)ans = -12.3693 288.4858/ -12.3693ans =-23.3227所以 =-23.322|)(|1Axnorm(X) ans =0.0074所以 2x|)(|1Ax(3)改变 ，取 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.003B2=a2+A; C2=B bb=32;23;33;31 rank(B2) ans =4 rank(C2)an
12、s =4rank(B2)= rank(C2)所以扰动矩阵有唯一解 x2=B2bx2 =1.06490.88931.02720.9859 x=Bb;x1=1;1;1;1;X=x-x1 %求 （设 X= ）xxX2 =-6.476110.4934-2.42921.4838norm(X2) %求 2x norm(X2)ans =12.655212.6552 就是 2x norm(X)/norm(x1) %求 2x norm(X2)/norm(x1) ans =6.3276所以 =6.32762x norm(a2)ans =0.0071 norm(a)/norm(A) %求 2A norm(a2)/n
13、orm(A)ans =2.3336e-004所以 =2.3336e-0042A norm(d)*norm(a2)*norm(x2)ans =9.0875 1-norm(d*(a2)ans =0.6943norm(X2)ans =12.65529.0875/0.6943ans =13.0887所以 = 13.0887， zhuiganfa %调用追赶法三对角方程组的解为：x(1)=0.83333333x(2)=0.66666667x(3)=0.50000000x(4)=0.33333333x(5)=0.16666667（2）边值解M-文件：sinit=bvpinit(0:pi,2;0)odefu
14、n=inline(-y(2);abs(y(2)-y(1)+exp(t)-3*sin(t),t,y);bcfun=inline(ya(1)+2;yb(1)-exp(pi)-3,ya,yb);sol=bvp4c(odefun,bcfun,sinit)t=linspace(0,pi,101);y=deval(sol,t);plot(t,y(2,:),sol.x,sol.y(2,:),o,sinit.x,sinit.y(2,:),s)legend(解曲线, 解点 ,粗略解)指令窗口:sinit = x: 0 1 2 3y: 2x4 doublesol = x: 1x17 doubley: 2x17 d
15、oubleyp: 2x17 doublesolver: bvp4c四、公元 1225 年，比萨的数学家 Leonardo 研究了方程 ，3210xx得到一个根 ，没有人知道他用什么方法得到这个值。对于这个*1.36807x方程，分别用下列方法：迭代法 ；迭代法1201kkxx；对的 Steffensen 加速方法；对的 Steffensen 加速2310kkxx方法；Newton 法。求方程的根（可取 ） ，计算到 Leonardo 所得到的准01x确度。解：(1)由题意编写 m 文件如下：functionx0,k,er,x=diedai(g,x0,wucha,max)%g是给定的迭代函数%x0是给定的初始值%wucha是规定的误差范围%max是所应许的最大迭代次数%k是迭代次数加1%x是不动点近似值%x（ x1，x2.，xn）X(1)=x
