MATLAB作为线性系统的一种分析和仿真工具，是理工科大学生应该掌握的技术工具，它作为一种编程语言和可视化工具，可解决工程、科学计算和数学学科中许多问题。
MATLAB建立在向量、数组和矩阵的基础上，使用方便，人机界面直观，输出结果可视化。其中，矩阵是MATLAB的核心！
下面将从基本规则和操作，编程和作图以及文件的操作等等多个方面，来讲解数学建模中MATLAB的一些常用方法。
一、
变量、函数与编程所形成的m文件是MATLAB操作的基本，在介绍它们的具体使用方法之前，先给出一些必须了解的基本规则。
1、 变量
MATLAB和其他编程工具一样，变量是必须的基本元素，它也是以字母开头，后接字母、数字或下划线的字符序列，用法也基本一样。其具体的命名规则是：
（1）变量名必须是不含空格的单个词；
（2）变量名区分大小写；
（3）变量名最多不超过19个字符；
（4）变量名必须以字母打头，之后可以是任意字母、数字或下划线，变量名中不允许使用标点符号。
特殊变量表
注意，自定义的变量不能和表中变量名相同。
2、数学运算符号及标点符号
（1）MATLAB的每条命令后，若为逗号或无标点符号，则显示命令的结果；若命令后为分号，则禁止显示结果；
（2）“%”后面所有文字为注释；
（3）“...”表示续行。
3、数学函数
基本上，常用的数学函数在MATLAB中都有相应的命令，部分如下：
4、M文件
MATLAB的内部函数是有限的，有时为了研究某一个函数的各种性态，需要为MATLAB定义新函数，为此必须编写函数文件. 函数文件是文件名后缀为M的文件，这类文件的第一行必须是一特殊字符function开始，格式为：
函数值的获得必须通过具体的运算实现，并赋给因变量。
M文件建立方法：
1. 在Matlab中，点:File->New->M-file
Matlab的应用程序也以M文件保存。
例：定义函数 f(x1,x2)=100(x2-x1 ) +(1-x1)
1.建立M文件：fun.m
function
f=100*(x(2)-x(1)^2)^2+(1-x(1))^2
2. 可以直接使用函数fun.m
例如：计算 f(1,2), 只需在Matlab命令窗口键入命令：
x=[1 2]
fun(x)
二、数组与矩阵
矩阵是MATLAB最基本的数据对象，MATLAB的大部分运算或命令都是在矩阵运算的意义下执行的。在MATLAB中，不需对矩阵的维数和类型进行说明，MATLAB会根据用户所输入的内容自动进行配置。
1、数
（1）创建简单的数组
x=[a b c d e f ]
x=first：last
x=first：increment：last
创建从first开始，加increment计数，last结束的行向量
x=linspace(first，last，n）
创建从first开始，到last结束，有n个元素的行向量
x=logspace(first，last，n）
创建从开始，到结束，有n个元素的对数分隔行向量。
（2）数组元素的访问
（i）访问一个元素：x(i)表示访问数组x的第i个元素；
（ii）访问一块元素：x(a ：b ：c)表示访问数组x的从第a个元素开始，以步长为b到第c个元素（但不超过c），b可以为负数，b缺损时为1；
（iii）直接使用元素编址序号：x([a
（3）数组的方向
产生列向量有两种方法：
（4）数组的运算
（i）标量-数组运算
设：a=[a1,a2,…,an], c=标量
则：a+c=[a1+c,a2+c,…,an+c]
（ii）数组-数组运算
设：a=[a1,a2,…,an], b=[b1,b2,…,bn]
则：a+b= [a1+b1,a2+b2,…,an+bn]
2、矩阵
（1）矩阵的建立
逗号或空格用于分隔某一行的元素，分号用于区分不同的行。除了分号，在输入矩阵时，按Enter键也表示开始一新行。输入矩阵时，严格要求所有行有相同的列。
例
特殊矩阵的建立：
a=[
b=zeros(m，n)
c=ones(m，n)
d=eye(m，n)
（2）矩阵中元素的操作
（a）矩阵A的第r行：A（r，：）
（b）矩阵A的第r列：A（：，r）
（c）依次提取矩阵A的每一列，将A拉伸为一个列向量：A（：）
（d）取矩阵A的第i1~i2行、第j1~j2列构成新矩阵:A(i1:i2, j1:j2)
（e）以逆序提取矩阵A的第i1~i2行，构成新矩阵:A(i2:-1：i1，：）
（f）以逆序提取矩阵A的第j1~j2列，构成新矩阵:A(:, j2:-1：j1 ）
（g）删除A的第i1~i2行，构成新矩阵:A(i1:i2，：)=[ ]
（h）删除A的第j1~j2列，构成新矩阵:A(：，j1:j2)=[ ]
（i）将矩阵A和B拼接成新矩阵：[A
（3）矩阵的运算
（i）标量-矩阵运算：同标量-数组运算；
（ii）矩阵-矩阵运算：
[2] 矩阵运算：
矩阵加法：A+B
矩阵乘法：A*B
方阵的行列式：det（A）
方阵的逆：inv（A）
方阵的特征值与特征向量：[V，D]=eig[A]。
三、 MATLAB编程
Matlab虽然提供了大量现成的函数供我们计算时使用，但当要真正解决一些实际问题时，这还远远不够，编程仍然是不可或缺的一步。Matlab也可以象C，FORTRAN等高级计算机语言一样，进行程序设计。下面简单介绍Matlab中一些重要的编程手段。
1、关系与逻辑运算
除了传统的数学运算，Matlab支持关系和逻辑运算。一个重要的应用是控制基于真/假命题的一系列Matlab命令的流程，或执行次序。
（1）关系操作符
（2）逻辑运算符
（3）控制流
MATLAB提供三种决策或控制流结构：for循环、while循环、if-else-end结构。
i、for循环：允许一组命令以固定的和预定的次数重复
ii、While循环
while
iii、If-Else-End结构
（a）有一个选择的一般形式是：
先建立M文件fun1.m定义函数f（x），再在Matlab命令窗口输入fun1(2),fun1(-1)即可。
（b） 有三个或更多的选择的一般形式是：
{commands}
end
end
end
……
先建立M文件fun2.m定义函数f（x），再在Matlab命令窗口输入fun2(2),fun2(0.5), fun2(-1)即可。
四、MATLAB作图
强大的图形功能是Matlab的优点之一，它能方便快速的出图，给我们的工作带来了巨大的便利。
1. 曲线图
Matlab作图是通过描点、连线来实现的，故在画一个曲线图形之前，必须先取得该图形上的一系列的点的坐标（即横坐标和纵坐标），然后将该点集的坐标传给Matlab函数画图。
命令为：
PLOT(X,Y,S)
其中，X,Y是向量，分别表示点集的横坐标和纵坐标，
S指定曲线的颜色、线形等：
y：黄色
.：点
x：x-符号
PLOT(X,Y)——画实线
PLOT(X,Y1,S1,X,Y2,S2,……,X,Yn,Sn)——将多条线画在一起
例 在[0,2*pi]用红线画sin(x),用绿圈画cos(x)。
解：x=linspace(0,2*pi,30);
y=sin(x);
z=cos(x);
plot(x,y,'r',x,z,’g0')
2.符号函数(显函数、隐函数和参数方程)画图
(1) ezplot
ezplot(‘f(x)’,[a,b])
ezplot(‘f(x,y)’,[xmin,xmax,ymin,ymax])
隐函数f(x,y)=0的函数图
ezplot(‘x(t)’,’y(t)’,[tmin,tmax])
例 在[0,pi]上画y=cos(x)的图形；
解：输入命令
例 在[0,2*pi]上画 ， 星形图；
解：输入命令
例 在[-2，0.5]，[0，2]上画隐函数 的图。
解：输入命令
ezplot('exp(x)+sin(x*y)',[-2,0.5,0,2])
(2) fplot
fplot(‘fun’,lims)
注意：
[1] fun必须是M文件的函数名或是独立变量为x的字符串；
[2] fplot函数不能画参数方程和隐函数图形，但在一个图上可以画多个图形。
例 在[-1，2]上画 的 图形
解：先建M文件myfun1.m
再输入命令：
fplot(‘myfun1’,[-1,2])
例 在[-2,2]范围内绘制函数tanh的图形
解：fplot(‘tanh’,[-2,2])
例 x、y的取值范围都在[- ， ]，画函数tanh(x),sin(x),cos(x)的图形
解：输入命令
3、空间曲线
(1)一条曲线：PLOT3(x,y,z,s)
其中，X，Y，Z是n维向量，分别表示曲线上点集的横坐标、纵坐标、函数值
S则指定颜色、线形等
例 在区间[0，10*pi]画出参数曲线x=sin(t)，y=cos(t)，z=t
解：t=0:pi/50:10*pi;
(2)多条曲线：PLOT3(x,y,z)
其中x，y，z是都是m*n矩阵，其对应的每一列表示一条曲线.
例 画多条曲线观察函数Z=(X+Y).^2.
解：x=-3:0.1:3;y=1:0.1:5;
（这里meshgrid(x,y)的作用是产生一个以向量x为行、向量y为列的矩阵）
4、空间曲面
(1)surf(x,y,z) 画出数据点（x，y，z）表示的曲面数据矩阵。分别表示数据点的横坐标、纵
坐标、函数值
例
解：x=-3:0.1:3;
(2)Mesh(x,y,z) 画网格曲面
例
解：x=-3:0.1:3;
(3)meshz(X,Y,Z) 在网格周围画一个curtain图(如,参考平面)
例
解：输入命令:
5、特殊二、三维图形
二维：
1、极坐标图：polar (theta,rho,s)
用角度theta（弧度表示）和极半径rho作极坐标图，用s指定线型。
2、散点图：scatter（X,Y,S,C）
3、平面等值线图：contour (x,y,z,n) 绘制n个等值线的二维等值线图
三维：
1、空间等值线图：contour 3(x,y,z,n) 其中n表示等值线数。
2、三维散点图：scatter3（X,Y,Z,S,C）
五、文件操作
Matlab提供了对数据文件建立、打开、读、写以及关闭等一系列函数，数据文件一般存放在磁盘等介质上，用文件名标识，系统对文件名没有特殊要求。文件数据格式有二种形式，一是二进制格式文件，二是ASCII文本文件，系统对这两类文件提供了不同的读写功能函数。
1、文件的打开与关闭
（1）fopen打开文件
在读写文件之前，必须先用fopen命令打开一个文件，并指定允许对该文件进行的操作。文件操作结束后，应及时关闭文件，以免数据的丢失或误修改。fopen函数格式为：
其中filename为文件名，permission为文件格式，可以是下列格式之一：
•
•
•
•
•
•
•
•
例如：打开一个名为std.dat的数据文件并进行读操作，其命令格式为：
Fid=fopen( ‘std.dat’, ’r’ )
上述打开格式均为二进制格式，如果想用ASCII文本格式，则必须在格式字符串中加上字符t，例如用’r t’表示以ASCII格式打开供读操作的数据文件。
（2）fclose关闭文件
文件在进行完读、写等操作后，应及时关闭，以保证文件的安全可靠。关闭文件命令格式为：
Sta=fclose(Fid)
Sta表示关闭文件操作的返回代码，若关闭成功，返回0，否则返回–1。
2、文件的读写操作
（1）二进制数据文件
fread 读二进制数据文件。格式为：
[A,COUNT]=fread(Fid,size,precision)
其中A为数据矩阵，COUNT返回所读取的数据元素个数。size为可选项，若不选用则读取整个文件内容，若选用它的值可以是下列值：
N
inf
[M,N]
precision用于控制所读数据的精度格式。
缺省格式为uchar，即无符号字符格式。例如：Fid=fopen(‘std.dat’, ‘r’);
以读数据方式打开数据文件std.dat，并按长整型数据格式读取文件的前100个数据放入向量A，然后关闭文件。
fwrite 函数以二进制格式向数据文件写数据，其格式为：
COUNT=fwrite (Fid, A, precision)
例如：Fid=fopen(‘magic5.bin’, ‘wb’);
fwrite(Fid, magic, ’int32’);
上述语句将矩阵magic中的数据写入文件magic5.bin中，数据格式为32位整型二进制格式。
例 建立一数据文件test.dat，用于存放矩阵A的数据。
已知
Fid=fopen(‘test.dat’, ‘w’)
cnt=fwrite(Fid, A, ‘float’)
fclose(Fid)
程序段将矩阵A的数据以二进制浮点数格式写入文件test.dat中。
Fid=fopen(‘test.dat’, ‘r’)
[B,cnt]=fread(Fid, [5,inf], ‘float’)
fclose(Fid)
读取文件test.dat的内容。
（2）文本文件
fscanf
[A,COUNT]= fscanf (Fid, format, size)
其中A为数据矩阵，用以存放读取的数据，COUNT返回所读取的数据元素个数。
format用以控制读取的数据格式，由%加上格式符组成，格式符为：
d, i, o, u, x, e, f, g, s, c与[. . .]
例如：
s=fscanf(fid, ‘%s’)
a=fscanf(fid, ‘]’)
b= fscanf(fid, ‘%6.2d’)
fprintf
其中A为要写入文件的数据矩阵，先按format格式化数据矩阵A，后写入到Fid所指定的文件。
例如：x = 0: 0.1: 1;
六、数据处理与多项式计算
求一组数据的最值，中值，累加和等等简单的数据处理和多项式求根等相关运算是我们最常遇见的一类问题，下面介绍如何使用MATLAB来解决这些问题。
1、最大值和最小值
MATLAB提供的求数据序列的最大值和最小值的函数分别为max和min，两个函数的调用格式和操作过程类似。
（1）求向量的最大值和最小值
求一个向量X的最大值的函数有两种调用格式，分别是：
(1) y=max(X)：返回向量X的最大值存入y，如果X中包含复数元素，则按模取最大值。
(2) [y,I]=max(X)：返回向量X的最大值存入y，最大值的序号存入I，如果X中包含复数元素，则按模取最大值。
求向量X的最小值的函数是min(X)，用法和max(X)完全相同。
例 求向量x的最大值。
命令如下：
x=[-43,72,9,16,23,47];
y=max(x)
[y,l]=max(x)
（2）求矩阵的最大值和最小值
求矩阵A的最大值的函数有3种调用格式，分别是：
(1) max(A)：返回一个行向量，向量的第i个元素是矩阵A的第i列上的最大值。
(2) [Y,U]=max(A)：返回行向量Y和U，Y向量记录A的每列的最大值，U向量记录每列最大值的行号。
(3) max(A,[],dim)：dim取1或2。dim取1时，该函数和max(A)完全相同；dim取2时，该函数返回一个列向量，其第i个元素是A矩阵的第i行上的最大值。
求最小值的函数是min，其用法和max完全相同。
（3）两个向量或矩阵对应元素的比较
函数max和min还能对两个同型的向量或矩阵进行比较，调用格式为：
(1) U=max(A,B)：A,B是两个同型的向量或矩阵，结果U是与A,B同型的向量或矩阵，U的
每个元素等于A,B对应元素的较大者。
(2) U=max(A,n)：n是一个标量，结果U是与A同型的向量或矩阵，U的每个元素等于A对应元素和n中的较大者。
min函数的用法和max完全相同。
2、求和与求积
数据序列求和与求积的函数是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行的各元素乘积。
3、平均值和中值
求数据序列平均值的函数是mean，求数据序列中值的函数是median。两个函数的调用格式为：
mean(X)：返回向量X的算术平均值。
median(X)：返回向量X的中值。
mean(A)：返回一个行向量，其第i个元素是A的第i列的算术平均值。
median(A)：返回一个行向量，其第i个元素是A的第i列的中值。
mean(A,dim)：当dim为1时，该函数等同于mean(A)；当dim为2时，返回一个列向量，其第i个元素是A的第i行的算术平均值。
median(A,dim)：当dim为1时，该函数等同于median(A)；当dim为2时，返回一个列向量，其第i个元素是A的第i行的中值。
[median] 是在一组数据中居于中间的数(特别注意的地方是：这组数据之前已经经过升序排列！！！)，即在这组数据中，有一半的数据比它大，有一半的数据比它小。如果这组数据包含偶数个数字，中值是位于中间的两个数的平均值。
4、累加和与累乘积
在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行的累乘积向量。
5、标准方差与相关系数
（1）求标准方差
在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。
（2）相关系数
MATLAB提供了corrcoef函数，可以求出数据的相关系数矩阵。corrcoef函数的调用格式为：
corrcoef(X)：返回从矩阵X形成的一个相关系数矩阵。此相关系数矩阵的大小与矩阵X一样。它把矩阵X的每列作为一个变量，然后求它们的相关系数。
corrcoef(X,Y)：在这里，X,Y是向量，它们与corrcoef([X,Y])的作用一样。
例 生成满足正态分布的10000×5随机矩阵，然后求各列元素的均值和标准方差，再求这5列随机数据的相关系数矩阵。
命令如下：
X=randn(10000,5);
M=mean(X)
D=std(X)
R=corrcoef(X)
6、排序
MATLAB中对向量X是排序函数是sort(X)，函数返回一个对X中的元素按升序排列的新向量。sort函数也可以对矩阵A的各列或各行重新排序，其调用格式为：
[Y,I]=sort(A,dim)：其中dim指明对A的列还是行进行排序。若dim=1，则按列排；若dim=2，则按行排。Y是排序后的矩阵，而I记录Y中的元素在A中位置。
7、多项式计算
MATLAB规定的多项式的一般形式为：
并用向量[ ]来表示。
（1）多项式的四则运算
（i）多项式加、减
对于次数相同的若干个多项式，可直接对多项式系数向量进行加、减的运算。如果多项式的次数不同，则应该把低次的多项式系数不足的高次项用零补足，使同式中的各多项式具有相同的次数。
（ii）多项式乘法
若A、B是由多项式系数组成的向量，则CONV函数将返回这两个多项式的乘积。调用它的命令格式为：
C=conv(A,B)
命令的结果C为一个向量，由它构成一个多项式。
例 求多项式 与 的乘积。
A=[1 8 0 0 -10];
B=[2 -1 3]
B =
C=conv(A,B)
C =
当A、B是由多项式系数组成的向量时，DECONV函数用来对两个多项式作除法运算。调用的命令格式为：
[Q,r]=deconv(A,B)
本命令的结果：多项式A除以多项式B获商多项式赋予Q（也为多项式系数向量）；获余项多项式赋予r（其系数向量的长度与被除多项式相同，通常高次项的系数为0）。
DECONV是CONV的逆函数，即有A=conv(B,Q)+r。
例 试求 与 相除。
A=[1 8 0 0 -10];
B=[2 -1 3];
[P,r]=deconv(A,B)
P =
r =
商多项式P为0.5x2+4.25x+1.375，
余项多项式r为 -11.375x-14.125。
（2）多项式求根
n次多项式具有n个根，当然这些根可能是实根，也可能含有若干对共轭复根。MATLAB提供的roots函数用于求多项式的全部根，其调用格式为：
x=roots(P)
其中P为多项式的系数向量，求得的根赋给向量x，即x(1),x(2),…,x(n)分别代表多项式的n个根。
例 求多项式 的根。
命令如下：
A=[1,8,0,0,-10];
x=roots(A)
若已知多项式的全部根，则可以用poly函数建立起该多项式，其调用格式为：
P=poly(x)
若x为具有n个元素的向量，则poly(x)建立以x为其根的多项式，且将该多项式的系数赋给向量P。
例6-22
(1) 计算f(x)=0 的全部根。
(2) 由方程f(x)=0的根构造一个多项式g(x)，并与f(x)进行对比。
命令如下：
P=[3,0,4,-5,-7.2,5];
X=roots(P)
G=poly(X)
（3）多项式的求值
MATLAB提供了两种求多项式值的函数：polyval与polyvalm，它们的输入参数均为多项式系数向量P和自变量x。两者的区别在于前者是代数多项式求值，而后者是矩阵多项式求值。
（i）代数多项式求值
polyval函数用来求代数多项式的值，其调用格式为：
Y=polyval(P,x)
若x为一数值，则求多项式在该点的值；若x为向量或矩阵，则对向量或矩阵中的每个元素求其多项式的值。
例 已知 ，分别取x=1.2和下面的矩阵的2×3个元素为自变量计算该多项式的值。
A=[1 8 0 0 -10];
x=1.2;
y1=polyval(A,x)
y1 =
x=[-1 1.2 -1.4;2 -1.8 1.6]
x =
（ii）矩阵多项式求值
polyvalm函数用来求矩阵多项式的值，其调用格式与polyval相同，但含义不同。polyvalm函数要求x为方阵，它以方阵为自变量求多项式的值。设A为方阵，P代表多项式 ，那么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))
（4）多项式的导函数
对多项式求导数的命令格式是：
p=polyder(P)：求多项式P的导函数
p=polyder(P,Q)：求P·Q的导函数
[p,q]=polyder(P,Q)：求P/Q的导函数，导函数的分子存入p，分母存入q。
上述函数中，参数P,Q是多项式的向量表示，结果p,q也是多项式的向量表示。
例 求有理分式 的导数。
命令如下：
P=[1];
Q=[1,0,5];
[p,q]=polyder(P,Q)
七、MATLAB解方程与函数极值
线性方程组与非线性方程的求解，以及函数极值的求取也是我们常常遇见，但不易解决的问题，这一部分将介绍如何使用MATLAB来解决这几种问题。
1、线性方程组求解
（1）直接解法
i．利用左除运算符的直接解法
对于线性方程组Ax=b，可以利用左除运算符“”求解：
例 用直接解法求解下列线性方程组。
命令如下：
A=[2,1,-5,1;1,-5,0,7;0,2,1,-1;1,6,-1,-4];
b=[13,-9,6,0]';
x=Ab
ii．利用矩阵的分解求解线性方程组
矩阵分解是指根据一定的原理用某种算法将一个矩阵分解成若干个矩阵的乘积。常见的矩阵分解有LU分解、QR分解、Cholesky分解，以及Schur分解、Hessenberg分解、奇异分解等。
a．LU分解
矩阵的LU分解就是将一个矩阵表示为一个交换下三角矩阵和一个上三角矩阵的乘积形式。线性代数中已经证明，只要方阵A是非奇异的，LU分解总是可以进行的。
MATLAB提供的lu函数用于对矩阵进行LU分解，其调用格式为：
[L,U]=lu(X)：产生一个上三角阵U和一个变换形式的下三角阵L(行交换)，使之满足X=LU。注意，这里的矩阵X必须是方阵。
[L,U,P]=lu(X)：产生一个上三角阵U和一个下三角阵L以及一个置换矩阵P，使之满足PX=LU。当然矩阵X同样必须是方阵。
实现LU分解后，线性方程组Ax=b的解x=U(Lb)或x=U(LPb)，这样可以大大提高运算速度。
例 用LU分解求解上例中的线性方程组。
命令如下：
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)
b．QR分解
对矩阵X进行QR分解，就是把X分解为一个正交矩阵Q和一个上三角矩阵R的乘积形式。QR分解只能对方阵进行。MATLAB的函数qr可用于对矩阵进行QR分解，其调用格式为：
[Q,R]=qr(X)：产生一个一个正交矩阵Q和一个上三角矩阵R，使之满足X=QR。
[Q,R,E]=qr(X)：产生一个一个正交矩阵Q、一个上三角矩阵R以及一个置换矩阵E，使之满足XE=QR。
实现QR分解后，线性方程组Ax=b的解x=R(Qb)或x=E(R(Qb))。
例 用QR分解求解上例中的线性方程组。
命令如下：
A=[2,1,-5,1;1,-5,0,7;0,2,1,-1;1,6,-1,-4];
b=[13,-9,6,0]';
[Q,R]=qr(A);
x=R(Qb)
或采用QR分解的第2种格式，命令如下：
[Q,R,E]=qr(A);
x=E*(R(Qb))
c．Cholesky分解
如果矩阵X是对称正定的，则Cholesky分解将矩阵X分解成一个下三角矩阵和上三角矩阵的乘积。设上三角矩阵为R，则下三角矩阵为其转置，即X=R'R。MATLAB函数chol(X)用于对矩阵X进行Cholesky分解，其调用格式为：
R=chol(X)：产生一个上三角阵R，使R'R=X。若X为非对称正定，则输出一个出错信息。
[R,p]=chol(X)：这个命令格式将不输出出错信息。当X为对称正定的，则p=0，R与上述格式得到的结果相同；否则p为一个正整数。如果X为满秩矩阵，则R为一个阶数为q=p-1的上三角阵，且满足R'R=X(1:q,1:q)。
实现Cholesky分解后，线性方程组Ax=b变成R‘Rx=b，所以x=R(R’b)。
例 用Cholesky分解求解上例中的线性方程组。
命令如下：
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为非正定矩阵。
（2）迭代解法
迭代解法非常适合求解大型系数矩阵的方程组。在数值分析中，迭代解法主要包括 Jacobi迭代法、Gauss-Serdel迭代法、超松弛迭代法和两步迭代法。
i．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
elseif nargin<3
end
D=diag(diag(A));
L=-tril(A,-1);
U=-triu(A,1);
B=D(L+U);
f=Db;
y=B*x0+f;
n=1;
while norm(y-x0)>=eps
end
例 用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)
ii．Gauss-Serdel迭代法
在Jacobi迭代过程中，计算时，已经得到，不必再用，即原来的迭代公式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
elseif nargin<3
end
D=diag(diag(A));
L=-tril(A,-1);
U=-triu(A,1);
G=(D-L)U;
f=(D-L)b;
y=G*x0+f;
n=1;
while norm(y-x0)>=eps
end
例 用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)
例 分别用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])
2、非线性方程数值求解
（1） 单变量非线性方程求解
其中fname是待求根的函数文件名，x0为搜索的起点。一个函数可能有多个根，但fzero函数只给出离x0最近的那个根。tol控制结果的相对精度，缺省时取tol=eps，trace指定迭代信息是否在运算中显示，为1时显示，为0时不显示，缺省时取trace=0。
例 求f(x)=x-10x+2=0在x0=0.5附近的根。
步骤如下：
(1) 建立函数文件funx.m。
(2) 调用fzero函数求根。
（2）非线性方程组的求解
其中X为返回的解，fun是用于定义需求解的非线性方程组的函数文件名，X0是求根过程的初值，option为最优化工具箱的选项设定。最优化工具箱提供了20多个选项，用户可以使用optimset命令将它们显示出来。如果想改变其中某个选项，则可以调用optimset()函数来完成。例如，Display选项决定函数调用时中间结果的显示方式，其中‘off’为不显示，‘iter’表示每步都显示，‘final’只显示最终结果。optimset(‘Display’,‘off’)将设定Display选项为‘off’。
例 求下列非线性方程组在(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 =
将求得的解代回原方程，可以检验结果是否正确，命令如下：
q=myfun(x)
q =1.0e-009 *
3、 函数极值
这两个函数的调用格式相似。其中fmin函数用于求单变量函数的最小值点。fname是被最小化的目标函数名，x1和x2限定自变量的取值范围。fmins函数用于求多变量函数的最小值点，x0是求解的初始值向量。
MATLAB没有专门提供求函数最大值的函数，但只要注意到-f(x)在区间(a,b)上的最小值就是f(x)在(a,b)的最大值，所以fmin(f,x1,x2)返回函数f(x)在区间(x1,x2)上的最大值。
例 求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=
八、插值与拟合
在工程中，常有这样的问题：给定一批数据点，需确定满足特定要求的曲线或曲面。如果要求所求曲线（面）通过所给所有数据点，这就是插值问题；在数据较少的情况下，这样做能取得较好的效果。但是，如果数据较多，那么，插值函数是一个次数很高的函数，比较复杂，同时，给定的数据一般是由观察测量所得，往往带有随机误差，因而，要求曲线（面）通过所有数据点就即不现实也不必要。如果不要求曲线（面）通过所有数据点，而是要求它反映对象整体的变化趋势，可得到更简单实用的近似函数，这就是数据拟合，又称曲线拟合或曲面拟合。函数插值与曲线拟合都是要根据一组数据构造一个函数作为近似，由于近似的要求不同，二者在数学方法上是完全不同的。
1、一维数据插值
在MATLAB中，实现这些插值的函数是interp1，其调用格式为：
Y1=interp1(X,Y,X1,'method')
函数根据X,Y的值，计算函数在X1处的值。X,Y是两个等长的已知向量，描述的是插值点，X1是一个向量或标量，描述被插值的点，Y1是一个与X1等长的插值结果。method是所要采用的插值方法，允许的取法有‘linear’ （线性插值）、‘nearest’（最邻近插值）、‘cubic’（立方插值）、‘spline’（三次样条插值）。
注意：所有的插值方法都要求x是单调的，并且xi不能够超过x的范围。
例 在1-12的11小时内，每隔1小时测量一次温度，测得的温度依次为：5，8，9，15，25，29，31，30，22，25，27，24。试估计每隔1/10小时的温度值。
hours=1:12;
temps=[5 8 9 15 25 29 31 30 22 25 27 24];
h=1:0.1:12;
t=interp1(hours,temps,h,'spline');
plot(hours,temps,'+',h,t,hours,temps,'r:')
xlabel('Hour'),ylabel('Degrees Celsius’)
2、二维数据插值
在MATLAB中，提供了解决二维插值问题的函数interp2，其调用格式为：
Z1=interp2(X,Y,Z,X1,Y1,'method')
其中X,Y是两个向量，分别描述两个参数的采样点，Z是与参数采样点对应的函数值，X1,Y1是两个向量或标量，描述被插值的点。Z1是根据相应的插值方法得到的插值结果。 method的取值与一维插值函数相同。X,Y,Z也可以是矩阵形式。
同样，要求x,y单调；x,y可取为矩阵，或x取行向量，y取为列向量，xi,yi的值分别不能超出x,y的范围。
3、散点数据的插值
在二维插值中，若插值点不单调，则应采用插值函数griddata，其调用格式为：
cz =griddata（X,Y,Z,X1,Y1,‘method’）
函数的用法与上相同。但要求X1取行向量，Y1取为列向量。
4、线性最小二乘拟合
在MATLAB中，用polyfit函数来求得最小二乘拟合多项式的系数，其调用格式为：
a=polyfit(x,y,m)
其中X，Y输入同长度的数组，m表示拟合多项式的次数。
例
| xi | 0.1 | 0.2 | 0.4 | 0.5 | 0.6 | 0.7 | 0.8 | 0.9 | 1 |
| yi | 1.978 | 3.28 | 6.16 | 7.34 | 7.66 | 9.58 | 9.48 | 9.30 | 11.2 |
输入以下命令：
x=0:0.1:1;
y=[-0.447 1.978 3.28 6.16 7.08 7.34 7.66 9.56 9.48 9.30 11.2];
A=polyfit(x,y,2)
z=polyval(A,x);
plot(x,y,'k+',x,z,'r')
5、非线性最小二乘拟合
Matlab提供了两个求非线性最小二乘拟合的函数：lsqcurvefit和lsqnonlin。两个命令都要先建立M-文件fun.m，在其中定义函数f(x)，但两者定义f(x)的方式是不同的。
（1）lsqcurvefit
已知数据点：
lsqcurvefit用以求含参量x（向量）的向量值函数F(x,xdata)=（F（x，xdata1），…，F（x，xdatan））T中的参变量x(向量) 。其调用格式为：
x = lsqcurvefit (‘fun’,x0,xdata,ydata,options);
其中fun是一个事先建立的定义函数F(x,xdata) 的M-文件, 自变量为x和xdata，x0为迭代初值，xdata,ydata为已知数据点，options选项见无约束优化
（2）lsqnonlin
已知数据点：xdata=（xdata1，xdata2，…，xdatan）
lsqnonlin用以求含参量x（向量）的向量值函数f(x)=(f1(x),f2(x),…,fn(x)) 中的参量x。其中 fi(x)=f(x，xdatai，ydatai）=F(x,xdatai)-ydatai。其调用格式为：
x= lsqnonlin （‘fun’，x0，options）
其中参变量含义同上。
例 用下面一组数据拟合 中的参数a，b，k
| t | 100 | 200 | 300 | 400 | 500 | 600 | 700 | 800 | 900 | 1000 |
| c | 4.54 | 4.99 | 5.35 | 5.65 | 5.90 | 6.10 | 6.26 | 6.39 | 6.50 | 6.59 |
解法1.
1）编写M-文件 curvefun1.m
function f=curvefun1(x,tdata)
f=x(1)+x(2)*exp(-0.02*x(3)*tdata)
2）输入命令
tdata=100:100:1000
cdata=1e-03*[4.54,4.99,5.35,5.65,5.90,6.10,6.26,6.39,
6.50,6.59];
x0=[0.2,0.05,0.05];
x=lsqcurvefit ('curvefun1',x0,tdata,cdata)
f= curvefun1(x,tdata)
3）运算结果为：
f =0.0043
x = 0.0063
4）结论：a=0.0063, b=-0.0034, k=0.2542
解法 2
1）编写M-文件 curvefun2.m
function f=curvefun2(x)
tdata=100:100:1000;
cdata=1e-03*[4.54,4.99,5.35,5.65,5.90,
f=x(1)+x(2)*exp(-0.02*x(3)*tdata)- cdata
2）输入命令:
x=lsqnonlin('curvefun2',x0)
f= curvefun2(x)
3）运算结果为
f =1.0e-003 *(0.2322
-0.1668
x =0.0063
4）结论：即拟合得a=0.0063
可以看出,两个命令的计算结果是相同的。
九、运筹与优化
优化问题，一般是指用“最好”的方式，使用或分配有限的资源，即劳动力、原材料、机器、资金等，使得费用最小或利润最大。这一部分主要讲述线性规划与非线性规划的Matlab方法。
1、 线性规划
模型：min z=cX
VLB≤X≤VUB
命令：[1] x=linprog（c，A，b，Aeq,beq, VLB，VUB）
注意：[1] 若没有等式约束, 则令Aeq=[ ], beq=[ ].
例
解:
2、 非线性规划
模型：min F(X)
G(X) 0
其中X为n维变元向量，G(X)与Ceq(X)均为非线性函数组成的向量，其它变量的含义与线性规划、二次规划中相同.用Matlab求解上述问题，基本步骤分三步：
1．首先建立M文件fun.m,定义目标函数F（X）:
function f=fun(X);
f=F(X);
2． 若约束条件中有非线性约束:G(X) 或Ceq(X)=0,则建立M文件nonlcon.m定义函数G(X)与Ceq(X):
function [G,Ceq]=nonlcon(X)
G=...
Ceq=...
3. 建立主程序.非线性规划求解的函数是fmincon,命令的基本格式:
x=fmincon(‘fun’,X0,A,b,Aeq,beq,VLB,VUB,’nonlcon’,options)
其中的参变量含义与线性规划中相同。
注意：
[1] fmincon函数提供了大型优化算法和中型优化算法。默认时，若在fun函数中提供了梯度（options参数的GradObj设置为’on’），并且只有上下界存在或只有等式约束，fmincon函数将选择大型算法。当既有等式约束又有梯度约束时，使用中型算法。
[2] fmincon函数的中型算法使用的是序列二次规划法。在每一步迭代中求解二次规划子问题，并用BFGS法更新拉格朗日Hessian矩阵。
[3] fmincon函数可能会给出局部最优解，这与初值X0的选取有关。
例 min
s.t.
1．先建立M文件 fun4.m,定义目标函数:
function f=fun4(x);
f=exp(x(1))*(4*x(1)^2+2*x(2)^2+4*x(1)*x(2)+2*x(2)+1);
2．再建立M文件mycon.m定义非线性约束：
function [g,ceq]=mycon(x)
g=[x(1)+x(2);1.5+x(1)*x(2)-x(1)-x(2);-x(1)*x(2)-10];
3．主程序youh3.m为:
x0=[-1;1];
A=[];b=[];
Aeq=[1 1];beq=[0];
vlb=[];vub=[];
[x,fval]=fmincon('fun4',x0,A,b,Aeq,beq,vlb,vub,'mycon')
4．运算结果为：
x = -1.2250
fval = 1.8951
十、多元回归
在实际问题中，我们常常需要知道变量与变量之间的关系。它们之间的关系一般来说可分为两类：一类叫做确定性关系，也叫函数关系，其特征是：一个变量随着其他变量的确定而确定。另一类叫相关关系，其特征是：变量之间的关系很难用一种精确的方法表示出来。但要注意这两种关系之间并不存在不可逾越的鸿沟。
回归分析就是处理变量之间的相关关系的一种数学方法，它是最常用的数理统计方法，能解决预测、控制、生产工艺优化等问题。
1、多元线性回归
（1）确定回归系数的点估计值：
b=regress( Y, X )
其中：
（2）求回归系数的点估计和区间估计、并检验回归模型：
其中bint为回归系数的区间估计，r为残差，rint为置信区间；stats用于检验回归模型的统计量，有三个数值：相关系数r2、F值、与F对应的概率p，相关系数r2越接近1，说明回归方程越显著；F > F1-α（k，n-k-1）时拒绝H0，F越大，说明回归方程越显著；与F对应的概率p 时拒绝H0，回归模型成立；alpha为显著性水平（缺省时为0.05）
（3）画出残差及其置信区间：rcoplot（r，rint）
例 测16名成年女子的身高与腿长所得数据如下：
| 身高 | 143 | 145 | 146 | 147 | 149 | 150 | 153 | 154 | 155 | 156 | 157 | 158 | 159 | 160 | 162 | 164 |
| 腿长 | 88 | 85 | 88 | 91 | 92 | 93 | 93 | 95 | 96 | 98 | 97 | 96 | 98 | 99 | 100 | 102 |
以身高x为横坐标，以腿长y为纵坐标，我们希望知道身高与腿长之间的关系。
解：1、输入数据：
2、回归分析及检验：
得结果：b =
即 ； 的置信区间为[-33.7017，1.5612], 的置信区间为[0.6047,0.834]; r2=0.9282, F=180.9531, p=0.0000
p<0.05, 可知回归模型 y=-16.073+0.7194x 成立。
3、
rcoplot(r,rint)
从残差图可以看出，除第二个数据外，其余数据的残差离零点均较近，且残差的置信区间均包含零点，这说明回归模型 y=-16.073+0.7194x能较好的符合原始数据，而第二个数据可视为异常点。
4、预测及作图：
z=b(1)+b(2)*x
plot(x,Y,'k+',x,z,'r')
2、多项式回归
rstool（x，y，’model’, alpha）
其中x 为n´m矩阵，y 为n维列向量，alpha为显著性水平（缺省时为0.05）；
model由下列4个模型中选择1个（用字符串输入，缺省时为线性模型）：
linear（线性）：
purequadratic（纯二次）：
interaction（交叉）：
quadratic（完全二次）：
例 设某商品的需求量与消费者的平均收入、商品价格的统计数据如下，建立回归模型，预测平均收入为1000、价格为6时的商品需求量。
| 需求量 | 100 | 75 | 80 | 70 | 50 | 65 | 90 | 100 | 110 | 60 |
| 收入 | 1000 | 600 | 1200 | 500 | 300 | 400 | 1300 | 1100 | 1300 | 300 |
| 价格 | 5 | 7 | 6 | 6 | 8 | 7 | 5 | 4 | 3 | 9 |
选择纯二次模型，即
输入命令：
x1=[1000 600 1200 500 300 400 1300 1100 1300 300];
x2=[5 7 6 6 8 7 5 4 3 9];
y=[100 75 80 70 50 65 90 100 110 60]';
x=[x1' x2'];
rstool(x,y,'purequadratic')
出现了图形框后，在左边图形下方的方框中输入1000，右边图形下方的方框中输入6。
则画面左边的“Predicted Y”下方的数据变为88.47981，即预测出平均收入为1000、价格为6时的商品需求量为88.4791。然后，在画面左下方的下拉式菜单中选”all”, 则beta、rmse和residuals都传送到Matlab工作区中。
在Matlab工作区中输入命令： beta, rmse
得结果：beta =
故回归模型为：
剩余标准差为4.5362, 说明此回归模型的显著性较好。
3、非线性回归
（1）回归：
（i）确定回归系数的命令：
beta为估计出的回归系数，r为残差，J为Jacobian矩阵；输入数据x、y分别为 矩阵和n维列向量，对一元非线性回归，x为n维列向量；model是事先用m-文件定义的非线性函数，beta0为回归系数的初值
（ii）非线性回归命令：nlintool（x，y，’model’, beta0，alpha）
（2）预测和预测误差估计：
[Y，DELTA]=nlpredci（’model’, x，beta，r，J）
求nlinfit 或nlintool所得的回归函数在x处的预测值Y及预测值的显著性为1-alpha的置信区间Y DELTA.
3745
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
