北京科技大学计算方法大作业

计算方法大作业机械电子工程系 老师：***注：本文本只有程序题，证明题全部在手写已交到理化楼204了。

2. 证明方程 310x x --=在[1,2]上有一实根*x ，并用二分法求这个根。

要求31||10k k x x -+-<。

请给出程序和运行结果。

证明：设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=x^3-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 = 10p = 1.3251953125000003. 用Newton 迭代法求方程 32210200x x x ++-=的一个正根，计算结果精确到7位有效数字. 要求给出程序和运行结果. 解： 取迭代初值01x = ,并设32()21020f x x x x =++-,则'2()3410f x x x =++． 牛顿迭代函数为32'2()21020()()3410f x x x x x x x f x x x ϕ++-=-=-++ 牛顿迭格式为：3212210203410k k k k k k k x x x x x x x +++-=-++ Matlab 程序如下： %定义zuoye3.m 文件function x=zuoye3(fname,dfname,x0,e,N) if nargin<5,N=500;end if nargin<4,e=1e-7;end x=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('x^3+2*x^2+10*x-20'),inline('3*x^2+4*x+10'),1,1e-7) 结果如下： 1.41176470588235 1.36933647058824 1.36880818861753 1.36880810782137 ans =1.368808107821374. 用牛顿迭代法求方程310x x --=在01x =附近的根. 要求给出程序和运行结果.解：令：3()1f x x x =--,则'2()31f x x =-． 牛顿迭代函数为3'2()1()()31f x x x x x x f x x ϕ--=-=-- 牛顿迭格式为：312131k k k k k x x x x x +--=-- Matalb 程序如下： %定义zuoye4.m 文件function x=zuoye4(fname,dfname,x0,e,N) if nargin<5,N=500;end if nargin<4,e=1e-7;end x=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('x^3-x-1'),inline('3*x^2-1'),1,1e-7)结果如下：1.500000000000001.347826086956521.325200398950911.324718173999051.324717957244791.32471795724475ans =1.324717957244756. 编写用全主元Gauss消去法解线性方程组的程序，并求解12345123451234512345123450.024*******4233433241634418x x x x x x x x x x x x x x x x x x x x x x x x x -+-+=⎧⎪-++++=⎪⎪+++-=⎨⎪-++++=⎪⎪+-++=⎩ 解：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 4] b=[11 14 4 16 18] function 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))>ss=abs(A(i,j));p=i;q=j; end end endif 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:1D(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 =2.94117647058824-3.823529411764721.000000000000000.941176470588245.941176470588247. 用追赶法解线性方程组12345 210001 121000 012100 001210 000120xxxxx-⎡⎤⎡⎤⎡⎤⎢⎥⎢⎥⎢⎥--⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥=--⎢⎥⎢⎥⎢⎥--⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥-⎣⎦⎣⎦⎣⎦要求给出程序和运行结果. 解：1000021000131000010021000221210024010000100121033001213500100001440001246001000055A LU -⎡⎤⎡⎤⎢⎥⎢⎥⎢⎥⎢⎥---⎡⎤⎢⎥⎢⎥⎢⎥--⎢⎥⎢⎥⎢⎥--⎢⎥⎢⎥⎢⎥===--⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥--⎢⎥--⎢⎥⎢⎥⎢⎥-⎣⎦⎢⎥⎢⎥⎢⎥⎢⎥-⎢⎥⎢⎥⎣⎦⎣⎦于是有求解Ax=b 即为求解{Ly bUx y ==，式中b=(1 0 0 0 0)T据12345100001100012020100033000104040015y y y y y ⎡⎤⎢⎥⎢⎥-⎡⎤⎡⎤⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥-⎢⎥⎢⎥⎢⎥=⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥-⎢⎥⎢⎥⎢⎥⎣⎦⎣⎦⎢⎥⎢⎥-⎢⎥⎣⎦ y=112131415⎡⎤⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦据21000301002400103500014600005-⎡⎤⎢⎥⎢⎥-⎢⎥⎢⎥-⎢⎥⎢⎥⎢⎥-⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦12345x x x x x ⎡⎤⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦=112131415⎡⎤⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦ x=5623121316⎡⎤⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦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);end x在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.833333333333330.666666666666670.500000000000000.333333333333330.166666666666679. 分别用Jacobi 迭代法和Gauss-Seidel 迭代法求解程组（编写程序）12312312310+3142-103-531014x x x x x x x x x +=⎧⎫⎪⎪+=⎨⎬⎪⎪++=⎩⎭取初值=TX （0）（0,0,0），精确到小数后面四位。

解：（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=0 while K<N for 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(‘已达到迭代次数上限’);end disp([‘k=’,num2str(k)]),x 运行结果如下： k=15 x=1.0000003279064230.9999998010069051.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奇异，所以此方程组无解。