matlab处理矩阵的原理,MALAB原理及编程

MATLAB学习摘要

数组操作和矩阵操作(Array Operations vs. Matrix Operations)

对矩阵的元素一个一个孤立进行的操作称作数组操作；而把矩阵视为一个整体进行的运算则成为矩阵操作。布朗数组操作(Boolean Array Operations)对矩阵的比较运算是数组操作，也就是说，是对每个元素孤立进行的。因此其结果就不是一个“真”或者“假”，而是一堆“真假”。这个结果就是布朗数组。

1.1矩阵的超越函数

(1)矩阵平方根sqrtm

sqrtm(A)计算矩阵A的平方根。

(2)矩阵对数logm

logm(A)计算矩阵A的自然对数。此函数输入参数的条件与输出结果间的关系和函数sqrtm(A)完全一样。

(3)矩阵指数expm、expm1、expm2、expm3

expm(A)、expm1(A)、expm2(A)、expm3(A)的功能都求矩阵指数eA。

(4)普通矩阵函数funm

funm(A,‘fun’)用来计算直接作用于矩阵A的由‘fun’指定的超越函数值。当fun取sqrt时，funm(A,‘sqrt’)可以计算矩阵A的平方根，与sqrtm(A)的计算结果一样。

1.2程序的暂停

暂停程序的执行可以使用pause函数，其调用格式为：

pause(延迟秒数)

如果省略延迟时间，直接使用pause，则将暂停程序，直到用户按任一键后程序继续执行。若要强行中止程序的运行可使用Ctrl+C命令。

例3-5输入一个字符，若为大写字母，则输出其对应的小写字母；若为小写字母，则输出其对应的大写字母；若为数字字符则输出其对应的数值，若为其他字符则原样输出。

c=input('请输入一个字符','s');

if c>='A' & c<='Z'

disp(setstr(abs(c)+abs('a')-abs('A')));

elseif c>='a'& c<='z'

disp(setstr(abs(c)- abs('a')+abs('A')));

elseif c>='0'& c<='9'

disp(abs(c)-abs('0'));

else

disp(c);

end

1.3 switch语句

switch语句根据表达式的取值不同，分别执行不同的语句，其语句格式为：

switch表达式

case表达式1

语句组1

case表达式2

语句组2

……

case表达式m

语句组m

otherwise

语句组n

end

当表达式的值等于表达式1的值时，执行语句组1，当表达式的值等于表达式2的值时，执行语句组2，…，当表达式的值等于表达式m的值时，执行语句组m，当表达式的值不等于case所列的表达式的值时，执行语句组n。当任意一个分支的语句执行完后，直接执行switch语句的下一句。

例3-6某商场对顾客所购买的商品实行打折销售，标准如下(商品价格用price来表示)：

price<200没有折扣

200≤price<5003%折扣

500≤price<10005%折扣

1000≤price<25008%折扣

2500≤price<500010%折扣

5000≤price14%折扣

输入所售商品的价格，求其实际销售价格。

程序如下：

price=input('请输入商品价格');

switch fix(price/100)

case {0,1}%价格小于200

rate=0;

case {2,3,4}%价格大于等于200但小于500

rate=3/100;

case num2cell(5:9)%价格大于等于500但小于1000

rate=5/100;

case num2cell(10:24)%价格大于等于1000但小于2500

rate=8/100;

case num2cell(25:49)%价格大于等于2500但小于5000

rate=10/100;

otherwise%价格大于等于5000

rate=14/100;

end

price=price*(1-rate)%输出商品实际销售价格

1.4 try语句

语句格式为：

try

语句组1

catch

语句组2

end

try语句先试探性执行语句组1，如果语句组1在执行过程中出现错误，则将错误信息赋给保留的lasterr变量，并转去执行语句组2。

1.5循环结构

1．for语句

for语句的格式为：

for循环变量=表达式1:表达式2:表达式3

循环体语句

end

其中表达式1的值为循环变量的初值，表达式2的值为步长，表达式3的值为循环变量的终值。步长为1时，表达式2可以省略。

例3-8一个三位整数各位数字的立方和等于该数本身则称该数为水仙花数。输出全部水仙花数。

程序如下：

for m=100:999

m1=fix(m/100);%求m的百位数字

m2=rem(fix(m/10),10);%求m的十位数字

m3=rem(m,10);%求m的个位数字

if m==m1*m1*m1+m2*m2*m2+m3*m3*m3

disp(m)

end

end

s=0;

a=[12,13,14;15,16,17;18,19,20;21,22,23];

for k=a

s=s+k;

end

disp(s');

39485766

例3-11从键盘输入若干个数，当输入0时结束输入，求这些数的平均值和它们之和。

程序如下：

sum=0;

cnt=0;

val=input('Enter a number (end in 0):');

while (val~=0)

sum=sum+val;

cnt=cnt+1;

end

if (cnt > 0)

sum

mean=sum/cnt

end

1.6 break语句和continue语句

与循环结构相关的语句还有break语句和continue语句。它们一般与if语句配合使用。

break语句用于终止循环的执行。当在循环体内执行到该语句时，程序将跳出循环，继续执行循环语句的下一语句。

continue语句控制跳过循环体中的某些语句。当在循环体内执行到该语句时，程序将跳过循环体中所有剩下的语句，继续下一次循环。

例3-12求[100，200]之间第一个能被21整除的整数。

程序如下：

for n=100:200

if rem(n,21)~=0

continue

end

break

end

n

1.7函数调用

函数调用的一般格式是：

[输出实参表]=函数名(输入实参表)

要注意的是，函数调用时各实参出现的顺序、个数，应与函数定义时形参的顺序、个数一致，否则会出错。函数调用时，先将实参传递给相应的形参，从而实现参数传递，然后再执行函数的功能。

(1)函数文件的基本结构

函数文件由function语句引导，其基本结构为：

function输出形参表=函数名(输入形参表)

注释说明部分

函数体语句

其中以function开头的一行为引导行，表示该M文件是一个函数文件。函数名的命名规则与变量名相同。输入形参为函数的输入参数，输出形参为函数的输出参数。当输出形参多于一个时，则应该用方括号括起来。

例3-14编写函数文件求半径为r的圆的面积和周长。

函数文件如下：

function [s,p]=fcircle(r)

%CIRCLEcalculate the area and perimeter of a circle of radii r

%r圆半径

%s圆面积

%p圆周长

%2004年7月30日编

s=pi*r*r;

p=2*pi*r;

例3-15利用函数文件，实现直角坐标(x,y)与极坐标(ρ,θ)之间的转换。

函数文件tran.m：

function [rho,theta]=tran(x,y)

rho=sqrt(x*x+y*y);

theta=atan(y/x);

调用tran.m的命令文件

x=input('Please input x=:');

y=input('Please input y=:');

[rho,the]=tran(x,y);

rho

the

在MATLAB中，函数可以嵌套调用，即一个函数可以调用别的函数，甚至调用它自身。一个函数调用它自身称为函数的递归调用。

例3-16利用函数的递归调用，求n！。

n!本身就是以递归的形式定义的：

显然，求n!需要求(n-1)!，这时可采用递归调用。递归调用函数文件factor.m如下：

function f=factor(n)

if n<=1

f=1;

else

f=factor(n-1)*n;%递归调用求(n-1)!

end

(2) 函数参数的可调性

在调用函数时，MATLAB用两个永久变量nargin和nargout分别记录调用该函数时的输入实参和输出实参的个数。只要在函数文件中包含这两个变量，就可以准确地知道该函数文件被调用时的输入输出参数个数，从而决定函数如何进行处理。

例3-17nargin用法示例。

函数文件examp.m：

function fout=charray(a,b,c)

if nargin==1

fout=a;

elseif nargin==2

fout=a+b;

elseif nargin==3

fout=(a*b*c)/2;

end

命令文件mydemo.m：

x=[1:3];

y=[1;2;3];

examp(x)

examp(x,y')

examp(x,y,3)

1.8程序举例

例3-19猜数游戏。首先由计算机产生[1,100]之间的随机整数，然后由用户猜测所产生的随机数。根据用户猜测的情况给出不同提示，如猜测的数大于产生的数，则显示“High”，小于则显示“Low”，等于则显示“You won”，同时退出游戏。用户最多可以猜7次。

例3-20用筛选法求某自然数范围内的全部素数。

素数是大于1，且除了1和它本身以外，不能被其他任何整数所整除的整数。用筛选法求素数的基本思想是：要找出2 - m之间的全部素数，首先在2- m中划去2的倍数(不包括2)，然后划去3的倍数(不包括3)，由于4已被划去，再找5的倍数(不包括5)，…，直到再划去不超过的数的倍数，剩下的数都是素数。

例3-21设，求s=。

求函数f(x)在[a,b]上的定积分，其几何意义就是求曲线y=f(x)与直线x=a，x=b，y=0所围成的曲边梯形的面积。为了求得曲边梯形面积，先将积分区间[a,b]分成n等分，每个区间的宽度为h=(b-a)/n，对应地将曲边梯形分成n等分，每个小部分即是一个小曲边梯形。近似求出每个小曲边梯形面积，然后将n个小曲边梯形的面积加起来，就得到总面积，即定积分的近似值。近似地求每个小曲边梯形的面积，常用的方法有：矩形法、梯形法以及辛普生法等。

2图形处理

2.1图形标注

有关图形标注函数的调用格式为：

title(图形名称)-标题

xlabel(x轴说明)

ylabel(y轴说明)

text(x,y,图形说明)-在指定位置添加说明。

legend(图例1,图例2,…)

2.2二维统计分析图

在MATLAB中，二维统计分析图形很多，常见的有条形图、阶梯图、杆图和填充图等，所采用的函数分别是：

bar(x,y,选项)-竖直条形图

barh:产生一个水平条形图

stairs(x,y,选项)

stem(x,y,选项)

fill(x1,y1,选项1,x2,y2,选项2,…)

例5-13分别以条形图、阶梯图、杆图和填充图形式绘制曲线y=2sin(x)。

程序如下：

x=0:pi/10:2*pi;

y=2*sin(x);

subplot(2,2,1);bar(x,y,'g');

title('bar(x,y,''g'')');axis([0,7,-2,2]);

subplot(2,2,2);stairs(x,y,'b');

title('stairs(x,y,''b'')');axis([0,7,-2,2]);

subplot(2,2,3);stem(x,y,'k');

title('stem(x,y,''k'')');axis([0,7,-2,2]);

subplot(2,2,4);fill(x,y,'y');

title('fill(x,y,''y'')');axis([0,7,-2,2]);

2.3三维曲面

1．产生三维数据

在MATLAB中，利用meshgrid函数产生平面区域内的网格坐标矩阵。其格式为：

x=a:d1:b; y=c:d2:d;

[X,Y]=meshgrid(x,y);

语句执行后，矩阵X的每一行都是向量x，行数等于向量y的元素的个数，矩阵Y的每一列都是向量y，列数等于向量x的元素的个数。

2．绘制三维曲面的函数

surf函数和mesh函数的调用格式为：

mesh(x,y,z,c)

surf(x,y,z,c)，

surf(z)，z为已知矩阵，则使用矩阵的行和列的索引作为x和y的坐标。

一般情况下，x,y,z是维数相同的矩阵。x,y是网格坐标矩阵，z是网格点上的高度矩阵，c用于指定在不同高度下的颜色范围。

例5-17绘制三维曲面图z=sin(x+sin(y))-x/10。

程序如下：

[x,y]=meshgrid(0:0.25:4*pi);

z=sin(x+sin(y))-x/10;

mesh(x,y,z);

axis([0 4*pi 0 4*pi -2.5 1]);

此外，还有带等高线的三维网格曲面函数meshc和带底座的三维网格曲面函数meshz。其用法与mesh类似，不同的是meshc还在xy平面上绘制曲面在z轴方向的等高线，meshz还在xy平面上绘制曲面的底座。

2.4标准三维曲面

sphere函数的调用格式为：

[x,y,z]=sphere(n)%生成三维数据网格(球面)

cylinder函数的调用格式为：

[x,y,z]= cylinder(R,n)

MATLAB还有一个peaks函数，称为多峰函数，常用于三维曲面的演示。(柱面)

2.5其他三维图形

在介绍二维图形时，曾提到条形图、杆图、饼图和填充图等特殊图形，它们还可以以三维形式出现，使用的函数分别是bar3、stem3、pie3和fill3。

bar3函数绘制三维条形图，常用格式为：

bar3(y)

bar3(x,y)

barsh():产生三位水平条形图

stem3函数绘制离散序列数据的三维杆图，常用格式为：

stem3(z)

stem3(x,y,z)

pie3函数绘制三维饼图，常用格式为：

pie3(x)

fill3函数等效于三维函数fill，可在三维空间内绘制出填充过的多边形，常用格式为：

fill3(x,y,z,c)

2.6色彩处理

1．颜色的向量表示

MATLAB除用字符表示颜色外，还可以用含有3个元素的向量表示颜色。向量元素在[0,1]范围取值，3个元素分别表示红、绿、蓝3种颜色的相对亮度，称为RGB三元组。

2．色图

色图(Color map)是MATLAB系统引入的概念。在MATLAB中，每个图形窗口只能有一个色图。色图是m×3的数值矩阵，它的每一行是RGB三元组。色图矩阵可以人为地生成，也可以调用MATLAB提供的函数来定义色图矩阵。

3．三维表面图形的着色

三维表面图实际上就是在网格图的每一个网格片上涂上颜色。surf函数用缺省的着色方式对网格片着色。除此之外，还可以用shading命令来改变着色方式。

shading faceted命令将每个网格片用其高度对应的颜色进行着色，但网格线仍保留着，其颜色是黑色。这是系统的缺省着色方式。

shading flat命令将每个网格片用同一个颜色进行着色，且网格线也用相应的颜色，从而使得图形表面显得更加光滑。

shading interp命令在网格片内采用颜色插值处理，得出的表面图显得最光滑。

例5-233种图形着色方式的效果展示。

程序如下：

[x,y,z]=sphere(20);

colormap(copper);

subplot(1,3,1);

surf(x,y,z);

axis equal

subplot(1,3,2);

surf(x,y,z);shading flat;

axis equal

subplot(1,3,3);

surf(x,y,z);shading interp;

axis equal

4光照处理

MATLAB提供了灯光设置的函数，其调用格式为：

light('Color',选项1,'Style',选项2,'Position',选项3)

例5-24光照处理后的球面。

程序如下：

[x,y,z]=sphere(20);

subplot(1,2,1);

surf(x,y,z);axis equal;

light('Posi',[0,1,1]);

shading interp;

hold on;

plot3(0,1,1,'p');text(0,1,1,' light');

subplot(1,2,2);

surf(x,y,z);axis equal;

light('Posi',[1,0,1]);

shading interp;

hold on;

plot3(1,0,1,'p');text(1,0,1,' light');

meshgrid函数

已知矩阵x,y,则[new_x,new_y]=meshgrid(x,y)形成新的矩阵，新矩阵的行数为x矩阵元素的个数，列数为y矩阵元素的个数。例如

a =

1235

2469

57119

b =

101214161820

[new_a,new_b]=meshgrid(a,b)

new_a =

125 2473611599

1252473611599

1252473611599

1252473611599

1252473611599

1252473611599

new_b =

101010101010101010101010

121212121212121212121212

141414141414141414141414

161616161616161616161616

181818181818181818181818

202020202020202020202020

5图形的裁剪处理

例5-25绘制三维曲面图，并进行插值着色处理，裁掉图中x和y都小于0部分。

程序如下：

[x,y]=meshgrid(-5:0.1:5);

z=cos(x).*cos(y).*exp(-sqrt(x.^2+y.^2)/4);

surf(x,y,z);shading interp;

pause%程序暂停

i=find(x<=0&y<=0);

z1=z;z1(i)=NaN;

surf(x,y,z1);shading interp;

为了展示裁剪效果，第一个曲面绘制完成后暂停，然后显示裁剪后的曲面。

2.7图像处理与动画制作

1．imread和imwrite函数

imread和imwrite函数分别用于将图像文件读入MATLAB工作空间，以及将图像数据和色图数据一起写入一定格式的图像文件。MATLAB支持多种图像文件格式，如.bmp、.jpg、.jpeg、.tif等。

2．image和imagesc函数

这两个函数用于图像显示。为了保证图像的显示效果，一般还应使用colormap函数设置图像色图。

例5-26有一图像文件flower.jpg，在图形窗口显示该图像。

程序如下：

[x,cmap]=imread('flower.jpg');%读取图像的数据阵和色图阵

image(x);colormap(cmap);

axis image off%保持宽高比并取消坐标轴

MATLAB提供getframe、moviein和movie函数进行动画制作。

1．getframe函数

getframe函数可截取一幅画面信息(称为动画中的一帧)，一幅画面信息形成一个很大的列向量。显然，保存n幅图面就需一个大矩阵。

2．moviein函数

moviein(n)函数用来建立一个足够大的n列矩阵。该矩阵用来保存n幅画面的数据，以备播放。之所以要事先建立一个大矩阵，是为了提高程序运行速度。

3．movie函数

movie(m,n)函数播放由矩阵m所定义的画面n次，缺省时播放一次。

例5-27绘制了peaks函数曲面并且将它绕z轴旋转。

程序如下

[X,Y,Z]=peaks(30);

surf(X,Y,Z)

axis([-3,3,-3,3,-10,10])

axis off;%取消坐标轴

shading interp;

colormap(hot);

m=moviein(20);%建立一个20列大矩阵

for i=1:20

view(-37.5+24*(i-1),30)%改变视点

m(:,i)=getframe;%将图形保存到m矩阵

end

movie(m,2);%播放画面2次

%%%%

disp(),fprint()

3数据分析与多项式的计算

(1)数据求和与求积

数据序列求和与求积的函数是sum和prod，其使用方法类似。设X是一个向量，A是一个矩阵，函数的调用格式为：

sum(X)：返回向量X各元素的和。

prod(X)：返回向量X各元素的乘积。

sum(A)：返回一个行向量，其第i个元素是A的第i列的元素和。

prod(A)：返回一个行向量，其第i个元素是A的第i列的元素乘积。

sum(A,dim)：当dim为1时，该函数等同于sum(A)；当dim为2时，返回一个列向量，其第i个元素是A的第i行的各元素之和。

prod(A,dim)：当dim为1时，该函数等同于prod(A)；当dim为2时，返回一个列向量，其第i个元素是A的第i行的各元素乘积。

(2)累加和与累乘积

在MATLAB中，使用cumsum和cumprod函数能方便地求得向量和矩阵元素的累加和与累乘积向量，函数的调用格式为：

cumsum(X)：返回向量X累加和向量。

cumprod(X)：返回向量X累乘积向量。

cumsum(A)：返回一个矩阵，其第i列是A的前i列的累加和向量。

cumprod(A)：返回一个矩阵，其第i列是A的前i列的累乘积向量。

cumsum(A,dim)：当dim为1时，该函数等同于cumsum(A)；当dim为2时，返回一个矩阵，其第i行是A的前i行的累加和向量。

cumprod(A,dim)：当dim为1时，该函数等同于cumprod(A)；当dim为2时，返回一个向量，其第i行是A的前i行的累乘积向量。

(3)求数据序列平均值的函数是mean，求数据序列中值的函数是median。排序函数：sort

(4)标准方差与相关系数

在MATLAB中，提供了计算数据序列的标准方差的函数std。对于向量X，std(X)返回一个标准方差。对于矩阵A，std(A)返回一个行向量，它的各个元素便是矩阵A各列或各行的标准方差。std函数的一般调用格式为：

Y=std(A,flag,dim)

其中dim取1或2。当dim=1时，求各列元素的标准方差；当dim=2时，则求各行元素的标准方差。flag取0或1，当flag=0时，按σ1所列公式计算标准方差，当flag=1时，按σ2所列公式计算标准方差。缺省flag=0，dim=1。

(5) 数据插值

①一维数据插值

在MATLAB中，实现这些插值的函数是interp1，其调用格式为：

Y1=interp1(X,Y,X1,'method')

函数根据X,Y的值，计算函数在X1处的值。X,Y是两个等长的已知向量，分别描述采样点和样本值，X1是一个向量或标量，描述欲插值的点，Y1是一个与X1等长的插值结果。method是插值方法，允许的取值有‘linear’、‘nearest’、‘cubic’、‘spline’。

MATLAB中有一个专门的3次样条插值函数Y1=spline(X,Y,X1)，其功能及使用方法与函数Y1=interp1(X,Y,X1,‘spline’)完全相同。

② 二维数据插值

在MATLAB中，提供了解决二维插值问题的函数interp2，其调用格式为：

Z1=interp2(X,Y,Z,X1,Y1,'method')

其中X,Y是两个向量，分别描述两个参数的采样点，Z是与参数采样点对应的函数值，X1,Y1是两个向量或标量，描述欲插值的点。Z1是根据相应的插值方法得到的插值结果。method的取值与一维插值函数相同。X,Y,Z也可以是矩阵形式。

同样，X1,Y1的取值范围不能超出X,Y的给定范围，否则，会给出“NaN”错误。

例6-13某实验对一根长10米的钢轨进行热源的温度传播测试。用x表示测量点0:2.5:10(米)，用h表示测量时间0:30:60(秒)，用T表示测试所得各点的温度(℃)。试用线性插值求出在一分钟内每隔20秒、钢轨每隔1米处的温度TI。

命令如下：

x=0:2.5:10;

h=[0:30:60]';

T=[95,14,0,0,0;88,48,32,12,6;67,64,54,48,41];

xi=[0:10];

hi=[0:20:60]';

TI=interp2(x,h,T,xi,hi)

(6) 曲线拟合

在MATLAB中，用polyfit函数来求得最小二乘拟合多项式的系数，再用polyval函数按所得的多项式计算所给出的点上的函数近似值。

polyfit函数的调用格式为：

[P,S]=polyfit(X,Y,m)

函数根据采样点X和采样点函数值Y，产生一个m次多项式P及其在采样点的误差向量S。其中X,Y是两个等长的向量，P是一个长度为m+1的向量，P的元素为多项式系数。

polyval函数的功能是按多项式的系数计算x点多项式的值。？？？

(7) 多项式运算

(1)乘conv(p1，p2)

除[Q，r]=deconv(p1，p2)

(2)多项式的导函数

对多项式求导数的函数是：

p=polyder(P)：求多项式P的导函数

p=polyder(P,Q)：求P·Q的导函数

[p,q]=polyder(P,Q)：求P/Q的导函数，导函数的分子存入p，分母存入q。上述函数中，参数P,Q是多项式的向量表示，结果p,q也是多项式的向量表示。

(3)多项式的求值

MATLAB提供了两种求多项式值的函数：polyval与polyvalm，它们的输入参数均为多项式系数向量P和自变量x。两者的区别在于前者是代数多项式求值，而后者是矩阵多项式求值。

① 代数多项式求值

polyval函数用来求代数多项式的值，其调用格式为：

Y=polyval(P,x)

若x为一数值，则求多项式在该点的值；若x为向量或矩阵，则对向量或矩阵中的每个元素求其多项式的值。

例6-19已知多项式x4+8x3-10，分别取x=1.2和一个2×3矩阵为自变量计算该多项式的值。

② 矩阵多项式求值

polyvalm函数用来求矩阵多项式的值，其调用格式与polyval相同，但含义不同。polyvalm函数要求x为方阵，它以方阵为自变量求多项式的值。设A为方阵，P代表多项式x3-5x2+8，那么polyvalm(P,A)的含义是：

A*A*A-5*A*A+8*eye(size(A))

而polyval(P,A)的含义是：

A.*A.*A-5*A.*A+8*ones(size(A))

例6-20仍以多项式x4+8x3-10为例，取一个2×2矩阵为自变量分别用polyval和polyvalm计算该多项式的值。

(4)多项式求根函数：x=roots(A)

A为多项式的系数矩阵(向量)。

若已知多项式的全部根，则可以用poly函数建立起该多项式，其调用格式为：

P=poly(x)

若x为具有n个元素的向量，则poly(x)建立以x为其根的多项式，且将该多项式的系数赋给向量P。？？

8．利用矩阵的分解求解线性方程组

矩阵分解是指根据一定的原理用某种算法将一个矩阵分解成若干个矩阵的乘积。常见的矩阵分解有LU分解、QR分解、Cholesky分解，以及Schur分解、Hessenberg分解、奇异分解等。

(1) LU分解

矩阵的LU分解就是将一个矩阵表示为一个交换下三角矩阵和一个上三角矩阵的乘积形式。线性代数中已经证明，只要方阵A是非奇异的，LU分解总是可以进行的。

MATLAB提供的lu函数用于对矩阵进行LU分解，其调用格式为：

[L,U]=lu(A)：产生一个上三角阵U和一个变换形式的下三角阵L(行交换)，使之满足A=LU。注意，这里的矩阵A必须是方阵。

[L,U,P]=lu(A)：产生一个上三角阵U和一个下三角阵L以及一个置换矩阵P，使之满足PA=LU。当然矩阵A同样必须是方阵。

实现LU分解后，线性方程组Ax=b的解x=U(Lb)或x=U(LPb)，这样可以大大提高运算速度。

例7-2用LU分解求解例7-1中的线性方程组。

命令如下：

A=[2,1,-5,1;1,-5,0,7;0,2,1,-1;1,6,-1,-4];

b=[13,-9,6,0]';

[L,U]=lu(A);

x=U(Lb)

或采用LU分解的第2种格式，命令如下：

[L,U ,P]=lu(A);

x=U(LP*b)

(2) QR分解

对矩阵X进行QR分解，就是把X分解为一个正交矩阵Q和一个上三角矩阵R的乘积形式。QR分解只能对方阵进行。MATLAB的函数qr可用于对矩阵进行QR分解，其调用格式为：

[Q,R]=qr(A)：产生一个一个正交矩阵Q和一个上三角矩阵R，使之满足X=QR。

[Q,R,E]=qr(A)：产生一个一个正交矩阵Q、一个上三角矩阵R以及一个置换矩阵E，使之满足XE=QR。

实现QR分解后，线性方程组Ax=b的解x=R(Qb)或x=E(R(Qb))。

(3) Cholesky分解

如果矩阵A是对称正定的，则Cholesky分解将矩阵A分解成一个下三角矩阵和上三角矩阵的乘积。设上三角矩阵为R，则下三角矩阵为其转置，即A=R'R。MATLAB函数chol(X)用于对矩阵A进行Cholesky分解，其调用格式为：

R=chol(A)：产生一个上三角阵R，使R'R=A。若X为非对称正定，则输出一个出错信息。

[R,p]=chol(A)：这个命令格式将不输出出错信息。当A为对称正定的，则p=0，R与上述格式得到的结果相同；否则p为一个正整数。如果A为满秩矩阵，则R为一个阶数为q=p-1的上三角阵，且满足R'R=X(1:q,1:q)。

实现Cholesky分解后，线性方程组Ax=b变成R’Rx=b，所以x=R(R’b)。

例7-4用Cholesky分解求解例7-1中的线性方程组。

命令如下：

A=[2,1,-5,1;1,-5,0,7;0,2,1,-1;1,6,-1,-4];

b=[13,-9,6,0]';

R=chol(A)

??? Error using ==> chol

Matrix must be positive definite

命令执行时，出现错误信息，说明A为非正定矩阵。

9迭代解法

迭代解法非常适合求解大型系数矩阵的方程组。在数值分析中，迭代解法主要包括Jacobi迭代法、Gauss-Serdel迭代法、超松弛迭代法和两步迭代法。

1．Jacobi迭代法

对于线性方程组Ax=b，如果A为非奇异方阵，即aii≠0(i=1,2,…,n)，则可将A分解为A=D-L-U，其中D为对角阵，其元素为A的对角元素，L与U为A的下三角阵和上三角阵，于是Ax=b化为：

x=D-1(L+U)x+D-1b

与之对应的迭代公式为：

x(k+1)=D-1(L+U)x(k)+D-1b

这就是Jacobi迭代公式。如果序列{x(k+1)}收敛于x，则x必是方程Ax=b的解。

Jacobi迭代法的MATLAB函数文件Jacobi.m如下：

function [y,n]=jacobi(A,b,x0,eps)

if nargin==3

eps=1.0e-6;

elseif nargin<3

error

return

end

D=diag(diag(A));%求A的对角矩阵

L=-tril(A,-1);%求A的下三角阵

U=-triu(A,1);%求A的上三角阵

B=D(L+U);

f=Db;

y=B*x0+f;

n=1;%迭代次数

while norm(y-x0)>=eps

x0=y;

y=B*x0+f;

n=n+1;

end

例7-5用Jacobi迭代法求解下列线性方程组。设迭代初值为0，迭代精度为10-6。

在命令中调用函数文件Jacobi.m，命令如下：

A=[10,-1,0;-1,10,-2;0,-2,10];

b=[9,7,6]';

[x,n]=jacobi(A,b,[0,0,0]',1.0e-6)

2．Gauss-Serdel迭代法

在Jacobi迭代过程中，计算时，x(k)已经得到，不必再用，即原来的迭代公式Dx(k+1)=(L+U)x(k)+b可以改进为Dx(k+1)=Lx(k+1)+Ux(k)+b，于是得到：

x(k+1)=(D-L)-1Ux(k)+(D-L)-1b

该式即为Gauss-Serdel迭代公式。和Jacobi迭代相比，Gauss-Serdel迭代用新分量代替旧分量，精度会高些。

Gauss-Serdel迭代法的MATLAB函数文件gauseidel.m如下：

function [y,n]=gauseidel(A,b,x0,eps)

if nargin==3

eps=1.0e-6;

elseif nargin<3

error

return

end

D=diag(diag(A));%求A的对角矩阵

L=-tril(A,-1);%求A的下三角阵

U=-triu(A,1);%求A的上三角阵

G=(D-L)U;

f=(D-L)b;

y=G*x0+f;

n=1;%迭代次数

while norm(y-x0)>=eps

x0=y;

y=G*x0+f;

n=n+1;

end

例7-6用Gauss-Serdel迭代法求解下列线性方程组。设迭代初值为0，迭代精度为10-6。

在命令中调用函数文件gauseidel.m，命令如下：

A=[10,-1,0;-1,10,-2;0,-2,10];

b=[9,7,6]';

[x,n]=gauseidel(A,b,[0,0,0]',1.0e-6)

例7-7分别用Jacobi迭代和Gauss-Serdel迭代法求解下列线性方程组，看是否收敛。

命令如下：

a=[1,2,-2;1,1,1;2,2,1];

b=[9;7;6];

[x,n]=jacobi(a,b,[0;0;0])

[x,n]=gauseidel(a,b,[0;0;0])

7.2非线性方程数值求解

7.2.1单变量非线性方程求解

在MATLAB中提供了一个fzero函数，可以用来求单变量非线性方程的根。该函数的调用格式为：

z=fzero('fname',x0,tol,trace)

其中fname是待求根的函数文件名，x0为搜索的起点。一个函数可能有多个根，但fzero函数只给出离x0最近的那个根。tol控制结果的相对精度，缺省时取tol=eps，trace指定迭代信息是否在运算中显示，为1时显示，为0时不显示，缺省时取trace=0。

例7-8求f(x)=x-10x+2=0在x0=0.5附近的根。

步骤如下：

(1)建立函数文件funx.m。

function fx=funx(x)

fx=x-10.^x+2;

(2)调用fzero函数求根。

z=fzero('funx',0.5)

z =

0.3758

7.2.2非线性方程组的求解

对于非线性方程组F(X)=0，用fsolve函数求其数值解。fsolve函数的调用格式为：

X=fsolve('fun',X0,option)

其中X为返回的解，fun是用于定义需求解的非线性方程组的函数文件名，X0是求根过程的初值，option为最优化工具箱的选项设定。最优化工具箱提供了20多个选项，用户可以使用optimset命令将它们显示出来。如果想改变其中某个选项，则可以调用optimset()函数来完成。例如，Display选项决定函数调用时中间结果的显示方式，其中‘off’为不显示，‘iter’表示每步都显示，‘final’只显示最终结果。optimset(‘Display’,‘off’)将设定Display选项为‘off’。

例7-9求下列非线性方程组在(0.5,0.5)附近的数值解。

(1)建立函数文件myfun.m。

function q=myfun(p)

x=p(1);

y=p(2);

q(1)=x-0.6*sin(x)-0.3*cos(y);

q(2)=y-0.6*cos(x)+0.3*sin(y);

(2)在给定的初值x0=0.5,y0=0.5下，调用fsolve函数求方程的根。

x=fsolve('myfun',[0.5,0.5]',optimset('Display','off'))

x =

0.6354

0.3734

将求得的解代回原方程，可以检验结果是否正确，命令如下：

q=myfun(x)

q =

1.0e-009 *

0.23750.2957

可见得到了较高精度的结果。

10.函数极值

MATLAB提供了基于单纯形算法求解函数极值的函数fmin和fmins，它们分别用于单变量函数和多变量函数的最小值，其调用格式为：

x=fmin('fname',x1,x2)

x=fmins('fname',x0)

这两个函数的调用格式相似。其中fmin函数用于求单变量函数的最小值点。fname是被最小化的目标函数名，x1和x2限定自变量的取值范围。fmins函数用于求多变量函数的最小值点，x0是求解的初始值向量。

例7-13求f(x)=x3-2x-5在[0,5]内的最小值点。

(1)建立函数文件mymin.m。

function fx=mymin(x)

fx=x.^3-2*x-5;

(2)调用fmin函数求最小值点。

x=fmin('mymin',0,5)

x=

0.8165

11函数积分问题

1)

2)被积函数由一个表格定义

在MATLAB中，对由表格形式定义的函数关系的求定积分问题用trapz(X,Y)函数。其中向量X,Y定义函数关系Y=f(X)。

例8-4用trapz函数计算定积分。

命令如下：

X=1:0.01:2.5;

Y=exp(-X);%生成函数关系数据向量

trapz(X,Y)

ans =

0.28579682416393

3) 二重定积分的数值求解

使用MATLAB提供的dblquad函数就可以直接求出上述二重定积分的数值解。该函数的调用格式为：

I=dblquad(f,a,b,c,d,tol,trace)

该函数求f(x,y)在[a,b]×[c,d]区域上的二重定积分。参数tol，trace的用法与函数quad完全相同。

例8-5计算二重定积分

(1)建立一个函数文件fxy.m：

function f=fxy(x,y)

global ki;

ki=ki+1;%ki用于统计被积函数的调用次数

f=exp(-x.^2/2).*sin(x.^2+y);

(2)调用dblquad函数求解。

global ki;ki=0;

I=dblquad('fxy',-2,2,-1,1)

ki

I =

1.57449318974494

ki =

1038

数之差分与差商？？常微分方程初值问题的数值解法??傅立叶变换？？

12符号运算

(1)．符号表达式的四则运算

符号表达式的加、减、乘、除运算可分别由函数symadd、symsub、symmul和symdiv来实现，幂运算可以由sympow来实现。

(2)．符号表达式的提取分子和分母运算

如果符号表达式是一个有理分式或可以展开为有理分式，可利用numden函数来提取符号表达式中的分子或分母。其一般调用格式为：

[n,d]=numden(s)

该函数提取符号表达式s的分子和分母，分别将它们存放在n与d中。

(3)．符号表达式的因式分解与展开

MATLAB提供了符号表达式的因式分解与展开的函数，函数的调用格式为：

factor(s)：对符号表达式s分解因式。

expand(s)：对符号表达式s进行展开。

collect(s)：对符号表达式s合并同类项。

collect(s,v)：对符号表达式s按变量v合并同类项。

(4)．符号表达式的化简

MATLAB提供的对符号表达式化简的函数有：

simplify(s)：应用函数规则对s进行化简。

simple(s)：调用MATLAB的其他函数对表达式进行综合化简，并显示化简过程。

(5)．符号表达式与数值表达式之间的转换

利用函数sym可以将数值表达式变换成它的符号表达式。

函数numeric或eval可以将符号表达式变换成数值表达式。

(6).符号表达式中变量的确定

MATLAB中的符号可以表示符号变量和符号常量。findsym可以帮助用户查找一个符号表达式中的的符号变量。该函数的调用格式为：

findsym(s,n)

函数返回符号表达式s中的n个符号变量，若没有指定n，则返回s中的全部符号变量。

(7)级数符号求和

求无穷级数的和需要符号表达式求和函数symsum，其调用格式为：

symsum(s,v,n,m)

其中s表示一个级数的通项，是一个符号表达式。v是求和变量，v省略时使用系统的默认变量。n和m是求和的开始项和末项。

(8)函数的泰勒级数

MATLAB提供了taylor函数将函数展开为幂级数，其调用格式为：

taylor(f,v,n,a)

该函数将函数f按变量v展开为泰勒级数，展开到第n项(即变量v的n-1次幂)为止，n的缺省值为6。v的缺省值与diff函数相同。参数a指定将函数f在自变量v=a处展开，a的缺省值是0。

(9)符号代数方程求解

在MATLAB中，求解用符号表达式表示的代数方程可由函数solve实现，其调用格式为：

solve(s)：求解符号表达式s的代数方程，求解变量为默认变量。

solve(s,v)：求解符号表达式s的代数方程，求解变量为v。

solve(s1,s2,…,sn,v1,v2,…,vn)：求解符号表达式s1,s2,…,sn组成的代数方程组，求解变量分别v1,v2,…,vn。

(10)符号常微分方程求解

在MATLAB中，用大写字母D表示导数。例如，Dy表示y'，D2y表示y''，Dy(0)=5表示y'(0)=5。D3y+D2y+Dy-x+5=0表示微分方程y'''+y''+y'-x+5=0。

符号常微分方程求解可以通过函数dsolve来实现，其调用格式为：

dsolve(e,c,v)该函数求解常微分方程e在初值条件c下的特解。参数v描述方程中的自变量，省略时按缺省原则处理，若没有给出初值条件c，则求方程的通解。

dsolve在求常微分方程组时的调用格式为：

dsolve(e1,e2,…,en,c1,…,cn,v1,…,vn)

该函数求解常微分方程组e1,…,en在初值条件c1,…,cn下的特解，若不给出初值条件，则求方程组的通解，v1,…,vn给出求解变量。

13．S函数的设计与应用

S函数称为系统函数(System Function)，它有固定的程序格式。用MATLAB语言可以编写S函数，此外还可以采用C、C++、FORTRAN和Ada等语言编写。

(1) 用MATLAB语言编写S函数

编写S函数有一套固定的规则，为此，Simulink提供了一个用M文件编写S函数的模板。该模板程序存放在toolboxsimulinkblocks目录下，文件名为sfuntmpl.m。用户可以从这个模板出发构建自己的S函数。

(2)主程序

S函数主程序的引导语句为：

function [sys,x0,str,ts]=fname(t, x, u, flag)(3)子程序

S函数M文件共有6个子程序，供Simulink在仿真的不同阶段调用。

MATLAB中的运算符和特殊字符说明

符号

符号用途说明

+

加

-

减

.＊

点乘详细说明help arith

＊

矩阵相乘

＾

矩阵求幂

.＾

点幂

左除 详细说明help slash

/

右除

.

点左除

./

点右除

kron

张量积 详细说明help kron

，

作分隔用，如把矩阵元素、向量参数、函数参数、几个表达式分隔开来

；

(a)写在一个表达式后面时，运算后命令窗口中不显示表达式的计算结果

(b)在创建矩阵的语句中指示一行元素的结束，例如m=[x y z;i j k]

：

(a)创建向量的表达式分隔符，如x=a:b:c (b)a(:,j)表示j列的所有行元素；a(i,:)表示i行的所有列元素；a(1:3,4)表示第四列的第1行至第3行元素

()

圆括号

［］

创建数组、向量、矩阵或字符串(字母型)

｛｝

创建单元矩阵(cell array)或结构(struct)

%

注释符,特别当编写自定义函数文件时，紧跟function后的注释语句，在你使用help函数名时会显示出来。

'

(a)定义字符串用(b)向量或矩阵的共轭转置符

.'

一般转置符

...

表示MATLAB表达式继续到下一行，增强代码可读性

＝

赋值符号

＝＝

等于关系运算符

＜,＞

小于，大于关系运算符 详细说明help relop

&

逻辑与

|

逻辑或

~

逻辑非

xor

逻辑异或

unique寻找集合中互异元素(去掉相同元素)

find查找非零、非NaN元素的索引值

union集合 并

intersect集合 交

setdiff集合 差

setxor集合 异或

小整理：MATLAB常用的基本数学函数

abs(x)：纯量的绝对值或向量的长度

angle(z)：复数z的相角(Phase angle)

sqrt(x)：开平方

real(z)：复数z的实部

imag(z)：复数z的虚部

conj(z)：复数z的共轭复数

round(x)：四舍五入至最近整数

fix(x)：无论正负，舍去小数至最近整数(返回沿0方向最接近的整数)

floor(x)：地板函数，即舍去正小数至最近整数(返回沿负无穷大方向最接近的整数)

ceil(x)：天花板函数，即加入正小数至最近整数(返回沿正无穷大方向最接近的整数)

rats(x)：将实数x化为分数表示

rat(x)：将实数x化为多项分数展开

sign(x)：符号函数(Signum function)。

当x<0时，sign(x)=-1；

当x=0时，sign(x)=0;

当x>0时，sign(x)=1。

rem(x,y)：求x除以y的余数

gcd (x,y)：整数x和y的最大公因数

lcm(x,y)：整数x和y的最小公倍数

exp(x)：自然指数

pow2(x)：2的指数

log(x)：以e为底的对数，即自然对数

log2(x)：以2为底的对数

log10(x)：以10为底的对数

小整理：MATLAB常用的三角函数

sin(x)：正弦函数cos(x)：馀弦函数

tan(x)：正切函数asin(x)：反正弦函数

acos(x)：反馀弦函数atan(x)：反正切函数

atan2(x,y)：四象限的反正切函数

sinh(x)：超越正弦函数cosh(x)：超越馀弦函数

tanh(x)：超越正切函数asinh(x)：反超越正弦函数

acosh(x)：反超越馀弦函数atanh(x)：反超越正切函数。

小整理：适用于向量的常用函数

min(x):向量(数组)x的元素的最小值

max(x):向量x的元素的最大值

mean(x):向量x的元素的平均值

median(x):向量x的元素的中位数

std(x):向量x的元素的标准差

diff(x):向量x的相邻元素的差

sort(x):对向量x的元素按升序进行排序

length(x):向量x的元素个数(确定一个数组的最大维)

size(x):确定数组的行数和列数

norm(x):向量x的欧氏(Euclidean)长度

sum(x):向量x的元素总和

prod(x):向量x的元素总乘积

cumsum(x):向量x的累计元素总和

cumprod(x):向量x的累计元素总乘积

dot(x, y):向量x和y的内积

cross(x, y):向量x和y的外积

下表即为MATLAB常用到的永久常数。

小整理：MATLAB的永久常数

i或j：基本虚数单位(即)

eps：系统的浮点(Floating-point)精确度

inf：无限大， 例如1/0

nan或NaN：非数值(Not a number)，例如0/0

pi：圆周率p(= 3.1415926...)

realmax：系统所能表示的最大数值

realmin：系统所能表示的最小数值

nargin:函数的输入引数个数

nargout:函数的输出引数个数

x=0:pi/10:10*pi;

y=sin(x);

h1=line(x,y,'Erasemode','Xor');

h2=line(x,y,'Erasemode','Xor');

h10=0.01;h11=0;

h20=-0.01;h21=0;

while 1

h11=h11+h10,y=sin(x+h11);

set(h1,'YData',y);

h21=h21+h20,y=sin(x+h21);

set(h2,'YData',y);drawnow;

end

axis:冻结当前坐标轴缩放比例，用于随后的绘图。

axis(x_min,x_max,y_min,y_max):指定坐标轴按给定范围进行缩放。

text(x_coordinate,y_coordinate,’string’):用于向图形添加一个文本框，该文本框放在指定的x和y坐标处，并包含指定的字符串。

loglog:产生一个x-y坐标图，两个坐标轴都采用对数坐标。(semilogxsemilogy)

polar:产生一个极坐标图

figure(n):打开一个名为figure(n)的图形窗口，在该窗口中绘制图形。

eps :可识别的最小差值

erf:误差函数