实验5 线性方程组的解法I.实验内容及要点一.注意以下函数的用法1.break:可以导致包含该命令的while、for循环终止。
也可以在if-end、switch-case、try-catch中导致中断2.continue:跳过位于其后的循环中的其他命令,执行循环的下一次迭代。
3.return:结束该命令所在函数的执行,把控制交给主调函数或命令窗口。
4.error(’message’):显示出错信息message,终止程序。
5.warning(‘message’):显示警告信息message,程序继续运行。
二.注意直接法中误差的判断(条件数的应用)和迭代法中收敛性的判断(见函数isshoulian)三.LU消元法的程序function x=xiaoyuan(a,b)[m n]=size(a); %可以讨论m n 的大小关系[l u]=lu(a);s=inv(l)*[a,b];x=ones(m,1);for i=m:-1:1h=s(i,m+1);for j=m:-1:1 %if j~=i% h=s(i,m+1)-h=h-x(j)*s(i,j);% (s(i,1:m) *xend % -s(i,i))endx(i)=h/s(i,i);end四.function y=isshoulian(a)s=size(a);if s(1)~=s(2),error('请输入方阵'),endn=s(1);for i=1:nm=0;for j=1:nif j~=im=m+abs(a(i,j));endendif abs(a(i,i))<my=0; %迭代不收敛returnendendy=1;%迭代收敛五.雅克比迭代function x=ydiedai(a,b,n)if isshoulian(a)==0warning('迭代不收敛')returnendl=length(b);t=b;b=zeros(l,1); %确保参与运算的是列向量for i=1:lb(i)=t(i);end[m m]=size(a);d=diag(diag(a));l=-tril(a,-1); %或l=-tril(a)+d;u=-triu(a,1); %或u=-triu(a)+d;b1=inv(d)*(l+u);f1=inv(d)*b;x=zeros(m,1);for i=1:n %常用while循环来设计带误差的终止条件x=b1*x+f1;end六.高斯——赛德尔迭代function x=gdiedai(A,b,x0,tol)l1=length(x0);h=zeros(l1,1); % x0=x0(:);for i=1:l1h(i)=x0(i);endl2=length(b);t=b;b=zeros(l2,1); %b=b(:);for i=1:l2b(i)=t(i);end[m n]=size(A); %.....d=diag(diag(A));l=-tril(A,-1); %或l=-tril(A)+d;u=-triu(A,1); %或u=-triu(A)+d;b2=inv(d-l)*u;f2=inv(d-l)*b;x1=h; %即x0x=b2*x1+f2;i=1;while abs(x-x1)>tol %常用范数来做判断x1=x;i=i+1;x=b2*x1+f2;endx;iII.课后作业一.略二.解:1.a=[3.0212.714 6.913;1.031 -4.273 1.121;5.084 -5.832 9.155];b=[12.648;-2.121;8.407];h=det(a) %判断a是否几近奇异,进而判断是否可能病态x=xiaoyuan(a,b)a(2,2)=-4.275;h1=det(a) %判断a是否几近奇异,进而判断是否可能病态x=xiaoyuan(a,b)2.解:3.a=hilb(10);x=ones(10,1);b=a*x;b=b.*(1+0.01);x1=xiaoyuan(a,b);x2=gdiedai(a,b,x,0.001);x3=ydiedai(a,b,3);[b x1 x2 x3]c=cond(a)p=max(abs(eig(a)))三.解:n=1000;b=[1:n]';a1=sparse(1:n,1:n,4,n,n);a2=sparse(2:n,1:n-1,1,n,n);a=a1+a2+a2';% 输出用稀疏矩阵求解的时间t1tic; x=a\b; t1=toc% 与满阵做比较aa=full(a);% 输出用满阵求解的时间tic; xx=aa\b; t2=toc% 为检验x与xx是否相同分别输出其分量之和y=sum(x)yy=sum(xx)四.解:1.% 本题可以转化为求解方程组% 例如a题,可转化为t1*sin(20*pi/180)=5;w+t2*sin(10*pi/180)=5;t1*cos(20*pi/180)-t1*cos(10*pi/180)=0 % 以下求解aa=[sin(20*pi/180) 0 0;0 sin(10*pi/180) 1;cos(20*pi/180) -cos(10*pi/180) 0];b=[5;5;0]x1=xiaoyuan(a,b)x2=gdiedai(a,b,[0 0 0],0.001)%看结果如何,若不行,就看范数x3=ydiedai(a,b,3) %是否小于1,即迭代是否收敛2.% 本题可以转化为求解方程组% 例如b题,可转化为l1*sin(20*pi/180)+l2*sin(10*pi/180)=d;l1*cos(20*pi/180)+l2*cos(10*pi/180)=h % 以下求解b题a=[sin(20*pi/180),sin(10*pi/180);cos(20*pi/180),cos(10*pi/180)];d=2;h=8;b=[d;h];L=xiaoyuan(a,b)3.% 本题可以转化为求解方程组% 例如c题,可转化为% t1*sin(40*pi/180)=5;% t2*sin(30*pi/180)+w1=5;% t1*cos(40*pi/180)-t2*cos(30*pi/180)=0;% t2*sin(30*pi/180)-t3*sin(20*pi/180)-w2=0;% t2*cos(30*pi/180)-t3*cos(20*pi/180)=0;% 以下求解c题a=[sin(40*pi/180) 0 0 0 0;0 sin(30*pi/180) 0 1 0;cos(40*pi/180) -cos(30*pi/180) 0 0 0; 0sin(30*pi/180) -sin(20*pi/180) 0 -1;0 cos(30*pi/180) -cos(20*pi/180) 0 0];b=[5;5;0;0;0];xiaoyuan(a,b)五.略六.略七.略八.略九.略十.略。