3 线性代数实验
1.1 基础训练
- 用左除运算符求解线性方程组
{ 7 x 1 + x 2 + 2 x 3 − x 4 = 6 , 2 x 1 + 6 x 2 + x 3 + 3 x 4 = 5 , x 1 − 4 x 2 + 10 x 3 + 3 x 4 = 4 , 2 x 1 − x 2 − 3 x 3 + 11 x 4 = 3. \left\{\begin{matrix}7x_1+x_2+2x_3-x_4=6,\\2x_1+6x_2+x_3+3x_4=5,\\x_1-4x_2+10x_3+3x_4=4,\\2x_1-x_2-3x_3+11x_4=3.\\\end{matrix}\right. ⎩ ⎨ ⎧7x1+x2+2x3−x4=6,2x1+6x2+x3+3x4=5,x1−4x2+10x3+3x4=4,2x1−x2−3x3+11x4=3.
返回解向量(1个返回参数)。注:要求本问题不使用符号计算函数solve解线性方程组。
首先我们根据方程定义系数矩阵A和常数矩阵B 设未知数所组成的列向量为X
由AX=B X=A^-1*B
A=[7 1 2 -1
2 6 1 3
1 -4 10 3
2 -1 -3 11];
B=[6;5;4;3];
X=A\B
X =
0.7299
0.3828
0.3953
0.2826
- 基础知识:假定现在要在平面上旋转一个物体,需要绘制物体的旋转过程。由于绘制过程需要知道物体上各点的坐标,而旋转过程中物体上点的坐标是不断变化的。因此,绘制过程需要跟踪物体坐标的变化. 问题可以归结为已知旋转角度 θ \theta θ,旋转前P点的坐标为 ( x , y ) (x,y) (x,y),计算旋转后得点 P ′ P^\prime P′,设其坐标为 ( x ′ , y ′ ) (x^\prime,y^\prime) (x′,y′). 这里P点代表任意待旋转的点. 旋转变换就是一种线性变换,看如何构造一个矩阵完成这种旋转变换,完成坐标的计算. 问题进一步化为构造矩阵 M = ( m 11 m 12 m 21 m 22 ) M=\left(\begin{matrix}m_{11}&m_{12}\\m_{21}&m_{22}\\\end{matrix}\right) M=(m11m21m12m22),满足
( x ′ y ′ ) = ( m 11 m 12 m 21 m 22 ) ( x y ) . \left(\begin{matrix}x^\prime\\y^\prime\\\end{matrix}\right)=\left(\begin{matrix}m_{11}&m_{12}\\m_{21}&m_{22}\\\end{matrix}\right)\left(\begin{matrix}x\\y\\\end{matrix}\right). (x′y′)=(m11m21m12m22)(xy).
通过建模可得旋转变换矩阵 M = ( c o s θ − s i n θ s i n θ c o s θ ) . M=\left(\begin{matrix}cos{\theta}&-sin{\theta}\\sin{\theta}&cos{\theta}\\\end{matrix}\right). M=(cosθsinθ−sinθcosθ).
请设计一个正六边形绕原点旋转一圈的动画程序. 正六边形的位置、大小自行设置。
N=6;
t=0:2*pi/N:2*pi;
x=2*sin(t)
y=2*cos(t)
z=[x;y];
plot(x,y)
hold on
for th=0:pi/5:2*pi
M=[cos(th) -sin(th)
sin(th) cos(th)];
hh=M*z;
a=hh(1,:);
b=hh(2,:);
plot(a,b)
pause(0.5)
end
![在这里插入图片描述](https://img-blog.csdnimg.cn/008c578a86aa49e58bd825589b30e52c.png#pic_center)
1.2实验任务
一. 实验任务
已知某种生物分3个年龄组, 0 2 0~2 0 2岁为第一组, 3 − 4 3-4 3−4岁为第二组, 5 6 5~6 5 6岁为第三组. 每两年为一个时段. 设 x j ( k ) x_j^{(k)} xj(k)为第j个年龄组在第k个时段的数量. 记 X ( k ) = ( x 1 ( k ) , x 2 ( k ) , x 3 ( k ) ) T X^{(k)}=(x_1^{(k)},x_2^{(k)},x_3^{(k)})^T X(k)=(x1(k),x2(k),x3(k))T. 根据生物繁殖率以及存活率,得到该种生物数量变化规律的数学模型 X ( k + 1 ) = ( 0 3 2 0.4 0 0 0 0.15 0 ) X ( k ) ( k = 0 , 1 , 2 , . . . ) X^{(k+1)}=\left(\begin{matrix}0&3&2\\0.4&0&0\\0&0.15&0\\\end{matrix}\right)X^{(k)} (k=0,1,2,...) X(k+1)= 00.40300.15200 X(k)(k=0,1,2,...),已知 X ( 0 ) = ( 10 , 20 , 15 ) T X^{(0)}=(10,20,15)^T X(0)=(10,20,15)T.
请编程计算2、4、6、8、10年后该种生物各年龄组数量,当k足够大时,生物各年龄组数量变化规律是什么,说明导致该规律的原因。
(2)如果使用某种杀虫剂后,生物数量变化规律由式子 X ( k + 1 ) = ( 0 3 r 2 r 0.4 0 0 0 0.15 0 ) X ( k ) X^{(k+1)}=\left(\begin{matrix}0&3r&2r\\0.4&0&0\\0&0.15&0\\\end{matrix}\right)X^{(k)} X(k+1)= 00.403r00.152r00 X(k)
决定, r r r反应使用杀虫剂的剂量 ( 0.1 ≤ r ≤ 0.9 ) (0.1\le r\le0.9) (0.1≤r≤0.9),请编写程序确定 r r r使得三个年龄组昆虫数量最后基本趋于稳定不变,并找出初始稳定值。
二. 实验目的
熟悉Matlab矩阵计算函数。
三. 实验过程
(1)
数学模型已经为我们构建完毕
设变化矩阵为A 根据递推关系得出关系式:
X ( k + 1 ) = A k X ( 0 ) X^{(k+1)}=A^kX^{(0)} X(k+1)=AkX(0)
计算2-10年的数量,实现程序如下
A=[0 3 2
0.4 0 0
0 0.15 0];
x=[10;20;15];
for i=2:2:10
X=[];
X=A^i*x;
fprintf('第%d年的数量分别为\n',i)
disp(vpa(X,4))
end
输出结果:
第2年的数量分别为
18.0
36.0
0.6
第4年的数量分别为
32.4
43.68
1.08
第6年的数量分别为
51.98
53.28
1.944
第8年的数量分别为
78.36
65.49
3.119
第10年的数量分别为
113.7
81.08
4.702
通过改变循环的次数 观察输出结果我们发现 三个年龄段的种群数量均随着时间的推移而增长 但是相较奇偶年数又有所不同 我们绘制三个年龄阶段的1-10年的数量
可以直观的看出变化规律:
年龄组1 偶数年份相对于前一年的奇数年有所下降 整体呈现上升的趋势
年龄组2偶数年份相对于前一年的奇数年有所上升 整体呈现上升的趋势
年龄组3 偶数年份相对于前一年的奇数年有所下降 整体呈现上升的趋势
下面我们来探讨原因
通过[x,y]=eig(A) 求得A的特征值和特征向量
x =
-0.9429 -0.9323 0.1393
-0.3302 0.3580 -0.5525
-0.0434 -0.0516 0.8218
y =
1.1424 0 0
0 1.0415 0
0 0 -0.1009
为了方便起见我们可以将特征值改为正值 同时改变特征向量的正负
设λ1=1.1424,λ2=-1.0415,λ3=-0.1009 α1,α2,α3为对应于λ1,λ2,λ3的特征向量,则三个特征向量线性无关。
因此x可由α1,α2,α3线性表示,即x=[x1;x2;x3]= c1α1+ c2α2+c3α3
syms c1 c2 c3;
C=[c1 c2 c3];
[c1,c2,c3]=solve(C*x(:,1)==10,C*x(:,2)==20,C*x(:,3)==15);
disp(vpa([c1,c2,c3],4))
其中[c1 c2 c3]=[ -17.36, 15.17, 31.4]
由于λ2=-1.0415,λ3=-0.1009对应的特征向量都有正有负 所以奇偶年份有所不同
$$
X(k)\ =A^kX(0)\ =\mathrm{\ }Ak(c1α1+ c2α2+c3α3) \
=\mathrm{\ }c1Akα1+ c2Akα2 + c3Akα3 \
=\mathrm{\ }c1λ1kα1+ c2λ2kα2 + c3λ3kα3 \
$$
当k趋于无穷时
X
→
−
17.36
∗
1.1
4
k
(
−
0.9429
−
0.3302
−
0.0434
)
+
15.71
∗
1.04
k
(
−
0.9323
0.3580
−
0.0516
)
X\rightarrow\mathrm{-17.36*1.1}\mathrm{4}^k \begin{pmatrix} -0.9429\\-0.3302\\-0.0434 \end{pmatrix} +15.71*1.04k \begin{pmatrix} -0.9323\\0.3580\\-0.0516 \end{pmatrix}
X→−17.36∗1.14k
−0.9429−0.3302−0.0434
+15.71∗1.04k
−0.93230.3580−0.0516
所以随着时间的推移 种群数量不断增加
(2)
通过(1)的推导 我们得知只有在主特征值等于1的时候才会在时间趋于无穷的时候出现均衡
即问题转化为当r为何值时特征值为1
根据特征值的定义|λI-A|=0 带入特征值1 便可以求出r的值
syms r
A=[0 3*r 2*r
0.4 0 0
0 0.15 0];
r=solve(det((eye(3)-A)))
解得r =25/33
下面计算稳定值 程序仿照(1)所示
syms r
A=[0 3*r 2*r
0.4 0 0
0 0.15 0];
r=solve(det((eye(3)-A)));
A=[0 3*r 2*r
0.4 0 0
0 0.15 0];
x=[10;20;15];
for i=2:1:100
X=[];
X=A^i*x;
fprintf('第%d年的数量分别为\n',i)
disp(vpa(X,4))
end
由于我们采用的是四位有效数字 当k>46的时候基本上趋于稳定:
三个种群的数量分别为:39.56 15.83 2.373
实验自评与改进方向
自评:关于第一问完成情况较好,但是第二问并没有通过编程找到r,而是借助第一问推导出的结果,通过特征值为1的条件去解得r值,理想的解决办法是通过数量相等去解得r值。
改进方向:
syms r;
A=[0 3*r 2*r
0.4 0 0
0 0.15 0];
x=[10;20;15];
X=[];
for i=100:100:200
X=[X,A^i*x];
X=vpa(X,2);
end
disp(X)
r=solve(X(1,1)==X(1,2),X(2,1)==X(2,2),X(3,1)==X(3,2))
但是由于方程次数太高还有精确度的问题,无法得出有效的结果。但是这样改进的方向是没有问题的。想用fix(100*X)/100来确定保留的小数位数,但是由于数据类型的原因无法实现。
实验体会,收获及建议
题目主要涉及的是矩阵的相关知识,但是也用到了许多之前学到的内容。这就需要及时的回顾一下前边的内容。同时也结合了上学期的矩阵相关知识,用matlab去解决实际问题。将知识真正的应用到了实际生活中。
更好的了解了上课所讲的增长模型,回顾了特征值和特征向量去分析解决迭代问题的方法。更为重要的是了解了用matlab解决问题的思路和方法。
在编程中遇到的一些意想不到的问题也让我意识到,要在平时多加练习,不能眼高手低。