Matlab線性代數(shù)實驗

Matlab線性代數(shù)實驗
8.1 實驗(Ⅰ)：用Matlab學線性代數(shù)
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;?? % 確定相關的參數(shù)
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的列向量指標。??????
? %有兩種不同的繪制方式，由參數(shù)l和p控制。
? %l=1,向量用紅色粗線條繪制；l=2,向量用藍色細線條繪制。p是控制箭頭的大小參數(shù)。?
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做線性代數(shù)題（矩陣的符號演算）
????????? ◆ 試題?? 設 ，求 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.高維線性離散動力系統(tǒng)???????
?參考程序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)???????? ％計算協(xié)方差矩陣，注意數(shù)據需按列向量存入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教程