北京科技大学计算方法大作业概要
计算方法大作业机械电子工程系老师：廖福成注：本文本只有程序题，证明题全部在手写已交到理化楼204了。2. 证明方程 在上有一实根，并用二分法求这个根。要求。请给出程序和运行结果。证明：设f(x)=x3-x-1则f(1)= -1， f(2)= 5，f(1)*f(2)= -5<0因此，方程在1,2上必有一实根。二分法求解程序：%预先定义homework2.m文件如下：function lc=homework2(x)lc=x3-x-1;在MALAB窗口运行：cleara=1;b=2;tol=10(-3);N=10000;k=0; fa=homework2(a); % f 需事先定义for k=1:Np=(a+b)/2; fp=homework2(p);if( fp=0 | (b-a)/2<tol)breakendif fa*fp<0 b=p; else a=p; endendk,p程序运行结果：k = 103. 用Newton迭代法求方程 的一个正根，计算结果精确到7位有效数字. 要求给出程序和运行结果. 解：取迭代初值 ,并设,则牛顿迭代函数为牛顿迭格式为：Matlab程序如下：%定义zuoye3.m文件function x=zuoye3(fname,dfname,x0,e,N)if nargin<5,N=500;endif nargin<4,e=1e-7;endx=x0;x0=x+2*e;k=0;while abs(x0-x)>e&k<N, k=k+1; x0=x;x=x0-feval(fname,x0)/feval(dfname,x0); disp(x)endif k=N,warning('已达上限次数');end在Matlab窗口中执行：zuoye3(inline('x3+2*x2+10*x-20'),inline('3*x2+4*x+10'),1,1e-7)结果如下：1.369336470588241.368808188617531.36880810782137ans = 1.368808107821374. 用牛顿迭代法求方程在旁边的根. 要求给出程序和运行结果. 解：令：,则牛顿迭代函数为牛顿迭格式为：Matalb程序如下：%定义zuoye4.m文件function x=zuoye4(fname,dfname,x0,e,N)if nargin<5,N=500;endif nargin<4,e=1e-7;endx=x0;x0=x+2*e;k=0;while abs(x0-x)>e&k<N, k=k+1; x0=x;x=x0-feval(fname,x0)/feval(dfname,x0); disp(x)endif k=N,warning('已达上限次数');end在Matlab窗口执行：zuoye4(inline('x3-x-1'),inline('3*x2-1'),1,1e-7)结果如下：1.500000000000001.347826086956521.325200398950911.324718173999051.324717957244791.32471795724475ans = 1.324717957244756. 编写用全主元Gauss消去法解线性方程组的程序，并求解解：Matlab程序如下：A=2 -1 4 -3 1;-1 1 2 1 3;4 2 3 3 -1;-3 1 3 2 4;1 3 -1 4 4b=11 14 4 16 18function x=zuoye6(A,b)n,v=size(b);D=A,b;eye(n),zeros(n,v),s1,m=size(D);for k=1:(n-1) s=abs(A(k,k);p=k;q=k; for i=k:n for j=k:n if abs(A(i,j)>s s=abs(A(i,j);p=i;q=j; end end end if p>k t=D(k,:); D(k,:)=D(p,:); D(p,:)=t; end if q>k t1=D(:,k); D(:,k)=D(:,q); D(:,q)=t1; end h=D(k+1:n,k)/D(k,k); D(k+1:n,k+1:m)=D(k+1:n,k+1:m)-h*D(k,k+1:m); D(k+1:n,k)=zeros(n-k,1);endfor k=n:-1:1 D(k,k:m)=D(k,k:m)/D(k,k); for r=1:k-1 D(r,:)=D(r,:)-D(r,k)*D(k,:); endendP=D(n+1:2*n,1:n);Q=D(1:n,n+1:m);x=P*Q在Matlab窗口中执行：A=0.02 -1 4 -3 1;-1 1 2 1 3;4 2 3 3 -1;-3 1 3 2 4;1 3 -1 4 4;b=11 14 4 16 18'zuoye6(A,b)运行结果如下：x = -3.82352941176472 1.000000000000007. 用追逐法解线性方程组 要求给出程序和运行结果. 解：于是有求解Ax=b即为求解，式中b=(1 0 0 0 0)T据 y=据= x=Matlab程序如下：%定义zuoye7.m文件function x=zuoye7(a,b,c,d)a1=0;a;n=length(b);q=zeros(n,1);p=zeros(n,1);%LU分解q(1)=b(1);for k=2:n,p(k)=a1(k)/q(k-1); q(k)=b(k)-p(k)*c(k-1); end%解Ly=dy=zeros(n,1);y(1)=d(1);for k=2:n, y(k)=d(k)-p(k)*y(k-1);end%解Ux=yx=zeros(n,1); x(n)=y(n)/q(n);for k=n-1:-1:1,x(k)=(y(k)-c(k)*x(k+1)/q(k);endx在Matlab窗口中执行： a=-1 -1 -1 -1' b=2 2 2 2 2' c=-1 -1 -1 -1' d=1 0 0 0 0' x=zuoye7(a,b,c,d)运行结果如下：x = 0.83333333333333 0.66666666666667 0.50000000000000 0.333333333333339. 分别用Jacobi迭代法和Gauss-Seidel迭代法求解程组（编写程序）取初值，精确到小数后面四位。解：（1） Jacobi迭代法的Matlab程序：% x0为初始向量，ep为精度，N为最大次数，x是近似解向量A=10 3 1;2 -10 3;1 3 10;b=14 -5 14;n=length(b);N=500;ep=1e-6;x0=zero(n,1);n=length(b);x0=zeros(n,1);x=zeros(n,1);k=0while K<Nfor i=1:nx(i)=(b(i)-A(I,1:i-1,i+1:n)*x0(1:i-1,i+1:n)/A(i,i);endif norm(x-x0,inf)<ep,break;end;x0=x;k=k+1;endif k=N,Warning(已达到迭代次数上限);enddisp(k=,num2str(k),x运行结果如下：k=15x= 1.000000327906423 0.999999801006905 1.000000327906432（2）Gauss-Seidel迭代法的Matlab程序：% x0为初始向量，ep为精度，N为最大次数，x是近似解向量Format long;clear;A=10 3 1;2 -10 3;1 3 10;b=14 -5 14;n=length(b);N=500;ep=1e-6;x0=zero(n,1);P=inf;%以下是Guass-Seidal迭代法程序D=diag(diag(A);U=-triu(A,1);L=-tril(A,-1);dD=det(D);if dD=0disp(请留意：因为对角矩阵D奇异，所以此方程组无解。)elsedisp(请留意：因为对角矩阵D非奇异，所以此方程组无解。)iD=inv(D-L);B_GS=ID*U;d_GS=iD*b;endk=0;while k<Nx=B_GS*x0+d_GS;if norm(x-x0,P)<ep,break;endx0=x;k=k+1;endif k=N,warning(已达到迭代次数上限)；enddisp(k=,num2str(k),x,%D,U,L,B_GS,d_GS,%jx=Ab,运行结果如下：请留意：因为对角矩阵D非奇异，所以此方程组有解。k=9x=1.000000043651052 1.000000031224555 0.99999998626752810编写幂法程序求矩阵 按模最大的特征值和对应的特征向量。解：幂法程序如下：% A 为n 阶方阵，u0 为初始向量，%epsilon 为上限，max1 为循环次数。lambda 返回按模最大的特征值，u 返回对应的特征向量。clear;format long;clc;max1=1000; epsilon=1e-6;A=1 1 0.5;1 1 0.25;0.5 0.25 2; u0=1;1;1; u=u0;v0=u0;lambda=0;k=0;err=1;R=;while (k<max1)&(err>epsilon)v=A*u; m,j=max(abs(v); m=v(j); u=v/m;err=abs(lambda-m);lambda=m; k=k+1;Temp=k;m;v;u;R=R Temp;endk, lambda, u, disp(R), xlswrite('d:1.xls', R,'sheet1', 'b1') ;%D,Lmd=eig(A) %D 之列为A 的特征向量，Lmd 之对角线上为A 的特征值运行结果如下：k = 23lambda = 2.536527202038736u = 0.748222175139906 0.649662222855908 1.000000000000000
