Matlab详细教程

HUASHI3483

4.8.1 For 回圈  

------------------------------------------------------------------------------
--  

for 回圈是用在须重复执行且执行次数有一定的算式，它的结构如下：   

for index = array   

command A   

end   

如果我们要计算一缆车离铁塔的速度 (v)，它的速度计算方式与且铁塔的距离 (d)有关，
假设以 10 公尺为判断值，则速度计算分为二个算式：   

   

假设有一个阵列 d 为缆车到铁塔的距离，则以下的for 回圈可计算速对应的速度   

＞> for k = 1:length(d)   

if d(k) <= 10   

velocity = 0.425 + 0.00175*d(k)^2;   

else   

velocity = 0.625 + 0.12*d - 0.00025*d(k)^2;   

end   

fprintf('d= %f velocity= %f\n',d(k),velocity)   

end  


另外几个例子   

＞> for n=1:10   

x(n)=sin(n*pi/10);   

end   

＞> disp(x)  


＞> for n=1:5   

for m=5:-1:1   

A(n,m)=n^2+m^2;   

end   

disp(n)   

end   

＞> disp(A)  


但是如果可以用阵列或是矩阵运算来取代以for 回圈计算，就应采用前者因为计算速度快
多了。上述的例子可改为   

＞> n=1:10;   

＞> x=sin(n*pi/10);  


使用 for 回圈的规则如下：   

上述的 for 回圈中的指标 (index) 须为是一变数。   
如果 array 代表阵列是空无一物，则回圈不会被执行，例如 k=1:0。   
如果 array 代表阵列是一纯量，则回圈会被执行一次，例如 k=1:1。   
如果 array 代表阵列是一向量，则回圈会被依序的执行，例如 k=1:b, b=[1 3 5]。   
如果 array 代表阵列是一矩阵，则回圈会被逐行依序的执行，例如 k=1:B, B=[1 2; 3  
4]。   
for 完整的语法为： for k = first:increment:last，其中的 first, increment, last
分别为初始值，增量，终止值。而回圈被执行的次数由以下的算式决定：   
floor((last-first)/increment)+1   

如果计算得到的值为负，则回圈不被执行。  



------------------------------------------------------------------------------
--  
   

HUASHI3483

5.1.1 基本矩阵运算元  

------------------------------------------------------------------------------
--  

我们在第二章已说明过 MATLAB 的运算是以阵列(array)及矩阵 (matrix) 方式在做运算
，而这二者在MATLAB的基本运算性质不同，阵列强调元素对元素的运算，而矩阵则采用线
性代数的运算方式。我们就来说明矩阵运算的特点。   


以下将阵列及矩阵的运算符号及其意义列出   

阵列运算符号 矩阵运算符号  功能   
+ + 加   
- - 减   
.* * 乘   
./ / 左除   
.\ \ 右除   
.^ ^ 次方   
.' ' 转置   



利用这些运算符号即可进行以下的矩阵运算。   

＞> A=[2 5 1; 7 3 8; 4 5 21; 16 13 0];   

＞> A' % A的转置矩阵   

A =   

2 7 4 16   

5 3 5 13   

1 8 21 0  


＞> A=[4 -1 3]; B=[-2 5 2];   

＞> dot_prod = sum(A.*B) % 二个阵列做内积   

dot_prod =   

-7   

＞> c=dot(A,B) % 以dot函数也可做内积运算   

c =   

-7  


＞> A=[4; -1; 3];   

＞> dot_prod = sum(A'.*B); % 如果A是行阵列则先做转置，再做内积   


＞> F=[2 5 -1]; G=[0 1 -3];   

＞> out_prod=F'*G; % 二矩阵做外积  


＞> A=[2,5,1; 0,3,-1];   

＞> B=[1,0,2; -1,4,-2; 5,2,1];   

＞> C=A*B % 矩阵相乘，注意二个矩阵的大小须相容   

C =   

2 22 -5   

-8 10 -7  


＞> A=[2 1; 4 3];   

＞> A^2 % 矩阵次方   

ans =   

4 1   

16 9  



------------------------------------------------------------------------------

HUASHI3483

5.1.2 矩阵多项式  

------------------------------------------------------------------------------
--  

函数polyvalm是以矩阵方式做多项式函数计算，有别于polyval是以阵列方式计算函数值
。它的语法为 polyvalm(a,X)，其中X为一矩阵而a则是一多项式。以下的例子可说明其用
法。   

＞> X=[1 1 1; 2 2 2; 3 3 3];   

＞> a=[1 1 1]; % 注意a=X*X+X+I   

＞> f=polyvalm(a,X)   

f =   

8 7 7   

14 15 14   

21 21 22   



------------------------------------------------------------------------------
--  
   

HUASHI3483

5.3.1 反矩阵、矩阵秩与行列式  

------------------------------------------------------------------------------
--  

一个正方矩阵A的反矩阵的定义是，所以此二矩阵相乘不论是或，结果皆为单位矩阵。但
是一矩阵如果是奇异(singular) 或是条件不足 (ill-conditioned)，其反矩阵并不存在
。条件不足的矩阵与一组联立方程组其中的方程式并不独立有关，而一矩阵的秩(rank)  
即是代表矩阵中独立方程式个数。如果一矩阵的秩数和其矩阵的列数相等，则此矩阵为非
奇异且其反矩阵存在。   


MATLAB的反矩阵函数和秩函数语法分别为inv(A), rank(A)，：例如：   

＞> A=[2 1; 4 3];   

＞> rank(A)   

2 % 表示A秩数为2且等于矩阵的列数   

＞> inv(A) % 反矩阵   

ans =   

1.5000 -0.5000   

-2.0000 1.0000   

＞> B=[2 1; 3 2; 4 5]; % B为奇异矩阵   

＞> rank(B)   

ans =   

2 % 表示B秩数为2，但是其列数为3   

＞> inv(B)   

??? Error using ==> inv   

Matrix must be square.   


相信大家都会计算矩阵行列式的值，但是如一矩阵大小超过4以上，行列式值的计算就会
繁复。MATLAB提供计算行列式的函数，其语法为det(A)，例如：   

＞> A=[1 3 0; -1 5 2; 1 2 1];   

＞> det(A) % 矩阵之行列式值   

ans =   

10  
沪p算就会繁复。MATLAB提供计算行列式的函数，其语法为det(A)，例如：   

＞> A=[1 3 0; -1 5 2; 1 2 1];   

＞> det(A) % 矩阵之行列式值   

ans =   

10  



------------------------------------------------------------------------------
--  
   

HUASHI3483

5.3.3 矩阵分解  

------------------------------------------------------------------------------
--  

矩阵分解 (decomposition, factorization)是多半将矩阵拆解为数个三角形矩阵
(triangular matrix)，依使用目的的不同，可分为三种矩阵分解法：1)三角分解法  
(Triangular Factorization)，2)QR 分解法 (QR Factorization)，3)奇异值分解法  
(Singular Value Decompostion)。  


(1) 三角分解法  


三角分解法是将原正方 (square) 矩阵分解成一个上三角形矩阵　或是排列(permuted)  
的上三角形矩阵　和一个下三角形矩阵，这样的分解法又称为LU分解法。它的用途主要在
简化一个大矩阵的行列式值的计算过程，求反矩阵，和求解联立方程组。不过要注意这种
分解法所得到的上下三角形矩阵并非唯一，还可找到数个不同的一对上下三角形矩阵，此
两三角形矩阵相乘也会得到原矩阵。   


我们举以下二个矩阵为例：   

   

利用三角分解法可将A和B二矩阵分别拆解为上下三角形矩阵   

   

注意B分解的矩阵得到的第一个矩阵[LB]是排列的下三角形矩阵，如果第二、三列互换，
则此变成完全的下三角形矩阵。   


以MATLAB函数计算上述的LU分解法，其语法为[L,U]=lu(A)，其中L代表下三角形矩阵U代
表上三角形矩阵。我们来看一个例子。   

＞> A = [1 2 -1, -2 -5 3; -1 -3 0]; B=[1 3 2; -2 -6 1; 2 5 7];   

＞> [L1,U1] = lu(A); [L2,U2] = lu(B);   

＞> L1; U1   

L1 = % 注意这个矩阵L1和之前的[LA]不相同   

-0.5 1 0   

1 0 0   

0.5 1 1   

U1 = % 注意这个矩阵U1和之前的[UA]不相同   

-2 -5 3   

0 -0.5 0.5   

0 0 -2   

＞> L2; U2   

L2 = % 注意这个矩阵L2和之前的[LB]不相同   

-0.5 0 1   

1 0 0   

-1 1 0   

U2 = % 注意这个矩阵U2和之前的[UB]不相同   

-2 -6 1   

0 -1 8   

0 0 2.5  


(2) QR分解法  


QR分解法是将矩阵分解成一个正规正交矩阵与上三角形矩阵。还记得先前我们介绍的正规
正交矩阵Q满足的条件吗！所以称为QR分解法与此正规正交矩阵的通用符号Q有关。   


MATLAB以qr函数来执行QR分解法，其语法为[Q,R]=qr(A)，其中Q代表正规正交矩阵，而R
代表上三角形矩阵。此外，原矩阵A不必为正方矩阵；如果矩阵A大小为，则矩阵Q大小为
，矩阵R大小为。   
(3) 奇异值分解法   

奇异值分解 (sigular value decomposition,SVD) 是另一种正交矩阵分解法；SVD是最可
靠的分解法，但是它比QR 分解法要花上近十倍的计算时间。[U,S,V]=svd(A)，其中U和V
代表二个相互正交矩阵，而S代表一对角矩阵。和QR分解法相同者，原矩阵A不必为正方矩
阵。   

使用SVD分解法的用途是解最小平方误差法和数据压缩。   



------------------------------------------------------------------------------

HUASHI3483

第六章 联立线性方程式  

------------------------------------------------------------------------------
--  
工程问题挑战：电动车性能  

许多工程问题的分析，例如电路分析及梁结构受静力的问题皆涉及解联立线性方程式。解
联立线性方程式有不同的方法，我们在这章仅介绍最直接的方法--矩阵运算方式求解。  



在这章中提到的工程问题挑战「电动车性能」，在工业界和学术界已发展了数十年，到目
前为止所遭遇的瓶颈还是在缺乏可提供高续航力的电源。请参考有关电动车的介绍及发展
现况。   

电动车资料库系统 http://www2.seeder.net.tw/evs/
EV Information Network http://www.radix.net/~futurev/



------------------------------------------------------------------------------
--  

第六章 联立线性方程式   
6.1 利用矩阵解法   
6.2 范例问题：电路分析   

------------------------------------------------------------------------------
--  

HUASHI3483

6.1 利用矩阵解法  

------------------------------------------------------------------------------
--  

假设一组联立线性方程式为   

   

我们习惯将上组方程式以矩阵方式表示如下   

AX=B   

其中 A 为等式左边各方程式的系数项，X 为欲求解的未知项，B 代表等式右边之已知项  


   

要解上述的联立方程式，我们可以利用在第六章介绍的矩阵左除 \ 做运算，即是 X=A\B
。   


如果将原方程式改写成 XA=B，且令 X, A 和 B 分别为   

   

注意上式的 X, B 已改写成列向量，A其实是前一个方程式中 A 的转置矩阵。上式的 X  
可以矩阵右除 / 求解，即是 X=B/A。   


若以反矩阵运算求解 AX=B, X=B，即是 X=inv(A)*B，或是改写成 XA=B, X=B，即是  
X=B*inv(A)。   


我们直接以下面的例子来说明这三个运算的用法：   

＞> A=[3 2 -1; -1 3 2; 1 -1 -1]; % 将等式的左边系数键入   

＞> B=[10 5 -1]'; % 将等式右边之已知项键入，B要做转置   

＞> X=A\B % 先以左除运算求解   

X = % 注意X为行向量   

-2   

5   

6   

＞> C=A*X % 验算解是否正确   

C = % C=B   

10   

5   

-1  


＞> A=A'; % 将A先做转置   

＞> B=[10 5 -1];   

＞> X=B/A % 以右除运算求解的结果亦同   

X = % 注意X为列向量   

10 5 -1   

＞> X=B*inv(A); % 也可以反矩阵运算求解  



------------------------------------------------------------------------------
--  
   

HUASHI3483

7.1.1 一维内插  

------------------------------------------------------------------------------
--  

线性内插是假设在二个已知数据中的变化为线性关系，因此可由已知二点的座标(a, b)去
计算通过这二点的斜线，公式如下：   

   

其中 a<b<c 在上式的 b 点即是代表要内插的点，f(b) 则是要计算的内插函数值。下图
即是一个以二种内插法的比较   

\pcxfile[12cm,5cm]{fig9_1.pcx}   

\caption{线性式与 spline 函数的曲线契合}  


线性内插是最简单的内插方法，但其适用范围很小；如果原来数据的函数f有极大的变化
，假设其数据点之间为线性变化并不合理。所以我们可以用二次、三次方程式或是另一种
称为spline函数来近似原来数据的函数。MATLAB的一维内插函数是interp1，其语法为
interp1(x,y,xi),interp1(x,y,xi,'method')；其中的x,y是原已知的数据的x,y值，而xi
则是要内插的数据点，另外method可以设定内插方法有 linear,cubic,spline，分别是一
次、三次方程式和spline函数，其中预设方法是linear。如果数据的变化较大，以  
spline函数内插所形成的曲线最平滑，所以效果最好。而三次方程式所得到的内插曲线平
滑度，则介于线性与spline函数之间。  


我们以下面的例子说明。假设有一个汽车引擎在定转速下，温度与时间（单位为sec）的
三次量测值如下  

time temp1 temp2 temp3   
0 0 0 0   
1 20 110 176   
2 60 180 220   
3 68 240 349   
4 77 310 450   
5 110 405 503   


其中温度的数据从 20oC变化到 503oC，如果要估计在t=2.6, 4.9 sec 的温度，可以下列
指令计算   

＞> x=[0 1 2 3 4 5]'; % 键入时间   

＞> y=[0 20 60 68 77 110]'; % 键入第一组时间   

＞> y1=interp1(x,y,2.6) % 要内插的数据点为 2.6   

y1 = % 对应 2.6 的函数值为 64.8   

64.8   

＞> y1=interp1(x,y,[2.6 4.9]) % 内插数据点为 2.6, 4.9，注意用[ ]将多个内插点放
在其中   

y1 =   

64.8   

106.7   

＞> y1=interp1(x,y,2.6,'cubic') % 以三次方程式对数据点 2.6 作内插   

y1 = % 对应 2.6 的函数值为 66.264   

66.264   

＞> y1=interp1(x,y,2.6,'spline') % 以spline函数对数据点 2.6 作内插   

y1 = % 对应 2.6 的函数值为 66.368   

66.368  


以下的例子还配合绘图功能，用以比较不同内插方法的差异。   

＞> h=1:12;   

＞> temp=[5 8 9 15 25 29 31 30 22 25 27 24]; % 这组温度数据变化较大   

＞> plot(h,temp,'--',h,temp,'+') % 将线性内插结果绘图   

＞> h_3=1:0.1:12 % 要每0.1小时估计一次温度值   

＞> t_3=interp1(h,temp,h_3,'cubic') % 以三次方程式做内插   

＞> t_s=interp1(h,temp,h_3,'spline') % 以spline函数做内插   

＞> hold on   

＞> subplot(1,2,1)   

＞> plot(h,temp,'--',h,temp,'+',h_3,t_3) % 将线性及三次方程式内插绘图   

＞> subplot(1,2,2)   

＞> plot(h,temp,'--',h,temp,'+',h_3,t_s) % 将线性方程式及spline内插绘图   

＞> hold off  



------------------------------------------------------------------------------
--  
   

HUASHI3483

7.1.3 Spline 内插  

------------------------------------------------------------------------------
--  

关于spline内插我们在7.1.1节已介绍过，它可以用interp1指定内插方式为spline来做。
另一种方式也可以用 spline(x,y,xi)来做，其中的x,y,xi的用法与interp1中的语法相同
。事实上这二种方法采用相同的spline 函数做运算，也就是当我们执行
interp1(x,y,xi,'spline')时，MATLAB即呼叫spline(x,y,xi)做运算，再将计算结果传回
interp1。   


我们还是以在 7.1.1 节相同的数据做 spline 内插示范，   

＞> x=[0 1 2 3 4 5]';   

＞> y=[0 20 60 68 77 110]';   

＞> y1=spline(x,y,2.6)   

y1 =   

67.3   

＞> y1=spline(x,y,[2.6,4.9])   

y1 =   

67.3 105.2  



------------------------------------------------------------------------------
--  

HUASHI3483

7.3.1 线性回归  

------------------------------------------------------------------------------
--  

我们以一简单数据组来说明什么是线性回归。假设有一组数据型态为 y=y(x)，其中   

x={0, 1, 2, 3, 4, 5}, y={0, 20, 60, 68, 77, 110}   

如果我们要以一个最简单的方程式来近似这组数据，则非一阶的线性方程式莫属。先将这
组数据绘图如下   



图中的斜线是我们随意假设一阶线性方程式 y=20x，用以代表这些数据的一个方程式。以
下将上述绘图的 MATLAB 指令列出，并计算这个线性方程式的 y 值与原数据 y 值间误差
平方的总合。   

＞> x=[0 1 2 3 4 5];   

＞> y=[0 20 60 68 77 110];   

＞> y1=20*x; % 一阶线性方程式的 y1 值   

＞> sum_sq = sum(y-y1).^2); % 误差平方总合为 573   

＞> axis([-1,6,-20,120])   

＞> plot(x,y1,x,y,'o'), title('Linear estimate'), grid  


如此任意的假设一个线性方程式并无根据，如果换成其它人来设定就可能采用不同的线性
方程式；所以我们须要有比较精确方式决定理想的线性方程式。我们可以要求误差平方的
总合为最小，做为决定理想的线性方程式的准则，这样的方法就称为最小平方误差
(least squares error)或是线性回归。MATLAB的polyfit函数提供了从一阶到高阶多项式
的回归法，其语法为polyfit(x,y,n)，其中x,y为输入数据组n为多项式的阶数，n=1就是
一阶的线性回归法。polyfit函数所建立的多项式可以写成   

   

从polyfit函数得到的输出值就是上述的各项系数，以一阶线性回归为例n=1，所以只有  
二个输出值。如果指令为coef=polyfit(x,y,n)，则coef(1)= , coef(2)=,...,coef(n+1)
= 。注意上式对n 阶的多项式会有 n+1 项的系数。我们来看以下的线性回归的示范：   

＞> x=[0 1 2 3 4 5];   

＞> y=[0 20 60 68 77 110];   

＞> coef=polyfit(x,y,1); % coef 代表线性回归的二个输出值   

＞> a0=coef(1); a1=coef(2);   

＞> ybest=a1*x+a0; % 由线性回归产生的一阶方程式   

＞> sum_sq=sum(y-ybest).^2); % 误差平方总合为 356.82   

＞> axis([-1,6,-20,120])   

＞> plot(x,ybest,x,y,'o'), title('Linear regression estimate'), grid  



------------------------------------------------------------------------------
--