当前位置：文档之家 › 大连理工矩阵上机作业

大连理工矩阵上机作业

  • 格式：docx
  • 大小：110.09 KB
  • 文档页数：12

下载文档原格式

  / 12

合集下载

大连理工大学矩阵与数值分析上机作业代码

大连理工大学矩阵与数值分析上机作业代码
则方程组有解 x = (1,1,⋯ ,1) 。 对 n = 10 和 n = 84 ， 分别用 Gauss 消去法和列主元消去法解
T
方程组，并比较计算结果。
1.1 程序
（1）高斯消元法 n=10 时， >> [A,b]=CreateMatrix(10) A= 6 8 0 0 0 0 0 0 0 0 1 6 8 0 0 0 0 0 0 0 0 1 6 8 0 0 0 0 0 0 0 0 1 6 8 0 0 0 0 0 0 0 0 1 6 8 0 0 0 0 0 0 0 0 1 6 8 0 0 0 0 0 0 0 0 1 6 8 0 0 0 0 0 0 0 0 1 6 8 0 0 0 0 0 0 0 0 1 6 8 0 0 0 0 0 0 0 0 1 6
1.3 M 文件
.m 1.3.1 CreateMatrix CreateMatrix.m function [A,b]=CreateMatrix(n) %用于存放习题1的题目信息，并建立构造题目中矩阵的函数 %对矩阵A赋值 A1=6*ones(1,n); A2=ones(1,n-1); A3=8*ones(1,n-1); A=diag(A1)+diag(A2,1)+diag(A3,-1); %对向量b赋值 b=15*ones(n,1); b(1)=7; b(n)=14;
10
，迭代次数上限取默认值
50，使用 Jacobi 法进行迭代。 >> test2 >> b=ones(20,1) >> x0=zeros(20,1) >> [x,n]=JacobiMethod(A,b,x0) x= 0.2766 0.2327 0.2159 0.2223 0.2227 0.2221 0.2222 0.2222 0.2222 0.2222 0.2222 0.2222 0.2222 0.2222 0.2221 0.2227 0.2223 0.2159 0.2327 0.2766

大连理工大学线性代数实验上机报告材料

大连理工大学线性代数实验上机报告材料

大连理工大学线性代数实验上机报告实验一首先随机生成五阶方阵AA=rand(5)A =0.8147 0.0975 0.1576 0.1419 0.6557 0.9058 0.2785 0.9706 0.4218 0.0357 0.1270 0.5469 0.9572 0.9157 0.8491 0.9134 0.9575 0.4854 0.7922 0.9340 0.6324 0.9649 0.8003 0.9595 0.6787>> B=rand(5)随机生成五阶方阵BB =0.7577 0.7060 0.8235 0.4387 0.4898 0.7431 0.0318 0.6948 0.3816 0.4456 0.3922 0.2769 0.3171 0.7655 0.6463 0.6555 0.0462 0.9502 0.7952 0.7094 0.1712 0.0971 0.0344 0.1869 0.7547>> b=rand(1,5)'随机生成列向量bb =0.27600.67970.65510.16260.1190计算A+B>> A+Bans =1.5725 0.8036 0.9811 0.5806 1.1455 1.6489 0.3103 1.6654 0.8033 0.48130.5192 0.8238 1.2743 1.6813 1.49541.5689 1.0037 1.4356 1.5874 1.6434 0.8035 1.0620 0.8347 1.1464 1.4334 计算A-B>> A-Bans =0.0570 -0.6085 -0.6658 -0.2969 0.1660 0.1627 0.2467 0.2758 0.0402 -0.4099 -0.2652 0.2700 0.6401 0.1502 0.2028 0.2579 0.9113 -0.4648 -0.0030 0.2246 0.4612 0.8678 0.7658 0.7726 -0.0760 计算A*B+B*A>> A*B+B*Aans =3.0288 2.3058 3.1439 2.7276 3.10342.9094 2.19673.0040 3.0737 3.25843.3422 2.1423 3.2104 3.5734 3.90494.1446 2.9794 4.3676 4.2354 4.9170 3.1350 1.7787 3.2289 3.1170 3.2815 求Ax=b的解>> x=A\bx =-0.98502.43963.3124-5.65151.7085验证克莱姆法则>> c=A(:,1)c =0.81470.90580.12700.91340.6324>> d=A(:,2)d =0.09750.5469 0.9575 0.9649>> e=A(:,3)e =0.1576 0.9706 0.9572 0.4854 0.8003>> f=A(:,4)f =0.1419 0.4218 0.91570.9595>> g=A(:,5)g =0.65570.03570.84910.93400.6787>> B1=[b';d';e';f';g']'B1 =0.2760 0.0975 0.1576 0.1419 0.6557 0.6797 0.2785 0.9706 0.4218 0.0357 0.6551 0.5469 0.9572 0.9157 0.8491 0.1626 0.9575 0.4854 0.7922 0.9340 0.1190 0.9649 0.8003 0.9595 0.6787>> B2=[c';b';e';f';g']'B2 =0.8147 0.2760 0.1576 0.1419 0.6557 0.9058 0.6797 0.9706 0.4218 0.0357 0.1270 0.6551 0.9572 0.9157 0.8491 0.9134 0.1626 0.4854 0.7922 0.9340 0.6324 0.1190 0.8003 0.9595 0.6787>> B3=[c';d';b';f';g']'B3 =0.8147 0.0975 0.2760 0.1419 0.6557 0.9058 0.2785 0.6797 0.4218 0.0357 0.1270 0.5469 0.6551 0.9157 0.8491 0.9134 0.9575 0.1626 0.7922 0.9340 0.6324 0.9649 0.1190 0.9595 0.6787>> B4=[c';d';e';b';g']'B4 =0.8147 0.0975 0.1576 0.2760 0.6557 0.9058 0.2785 0.9706 0.6797 0.0357 0.1270 0.5469 0.9572 0.6551 0.8491 0.9134 0.9575 0.4854 0.1626 0.9340 0.6324 0.9649 0.8003 0.1190 0.6787>> B5=[c';d';e';f';b']'B5 =0.8147 0.0975 0.1576 0.1419 0.2760 0.9058 0.2785 0.9706 0.4218 0.6797 0.1270 0.5469 0.9572 0.9157 0.6551 0.9134 0.9575 0.4854 0.7922 0.1626 0.6324 0.9649 0.8003 0.9595 0.1190>> x1=det(B1)/det(A)x1 =-0.9850>> x2=det(B2)/det(A) x2 =2.4396>> x3=det(B3)/det(A) x3 =3.3124>> x4=det(B4)/det(A) x4 =-5.6515>> x5=det(B5)/det(A)x5 =1.7085计算A的行列式>> det(A)ans =-0.0250计算B的行列式>> det(B)ans =0.0647求A的逆>> inv(A)ans =3.1375 -0.8078 -1.8788 -4.21945.1680-8.6076 3.5314 2.8907 13.7204 -14.3665 -6.2824 3.7220 3.6132 10.0084 -12.4190 13.6173 -6.8822 -6.3938 -23.5288 27.5825 -2.5292 1.0729 2.4193 5.8870 -7.2671 求B的逆>> inv(B)ans =-0.4430 3.4997 1.3255 -2.6005 -0.4697 1.4047 -1.1626 0.2422 -0.4475 -0.0119 0.7210 -1.8189 -2.0635 2.4434 0.0765 -0.6122 -0.1837 2.0165 0.0375 -1.2564 0.0384 -0.5157 -0.7370 0.5267 1.7407 求A的秩>> rank(A)ans =5求B的秩>> rank(B)ans =5求A*B的行列式>> det(A*B)ans =-0.0016求A*B的逆>> inv(A*B)ans =-74.0649 35.0433 31.2288 121.5740 -137.3442 6.8291 -1.2718 -2.2922 -8.9951 8.6972 63.9620 -31.4202 -29.5061 -105.6918 122.3246 -9.3196 5.7452 4.6259 11.9660 -15.4028 11.9582 -6.3521 -3.3817 -16.7574 18.6360>> rank(A*B)ans =5>> det(A)*det(B)ans =-0.0016验证 （1）>> (A*B)'ans =0.9569 1.5566 1.6237 2.2732 2.25520.6922 0.9401 0.4969 0.9371 0.80900.9461 1.6492 1.6875 2.3563 2.38000.7507 1.5887 1.8840 1.9421 2.14811.1399 1.52122.2149 2.4545 2.4497()()111,,T T T AB B A AB B A AB BA---==≠>> B'*A'ans =0.9569 1.5566 1.6237 2.2732 2.2552 0.6922 0.9401 0.4969 0.9371 0.8090 0.9461 1.6492 1.6875 2.3563 2.38000.7507 1.5887 1.8840 1.9421 2.14811.1399 1.52122.2149 2.4545 2.4497 （2）>> inv(B)*inv(A)ans =-74.0649 35.0433 31.2288 121.5740 -137.3442 6.8291 -1.2718 -2.2922 -8.9951 8.6972 63.9620 -31.4202 -29.5061 -105.6918 122.3246 -9.3196 5.7452 4.6259 11.9660 -15.4028 11.9582 -6.3521 -3.3817 -16.7574 18.6360 （3）>> A*Bans =0.9569 0.6922 0.9461 0.7507 1.13991.5566 0.9401 1.6492 1.5887 1.52121.6237 0.4969 1.6875 1.88402.21492.2732 0.9371 2.3563 1.9421 2.4545 2.2552 0.8090 2.3800 2.1481 2.4497>> B*Aans =2.0719 1.6135 2.1978 1.9769 1.9635 1.3528 1.2566 1.3549 1.4850 1.7372 1.7186 1.6454 1.5229 1.6894 1.6900 1.8714 2.0423 2.0113 2.2932 2.4625 0.8797 0.9697 0.8489 0.9690 0.8317 求矩阵X使得AXB=C首先随机生成五阶方阵C>> C=rand(5)C =0.4984 0.7513 0.9593 0.8407 0.35000.9597 0.2551 0.5472 0.2543 0.19660.3404 0.5060 0.1386 0.8143 0.25110.5853 0.6991 0.1493 0.2435 0.61600.2238 0.8909 0.2575 0.9293 0.4733X=A 的逆*B 的逆>> X=inv(A)*C*inv(B)X =3.8432 -13.8858 2.1418 9.4404 -4.5871-9.3312 41.9602 -7.9101 -28.4683 14.8942-7.8738 35.1218 -5.4107 -22.8861 10.158116.7545 -75.6079 14.6784 49.3951 -24.7450-3.5568 17.0848 -2.9018 -11.2670 5.4559实验二1. 验证：对于一般的方阵A,B,C,D ，首先随机生成方阵A,B,C,D A=rand(5)A BA DB CC D ≠-A =0.8258 0.1067 0.8687 0.4314 0.1361 0.5383 0.9619 0.0844 0.9106 0.8693 0.9961 0.0046 0.3998 0.1818 0.5797 0.0782 0.7749 0.2599 0.2638 0.5499 0.4427 0.8173 0.8001 0.1455 0.1450>> B=rand(5)B =0.8530 0.0760 0.4173 0.4893 0.7803 0.6221 0.2399 0.0497 0.3377 0.3897 0.3510 0.1233 0.9027 0.9001 0.2417 0.5132 0.1839 0.9448 0.3692 0.4039 0.4018 0.2400 0.4909 0.1112 0.0965>> C=rand(5)C =0.1320 0.2348 0.1690 0.5470 0.1835 0.9421 0.3532 0.6491 0.2963 0.3685 0.9561 0.8212 0.7317 0.7447 0.6256 0.5752 0.0154 0.6477 0.1890 0.7802 0.0598 0.0430 0.4509 0.6868 0.0811>> D=rand(5)D =0.9294 0.3063 0.6443 0.9390 0.2077 0.7757 0.5085 0.3786 0.8759 0.3012 0.4868 0.5108 0.8116 0.5502 0.4709 0.4359 0.8176 0.5328 0.6225 0.2305 0.4468 0.7948 0.3507 0.5870 0.8443>> Z=[A,B;C,D]Z =0.8258 0.1067 0.8687 0.4314 0.13610.8530 0.0760 0.4173 0.4893 0.78030.5383 0.9619 0.0844 0.9106 0.86930.6221 0.2399 0.0497 0.3377 0.38970.9961 0.0046 0.3998 0.1818 0.57970.3510 0.1233 0.9027 0.9001 0.24170.0782 0.7749 0.2599 0.2638 0.54990.5132 0.1839 0.9448 0.3692 0.40390.4427 0.8173 0.8001 0.1455 0.14500.4018 0.2400 0.4909 0.1112 0.09650.1320 0.2348 0.1690 0.5470 0.18350.9294 0.3063 0.6443 0.9390 0.20770.9421 0.3532 0.6491 0.2963 0.36850.7757 0.5085 0.3786 0.8759 0.30120.9561 0.8212 0.7317 0.7447 0.62560.4868 0.5108 0.8116 0.5502 0.47090.5752 0.0154 0.6477 0.1890 0.78020.4359 0.8176 0.5328 0.6225 0.23050.0598 0.0430 0.4509 0.6868 0.08110.4468 0.7948 0.3507 0.5870 0.8443求Z的行列式>> det(Z)ans =-0.0295求det(A)*det(D)-det(B)*det(C)>> det(A)*det(D)-det(B)*det(C)ans =1.8656e-004随机生成对角矩阵A>> A=diag([rand rand rand rand rand])A =0.1948 0 0 0 0 0 0.2259 0 0 0 0 0 0.1707 0 0 0 0 0 0.2277 0 0 0 0 0 0.4357 随机生成对角矩阵B>> B=diag([rand rand rand rand rand])B =0.3111 0 0 0 0 0 0.9234 0 0 0 0 0 0.4302 0 0 0 0 0 0.1848 0 0 0 0 0 0.9049 随机生成对角矩阵C>> C=diag([rand rand rand rand rand])C =0.9797 0 0 0 0 0 0.4389 0 0 0 0 0 0.1111 0 0 0 0 0 0.2581 0 0 0 0 0 0.4087 随机生成对角矩阵D>> D=diag([rand rand rand rand rand])D =0.5949 0 0 0 00 0.2622 0 0 00 0 0.6028 0 00 0 0 0.7112 00 0 0 0 0.2217>> Z=[A,B;C,D]Z =0.1948 0 0 0 00.3111 0 0 0 00 0.2259 0 0 00 0.9234 0 0 00 0 0.1707 0 00 0 0.4302 0 00 0 0 0.2277 0 0 0 0 0.1848 00 0 0 0 0.4357 0 0 0 0 0.90490.9797 0 0 0 00.5949 0 0 0 00 0.4389 0 0 00 0.2622 0 0 00 0 0.1111 0 00 0 0.6028 0 00 0 0 0.2581 00 0 0 0.7112 00 0 0 0 0.40870 0 0 0 0.2217计算Z的行列式>> det(Z)ans =-1.1243e-004计算det(A)*det(D)-det(B)*det(C)>> det(A)*det(D)-det(B)*det(C)ans =-9.3107e-005计算A*D-B*C的行列式>> det(A*D-B*C)ans =-1.1243e-004实验三求A列向量组的一个最大无关组，并把不属于极大无关组的向量利用极大无关组表示.N= 200865083;a=83;b=86;c=50;d=88;e=28;f=63;g=83;h=60;>>A=[a,b,c,d,3,4;1,2,3,4,4,3;12,15,22,17,5,7;e,f,g,h, 8,0];>> B=rref(A)B =1.0000 0 0 0 -0.3548 0.46560 1.0000 0 0 -1.4905 -2.00200 0 1.0000 0 0.0473 0.39500 0 0 1.0000 1.79841.3383所以a1,a2,a3,a4是一个极大无关组。

大连理工矩阵上机作业

大连理工矩阵上机作业

第一题Lagrange插值函数function y=lagrange(x0,y0,x);n=length(x0);m=length(x);for i=1:mz=x(i);s=0.0;for k=1:np=1.0;for j=1:nif j~=kp=p*(z-x0(j))/(x0(k)-x0(j));endends=p*y0(k)+s;endy(i)=s;endx0=[1:10];y0=[67.052,68.008,69.803,72.024,73.400,72.063,74.669,74.487,74.065,76 .777];lagrange(x0,y0,17)ans=-1.9516e+12x=[1:0.1:10];x=x';plot(x0,y0,'r');hold onplot(x,y,'k');legend('原函数','拟合函数')拟合图像如下拟合函数出现了龙格现象，运用多项式进行插值拟合时，效果并不好，高次多项式会因为误差的不断累积，导致龙格现象的发生。

第二题function fun =nihe(n)m=[67.052*10^6,68.008*10^6,69.803*10^6,72.024*10^6,73.400*10^6,72.063 *10^6,74.669*10^6,74.487*10^6,74.065*10^6,76.777*10^6];w=[1,2,3,4,5,6,7,8,9,10];d1=0;d2=0;d3=0;y1=polyfit(m,w,1);y2=polyfit(m,w,2);y3=polyfit(m,w,3);y2=poly2sym(s2);y3=poly2sym(s3);y4=poly2sym(s4);f1=subs(y1,17);f2=subs(y2,17);f3=subs(y3,17);for p=1:10;d1=d1+(subs(y1,w(p))-m(p))^2;d2=d2+(subs(y2,w(p))-m(p))^2;d3=d3+(subs(y3,w(p))-m(p))^2;endd1=sqrt(d1);d2=sqrt(d2);d3=sqrt(d3);fun=[f1 f2 f3;d2 d3 d4];return;结果三次函数的均方误差最小，拟合的最好。

大连理工大学优化方法上机作业

大连理工大学优化方法上机作业

大连理工大学优化方法上机作业本页仅作为文档页封面，使用时可以删除This document is for reference only-rar21year.March优化方法上机大作业学院：电子信息与电气工程学部姓名：学号：指导老师：上机大作业（一）%目标函数function f=fun(x)f=100*(x(2)-x(1)^2)^2+(1-x(1))^2;end%目标函数梯度function gf=gfun(x)gf=[-400*x(1)*(x(2)-x(1)^2)-2*(1-x(1));200*(x(2)-x(1)^2)]; End%目标函数Hess矩阵function He=Hess(x)He=[1200*x(1)^2-400*x(2)+2,-400*x(1);-400*x(1), 200;];end%线搜索步长function mk=armijo(xk,dk)beta=0.5; sigma=0.2;m=0; maxm=20;while (m<=maxm)if(fun(xk+beta^m*dk)<=fun(xk)+sigma*beta^m*gfun(xk)'*dk) mk=m; break;endm=m+1;endalpha=beta^mknewxk=xk+alpha*dkfk=fun(xk)newfk=fun(newxk)%最速下降法function [k,x,val]=grad(fun,gfun,x0,epsilon)%功能：梯度法求解无约束优化问题：minf(x)%输入：fun，gfun分别是目标函数及其梯度，x0是初始点，% epsilon为容许误差%输出：k是迭代次数，x，val分别是近似最优点和最优值maxk=5000; %最大迭代次数beta=0.5; sigma=0.4;k=0;while(k<maxk)gk=feval(gfun,x0); %计算梯度dk=-gk; %计算搜索方向if(norm(gk)<epsilon), break;end%检验终止准则m=0;mk=0;while(m<20) %用Armijo搜索步长if(feval(fun,x0+beta^m*dk)<=feval(fun,x0)+sigma*beta^m*gk'*dk) mk=m;break;endm=m+1;endx0=x0+beta^mk*dk;k=k+1;endx=x0;val=feval(fun,x0);>> x0=[0;0];>> [k,x,val]=grad('fun','gfun',x0,1e-4)迭代次数：k =1033x =0.99990.9998val =1.2390e-008%牛顿法x0=[0;0];ep=1e-4;maxk=10;k=0;while(k<maxk)gk=gfun(x0);if(norm(gk)<ep)x=x0miny=fun(x)k0=kbreak;elseH=inv(Hess(x0));x0=x0-H*gk;k=k+1;endendx =1.00001.0000miny =4.9304e-030迭代次数k0 =2%BFGS方法function [k,x,val]=bfgs(fun,gfun,x0,varargin) %功能：梯度法求解无约束优化问题：minf(x)%输入：fun，gfun分别是目标函数及其梯度，x0是初始点，% epsilon为容许误差%输出：k是迭代次数，x，val分别是近似最优点和最优值N=1000;epsilon=1e-4;beta=0.55;sigma=0.4;n=length(x0);Bk=eye(n);k=0;while(k<N)gk=feval(gfun,x0,varargin{:});if(norm(gk)<epsilon), break;enddk=-Bk\gk;m=0;mk=0;while(m<20)newf=feval(fun,x0+beta^m*dk,varargin{:});oldf=feval(fun,x0,varargin{:});if(newf<=oldf+sigma*beta^m*gk'*dk)mk=m;break;endm=m+1;endx=x0+beta^mk*dk;sk=x-x0;yk=feval(gfun,x,varargin{:})-gk;if(yk'*sk>0)Bk=Bk-(Bk*sk*sk'*Bk)/(sk'*Bk*sk)+(yk*yk')/(yk'*sk);endk=k+1;x0=x;endval=feval(fun,x0,varargin{:});>> x0=[0;0];>> [k,x,val]=bfgs('fun','gfun',x0)k =20x =1.00001.0000val =2.2005e-011%共轭梯度法function [k,x,val]=frcg(fun,gfun,x0,epsilon,N)if nargin<5,N=1000;endif nargin<4, epsilon=1e-4;endbeta=0.6;sigma=0.4;n=length(x0);k=0;while(k<N)gk=feval(gfun,x0);itern=k-(n+1)*floor(k/(n+1));itern=itern+1;if(itern==1)dk=-gk;elsebetak=(gk'*gk)/(g0'*g0);dk=-gk+betak*d0; gd=gk'*dk;if(gd>=0),dk=-gk;endendif(norm(gk)<epsilon),break;endm=0;mk=0;while(m<20)if(feval(fun,x0+beta^m*dk)<=feval(fun,x0)+sigma*beta^m*gk'*dk) mk=m;break;endm=m+1;endx=x0+beta^m*dk;g0=gk; d0=dk;x0=x;k=k+1;endval=feval(fun,x);>> x0=[0;0];[k,x,val]=frcg('fun','gfun',x0,1e-4,1000)k =122x =1.00011.0002val =7.2372e-009上机大作业（二）%目标函数function f_x=fun(x)f_x=4*x(1)-x(2)^2-12;%等式约束条件function he=hf(x)he=25-x(1)^2-x(2)^2;end%不等式约束条件function gi_x=gi(x,i)switch icase 1gi_x=10*x(1)-x(1)^2+10*x(2)-x(2)^2-34;case 2gi_x=x(1);case 3gi_x=x(2);otherwiseend%求目标函数的梯度function L_grad=grad(x,lambda,cigma)d_f=[4;2*x(2)];d_g(:,1)=[-2*x(1);-2*x(2)];d_g(:,2)=[10-2*x(1);10-2*x(2)];d_g(:,3)=[1;0];d_g(:,4)=[0;1];L_grad=d_f+(lambda(1)+cigma*hf(x))*d_g(:,1);for i=1:3if lambda(i+1)+cigma*gi(x,i)<0L_grad=L_grad+(lambda(i+1)+cigma*gi(x,i))*d_g(:,i+1);continueendend%增广拉格朗日函数function LA=lag(x,lambda,cee)LA=fun(x)+lambda(1)*hf(x)+0.5*cee*hf(x)^2;for i=1:3LA=LA+1/(2*cee)*(min(0,lambda(i+1)+cee*gi(x,i))^2-lambda(i+1)^2); endfunction xk=BFGS(x0,eps,lambda,cigma)gk=grad(x0,lambda,cigma);res_B=norm(gk);k_B=0;a_=1e-4;rho=0.5;c=1e-4;length_x=length(x0);I=eye(length_x);Hk=I;while res_B>eps&&k_B<=10000dk=-Hk*gk;m=0;while m<=5000if lag(x0+a_*rho^m*dk,lambda,cigma)-lag(x0,lambda,cigma)<=c*a_*rho^m*gk'*dkmk=m;break;endm=m+1;endak=a_*rho^mk;xk=x0+ak*dk;delta=xk-x0;y=grad(xk,lambda,cigma)-gk;Hk=(I-(delta*y')/(delta'*y))*Hk*(I-(y*delta')/(delta'*y))+(delta*delta')/(delta'*y);k_B=k_B+1;x0=xk;gk=y+gk;res_B=norm(gk);end%增广拉格朗日法function val_min=ALM(x0,eps)lambda=zeros(4,1);cigma=5;alpha=10;k=1;res=[abs(hf(x0)),0,0,0];for i=1:3res(1,i+1)=norm(min(gi(x0,i),-lambda(i+1)/cigma)); endres=max(res);while res>eps&&k<1000xk=BFGS(x0,eps,lambda,cigma);lambda(1)=lambda(1)+cigma*hf(xk);for i=1:3lambda(i+1)=lambda(i+1)+min(0,lambda(i+1)+gi(x0,1)); endk=k+1;cigma=alpha*cigma;x0=xk;res=[norm(hf(x0)),0,0,0];for i=1:3res(1,i+1)=norm(min(gi(x0,i),-lambda(i+1)/cigma)); endres=max(res);endval_min=fun(xk);fprintf('k=%d\n',k);fprintf('fmin=%.4f\n',val_min);fprintf('x=[%.4f;%.4f]\n',xk(1),xk(2));>> x0=[0;0];>> val_min=ALM(x0,1e-4)k=10fmin=-31.4003x=[1.0984;4.8779]val_min =-31.4003上机大作业（三）A=[1 1;-1 0;0 -1];n=2;b=[1;0;0];G=[0.5 0;0 2];c=[2 4];cvx_solver sdpt3cvx_beginvariable x(n)minimize (x'*G*x-c*x)subject toA*x<=bcvx_enddisp(x)Status: SolvedOptimal value (cvx_optval): -2.40.40000.6000A=[2 1 1;1 2 3;2 2 1;-1 0 0;0 -1 0;0 0 -1]; n=3;b=[2;5;6;0;0;0];C=[-3 -1 -3];cvx_solver sdpt3cvx_beginvariable x(n)minimize (C*x)subject toA*x<=bcvx_enddisp(x)Status: SolvedOptimal value (cvx_optval): -5.40.20000.00001.600011。

大连理工大学矩阵与数值分析上机作业

大连理工大学矩阵与数值分析上机作业
s=s+abs(x(i));
end
case2%2-范数
fori=1:n
s=s+x(i)^2;
end
s=sqrt(s);
caseinf%无穷-范数
s=max(abs(x));
end
计算向量x，y的范数
Test1.m
clearall;
clc;
n1=10;n2=100;n3=1000;
x1=1./[1:n1]';x2=1./[1:n2]';x3=1./[1:n3]';
xlabel('x');ylabel('p(x)');
运行结果：
x=2的邻域：
x =
1.6000 1.8000 2.0000 2.2000 2.4000
相应多项式p值：
p =
1.0e-003 *
-0.2621 -0.0005 0 0.0005 0.2621
p(x)在 [1.95,20.5]上的图像
程序：
[L,U]=LUDe.(A);%LU分解
xLU=U\(L\b)
disp('利用PLU分解方程组的解:');
[P,L,U] =PLUDe.(A);%PLU分解
xPLU=U\(L\(P\b))
%求解A的逆矩阵
disp('A的准确逆矩阵:');
InvA=inv(A)
InvAL=zeros(n);%利用LU分解求A的逆矩阵
0 0 0.5000 -0.2500 -0.1250 -0.0625 -0.0625
0 0 0 0.5000 -0.2500 -0.1250 -0.1250
0 0 0 0 0.5000 -0.2500 -0.2500

大连理工大学矩阵分析matlab上机作业

大连理工大学矩阵分析matlab上机作业
x=zeros(n,1); %为列向量 x 预分配存储空间 y=1:n; %定义行向量 y y=y'; %把行向量 y 改成列向量 for i=1:n
x(i)=1/i; %按要求给向量 x 赋值，其值递减 end normx1=norm(x,1); %求解向量 x 的 1 范数 normx1 normx2=norm(x,2); %求解向量 x 的 2 范数 normx2 normxinf=norm(x,inf); %求解向量 x 的无穷范数 normxinf normy1=norm(y,1); %求解向量 y 的 1 范数 normy1 normy2=norm(y,2); %求解向量 y 的 2 范数 normy2 normyinf=norm(y,inf); %求解向量 y 的无穷范数 normyinf z1=[normx1,normx2,normxinf]; z2=[normy1,normy2,normyinf]; end
for i=2:n
for j=i:n U(i,j)=A(i,j)-L(i,1:i-1)*U(1:i-1,j);
式
%Doolittle 分解计算上三角矩阵的公
L(j,i)=(A(j,i)-L(j,1:i-1)*U(1:i-1,i))/U(i,i); %Doolittle 分解计算下三角矩 阵的公式
end
1 1 1 ������ x = (1, 2 , 3 , … , ������) ,
������ = (1,2, … , ������)������.
对n = 10，100，1000甚至更大的n计算其范数，你会发现什么结果?你能否修改
你的程序使得计算结果相对精确呢？
1.1 源代码
function [z1,z2]=norm_vector(n) %向量 z1 的值为向量 x 的是三种范数，向量 z2 的值为向量 y 的三 种范数，n 为输入参数

大连理工_2012矩阵与数值分析大作业

矩阵与数值分析学生：学号：任课老师：金光日教学班号：（2）班院系：电子信息与电气工程学部《矩阵与数值分析》课程数值实验题目1.给定n 阶方程组A x b =，其中6186186186A ⎛⎫ ⎪ ⎪⎪= ⎪ ⎪ ⎪⎝⎭，7151514b ⎛⎫ ⎪⎪ ⎪= ⎪ ⎪⎪⎝⎭则方程组有解(1,1,,1)T x = 。

对10n =和84n =，分别用Gauss 消去法和列主元消去法解方程组，并比较计算结果。

1答： 程序1. Gauss 消元法function x=DelGauss(A,b) % Gauss 消去法 [n,m]=size(A); det=1; %存储行列式值 x=zeros(n,1); for k=1:n-1 for i=k+1:n if A(k,k)==0 return endm=A(i,k)/A(k,k); for j=k+1:nA(i,j)=A(i,j)-m*A(k,j); endb(i)=b(i)-m*b(k); enddet=det*A(k,k); %计算行列式enddet=det*A(n,n);for k=n:-1:1 %回代求解for j=k+1:nb(k)=b(k)-A(k,j)*x(j);endx(k)=b(k)/A(k,k);end2. 列主元Gauss消去法：function x=detGauss(A,b)% Gauss列主元消去法[n,m]=size(A);nb=length(b);det=1; %存储行列式值x=zeros(n,1);for k=1:n-1amax=0; %选主元for i=k:nif abs(A(i,k))>amaxamax=abs(A(i,k));r=i;endendif amax<1e-10return;endif r>k %交换两行for j=k:nz=A(k,j);A(k,j)=A(r,j);A(r,j)=z;endz=b(k);b(k)=b(r);b(r)=z;det=-det;endfor i=k+1:n %进行消元m=A(i,k)/A(k,k);for j=k+1:nA(i,j)=A(i,j)-m*A(k,j);endb(i)=b(i)-m*b(k);enddet=det*A(k,k);enddet=det*A(n,n);for k=n:-1:1 %回代求解for j=k+1:nb(k)=b(k)-A(k,j)*x(j);endx(k)=b(k)/A(k,k);end矩阵A和b的构造clc;clear;n=10;%n=84;A=eye(n)*6+diag(ones(1,n-1)*8,-1)+diag(ones(1,n-1),1); b=[7,15*ones(1,n-2),14]';计算结果：(1)n=10时Gauss消元法>>x=DelGauss(A,b)x =1.00001.00001.00001.00001.00001.00001.00001.00001.00001.0000列主元Gauss消去法>>x=detGauss(A,b)x =1111111111(2) n=84时Gauss消元法>>x=DelGauss(A,b) x =1.0e+008 *0.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 -0.0000 0.0000 -0.0000 0.0000 -0.0000 0.0000 -0.0000 0.0000 -0.0000 0.0000 -0.0000 0.0000 -0.0001 0.0002 -0.0003 0.0007 -0.0013 0.0026 -0.0052 0.0105 -0.0209 0.0419 -0.0836 0.16650.6501-1.25822.3487-4.02635.3684列主元Gauss消去法>>x=detGauss(A,b) x =1.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.00001.0000 1.0000 1.00001.00001.0000 1.0000结果分析由上述实验结果可知，对于n=10采用Gauss 消去法和Gauss 列主元消去法得到的实验结果是相同的，而对于n=84，Gauss 消去法所得到的结果是错误的，Gauss 列主元消去法得到的结果是正确的。

大工矩阵上机题答案

上机题一．（1）s=0;for j=2:100;s=s+1/(j^2-1); endss =0.7400s=0;for j=100:-1:2;s=s+1/(j^2-1); endss =0.7400（2）s=0;for j=2：10000;s=s+1/(j^2-1); endss =0.7499for j=10000:-1:2;s=s+1/(j^2-1); endss =0.7499（3）s=0;for j=2：1000000;s=s+1/(j^2-1); endss =0.7500s=0;for j=1000000:-1:2; s=s+1/(j^2-1); endss =0.7500二1、Jacobi 迭代法算法：对于线性方程组Ax=b ，如果A 为非奇异方程，则可将A 分解为：A=D-L-U 其中D 为对角阵，其元素为A 的对角元素，L 与U 为A 的下三角阵和上三角阵。

于是Ax=b 化为：111()k k x D L U x D b --+=++，其中1()J B D L U -=+，1f D b -=。

程序：function x=jacobi(A,b,x0)A=[-2 1 0 0;1 -2 1 0;0 1 -2 1;0 0 1 -2];b=[-1 ;0; 0; 0];x0=[0;0 ;0 ;0];D=diag(diag(A));U=-triu(A,1);L=-tril(A,-1);B=D\(L+U);f=D\b;x=B*x0+f;n=1;while norm(x-x0,2)>=1.0e-6x0=x;x=B*x0+f;n=n+1;endfprintf('迭代次数为：')nfprintf('方程组的解为：')计算结果：迭代次数为：n =60方程组的解为：ans =0.80000.60000.40000.2000.Gauss-seidel 迭代法function x=Gaussseidel(A,b,x0)A=[-2 1 0 0;1 -2 1 0;0 1 -2 1;0 0 1 -2];b=[-1 ;0; 0; 0];x0=[0;0 ;0 ;0];D=diag(diag(A));U=-triu(A,1);L=-tril(A,-1);B=(D-L)\U;f=(D-L)\b;x=B*x0+f;n=1;while norm(x-x0,2)>=1.0e-6x0=x;x=B*x0+f;n=n+1;endfprintf('迭代次数')nfprintf('方程组的解为')迭代次数为：n =31方程组的解为：ans =0.80000.60000.40000.20002.用Gauss列主元消去法算法：Gauss列主元消去法是在Gauss消去法中增加选主元的过程，即在第k步（k=1,2,3,…）消元时，首先在第k列主对角元以下（含对角元）元素中挑选绝对值最大的数（即为列主元），并通过初等行变换，使得该数位于主对角线上，然后再继续消元。

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

（1）从大到小的顺序的计算程序：function y=snd(n)format longy=0;if n<2disp('请输入大于1的数！')elses=0;i=2;while i<=ns=single(s+(1/(i^2-1))); i=i+1;endy=s;end（2）从小到大的顺序的计算程序：function y=snx(n)format longy=0;if n<2disp('请输入大于1的数！')elses=0;i=n;while 1s=single(s+(1/(i^2-1)));i=i-1;if i==1breakendendy=s;end（3）按两种顺序分别计算并指出有效位数（编制程序时用单精度）S的计算结果：①210S的计算结果：②410S的计算结果：③610计算时的有效位数为七位数。

① 秦九昭算法计算程序：function y=qjz(a,x) j=3; i=size(a,2); switch i case 1 y=a(1); case 2y=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; end② 计算在点23的值。

计算结果如下：当23=x 时()86652=x f 。

①Gauss法计算程序和结果：程序：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;-7;7;10]; [a,b]=gauss(A,B);j=size(a,2);while j>=1k=j+1;s=b(j);while k<=9s=s-x(k)*a(j,k);k=k+1;endx(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_ir=a(i,j)/a(j,j);a(i,:)=a(i,:)-r*a(j,:);b(i,:)=b(i,:)-r*b(j,:);i=i+1;endj=j+1;endx=a;y=b;运行的结果为：()T 289.0-345.0.0=。

大连理工大学线性代数实验上机报告材料

大连理工大学线性代数实验上机报告实验一首先随机生成五阶方阵AA=rand(5)A =0.8147 0.0975 0.1576 0.1419 0.6557 0.9058 0.2785 0.9706 0.4218 0.0357 0.1270 0.5469 0.9572 0.9157 0.8491 0.9134 0.9575 0.4854 0.7922 0.9340 0.6324 0.9649 0.8003 0.9595 0.6787>> B=rand(5)随机生成五阶方阵BB =0.7577 0.7060 0.8235 0.4387 0.4898 0.7431 0.0318 0.6948 0.3816 0.4456 0.3922 0.2769 0.3171 0.7655 0.6463 0.6555 0.0462 0.9502 0.7952 0.7094 0.1712 0.0971 0.0344 0.1869 0.7547>> b=rand(1,5)'随机生成列向量bb =0.27600.67970.65510.16260.1190计算A+B>> A+Bans =1.5725 0.8036 0.9811 0.5806 1.1455 1.6489 0.3103 1.6654 0.8033 0.48130.5192 0.8238 1.2743 1.6813 1.49541.5689 1.0037 1.4356 1.5874 1.6434 0.8035 1.0620 0.8347 1.1464 1.4334 计算A-B>> A-Bans =0.0570 -0.6085 -0.6658 -0.2969 0.1660 0.1627 0.2467 0.2758 0.0402 -0.4099 -0.2652 0.2700 0.6401 0.1502 0.2028 0.2579 0.9113 -0.4648 -0.0030 0.2246 0.4612 0.8678 0.7658 0.7726 -0.0760 计算A*B+B*A>> A*B+B*Aans =3.0288 2.3058 3.1439 2.7276 3.10342.9094 2.19673.0040 3.0737 3.25843.3422 2.1423 3.2104 3.5734 3.90494.1446 2.9794 4.3676 4.2354 4.9170 3.1350 1.7787 3.2289 3.1170 3.2815 求Ax=b的解>> x=A\bx =-0.98502.43963.3124-5.65151.7085验证克莱姆法则>> c=A(:,1)c =0.81470.90580.12700.91340.6324>> d=A(:,2)d =0.09750.5469 0.9575 0.9649>> e=A(:,3)e =0.1576 0.9706 0.9572 0.4854 0.8003>> f=A(:,4)f =0.1419 0.4218 0.91570.9595>> g=A(:,5)g =0.65570.03570.84910.93400.6787>> B1=[b';d';e';f';g']'B1 =0.2760 0.0975 0.1576 0.1419 0.6557 0.6797 0.2785 0.9706 0.4218 0.0357 0.6551 0.5469 0.9572 0.9157 0.8491 0.1626 0.9575 0.4854 0.7922 0.9340 0.1190 0.9649 0.8003 0.9595 0.6787>> B2=[c';b';e';f';g']'B2 =0.8147 0.2760 0.1576 0.1419 0.6557 0.9058 0.6797 0.9706 0.4218 0.0357 0.1270 0.6551 0.9572 0.9157 0.8491 0.9134 0.1626 0.4854 0.7922 0.9340 0.6324 0.1190 0.8003 0.9595 0.6787>> B3=[c';d';b';f';g']'B3 =0.8147 0.0975 0.2760 0.1419 0.6557 0.9058 0.2785 0.6797 0.4218 0.0357 0.1270 0.5469 0.6551 0.9157 0.8491 0.9134 0.9575 0.1626 0.7922 0.9340 0.6324 0.9649 0.1190 0.9595 0.6787>> B4=[c';d';e';b';g']'B4 =0.8147 0.0975 0.1576 0.2760 0.6557 0.9058 0.2785 0.9706 0.6797 0.0357 0.1270 0.5469 0.9572 0.6551 0.8491 0.9134 0.9575 0.4854 0.1626 0.9340 0.6324 0.9649 0.8003 0.1190 0.6787>> B5=[c';d';e';f';b']'B5 =0.8147 0.0975 0.1576 0.1419 0.2760 0.9058 0.2785 0.9706 0.4218 0.6797 0.1270 0.5469 0.9572 0.9157 0.6551 0.9134 0.9575 0.4854 0.7922 0.1626 0.6324 0.9649 0.8003 0.9595 0.1190>> x1=det(B1)/det(A)x1 =-0.9850>> x2=det(B2)/det(A) x2 =2.4396>> x3=det(B3)/det(A) x3 =3.3124>> x4=det(B4)/det(A) x4 =-5.6515>> x5=det(B5)/det(A)x5 =1.7085计算A的行列式>> det(A)ans =-0.0250计算B的行列式>> det(B)ans =0.0647求A的逆>> inv(A)ans =3.1375 -0.8078 -1.8788 -4.21945.1680-8.6076 3.5314 2.8907 13.7204 -14.3665 -6.2824 3.7220 3.6132 10.0084 -12.4190 13.6173 -6.8822 -6.3938 -23.5288 27.5825 -2.5292 1.0729 2.4193 5.8870 -7.2671 求B的逆>> inv(B)ans =-0.4430 3.4997 1.3255 -2.6005 -0.4697 1.4047 -1.1626 0.2422 -0.4475 -0.0119 0.7210 -1.8189 -2.0635 2.4434 0.0765 -0.6122 -0.1837 2.0165 0.0375 -1.2564 0.0384 -0.5157 -0.7370 0.5267 1.7407 求A的秩>> rank(A)ans =5求B的秩>> rank(B)ans =5求A*B的行列式>> det(A*B)ans =-0.0016求A*B的逆>> inv(A*B)ans =-74.0649 35.0433 31.2288 121.5740 -137.3442 6.8291 -1.2718 -2.2922 -8.9951 8.6972 63.9620 -31.4202 -29.5061 -105.6918 122.3246 -9.3196 5.7452 4.6259 11.9660 -15.4028 11.9582 -6.3521 -3.3817 -16.7574 18.6360>> rank(A*B)ans =5>> det(A)*det(B)ans =-0.0016验证 （1）>> (A*B)'ans =0.9569 1.5566 1.6237 2.2732 2.25520.6922 0.9401 0.4969 0.9371 0.80900.9461 1.6492 1.6875 2.3563 2.38000.7507 1.5887 1.8840 1.9421 2.14811.1399 1.52122.2149 2.4545 2.4497()()111,,T T T AB B A AB B A AB BA---==≠>> B'*A'ans =0.9569 1.5566 1.6237 2.2732 2.2552 0.6922 0.9401 0.4969 0.9371 0.8090 0.9461 1.6492 1.6875 2.3563 2.38000.7507 1.5887 1.8840 1.9421 2.14811.1399 1.52122.2149 2.4545 2.4497 （2）>> inv(B)*inv(A)ans =-74.0649 35.0433 31.2288 121.5740 -137.3442 6.8291 -1.2718 -2.2922 -8.9951 8.6972 63.9620 -31.4202 -29.5061 -105.6918 122.3246 -9.3196 5.7452 4.6259 11.9660 -15.4028 11.9582 -6.3521 -3.3817 -16.7574 18.6360 （3）>> A*Bans =0.9569 0.6922 0.9461 0.7507 1.13991.5566 0.9401 1.6492 1.5887 1.52121.6237 0.4969 1.6875 1.88402.21492.2732 0.9371 2.3563 1.9421 2.4545 2.2552 0.8090 2.3800 2.1481 2.4497>> B*Aans =2.0719 1.6135 2.1978 1.9769 1.9635 1.3528 1.2566 1.3549 1.4850 1.7372 1.7186 1.6454 1.5229 1.6894 1.6900 1.8714 2.0423 2.0113 2.2932 2.4625 0.8797 0.9697 0.8489 0.9690 0.8317 求矩阵X使得AXB=C首先随机生成五阶方阵C>> C=rand(5)C =0.4984 0.7513 0.9593 0.8407 0.35000.9597 0.2551 0.5472 0.2543 0.19660.3404 0.5060 0.1386 0.8143 0.25110.5853 0.6991 0.1493 0.2435 0.61600.2238 0.8909 0.2575 0.9293 0.4733X=A 的逆*B 的逆>> X=inv(A)*C*inv(B)X =3.8432 -13.8858 2.1418 9.4404 -4.5871-9.3312 41.9602 -7.9101 -28.4683 14.8942-7.8738 35.1218 -5.4107 -22.8861 10.158116.7545 -75.6079 14.6784 49.3951 -24.7450-3.5568 17.0848 -2.9018 -11.2670 5.4559实验二1. 验证：对于一般的方阵A,B,C,D ，首先随机生成方阵A,B,C,D A=rand(5)A BA DB CC D ≠-A =0.8258 0.1067 0.8687 0.4314 0.1361 0.5383 0.9619 0.0844 0.9106 0.8693 0.9961 0.0046 0.3998 0.1818 0.5797 0.0782 0.7749 0.2599 0.2638 0.5499 0.4427 0.8173 0.8001 0.1455 0.1450>> B=rand(5)B =0.8530 0.0760 0.4173 0.4893 0.7803 0.6221 0.2399 0.0497 0.3377 0.3897 0.3510 0.1233 0.9027 0.9001 0.2417 0.5132 0.1839 0.9448 0.3692 0.4039 0.4018 0.2400 0.4909 0.1112 0.0965>> C=rand(5)C =0.1320 0.2348 0.1690 0.5470 0.1835 0.9421 0.3532 0.6491 0.2963 0.3685 0.9561 0.8212 0.7317 0.7447 0.6256 0.5752 0.0154 0.6477 0.1890 0.7802 0.0598 0.0430 0.4509 0.6868 0.0811>> D=rand(5)D =0.9294 0.3063 0.6443 0.9390 0.2077 0.7757 0.5085 0.3786 0.8759 0.3012 0.4868 0.5108 0.8116 0.5502 0.4709 0.4359 0.8176 0.5328 0.6225 0.2305 0.4468 0.7948 0.3507 0.5870 0.8443>> Z=[A,B;C,D]Z =0.8258 0.1067 0.8687 0.4314 0.13610.8530 0.0760 0.4173 0.4893 0.78030.5383 0.9619 0.0844 0.9106 0.86930.6221 0.2399 0.0497 0.3377 0.38970.9961 0.0046 0.3998 0.1818 0.57970.3510 0.1233 0.9027 0.9001 0.24170.0782 0.7749 0.2599 0.2638 0.54990.5132 0.1839 0.9448 0.3692 0.40390.4427 0.8173 0.8001 0.1455 0.14500.4018 0.2400 0.4909 0.1112 0.09650.1320 0.2348 0.1690 0.5470 0.18350.9294 0.3063 0.6443 0.9390 0.20770.9421 0.3532 0.6491 0.2963 0.36850.7757 0.5085 0.3786 0.8759 0.30120.9561 0.8212 0.7317 0.7447 0.62560.4868 0.5108 0.8116 0.5502 0.47090.5752 0.0154 0.6477 0.1890 0.78020.4359 0.8176 0.5328 0.6225 0.23050.0598 0.0430 0.4509 0.6868 0.08110.4468 0.7948 0.3507 0.5870 0.8443求Z的行列式>> det(Z)ans =-0.0295求det(A)*det(D)-det(B)*det(C)>> det(A)*det(D)-det(B)*det(C)ans =1.8656e-004随机生成对角矩阵A>> A=diag([rand rand rand rand rand])A =0.1948 0 0 0 0 0 0.2259 0 0 0 0 0 0.1707 0 0 0 0 0 0.2277 0 0 0 0 0 0.4357 随机生成对角矩阵B>> B=diag([rand rand rand rand rand])B =0.3111 0 0 0 0 0 0.9234 0 0 0 0 0 0.4302 0 0 0 0 0 0.1848 0 0 0 0 0 0.9049 随机生成对角矩阵C>> C=diag([rand rand rand rand rand])C =0.9797 0 0 0 0 0 0.4389 0 0 0 0 0 0.1111 0 0 0 0 0 0.2581 0 0 0 0 0 0.4087 随机生成对角矩阵D>> D=diag([rand rand rand rand rand])D =0.5949 0 0 0 00 0.2622 0 0 00 0 0.6028 0 00 0 0 0.7112 00 0 0 0 0.2217>> Z=[A,B;C,D]Z =0.1948 0 0 0 00.3111 0 0 0 00 0.2259 0 0 00 0.9234 0 0 00 0 0.1707 0 00 0 0.4302 0 00 0 0 0.2277 0 0 0 0 0.1848 00 0 0 0 0.4357 0 0 0 0 0.90490.9797 0 0 0 00.5949 0 0 0 00 0.4389 0 0 00 0.2622 0 0 00 0 0.1111 0 00 0 0.6028 0 00 0 0 0.2581 00 0 0 0.7112 00 0 0 0 0.40870 0 0 0 0.2217计算Z的行列式>> det(Z)ans =-1.1243e-004计算det(A)*det(D)-det(B)*det(C)>> det(A)*det(D)-det(B)*det(C)ans =-9.3107e-005计算A*D-B*C的行列式>> det(A*D-B*C)ans =-1.1243e-004实验三求A列向量组的一个最大无关组，并把不属于极大无关组的向量利用极大无关组表示.N= 200865083;a=83;b=86;c=50;d=88;e=28;f=63;g=83;h=60;>>A=[a,b,c,d,3,4;1,2,3,4,4,3;12,15,22,17,5,7;e,f,g,h, 8,0];>> B=rref(A)B =1.0000 0 0 0 -0.3548 0.46560 1.0000 0 0 -1.4905 -2.00200 0 1.0000 0 0.0473 0.39500 0 0 1.0000 1.79841.3383所以a1,a2,a3,a4是一个极大无关组。

  1. 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
  2. 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
  3. 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。

第一题Lagrange插值函数function y=lagrange(x0,y0,x);n=length(x0);m=length(x);for i=1:mz=x(i);s=0.0;for k=1:np=1.0;for j=1:nif j~=kp=p*(z-x0(j))/(x0(k)-x0(j));endends=p*y0(k)+s;endy(i)=s;endx0=[1:10];y0=[67.052,68.008,69.803,72.024,73.400,72.063,74.669,74.487,74.065,76 .777];lagrange(x0,y0,17)ans= -1.9516e+12x=[1:0.1:10];x=x';plot(x0,y0,'r');hold onplot(x,y,'k');legend('原函数','拟合函数')拟合图像如下拟合函数出现了龙格现象，运用多项式进行插值拟合时，效果并不好，高次多项式会因为误差的不断累积，导致龙格现象的发生。

第二题function fun =nihe(n)m=[67.052*10^6,68.008*10^6,69.803*10^6,72.024*10^6,73.400*10^6,72.063 *10^6,74.669*10^6,74.487*10^6,74.065*10^6,76.777*10^6];w=[1,2,3,4,5,6,7,8,9,10];d1=0;d2=0;d3=0;y1=polyfit(m,w,1);y2=polyfit(m,w,2);y3=polyfit(m,w,3);y2=poly2sym(s2);y3=poly2sym(s3);y4=poly2sym(s4);f1=subs(y1,17);f2=subs(y2,17);f3=subs(y3,17);for p=1:10;d1=d1+(subs(y1,w(p))-m(p))^2;d2=d2+(subs(y2,w(p))-m(p))^2;d3=d3+(subs(y3,w(p))-m(p))^2;endd1=sqrt(d1);d2=sqrt(d2);d3=sqrt(d3);fun=[f1 f2 f3;d2 d3 d4];return;结果三次函数的均方误差最小，拟合的最好。

函数function f=fun(x)syms aa=x;f=a*a*a+a*a+a-3;梯度函数function df=dfun(x)df=3*x*x+2*x+1;Newton法function result=didainewton(x0)k=0;xk=x0;xi=1;e0=abs(x0-xi);ek=e0;m=zeros(7,1);n=zeros(7,1);p=zeros(7,1);result=zeros(7,3);while k<7ak=feval('fun',xk);bk=feval('dfun',xk);xk=xk-ak/bk;e0=ek;k=k+1;m(k)=xk;ek=abs(m(k)-xi);jingdu=ek/(e0*e0);n(k)=ek;p(k)=jingdu;endresult=[m,n,p];return;计算结果Newton迭代法的收敛速度快，至少是平方收敛的，GAUSS消去法function x=DelGauss(N)%GaussÏûÈ¥·¨syms M;M=N;a=zeros(M);b=ones(M,1);for i=1:Mfor j=1:Ma(i,j)=1/(i+j-1);endend[n,m]=size(a);nb=length(b); det=1;x=zeros(n,1);for k=1:n-1for i=k+1:nif a(k,k)==0returnendm=a(i,k)/a(k,k);for j=k+1:na(i,j)=a(i,j)-m*a(k,j);endb(i)=b(i)-m*b(k);enddet=det*a(k,k);enddet=det*a(n,n);for k=n:-1:1for j=k+1:nb(k)=b(k)-a(k,j)*x(j);endx(k)=b(k)/a(k,k)endx=x'N=4X=[-4.0000 60.0000 -180.0000 140.0000]TN=8X=[ -8.0000001, 504.00001, -7560.0, 46200.0, -138600.0, 216216.0, -168168.0, 51480.0]TJacobbifunction x=jacobi(N)syms MM=N;A=zeros(M);b=ones(M,1);for i=1:Mfor j=1:MA(i,j)=1/(i+j-1);endendp=zeros(M,1);q=zeros(M,1);for i=1:Mq(i)=A(i,i);endD=diag(q);L=-tril(A,-1);U=-triu(A,1);B=(inv(D))*(L+U);f=(inv(D))*b;x=B*p+f;n=1;sigmal=1e-8;while norm((p-x))>=sigmalp=x;x=B*p+f;n=n+1;endeval('x');N=4时R（B）=2.5821>1，迭代法不收敛N=8时R（B）= 6.0422>1,迭代法不收敛Gauss-Seidelfunction x=gaussseidel(N)syms MM=N;A=zeros(M);b=ones(M,1);for i=1:Mfor j=1:MA(i,j)=1/(i+j-1);endendp=zeros(M,1);L=-tril(A,-1);D=diag(diag(A));U=-triu(A,1);G=(inv(D-L))*U;f=(inv(D-L))*b;x=G*p+f;n=1;detal=1e-4;while abs(x-p)>=detalp=x;x=G*p+f;n=n+1;endeval('x');结果N=4X=[2.8527 -14.7133 -3.5388 26.7350]TN=8X=[-3.9202 28.0256 -9.4820 -41.3197 -34.0806 -6.4475 28.8462 65.3427]T共轭梯度法function result=cg(N)syms MM=N;x0=zeros(M,1);A=zeros(M);b=ones(M,1);for i=1:Mfor j=1:MA(i,j)=1/(i+j-1);endendk=0;r0=b-A*x0;p0=r0;sigmal=1e-4;rk=r0;pk=p0;while dot(r0,r0)>=sigmalak=(dot(r0,r0))/(dot(pk,A*pk));ri=r0;xk=x0+ak*pk;rk=r0-ak*A*pk;bk=(dot(rk,rk))/(dot(ri,ri));pk=rk+bk*p0;p0=pk;r0=rk;x0=xk;k=k+1;endresult=x0;return;结果N=4X=[-4.0000 60.0000 -180.0000 140.0000]TN=8X=[ 4.5773 -72.4211 199.8408 -9.9954 -174.8581 -171.8511 -10.8652 269.4734]T第五题function f=san(x1,y)n=length(x1);h=zeros(1,n-1);for p=1:n-1h(p)=x1(p+1)-x1(p);endfor p=1:n-2r(p)=h(p+1)/(h(p+1)+h(p));u(p+1)=h(p)/(h(p+1)+h(p));endu(1)=1;r(n-1)=1;c=2*ones(1,n);A=diag(r,-1)+diag(u,1)+diag(c);g=zeros(n,1);g(1)=3*(y(2)-y(1))/h(1);g(n)=3*(y(n)-y(n-1))/h(n-1);for p=1;n-2;g(p+1)=3*((y(p+2)-y(p+1))/h(p+1)*u(p+1)+r(p)*(y(p+1)-y(p))/h(p)); endm=A\g;m=m';f=m;returnx1=[0:4];y=[1,3,3,4,2];san(x1,y)ans =2.5000 1.0000 -0.5000 1.0000 -3.5000利用教科书p179的公式（5-35）function fun=sanci(m,h,xk,y)syms xn=length(m);l=length(xk);for i=1:n-1fun=((h(i)+2*(x-xk(i)))*y(i)/h(i)+(x-xk(i))*m(i))*(x-xk(i+1))^2/h(i)^2+((h(i)-2*(x-xk(i+1)))*y(i+1)/h(i)+(x-xk(i+1))*m(i+1))*(x-xk(i))^2/h(i)^2endend0<=x<1((9*x)/2 + 1)*(x - 1)^2 - x^2*(5*x - 8)1<=x<2(7*x - 4)*(x - 2)^2 - ((13*x)/2 - 16)*(x - 1)^22<=x<3((11*x)/2 - 8)*(x - 3)^2 - (7*x - 25)*(x - 2)^23<=x<4(9*x - 23)*(x - 4)^2 - ((15*x)/2 - 32)*(x - 3)^2第六题函数function f=fun(x)syms a;a=x;f=a*a*a+2*a*a+10*a-100;end梯度函数function grad=gfun(x)syms aa=x;grad=3*a*a+4*a+10;迭代函数f unction result=xianjiefa(x0)sigmal=1e-6;maxi=500;m=zeros(500,1);d0=fun(x0)/gfun(x0);x1=x0-d0;k=1;xk=x1;m(1)=x1;while k<maxidk=fun(xk)*(xk-x0)/(fun(xk)-fun(x0)); x0=xk;xk=xk-dk;k=k+1;m(k)=xk;if abs(xk-x0)<sigmalbreak;endendn=zeros(k,1);for i=1:kn(i)=m(i);endresult=(n);第七题function fun=intgra(n)syms x;syms m;sum1=0;sum2=0;m=n;f=(exp(3*x))*cos(pi*x);x0=0:2*pi/m:2*pi;a=0;b=2*pi;h=(b-a)/m;for i=1:m-1xk=a+i*h;sum1=sum1+subs(f,xk);endfor i=0:n-1xk=a+(i+1/2)*h;sum2=sum2+subs(f,xk);endt=((b-a)/(2*n))*(subs(f,a)+2*sum1+subs(f,b));s=((b-a)/(6*n))*(subs(f,a)+2*sum1+4*sum2+subs(f,b)); z=int(f,x,0,2*pi);T=vpa(t,2);S=vpa(s,2);Z=vpa(z,2)Dt=abs(T-Z);Ds=abs(S-Z);fun=[T;S;Dt;Ds];return;计算结果第八题Gauss-Chebyshev函数function result=chebyshev(f,n)x=zeros(1,n); A=zeros(1,n);m=n-1;for i=0:mx(i+1)=cos((2*i+1)*pi/(2*n));A(i+1)=pi/n;endresult=double(sum(A(1:n).*subs(f,x(1:n))));returnGauss-Legendre函数function result=lanendre(f,a,x)t=pi/4+pi/4*x;result=double(pi/4*sum(a.*subs(f,t)));return;计算过程syms xf1=x^2f2=sin(x)/xx2=[-0.57735,0.57735]a2=[1,1];x3=[-0.77460,0,0.77460]a3=[0.55556,0.88889,0.55556]x5=[-0.90618,-0.53847,0,0.53847,-0.90618]a5=[0.23693,0.47863,0.56889,0.47863,0.23693]fun1=[chebyshev(f1,2),chebyshev(f1,3),chebyshev(f1,5)]fun2=[lanendre(f2,a2,x2),lanendre(f2,a3,x3),lanendre(f2,a5,x5)] real1=double(int(f1/sqrt(1-x^2),-1,1))real2=double(int(f2,0,pi/2))结果第九题Euler法function [t,x]=Euler(f,x0,h)t=0:h:1;m=length(t);x=zeros(1,m);x(1)=x0;for i=1:m-1x(i+1)=x(i)+h*subs(f,{'t','x'},{t(i),x(i)});endend改进Euler法function [t,x]=ImproveEuler(f,x0,h)t=0:h:1;m=length(t);x=zeros(1,m);x(1)=x0;for i=1:m-1k1=subs(f,{'t','x'},{t(i),x(i)});x1=x(i)+h*k1;k2=subs(f,{'t','x'},{t(i+1),x1});x(i+1)=x(i)+h/2*(k1+k2);endendRunge_Kutta法function [t,x]=runge_kutta(f,x0,h)t=0:h:1;m=length(t);x=zeros(1,m);x(1)=x0;for i=1:m-1k1=h*subs(f,{'t','x'},{t(i),x(i)});k2=h*subs(f,{'t','x'},{t(i)+h/2,x(i)+k1/2}); k3=h*subs(f,{'t','x'},{t(i)+h/2,x(i)+k2/2}); k4=h*subs(f,{'t','x'},{t(i)+h,x(i)+k3});x(i+1)=x(i)+(k1+2*k2+2*k3+k4)/6;endend调用函数进行计算syms t xf=x*cos(t);x0=1;h=[0.1,0.01,0.001];for i=1:3[t,x1]=Euler(f,x0,h(i));plot(t,x1);hold on;end>>clcfor i=1:3[t,x2]=Improved_Euler(f,x0,h(i));plot(t,x2);hold on;end>>clcfor i=1:3[t,x3]=Runge_Kutta4(f,x0,h(i)); plot(t,x3);hold on;end。