大连理工大学-矩阵大作业

1、（1）从大到小的顺序的计算程序：function y=snd(n)format long y=0;if n<2 disp('请输入大于1的数！')else s=0; i=2;while i<=n s=single(s+(1/(i2-1); i=i+1; end y=s;end（2）从小到大的顺序的计算程序：function y=snx(n)format long y=0;if n<2 disp('请输入大于1的数！')else s=0; i=n; while 1 s=single(s+(1/(i2-1); i=i-1; if i=1 break

2、 end end y=s;end（3）按两种顺序分别计算并指出有效位数（编制程序时用单精度）1 的计算结果：2 的计算结果：3 的计算结果：计算时的有效位数为七位数。1 秦九昭算法计算程序：function y=qjz(a,x)j=3;i=size(a,2);switch i case 1 y=a(1); case 2 y=a(1)*x+a(2); otherwise p=a(1)*x+a(2); while j<=i p=p*x+a(j); j=j+1; end y=p;end2 计算在点23的值。计算结果如下：当时。1 Gauss法计算程序和结果：程序：A(1,:)=31,-13,0

3、,0,0,-10,0,0,0;A(2,:)=-13,35,-9,0,-11,0,0,0,0;A(3,:)=0,-9,31,-10,0,0,0,0,0;A(4,:)=0,0,-10,79,-30,0,0,0,-9;A(5,:)=0,0,0,-30,57,-7,0,-5,0;A(6,:)=0,0,0,0,-7,47,-30,0,0;A(7,:)=0,0,0,0,0,-30,41,0,0;A(8,:)=0,0,0,0,-5,0,0,27,-2;A(9,:)=0,0,0,-9,0,0,0,-2,29;B=-15;27;-23;0;-20;12;-7;7;10;a,b=gauss(A,B);j=size

4、(a,2);while j>=1 k=j+1; s=b(j); while k<=9 s=s-x(k)*a(j,k); k=k+1; end x(j)=s/a(j,j); j=j-1;enddisp(x)function x,y=gauss(a,b)num_i=size(a,1);j=1;while j<=(num_i-1) i=j+1; while i<=num_i r=a(i,j)/a(j,j); a(i,:)=a(i,:)-r*a(j,:); b(i,:)=b(i,:)-r*b(j,:); i=i+1; end j=j+1;endx=a;y=b;运行的结果为：。2

5、 列主元消去法计算程序和结果：A(1,:)=31,-13,0,0,0,-10,0,0,0;A(2,:)=-13,35,-9,0,-11,0,0,0,0;A(3,:)=0,-9,31,-10,0,0,0,0,0;A(4,:)=0,0,-10,79,-30,0,0,0,-9;A(5,:)=0,0,0,-30,57,-7,0,-5,0;A(6,:)=0,0,0,0,-7,47,-30,0,0;A(7,:)=0,0,0,0,0,-30,41,0,0;A(8,:)=0,0,0,0,-5,0,0,27,-2;A(9,:)=0,0,0,-9,0,0,0,-2,29;B=-15;27;-23;0;-20;12

6、;-7;7;10;a,b=lzy(A,B);j=size(a,2);while j>=1 k=j+1; s=b(j); while k<=9 s=s-x(k)*a(j,k); k=k+1; end x(j)=s/a(j,j); j=j-1;enddisp(x)function a1,b1=lzy(a,b)num_i,num_j=size(a);ab=zeros(num_i,num_j+1);for k=1:num_j ab(:,k)=a(:,k);endab(:,num_j+1)=b(:,1);j=1;while j<num_j max,max_i=searmax(j,j,a

7、b); i_nu=ab(j,:); ab(j,:)=ab(max_i,:); ab(max_i,:)=i_nu; m=j+1; while m<=num_i for n=j:num_j+1 ab(m,n)=ab(m,n)-(ab(m,j)/max)*ab(j,n); end m=m+1; end j=j+1;enda1=zeros(num_i,num_j);for k=1:num_i a1(:,k)=ab(:,k); endb1=ab(:,num_i+1);function b,c=searmax(i,j,a)num_i=size(a,1);k=i;m=abs(a(k,j);c=k;wh

8、ile k<num_i k=k+1; if m>=abs(a(k,j) continue else m=abs(a(k,j); c=k; endendb=a(c,j);运行的结果为：（1） LU分解的计算程序及结果：function l,u=lufz(a)num_i,num_j=size(a);l=eye(num_i);u=eye(num_i);for j=1:num_j u(1,j)=a(1,j); l(j,1)=a(j,1)/u(1,1);endi=2;while i<=num_i j=i; while j<num_i s=0; for k=1:i-1 s=s+l(

9、i,k)*u(k,j); end u(i,j)=a(i,j)-s; s=0; for k=1:i-1 s=s+l(j+1,k)*u(k,i); end l(j+1,i)=(a(j+1,i)-s)/u(i,i); j=j+1; end s=0; for k=1:i-1 s=s+l(i,k)*u(k,num_i); end u(i,num_i)=a(i,num_i)-s; i=i+1;end输入要求解的矩阵后运行的结果为：（2） 带列主元的LU分解计算程序及结果由于Matlab中自带带列主元的LU分解函数，故计算程序如下：A(1,:)=31,-13,0,0,0,-10,0,0,0;A(2,:)=-

10、13,35,-9,0,-11,0,0,0,0;A(3,:)=0,-9,31,-10,0,0,0,0,0;A(4,:)=0,0,-10,79,-30,0,0,0,-9;A(5,:)=0,0,0,-30,57,-7,0,-5,0;A(6,:)=0,0,0,0,-7,47,-30,0,0;A(7,:)=0,0,0,0,0,-30,41,0,0;A(8,:)=0,0,0,0,-5,0,0,27,-2;A(9,:)=0,0,0,-9,0,0,0,-2,29;L,U,P=lu(A);运行结果如下：为单位阵。（3） 求A的逆矩阵的计算程序及结果：A(1,:)=31,-13,0,0,0,-10,0,0,0;A

11、(2,:)=-13,35,-9,0,-11,0,0,0,0;A(3,:)=0,-9,31,-10,0,0,0,0,0;A(4,:)=0,0,-10,79,-30,0,0,0,-9;A(5,:)=0,0,0,-30,57,-7,0,-5,0;A(6,:)=0,0,0,0,-7,47,-30,0,0;A(7,:)=0,0,0,0,0,-30,41,0,0;A(8,:)=0,0,0,0,-5,0,0,27,-2;A(9,:)=0,0,0,-9,0,0,0,-2,29;num_i,num_j=size(A);A_n=eye(num_i);E=eye(num_i);L,U=lufz(A);for i=1

12、:num_i y=hd(1,L,E(:,i); A_n(:,i)=hd(0,U,y);enddisp(A_n)function l,u=lufz(a)num_i,num_j=size(a);l=eye(num_i);u=eye(num_i);for j=1:num_j u(1,j)=a(1,j); l(j,1)=a(j,1)/u(1,1);endi=2;while i<=num_i j=i; while j<num_i s=0; for k=1:i-1 s=s+l(i,k)*u(k,j); end u(i,j)=a(i,j)-s; s=0; for k=1:i-1 s=s+l(j+

13、1,k)*u(k,i); end l(j+1,i)=(a(j+1,i)-s)/u(i,i); j=j+1; end s=0; for k=1:i-1 s=s+l(i,k)*u(k,num_i); end u(i,num_i)=a(i,num_i)-s; i=i+1;endfunction x=hd(f,a,b)num_i,num_j=size(a);x=zeros(num_i,1);switch f case 0 i=num_i-1; x(num_i)=b(num_i)/a(num_i,num_j); while i>=1 s=0; for k=i+1:num_i s=s+a(i,k)*

14、x(k); end x(i)=(b(i)-s)/a(i,i); i=i-1; end case 1 x(1)=b(1)/a(1,1); j=2; while j<=num_i s=0; for k=1:j-1 s=s+a(j,k)*x(k); end x(j)=(b(j)-s)/a(j,j); j=j+1; end otherwise disp('error！请输入正确的参数！')end运行结果：（4） 求A的行列式的计算程序及结果：A(1,:)=31,-13,0,0,0,-10,0,0,0;A(2,:)=-13,35,-9,0,-11,0,0,0,0;A(3,:)=0,

15、-9,31,-10,0,0,0,0,0;A(4,:)=0,0,-10,79,-30,0,0,0,-9;A(5,:)=0,0,0,-30,57,-7,0,-5,0;A(6,:)=0,0,0,0,-7,47,-30,0,0;A(7,:)=0,0,0,0,0,-30,41,0,0;A(8,:)=0,0,0,0,-5,0,0,27,-2;A(9,:)=0,0,0,-9,0,0,0,-2,29;num_i,num_j=size(A);L,U=lufz(A);s=1;for i=1:num_i s=s*U(i,i);enddisp('矩阵的行列式值为：',num2str(s)运行程序后结果

16、与调用matlab中det（）函数结果如下：（1） 求cholesky分解程序及结果：function l=chole(a)num_i,num_j=size(a);l=zeros(num_i);l(1,1)=sqrt(a(1,1);for k=2:num_i l(k,1)=a(k,1)/l(1,1); endj=2;while j<=num_j s1=0; for k=1:j-1 s1=s1+l(j,k)2; end l(j,j)=sqrt(a(j,j)-s1); i=j+1; while i<=num_i s2=0; for k=1:j-1 s2=s2+l(i,k)*l(j,k)

17、; end l(i,j)=(a(i,j)-s2)/l(j,j); i=i+1; end j=j+1;end运行程序后的结果为：（2） 求解方程组程序及结果：a=2,1,-1,1;1,5,2,7;-1,2,10,1;1,7,1,11;b=13;-9;6;0;L=chole(a);y=hd(1,L,b);x=hd(0,L',y);disp(x)运行后的结果为:程序和结果如下：A=4,2,0,0;3,-2,1,0;0,2,5,3;0,0,-1,6;B=6;2;10;5;L,U=sjfj(A);%y=hd(1,L,B);x=hd(0,U,y);disp('方程组的解为：')di

18、sp(x)function l,u=sjfj(a)num_i=size(a,1);i=2;while i<=num_i a(i,i-1)=a(i,i-1)/a(i-1,i-1); a(i,i)=a(i,i)-a(i,i-1)*a(i-1,i); i=i+1;endl=tril(a);u=triu(a);for j=1:num_i; l(j,j)=1;end运行程序后结果为：编程求解矩阵A的QR分解：（1） QR分解计算程序：function Q,R=hqr(A)n,n=size(A);Q=eye(n);for i=1:(n-1) E=eye(n-i+1); e1=E(:,1); a=ze

19、ros(1,n-i+1)' for j=1:(n-i+1) a(j)=A(j+i-1,i); end av=sqrt(a'*a); w=a-av*e1; h=E-(2/(w'*w)*(w*w'); q=eye(n); for k=i:n for j=i:n q(k,j)=h(k-i+1,j-i+1); end end A=q*A; Q=q*Q;endR=A;Q=Q'（2） 输入矩阵A后的计算结果：（1） Jacobi迭代法求解程序与结果：a=0;b=0;c=0;while 1 a0=0.25*(7+b-c); b0=(1/8)*(-21-4*a-c);

20、c0=(1/5)*(15+2*a-b); A=abs(a-a0),abs(b-b0),abs(c-c0); f=max(A); if f<=0.001 x=a0,b0,c0; break else a=a0; b=b0; c=c0; endenddisp('方程组的解为x')disp(x)运行后的结果为：（2） Gauss-Seidel迭代法求解的程序与结果：a=0;b=0;c=0;while 1 a0=a; b0=b; c0=c; a=(1/4)*(7+b-c); b=(1/8)*(-21-4*a-c); c=(1/5)*(15+2*a-b); A=abs(a-a0),

21、abs(b-b0),abs(c-c0); f=max(A); if f<=0.001 x=a,b,c; break endenddisp('方程组的解为：')disp(x)运行后的结果为：（1） 计算在上的根的程序：syms xf=2*x3-5*x-1;df=diff(f,x);g=x-(f/df);x0=1;while 1 x1=subs(g,x0); dx=abs(x1-x0); if dx<=0.001 disp('Newton法的解为：x=',num2str(x1) break else x0=x1; endendx0=1;x1=1.1;wh

22、ile 1 x2=x1-(subs(f,x1)/(subs(f,x1)-subs(f,x0)*(x1-x0); dx=abs(x2-x1); if dx<=0.001 disp('割线法的解为：x=',num2str(x2) break else x0=x1; x1=x2; endend运行后结果为：Newton法，割线法。（2） 计算在上的根的程序：syms xf=exp(x)*sin(x);df=diff(f,x);g=x-(f/df);x0=-2;while 1 x1=subs(g,x0); dx=abs(x1-x0); if dx<=0.001 disp(&

23、#39;Newton法的解为：x=',num2str(x1) break else x0=x1; endendx0=-2;x1=-2.1;while 1 x2=x1-(subs(f,x1)/(subs(f,x1)-subs(f,x0)*(x1-x0); dx=abs(x2-x1); if dx<=0.001 disp('割线法的解为：x=',num2str(x2) break else x0=x1; x1=x2; endend运行后结果为：Newton法，割线法。syms xf=x*cos(x)-2;x0=-4;x2=-2;while 1x1=0.5*(x0+x2

24、);f0=subs(f,x0);f1=subs(f,x1);c=f0*f1;if c<0 x2=x1;else x0=x1;endl=abs(x2-x0);if l<=0.001 disp('二分法的解为x=',num2str(x0) breakendend运行后的结果为。（1）Lagrange插值程序：function yv=lagran(xd,yd,xv)syms xnum=length(xd);%计算节点的个数f=0;for k=1:num%创建lagrange插值函数 index=setdiff(1:num,k); f=f+yd(k)*prod(x-xd(i

25、ndex)./(xd(k)-xd(index);endyv=subs(f,xv);%计算xv点插值函数值（2） 绘制函数图程序：syms xf=1./(1+x.2);xd2=-5:2:5;yd2=subs(f,xd2);x0=-5:0.1:5;y=subs(f,x0);plot(x0,y)hold ony2=lagran(xd2,yd2,x0);plot(x0,y2)（3） 运行程序后得到的步长分别为2，1，的插值函数图与原图形的比较图如下：（1）三次样条插值程序：function yv=csi(xd,yd,xv,h,mf,ml)num_xd=length(xd);n=num_xd-1;a_m

26、=1:num_xd-2;for i=1:num_xd-2 a_m(i)=2;enda_np=1:num_xd-3;for i=1:num_xd-3 a_np(i)=0.5;endA=diag(a_m)+diag(a_np,-1)+diag(a_np,1);g=1:n-1;for k=1:n-1 g(k)=3*(0.5/h)*(yd(k+2)-yd(k+1)+(0.5/h)*(yd(k+1)-yd(k);endb=1:n-1;b(1)=g(1)-0.5*mf;b(n-1)=g(n-1)-0.5*ml;for k=2:n-2 b(k)=g(k); endb=b.'M=Ab;m=1:n+1;

27、m(1)=mf;m(n+1)=ml;for i=2:n m(i)=M(i-1); endnum_xv=length(xv);yv=1:num_xv;for i=1:num_xv for j=1:num_xd if xv(i)<xd(j) k1=j-1; k2=k1+1; break end end yv1=(h+2*(xv(i)-xd(k1)*(xv(i)-xd(k2)2*(yd(k1)/h3); yv2=(h-2*(xv(i)-xd(k2)*(xv(i)-xd(k1)2*(yd(k2)/h3); yv3=(xv(i)-xd(k1)*(xv(i)-xd(k2)2*(m(k1)/h2);

28、yv4=(xv(i)-xd(k2)*(xv(i)-xd(k1)2*(m(k2)/h2); yv(i)=yv1+yv2+yv3+yv4;end（2） 绘图程序（改变程序中的h值即可改变步长）：clear allclch=0.5;xd=-5:h:5;yd=1./(1+xd.2);xv=-5:0.1:5;y=1./(1+xv.2);mf=0.014793;ml=-0.014793;yv=csi(xd,yd,xv,h,mf,ml);plot(xv,y)hold onplot(xv,yv)（3） 步长为2，1，时原函数图与插值函数图的比较：（1） 复化梯形公式计算积分程序：function t=trap

29、e(f,a,b,n)h=(b-a)/n;index=(a+h):h:(b-h);t1=sum(subs(f,index);t=(b-a)/(2*n)*(subs(f,a)+2*t1+subs(f,b);（2） 复化Simpson公式计算程序：function s=simpson(f,a,b,n)h=(b-a)/n;index1=(a+0.5*h):h:(b-0.5*h);index2=(a+h):h:(b-h);s1=sum(subs(f,index1);s2=sum(subs(f,index2);s=(b-a)/(6*n)*(subs(f,a)+4*s1+2*s2+subs(f,b);（3） 区间数不同时的计算结果：（其中）（1） Gauss积分法计算积分程序：function s=gauss(f,a,b,n)switch n%Gauss case 2 x=-1/sqrt(3),1/sqrt(3); w=1,1; case 3 x=-0.7745966692,0,0.7745966692; w=0.5555555556,0.8888888889,0.5555555556; case 4 x=-0.8611363