没有合适的资源?快使用搜索试试~ 我知道了~
MATLAB的线性代数计算.doc
1.该资源内容由用户上传,如若侵权请联系客服进行举报
2.虚拟产品一经售出概不退款(资源遇到问题,请及时私信上传者)
2.虚拟产品一经售出概不退款(资源遇到问题,请及时私信上传者)
版权申诉
5星 · 超过95%的资源 1 下载量 118 浏览量
2022-07-06
03:08:58
上传
评论
收藏 267KB DOC 举报
温馨提示
试读
40页
MATLAB的线性代数计算
资源推荐
资源详情
资源评论
第二章 MATLAB的線性代數計算
本章先介紹用MATLAB解線性方程組的方法, 應用此方法, 說明線性代數中有關線
性組合,線性相依,線性獨立的概念與判斷. 另外也討論並估計線性方程組近似值解之正
確度. 最後說明LU-Factorization與 Choleski-Decomposition 及其應用.
(一) 解線性方程組 Ax= b
(1) 矩陣 A 是一個 upper triangular matrix, 主對角線上的元素不為零
A=[4 -1 3;0 2 5;0 0 8];
b=[1 0 2]';
n=3; X=zeros(n,1); % 給初始值 X=[0 0 0]'
for j=n:-1:1 % 利用loop來執行Backward Subsitution
X(j)=(b(j)-A(j,: )*X)/A(j,j);
end, X
X =
-0.0938
-0.6250
0.2500
(2) 矩陣 A 是一般矩陣, 而且是 nonsingular matrix 則利用 Gaussian
Elimination
Algorithm採用 maximum column pivot 將其化為 triangular matrix,
以求解
A=[2 2 -3;3 1 -2;6 8 0];
b=[2 2 30]';
w=[A b]; % 建一擴增矩陣(augmented matrix)
p=[1 2 3]'; % 初始的 pivot vector
pivot=w(3,1); % 選定第一個 pivot element
p=[3 2 1]; % 更新後的 pivot vector
w(1,:)=(-w(1,1)/pivot)*w(3,:)+w(1,:) % 使(1,1)entry為0
w =
0 -0.6667 -3.0000 -8.0000
3.0000 1.0000 -2.0000 2.0000
6.0000 8.0000 0 30.0000
w(2,:)=(-w(2,1)/pivot)*w(3,:)+w(2,:) % 使(2,1)entry為0
w =
0 -0.6667 -3.0000 -8.0000
0 -3.0000 -2.0000 -13.0000
6.0000 8.0000 0 30.0000
第二章 第1頁
pivot=w(2,2); % 選定第二個 pivot element
w(1,:)=(-w(1,2)/pivot)*w(2,:)+w(1,:) % 使(1,2)entry為0
p=[3 2 1]'; % 更新後的 pivot vector
w =
0 0 -2.5556 -5.1111
0 -3.0000 -2.0000 -13.0000
6.0000 8.0000 0 30.0000
C=w(:,1:3); B=w(:,4); % separate coefficient matrix and
% right-hand side
x=zeros(3,1); % 給x的初始值x=[0 0 0]'
for j=3:-1:1 % 利用loop來執行類似Backward Subsitution
% 以求解
x(j)=(B(p(j))-C(p(j),:)*x)/C(p(j),j);
end, x
x =
1
3
2
xx=A\b % 用MATLAB內部的除法"\"來解
abs(x-xx) % 估計用上述Gaussian Elimination和
% A\b所求得的解之間的誤差
xx =
1.0000
3.0000
2.0000
ans =
1.0e-015 *
0
0
0.4441
(3) Row Reduced Echelon Form (rref)
A=[2 -1 0 2 3;-1 -4 3 1 2;5 1 -3 5 7];
b=[5 1 10]';
aug=[A b]; %建一個擴增矩陣(augmented matrix)
c=rref(aug)
c =
1.0000 0 0 0 0 2.0000
0 1.0000 0 -2.0000 -3.0000 -1.0000
0 0 1.0000 -2.3333 -3.3333 -0.3333
【說明】
第二章 第2頁
由 矩陣 c 的 rank=3 我們得知此方程組有解, 因為 rank(A) = rank([A
b])。
而且自由度(degree of freedom)為2, 其解形式為:
x1 = 2.0000 ;
x2 = - 1.0000 + 3.0000 t +2.0000 s ;
x3 = - 0.3333 + 3.3333 t +2.3333 s ;
x4 = s ;
x5 = t ; s, t 為任意數。
(4) 利用 A\b 來解 Nonsingular Linear Systems
例 1.
A=[1 -4 2 1;-4 0 3 -1;2 3 0 5;1 -1 5 -5];
b=[8 -1 17 -4]';
x=A\b
x =
1
0
2
3
【說明】
(1) 我們利用A\b 的方法解出 x=[1 0 2 3]' 為此方程式的唯一解。
(2) Matlab中的A\b即是解Ax=b by Pivoting Gaussian Elimination。
註:若
A
為
Singular matrix
,則有下列兩種情形出現:
例 2.
A=[1 1 1;2 2 2;1 2 1];
b=[1 0 1]';
x=A\b
Warning: Matrix is singular to working precision.
x =
Inf
Inf
Inf
【說明】
因為矩陣 A 為 Singular Matrix ,我們無法找到唯一的 x ,
使得 Ax=b,故此方程式無解。
例 3.
A=[1 2 3;4 5 6;7 8 9];
b=[1 0 0]';
xcomp=A\b
第二章 第3頁
Warning: Matrix is close to singular or badly scaled.
Results may be inaccurate. RCOND = 2.055969e-018.
xcomp =
1.0e+015 *
-4.5036
9.0072
-4.5036
【說明】
當我們用A\ b來解此題時,Matlab警告其RCOND (reciprocal of
condition number)趨近於 eps.,表示矩陣A近似於一個
Singular matrix
,故我們可得知其xcomp無意義。
(5) Cramer's Rule Algorithm
A=[5 2 1;-3 6 4;0 2 1];
b=[1 -1 1]';
x=zeros(3,1);
for j=1:3
Aj=A;
Aj(:,j)=b; % A的第j行由b取代
x(j)=det(Aj)/det(A);
end, x
x =
0
2.5000
-4.0000
【說明】
雖然我們得到一組 x=[0 2.5 -4]' 使得Ax=b, 但是因為計算det(A)費時,
所以較少使用。
(6) 有關的一些應用問題
例1. 求一 通過 (1,6) ,(-1,8),(2,11) 三點的二次方程式 ,並畫出此方程
式。
【解】
設此方程式為 p(t)= a t
2
+ b t + c
則 p(1)=6,p(-1)=8,and p(2)=11;
i.e. 解下列線性方程式的解
a + b + c = 6
a - b + c = 8
4 a + 2 b + c = 11
Let
第二章 第4頁
A=[1 1 1;1 -1 1;4 2 1];
b=[6 8 11]';
x=A\b
x =
2
-1
5
故 p(t)= 2 t
2
- t +5
為通過 (1,6), (-1,8), (2,11)三點的二次方
程式。
畫 p(t) 的圖形:
v=[2 -1 5]; %建一個代表p(t)係數的向量
t=-4:0.05:4;
y=polyval(v,t);
plot(t,y)
-4 -3 -2 -1 0 1 2 3 4
0
5
10
15
20
25
30
35
40
45
【說明】
通過 n+1個點之n次方程式 P(t)=a
n
t
n
+a
n-1
t
n-1
+…+a
1
t + a
0
則 Ay=c;
其中
A=[ t
1
n
t
1
n-1
… t
1
1
. . .
. . .
. . .
t
n+1
n
t
n+1
n-1
… t
n+1
1 ]
y= [a
n
a
n-1
… a
1
a
0
]';
c= [ c
n
c
n-1
… c
1
c
0
]';
第二章 第5頁
剩余39页未读,继续阅读
资源评论
- yhiuhn2023-09-01超赞的资源,感谢资源主分享,大家一起进步!
老帽爬新坡
- 粉丝: 83
- 资源: 2万+
下载权益
C知道特权
VIP文章
课程特权
开通VIP
上传资源 快速赚钱
- 我的内容管理 展开
- 我的资源 快来上传第一个资源
- 我的收益 登录查看自己的收益
- 我的积分 登录查看自己的积分
- 我的C币 登录后查看C币余额
- 我的收藏
- 我的下载
- 下载帮助
安全验证
文档复制为VIP权益,开通VIP直接复制
信息提交成功