矩阵和数值分析实验作业.docx

1、矩阵和数值分析实验作业一 求累加值原题：设，分别编制从小到大和从大到小的顺序程序计算,并指出有效位数。解：1.从小到大a. 程序代码：lejia1.mformat longN=input(请输入N=); s=0;for j=2:N s=s+1/(j*j-1);end%求近似值s1=s*10(5);s2=fix(s1);%取五位数字S=10(-5)*s2;a=0.5*(1.5-1/N-1/(N+1);y=abs(S-a);for i=0:100 if y=1/2*10(-i) break end%判断有效数字的位数enddisp(结果是:);Sdisp(有效数字的位数是:);i-1b. 运行结果

2、：2从大到小a.程序代码：leijia2.mformat longN=input(please input N:); s=0;for j=2:N s=s+1/(N-j+2)2-1;end%求近似值s1=s*10(5);s2=fix(s1);S=10(-5)*s2;%取五位数字a=0.5*(1.5-1/N-1/(N+1);y=abs(S-a);for i=0:100 if y=1/2*10(-i) break end%判断有效数字位数enddisp(结果是:);Sdisp(有效数字的位数是:);i-1b.运行结果：二 解线性方程组1. 原题：分别Jacobi迭代法和Gauss-Seidel迭代法

3、求解线性方程组迭代法计算停止条件为：。解：1.Jacobi迭代法：a. 程序代码：Jacobi.mA=input(请输入矩阵 A=);b=input(请输入列向量b=);X=input(请输入初始值X=);S=input(请输入终止条件S=);max1=100;%最大循环次数 B=diag(A); G=diag(B);%求对角阵 D=inv(G);%求对角阵的逆矩阵 P=G-A; for k=1:max1 Y=D*P*X+D*b; Z=abs(Y-X); t=max(Z); X=Y; If t.G-S迭代法：a. 程序代码：GS.mA=input(请输入矩阵 A=);b=input(请输入列向

4、量b=);X=input(请输入初始值X=);S=input(请输入终止条件S=);max1=100;%最大循环次数 C=diag(A); P=diag(C);%求A的对角阵 U=-triu(A)+P;%求A的下三角阵 G=A+U; D=inv(G);for k=1:max1 Y=D*U*X+D*b; Z=abs(Y-X); t=max(Z); X=Y; if(t.Gauss列主元消去法：a. 程序代码：GAUSS.mA=input(请输入系数矩阵 A=);b=input(请输入列向量b=);G=A,b;B=;%定义一个空矩阵t=size(b);%求系数矩阵的维数k=t(1,1);%求系数矩阵

5、的维数for i=1:k C=0; for e=i:k C(e)=abs(G(e,i); end l=max(C);%挑选最大值 l2=find(C=l); l1=l2(1,1); T=G(l1,:); G(l1,:)=G(i,:); G(i,:)=T;%选主元，交换最大值所在的行 disp(主元调整后的矩阵：); G for j=i+1:k P=-G(j,i)/G(i,i); B(j,:)=P*G(i,:)+G(j,:); G(j,:)=B(j,:); end%进行消元 disp(消元后的矩阵：); Gendb=G(:,k+1);G(:,k+1)=;X=;for m=1:k X(k-m+1)

6、=b(k-m+1)/G(k-m+1,k-m+1); b=b-X(k-m+1)*G(:,k-m+1);end%求方程组的解disp(方程组的解：);Xb. 运行结果：2.QR法：a. 程序代码：QR.mA=input(请输入系数矩阵 A=);b=input(请输入列向量b=);m=size(A);k=m(1,1);%计算A的维数B=A;M=eye(size(A);%定义一个与A同维数的单位矩阵Q=M;for r=1:k-1; I=eye(size(B); e=I(:,1); a=B(:,1);%取出B的第一列 d=norm(a); w=a-d*e; c=rot90(w); H=I-2*w*c/(

7、c*w);%求变换矩阵 P=H*B; B=P; B(1,:)=; B(:,1)=; for i=1:k-r+1 for j=1:k-r+1 M(r+i-1,r+j-1)=H(i,j); end end %将变换矩阵H变成k维 Z=M; M=eye(size(A); Q=Z*Q;endQ R=Q*A; RX=;b=Q*b;for u=1:k X(k-u+1)=b(k-u+1)/R(k-u+1,k-u+1); b=b-X(k-u+1)*R(:,k-u+1);end%计算方程组的解disp(方程组的解：);Xb. 运行结果：三 非线性方程的迭代解法1.原题：用Newton迭代法求解方程的根，计算停止

8、的条件为：。解：a.程序代码：newton.msyms x;f=input(请以x为变量输入函数f=);s=input(请输入终止条件s=);X=input(请输入初始值X=);max1=100;%最大循环次数 p=diff(f);%求f的导函数for k=1:max1 y=X-subs(f,X)/subs(p,X);%迭代公式 m=abs(y-X); X=y; if(ms)%终止条件 break endenddisp(方程的根：);Xa. 运行结果：2原题：利用Newton迭代法求多项式的所有实零点，注意重根的问题。解：a. 程序代码：new.mfunction X=newton1(X,s)

9、syms x;f=x4-5.4*x3+10.56*x2-8.954*x+2.7951;max1=20;p=diff(f);for k=1:max1 y=X-subs(f,X)/subs(p,X); m=abs(y-X); X=y; if(ms) break endend%定义一个newton迭代法求线性方程组的的函数syms x;f=input(请以x为变量输入多项式f=);c=input(请输入多项式的次数c=);a=input(请输入根所在的估计的最大区间的左端点a=);b=input(请输入根所在的估计的最大区间的右端点b=);n=input(请输入将a,b区间分成的个数n=);s=in

10、put(请输入newton搜索时的终止条件s=);A=a:(b-a)/n:b;%将区间n等分disp(n等分后的区间为:);Afor i=1:n Y(i)=subs(f,A(i); Y(i+1)=subs(f,A(i+1); if Y(i)*Y(i+1)=100*s break end%判断根的重数 end if e1 disp(根的重数：); X(j) disp(的重数：); disp(e); p=p+e-1;%对根重数大于1的重数进行累加 end if p=k-c break %判断重数是否已达到总重数 endendb. 运行结果：四 数值积分原题：分别用复化梯形公式和Romberg公式计

11、算积分要求误差不超过10-5，并给出节点个数。解：1.复化梯形公式：a. 程序代码：fuhuatixing.mformat long;syms x; f=input(请以x为自变量输入函数f=);a=input(请输入积分区域左端点a=);b=input(请输入积分区域右端点b=);s=input(请输入终止条件s=);n=100;X=a:(b-a)/n:b;t1=0;fa=subs(f,a);fb=subs(f,b);for i=2:n t1=t1+subs(f,X(i); T(i)=(b-a)*(fa+fb+2*t1)/(2*n);%用复化梯形公式计算积分结果 if (T(i)-T(i-1

12、).Romberg公式：a. 程序代码：romberg.msyms x; f=input(请以x为自变量输入函数f=);a=input(请输入积分区域左端点a=);b=input(请输入积分区域右端点b=);s=input(请输入终止条件s=);T(1,1)=1/2*(subs(f,a)+subs(f,b);n=100;for i=2:n q1=0; for e=1:(2(i-2) q1=q1+subs(f,(a+(2*(e-1)+1)*(b-a)/2(i-1); end q=(b-a)*q1/2(i-1); T(i,1)=1/2*T(i-1,1)+q;%计算复化梯形得出的值 for j=2:i T(i,j)=(4(j-1)*T(i,j-1)-T(i-1,j-1)/(4(j-1)-1);%计算n大于2时复化newton-cotes的值 end if abs(T(i,i)-T(i-1,i-1)s%判断条件 break endend disp(T数表:); Tdisp(积分结果:); T(i,i)b. 运行结果：