Transcript Document
§2.4 求解线性方程组
一、求解齐次线性方程组AX=0
方法:null(A) 或 null(A, ’r’)
注:
• 此命令给出齐次线性方程组AX=0的基
础解系,即解空间的一组基。
• 前者是数值解,后者是有理形式的基。
• 从而可写出通解:基础解系的线性组
合。
例、 x x 2 x x 0,
1
2
3
4
2 x1 x2 x3 x4 0,
2 x 2 x x 2 x 0.
2
3
4
1
解:>> A=[1 1 2 -1;2 1 1 -1;2 2 1 2];
>> null(A)
ans =
0.3621
-0.8148
0.3621
0.2716
解:>>A=[1 1 2 -1;2 1 1 -1;2 2 1 2];
>>null(A,’r’)
ans =
4/3
-3
4/3
1
通解为X=k(4/3 -3 4/3 1)^T,
k为任意常数.
注:其中的‘r’表示格式以有理数的形
式显示,即是将系数矩阵化成简化行简
阶梯形求得的基础解系。
若使用null(A)命令,得到的是解空
间的一组标准正交基,即X’*X=E.此
时的X是近似值,待入有 AX不等于
零,近似为零。
如上例 >>x= null(A)
x=
0.3621
-0.8148
0.3621
0.2716
>> x'*x
ans =
1.0000
>> A*x
ans =
1.0e-15 *
0
-0.3886
-0.6661
例、 x 1 2x 2 x 3 x 4 0,
3x 1 6x 2 x 3 3x 4 0,
5x 10x x 5x 0.
2
3
4
1
解:>> A=[1 2 1 -1;3 6 -1 -3;5 10 1 -5];
>> null(A,'r')
ans =
-2 1
1 0
0 0
0 1
B=null(A)
B=
-0.3229
-0.2891
-0.0000
-0.9012
0.8538
-0.4997
-0.0000
-0.1456
>> B'*B
ans =
1.0000 -0.0000
-0.0000 1.0000
二、求解非齐次线性方程组
Amn X b
1、当AX=b有唯一解时,X=A\b给出唯
一解;特别当A是方阵时也可用
X=inv(A)*b
例、
x1 2 x2 4 x4 3,
x x 4 x 9 x 22,
1 2
3
4
2
x
3
x
x
5
x
3
,
1
2
3
4
3 x1 2 x2 5 x3 x4 3
解:>>A =[1 2 0 -4;1 -1 -4 9;2 -3 1 5;3 -2 -5 1];
>> b=[-3;22;-3;3];
>> rank(A),rank([A,b])
ans =
>> X=A\b(或inv(A)*b)
4
X=
ans =
-1
4
3
-2
2
例、 >> A=[1 2 3;4 5 6;7 8 0;2 5 8];
>> b=[6 15 15 15]';
>> rank(A),rank([A,b])
ans =
3
>> A\b
ans=
ans =
3
1.0000
1.0000
1.0000
2、当AX=b有无穷多解时,用
pinv(A)*b给出一个特解.
例、 x1 2 x2 4 x3 3 x4 1,
3 x 5 x 6 x 4 x 2,
1
2
3
4
4 x1 5 x2 2 x3 3 x4 1,
3 x1 8 x2 24x3 19x4 5.
解:
>> b=[1;2;1;5];
>> A=[1 2 4 -3;3 5 6 -4;4 5 -2 3;3 8 24 -19];
>> rank(A),rank([A,b])
ans =
>> pinv(A)*b
2
ans =
ans=
0.1173
2
0.1732
0.1006
-0.0447
注:这是它的一个特解,可以检验:
A*(pinv(A)*b)=b,用null命令可求到对应
的齐次线性方程组的基础解系,从而可
求到通解.
8
7
>> null(A,'r')
6
5
ans =
1 , 2
1
0
8 -7
0
1
-6 5
1 0
0
1
即为基础解系。
若想用线代教材中的高斯消元法求解,
我们可以把增广矩阵化成简化行阶梯形,
从而得到通解。
如上例
>> B=[A,b];
>> rref(B)
ans =
1 0 -8 7 -1
0 1 6 -5 1
0 0 0 0 0
0 0 0 0 0
对应的齐次的基础解系为:
8
7
6
5
1 , 2
1
0
0
1
非齐次的特解为:
1
1
0
0
0
故通解为:
X k11 k 2 2 0
7 1
8
5 1
6
k1 k 2 ,
0
0
1
1 0
0
k1 , k 2为任意常数.
例、求解下面的非齐次线性方程组
2x 1 3x 2 x 3
x 1 2x 2 4x 3
x
2
x
8
x
3
3
2
1
4x 1 x 2 9x 3
4
5
13
6
解: >> A=[2 3 1;1 -2 4;3 8 -2;4 -1 9];
>> b=[4;-5;13;-6];
>> rref([A,b])
ans =
1
0
0
0
0
1
0
0
2
-1
0
0
-1
2
0
0
由此可得通解为:
2 1
X k 1 2 , k R .
1 0
例、 2 x1 x2 x3 x4 1,
4 x1 2 x2 2 x3 x4 2,
2 x x x x 1.
1 2 3 4
解: >> A=[2 1 -1 1;4 2 -2 1;2 1 -1 -1];
>> b=[1;2;1];
>> rank(A),rank([A,b])
ans = 2
ans = 2
>> pinv(A)*b
ans = 0.3333
0.1667
-0.1667
-0.0000
>> A*pinv(A)*b
ans = 1.0000
2.0000
1.0000
注:
pinv(A)*b
给出它的
一个特解.
>> x=A\b
Warning: Rank deficient, rank = 2, tol =
1.884111e-15.
x = 0.5000
0
0
注:A\b也
0.0000
给出它的
>> A*x
一个特解.
ans = 1.0000
2.0000
1.0000
>> null(A,'r')
ans =
-0.5000 0.5000
1.0000
0
0
1.0000
0
0
给出对应的齐次基
础解系。从而得通
解
3、当求解非齐次线性方程组无解时
可用pinv(A)*b给出一个最小二
乘解,即使误差向量AX-b的平
方和最小化的解。
例、AX=b
>> b=[1;2;1];
>> A=[1 2 3;4 5 6;2 4 6;];
解:
>> B=[A,b];
>> rank(A),rank(B)
ans = 2
ans= 3
>> pinv(A)*b
ans =
0.3222
0.1556
-0.0111
最小二乘解!
例、
>> A=[1 2 3;4 5 6;7 8 0;2 5 8];
>> b=[366 804 351 514]';
>> rank(A),rank([A,b])
ans =
3
ans=
4
注:可知AX=b无解。
>> pinv(A)*b
ans =
247.9818
-173.1091
114.9273
最小二乘解!
作业
1、至少用三种方法解方程组AX=b,
如矩阵除法、求逆矩阵、初等变
换等,其中
0
8
3 5
0
1 8 2 1
2
A
b
0 5 9
3
1
7 0 4 5
6
2、判断AX=b的解,并用适当命令求解,
其中
2 3
1
1 2 4
(1)A
3
8
2
4 1 9
4
5
b
13
6
x 1 2x 2 3x 3 1,
(2)2x 1 2x 2 5x 3 2,
3x 5x x 3.
2
3
1
4x 1 2x 2 x 3 2,
(3)3x 1 x 2 2x 3 10,
11x 3x 8.
2
1
x 1
x 1
(4)x 1
x
1
x 1
3x 2 5x 3 4x 4 1,
3x 2 5x 3 2x 4 x 5 1,
2x 2 x 3 x 4 x 5 3,
4x 2 x 3 x 4 x 5 3,
2x 2 x 3 x 4 x 5 1;
3、解方程组,并验证得出的解满
足原方程.
3 2 13
7
4
2 1 2
9
7
(1)
X
1
2 15 3 2
1 2 11 5
0
5 7 6 5 1
24 96
7 10 8 7 2
34 136
(2) 6 8 10 9 3 X 36 144
5 7 9 10 4
35 140
1 2 3 4 5
15 60
§2.5、特征值、特征向量的计算
一、如何理解矩阵A的特征值、
特征向量
?
在Matlab中可用命令eigshow来显示
X和AX在单位圆周上的运动效果
例、
(1)
5
A 4
0
0
3
4
>> A=[5/4 0;0 3/4];
>> eigshow(A)
若X为单位特征向量,则X与AX共线,
且AX的长就是相应的特征值。
(2)
(3)
1 0
A
0 1
0 1
A
1 0
二、二维和三维图形
1、命令:plot(x,y)
例、y=sin(x)
>> x=0:pi/100:2*pi;
>> y=sin(x);
>> plot(x,y)
注:可增加注记。
>> plot(x,y),...
xlabel('x'),ylabel('sin(x)'),
title('plot of the Sine function')
或
>> plot(x,y),...
xlabel('x'),...
ylabel('sin(x)'),...
title('plot of Sine function')
注:可在同一坐标系下画出不同函
数图形
例、y=sin(x),y=cos(x)
>> x=0:pi/100:2*pi;
y=sin(x);
plot(x,y)
>> hold on %添加到同一张图中
>> y2=cos(x);
>> plot(x,y2,'r:'),...
legend('sin','cos')
2、三维图形
命令:
[x,y]=meshgrid(初值:步长:终值);
z=f(x,y);
surf(x,y,z)
例、
z xe
x 2 y 2
>> [x,y]=meshgrid(-2:.2:2);
>> z=x.*exp(-x.^2-y.^2);
>> surf(x,y,z)
三、eig(A): 表示求矩阵A的特征值。
例、 >> A=[0 1;-1 0];
1. >> eig(A)
ans =
0.0000 + 1.0000i
0.0000 - 1.0000i
注:我们知道矩阵A的特征值是它的
特征多项式的根,所有我们也可求出
它的特征多项式。
2 . >> p=poly(A)
p=
1 0 1
注:输出的这些数是特征多项式中
各项的系数,按降幂排列。
>> x=-8:0.1:8;
>> y=polyval(p,x);
>> plot(x,y)
>> syms a
>> f=a^2+1;
>> solve(f==0,a)
ans =
i
-i
例、
1 1 0
A 4 3 0
1 0 2
>> A=[-1 1 0;-4 3 0;1 0 2];
>> poly(A)
ans =
1 -4
5 -2
>> x=-4:1:4;
>> y=polyval(p,x);
>> plot(x,y)
>> syms a;
>> f=a^3-4*a^2+5*a-2;
>> solve(f==0,a)
ans =
1
1
2
二、[变量名,矩阵名]=eig(A): 表示求矩阵
A的特征值和特征向量。
例、 >> [X,D]=eig(A)
注:结果中D的主对角线上 的元
素是 特征值,X的列是相应的特
征向量
X=
0.7071 + 0.0000i 0.7071 + 0.0000i
0.0000 + 0.7071i 0.0000 - 0.7071i
D=
0.0000 + 1.0000i 0.0000 + 0.0000i
0.0000 + 0.0000i 0.0000 - 1.0000i
注:在Matlab中,如没有特殊声明,
所以的运算都是进行数值计算,结果
会自动以近似值的形式给出。若需要
得到精确的表达式,可声明为
symbolic对象,(符号对象,需要有
symbolic math toolbox工具箱的支持)。
命令为sym
例、求矩阵A的特征值与特征向量
1 2 2 2
2 1 2 2
A
2 2 1 2
2 2 2 1
>> A=[1 -2 -2 -2;-2 1 -2 -2;-2 -2 1 -2;-2 -2 -2 1];
>> [V,D]=eig(A)
V=
0.5000 -0.2113 0.7887 0.2887
0.5000 0.7887 -0.2113 0.2887
0.5000 -0.5774 -0.5774 0.2887
0.5000
0
0
-0.8660
D=
-5
0
0
0
0
3
0
0
0
0
3
0
0
0
0
3
>> A=sym([1 -2 -2 -2;-2 1 -2 -2;-2 -2 1 -2;2 -2 -2 1]);
>> [V,D]=eig(A)
V=
[ 1, -1, -1, -1]
[ 1, 1, 0, 0]
[ 1, 0, 1, 0]
[ 1, 0, 0, 1]
D=
-5
0
0
0
0
3
0
0
0
0
3
0
0
0
0
3
注:对于重特征值的,可由给出的特征
向量判断该矩阵是否可对角化。
例、 >> A=[5 6 -3;-1 0 1;1 2 1];
>> [V,D]=eig(A)
D=
2.0000
0
0
0 2.0000
0
0
0 2.0000
V=
0.9045 -0.9045 -0.8393
-0.3015 0.3015 0.5116
0.3015 -0.3015 0.1840
注:显然第一列、第二列相同,线性相
关,所以它没有三个线性无关特征
向量,故矩阵A不可对角化。
作业
1、考虑矩阵
5 3
5 3
,B
A
3 5
3 5
(1)使用Matlab 计算A,B的特征
值,并判断是否是相同 类型特征
值。
(2)使用命令分别形成它们的特
征多项式,并画出它们的图形,
x=-8:0.1:8.
(3)从图形看,关于B的特征值你能
得到什么结论?
2、求出矩阵A和B的特征值和特征
向量,并判断是否可对角化.其中
1 2 3 4
5 6 7 8
A
9 10 11 12
13 14 15 16
5 7 6 5
7 10 8 7
B
6 8 10 9
5 7 9 10
3、设矩阵A
0 2 1
A 2 5 2
4 8 3
(1)求矩阵A的特征值;
(2)问A能否相似于对角阵,
1
若相似,求可逆矩阵P,使得 P AP
为对角阵;
(3)求 A 100
4、按如下方式构造一个对称矩阵:
A=round(5*rand(6));
A=A+A'
并完成如下操作:
(1)用合适的命令求A的特征
值、特征值的和、特征值
的乘积、A的行列式及A的
迹,并比较这些结果,哪
些是相等的;
(2)用命令[X,D]求A的特征
向量,并计算X^(-1)AX并比
较此结果和D的关系;
计算A^(-1)和XD^(-1)X^(-1)
比较这两个结果。
5、令
A=ones(10)+eye(10),
求A的特征值、A-E的秩及A的
迹,并从理论上解释为什么1
是它的九重特征值,而另一
个特征值是11,写出你的理
由。