Matlab线性代数实验

Matlab线性代数实验
8.1 实验(Ⅰ)：用Matlab学线性代数
8.1.1实验与观察：向量组的线性关系和解线性方程组
1.  用线性组合的方式产生向量组
【     clear  n=3; m=2; a=-10; b=10; 
  rand('seed',32), A = unifrnd(a,b,[n,m])       】
【    x = unifrnd(-1,1,[1,m]),   A(:,3)=x(1)*A(:,1)+x(2)*A(:,2)        】
2.Gauss消元法和向量组的线性关系的判定
【   B=rref(A)   】
【     r=2;m=2;s=5;
   X=[-B(1:2,r+1:m+s) ;eye(m+s-r)]   %基础解系      】
【     r1=rank(A(:,1:2)), r2=rank(A(:,1:3))      】
【      B=rref(A(:,1:3))   】    
【     B=rref(A);C=B(1:2,:)  】                          
3 . 观察程序    
zxy8_1.m    （主程序，产生向量并画图）
【   clf,n=3;m=2;s=5;   % 确定相关的参数
a=-10;b=10;
rand('seed',32),A = unifrnd(a,b,[n,m]),  %产生m个n维向量（生成向量）
r=[1:m]; l=1;p=1.2;
zxy8_1plot(A,r,l,p)          %将向量画出
for i=m+1:m+s                %产生组合向量
    x = unifrnd(-1,1,[1,m]);
    A(:,i)=zeros(n,1);
   for j=1:m
    A(:,i)=A(:,i)+x(j)*A(:,j);
   end
end
hold on,r=[m:m+s];l=2;p=0.5;
zxy8_1plot(A,r,l,p)               %用另一种方式画组合向量   】
其中调用了画图子程序zxy8_1plot
【  function out=zxy8_1plot(A,r,l,p)
  %A为若干3维向量拼成的矩阵, 绘制这些向量；r为向量，指明要绘制A的列向量指标。      
  %有两种不同的绘制方式，由参数l和p控制。
  %l=1,向量用红色粗线条绘制；l=2,向量用蓝色细线条绘制。p是控制箭头的大小参数。 
a=-10;b=10;k=50;
x=linspace(a,b,k);y=x;
[X,Y]=meshgrid(x,y);axis([a b a b a b]);
xlabel('x1','fontsize',14),ylabel('x2','fontsize',14),zlabel('x3','fontsize',14)
if l==1             %第一种绘图方式
  for i=r
    h1=plot3([0 A(1,i)],[0 A(2,i)],[0 A(3,i)],'rs','linewidth',3);        %画线段
    u(1)=A(1,i);v(1)=A(2,i);w(1)=A(3,i);u(2)=eps; v(2)=eps; w(2)=eps;
    h=quiver3([0 A(1,i)],[0 A(2,i)],[0,A(3,i)],u,v,w,p);                    %画箭头
    set(h,'linewidth',1,'color','red'),   axis([a b a b a b]),hold on,
    text(A(1,i),A(2,i),A(3,i),['\leftarrowA',int2str(i)],'fontsize',14);  %文字标注
  end
  Z=zxyplane(X,Y,A(:,1)',A(:,2)',0,0,0);              %计算平面的z坐标
  mesh(X,Y,Z), axis([a b a b a b]),hidden off,       %画平面
elseif l==2        %第二种绘图方式，结构同上
   for i=r
     h1=plot3([0 A(1,i)],[0 A(2,i)],[0 A(3,i)],'o--');
     set(h1,'linewidth',2,'color','blue')
     u(1)=A(1,i);v(1)=A(2,i);w(1)=A(3,i);u(2)=eps;v(2)=eps; w(2)=eps;
     h=quiver3([0 A(1,i)],[0 A(2,i)],[0,A(3,i)],u,v,w,p);
     set(h,'linewidth',1,'color','blue'),     axis([a b a b a b]),hold on,  
     text(A(1,i),A(2,i),A(3,i),['\leftarrowA',int2str(i)],'fontsize',14);
   end
 end             】
zxyplane.m      （计算平面的z坐标）
【   function z=zxyplane(x,y,a,b,x0,y0,z0)
 %向量a,b叉积构成平面的法矢t,平面过（x0,y0,z0）点
t(1)=det([a(2) a(3);b(2) b(3)]);
t(2)=-det([a(1) a(3);b(1) b(3)]);
t(3)=det([a(1) a(2);b(1) b(2)]);
if t(3)~=0
   z=z0-(t(1)/t(3))*(x-x0)-(t(2)/t(3))*(y-y0);
   k=find(z<=-8&z>=8);z(k)=NaN;
elseif t(3)==0
   disp('t(3)=0,this programme can not finish the work you want to do ')
end                            】
8.1.2应用、思考与练习
1. 观察极大线性无关组的意义
【     B=rref(A), C=B(1:2,:),     A(:,1:2)*C,  A               】                                             
2. 平面四连杆机构的设计      
     
3. 用Matlab做线性代数题（矩阵的符号演算）
          ◆ 试题   设 ，求 a的值使得  1. a4不能由a1、a2、a3线性表示。  2. a4可以由a1、a2、a3线性表示，并写出表达式。
      【   syms a
a1=[1;4;0;2]; a2=[2;7;1;3];a3=[0;1;-1;1];a4=[3;10;a;4];
A=[a1,a2,a3,a4]
for i=2:4             %行初等变换
  A(i,:)=A(i,:)-A(1,:)*A(i,1);
end
  A(2,:)=A(2,:)/A(2,2);
for i=3:4
  A(i,:)=A(i,:)-A(2,:)*A(i,2);
end
A   】
 
8.2 实验(Ⅱ)：矩阵的相似化简
8.2.1 实验与观察：矩阵的特征―相似标准形的作用
 
1. 逼近直线的迭代点列
2. 估计直线－特征值、特征向量
【  [P,D] = eig(A)     】
【  x=linspace(a,b,30);  [pc,lamda]=eig(A),pc=-pc;   z1=pc(2,1)/pc(1,1)*x;
     plot(x,z1,'linewidth',3)              】
3.特征值和特征向量决定迭代性质?
 【  x0=[1 1]'; A=[1/5,99/100;1,0]; 
       [P,D]=eig(A);
       y0=inv(P)*x0;y=y0;
for i=1:50
 y=[D^i*y0,y];
end
x=P*y;
plot(y(1,:),y(2,:),'o',x(1,:),x(2,:),'*'),legend('Y','X=PY')       】
4. 观测程序说明
 zxy8_2.m的源代码如下：
【    clear，clf
a=-20*100;b=-a;c=a;d=b;p=0.1;  %设定画图范围
n=100;
A=[1/5,99/100;1,0];            %设定矩阵
axis([a b c d]),grid on,hold on
button=1
while button==1
[xi,yi,button]=ginput(1);           %用鼠标选初始点
plot(xi,yi,'o'),hold on,           
X0=[xi;yi];X=X0;
  for i=1:n
     X=[A*X,X0];                    %用这种方式作迭代，并画图
     h=plot(X(1,1),X(2,1),'R.',X(1,1:2),X(2,1:2),'r:');hold on
     quiver([X(1,2),1]',[X(2,2),1]',[X(1,1)-X(1,2),0]',[X(2,1)-X(2,2),0]',p)
    set(h,'MarkerSize',6), grid,
end
end
pause
x=linspace(a,b,30);              %画最大特征值所对应的特征向量所决定的直线[pc,lamda]=eig(A),pc=-pc;
z1=pc(2,1)/pc(1,1)*x;
z2=pc(2,2)/pc(1,2)*x;
h=plot(x,z1),set(h,'linewidth',2)                                】
8.2.2 应用、思考与练习
1. 植物基因的分布、杂交育种问题
3.高维线性离散动力系统       
 参考程序zxy8_3.m . (计算图8.6和8.7)
【   clear,clf
a=-2;b=-a;c=a;d=b;p=0.0001*abs(a); %设定画图范围
n=30;
A=[2,0;0,1/2]             %设定矩阵
axis([a b c d]),hold on
button=1                 
while button==1
[xi,yi,button]=ginput(1);   %用鼠标选初始点
plot(xi,yi,'o'),hold on,
X0=[xi;yi];X=X0;
  for i=1:n
     X=[A*X,X0];               %用这种方式作迭代，并画图
     h=plot(X(1,1),X(2,1),'R.',X(1,1:2),X(2,1:2),'r:');hold on
     quiver([X(1,2),1]',[X(2,2),1]',[X(1,1)-X(1,2),0]',[X(2,1)-X(2,2),0]',p)
   set(h,'MarkerSize',6), grid,
end
end             】
3. 主成分分析和线性变换
 下面是相应的程序：
【  x1=[-0.931 -3.931 20.069 -6.931 -4.931  20.931 -3.931 -12.931 28.069 0.069 -16.931 16.069 -11.931 -8.931 -7.931  0.069 0.069 -16.931 -8.931 -3.931 14.069 11.069 -15.931 21.069 -10.934 7.069 15.069 8.069 16.069];
     x2=[-0.379 -0.379 12.621 -4.379 -3.379 -11.379 -1.379  -3.379  6.621 2.621  -9.379  3.621  -7.379 -2.379 -3.379  3.621 0.379  -6.379 -3.379 -0.379  2.621  7.621  -7.379  4.621  -3.379 4.621  9.621 4.621  5.621]; 
       X=[x1',x2'];C=cov(X)         ％计算协方差矩阵，注意数据需按列向量存入X中。
       [P,latent,explained] = pcacov(C)       ％主成分分析，P即为所需的变换矩阵。    】
【     clf,  a=-20;b=30;c=-20;d=30;
  t=linspace(a,b,20);  z1= P(2,1)/P(1,1)*t;  z2=P(2,2)/P(1,2)*t;hold on,
  plot(x1,x2,'o',[a b],[0 0],':',[0 0],[c d],'b:',t,z1,t,z2),
   xlabel('x1'),ylabel('x2') ,   axis([a b c d]),axis('equal')                    】
  zxy8_5.m
【    clear, clf,a=-20;b=30,n=8，x=linspace(a,b,n);y=x;
[x1,x2]=meshgrid(x);[x4,x3]=meshgrid(x);      %生成网格线
A=[0 ,2.2;1 0.2];   a=pi/6;
for k=1:n                             %作网格线的线性变换
 for l=1:n
    z=A*[x1(k,l),x2(k,l)]';       %垂直线的的变换
    y1(k,l)=z(1);y2(k,l)=z(2);  %将变换后的坐标适当排列以便作图
    z=A*[x3(k,l),x4(k,l)]';       %水平线的的变换
    y3(k,l)=z(1); y4(k,l)=z(2);
 end
end
t=linspace(a,b,30);z1=5*cos(t);z2=15+5*sin(t);      %产生一个椭圆
z3=A*[z1;z2];                 %对椭圆的变换
figure(1)                        %在第一个图形窗口画x1-x2平面的内容plot(x1,x2,x3,x4,z1,z2),axis('equal')
figure(2)                        %在第二个图形窗口画y1-y2平面的内容

200多个MATLAB经典教程和MATLAB论文请查看：matlab教程